R devel - Apr 2008 - Bug in ci.plot(HH Package) (PR#11163)

Full_Name: Yasuhiro Nakajima
Version: 2.6.1
OS: WinXP SP2
Submission from: (NULL) (202.237.255.13)


Dear all,

I noticed the following behaviour of ci.plot in HH Package(ver.2.1-9):
> library(HH)
> data(women, package="datasets")
> attach(women)
> ft <- lm(height~weight)
> windows()
> ci.plot(ft,conf.level=0.95)
> windows()
> ci.plot(ft,conf.level=0.999)
I tried to change the confidence interval, but I couldn't.
CI was "always" 0.95. 

I have found wrong argument in predict function in ci.plot.
>>   predict(....., conf.level=conf.level)
I think the correct is "level=conf.level".

Is that right?


--- ci.plot source code ---
"ci.plot" <-
function(lm.object, ...)
  UseMethod("ci.plot")

"ci.plot.lm" <-
function(lm.object,
         xlim=range(data[, x.name]),
         newdata,
         conf.level=.95,
         data=model.frame(lm.object),
         newfit,
         ylim=range(newfit$pi.fit),
         pch=16,
         main.cex=1,
         main=list(paste(100*conf.level,
           "% confidence and prediction intervals for ",
           substitute(lm.object), sep=""), cex=main.cex), ...
         ) {
  formula.lm <- formula(lm.object)
  x.name <- as.character(formula.lm[[3]])
  missing.xlim <- missing(xlim)       ## R needs this
  missing.newdata <- missing(newdata) ## R needs this
  if.R(s={
    ## Save a copy of the data.frame in frame=0 to put it where
    ## model.frame.lm needs to find it when the example data is
    ## run through Splus CMD check.
    my.data.name <- as.character(lm.object$call$data)
    if (length(my.data.name)==0)
      stop("Please provide an lm.object calculated with an explicit
'data=my.data.frame' argument.")
    undo.it <- (!is.na(match(my.data.name, objects(0))))
    if (undo.it) old.contents <- get(my.data.name, frame=0)
    my.data <- try(get(my.data.name))
    if (class(my.data)=="Error")
      my.data <- try(get(my.data.name, frame=sys.parent()))
    if (class(my.data)=="Error")
      stop("Please send me an email with a reproducible situation that got
you
here. (rmh at temple.edu)")
    assign(my.data.name, my.data, frame=0)
  },r={})
  default.newdata <- data.frame(seq(xlim[1], xlim[2], length=51))
  names(default.newdata) <- x.name
  if (missing.xlim) xlim <- xlim + diff(xlim)*c(-.02,.02) ## needed
  if (missing.newdata) {
    newdata <- default.newdata
    newdata.x <- numeric()
  }
  else {
    if (is.na(match(x.name, names(newdata))))
      stop(paste("'newdata' must be a data.frame containing a
column named '",
                 x.name, "'", sep=""))
    if (missing.xlim)
      xlim=range(xlim, newdata[[x.name]])
    newdata.x <- as.data.frame(newdata)[,x.name]
    newdata <- rbind(as.data.frame(newdata)[,x.name, drop=FALSE],
                     default.newdata)
    newdata <- newdata[order(newdata[,x.name]), , drop=FALSE]
  }
  if (missing.xlim) xlim <- xlim + diff(xlim)*c(-.02,.02) ## repeat is needed
  if (missing(newfit)) newfit <-
    if.R(s={
      
      prediction <-
        predict(lm.object, newdata=newdata,
                se.fit=TRUE, ci.fit=TRUE, pi.fi=TRUE,
                conf.level=conf.level)
      {
        ## restore frame=0
        if (undo.it) assign(my.data.name, old.contents, frame=0)
        else remove(my.data.name, frame=0)
      }
      prediction
    }
         ,r={
           new.p <-
             predict(lm.object, newdata=newdata,
                     se.fit=TRUE, conf.level=conf.level,
                     interval = "prediction")
           new.c <-
             predict(lm.object, newdata=newdata,
                     se.fit=TRUE, conf.level=conf.level,
                     interval = "confidence")
           tmp <- new.p
           tmp$ci.fit <- new.c$fit[,c("lwr","upr"),
drop=FALSE]
           dimnames(tmp$ci.fit)[[2]] <-
c("lower","upper")
           attr(tmp$ci.fit,"conf.level") <- conf.level
           tmp$pi.fit <- new.p$fit[,c("lwr","upr"),
drop=FALSE]
           dimnames(tmp$pi.fit)[[2]] <-
c("lower","upper")
           attr(tmp$pi.fit,"conf.level") <- conf.level
           tmp$fit <- tmp$fit[,"fit", drop=FALSE]
           tmp
         })
  tpgsl <- trellis.par.get("superpose.line")
  tpgsl <- Rows(tpgsl, 1:4)
  tpgsl$col[1] <- 0
  xyplot(formula.lm, data=data, newdata=newdata, newfit=newfit,
         newdata.x=newdata.x,
         xlim=xlim, ylim=ylim, pch=pch,
         panel=function(..., newdata.x) {
           panel.ci.plot(...)
           if (length(newdata.x) > 0)
             panel.rug(x=newdata.x)
         },
         main=main,
         key=list(border=TRUE,
           space="right",
           text=list(c("observed","fit","conf
int","pred int")),
           points=list(
             pch=c(pch,32,32,32),
             col=c(trellis.par.get("plot.symbol")$col, tpgsl$col[2:4])
             ),
           lines=tpgsl),
         ...)
}

Please report bugs in contributes packages to the package maintainer 
(CCing), not to R-bugs.

Best,
Uwe Ligges


nakajima at hokkaido-iri.go.jp wrote:
> Full_Name: Yasuhiro Nakajima
> Version: 2.6.1
> OS: WinXP SP2
> Submission from: (NULL) (202.237.255.13)
> 
> 
> Dear all,
> 
> I noticed the following behaviour of ci.plot in HH Package(ver.2.1-9):
> 
>> library(HH)
>> data(women, package="datasets")
>> attach(women)
>> ft <- lm(height~weight)
>> windows()
>> ci.plot(ft,conf.level=0.95)
>> windows()
>> ci.plot(ft,conf.level=0.999)
> 
> I tried to change the confidence interval, but I couldn't.
> CI was "always" 0.95. 
> 
> I have found wrong argument in predict function in ci.plot.
> 
>>>   predict(....., conf.level=conf.level)
> 
> I think the correct is "level=conf.level".
> 
> Is that right?
> 
> 
> --- ci.plot source code ---
> "ci.plot" <-
> function(lm.object, ...)
>   UseMethod("ci.plot")
> 
> "ci.plot.lm" <-
> function(lm.object,
>          xlim=range(data[, x.name]),
>          newdata,
>          conf.level=.95,
>          data=model.frame(lm.object),
>          newfit,
>          ylim=range(newfit$pi.fit),
>          pch=16,
>          main.cex=1,
>          main=list(paste(100*conf.level,
>            "% confidence and prediction intervals for ",
>            substitute(lm.object), sep=""), cex=main.cex), ...
>          ) {
>   formula.lm <- formula(lm.object)
>   x.name <- as.character(formula.lm[[3]])
>   missing.xlim <- missing(xlim)       ## R needs this
>   missing.newdata <- missing(newdata) ## R needs this
>   if.R(s={
>     ## Save a copy of the data.frame in frame=0 to put it where
>     ## model.frame.lm needs to find it when the example data is
>     ## run through Splus CMD check.
>     my.data.name <- as.character(lm.object$call$data)
>     if (length(my.data.name)==0)
>       stop("Please provide an lm.object calculated with an explicit
> 'data=my.data.frame' argument.")
>     undo.it <- (!is.na(match(my.data.name, objects(0))))
>     if (undo.it) old.contents <- get(my.data.name, frame=0)
>     my.data <- try(get(my.data.name))
>     if (class(my.data)=="Error")
>       my.data <- try(get(my.data.name, frame=sys.parent()))
>     if (class(my.data)=="Error")
>       stop("Please send me an email with a reproducible situation that
got you
> here. (rmh at temple.edu)")
>     assign(my.data.name, my.data, frame=0)
>   },r={})
>   default.newdata <- data.frame(seq(xlim[1], xlim[2], length=51))
>   names(default.newdata) <- x.name
>   if (missing.xlim) xlim <- xlim + diff(xlim)*c(-.02,.02) ## needed
>   if (missing.newdata) {
>     newdata <- default.newdata
>     newdata.x <- numeric()
>   }
>   else {
>     if (is.na(match(x.name, names(newdata))))
>       stop(paste("'newdata' must be a data.frame containing a
column named '",
>                  x.name, "'", sep=""))
>     if (missing.xlim)
>       xlim=range(xlim, newdata[[x.name]])
>     newdata.x <- as.data.frame(newdata)[,x.name]
>     newdata <- rbind(as.data.frame(newdata)[,x.name, drop=FALSE],
>                      default.newdata)
>     newdata <- newdata[order(newdata[,x.name]), , drop=FALSE]
>   }
>   if (missing.xlim) xlim <- xlim + diff(xlim)*c(-.02,.02) ## repeat is
needed
>   if (missing(newfit)) newfit <-
>     if.R(s={
>       
>       prediction <-
>         predict(lm.object, newdata=newdata,
>                 se.fit=TRUE, ci.fit=TRUE, pi.fi=TRUE,
>                 conf.level=conf.level)
>       {
>         ## restore frame=0
>         if (undo.it) assign(my.data.name, old.contents, frame=0)
>         else remove(my.data.name, frame=0)
>       }
>       prediction
>     }
>          ,r={
>            new.p <-
>              predict(lm.object, newdata=newdata,
>                      se.fit=TRUE, conf.level=conf.level,
>                      interval = "prediction")
>            new.c <-
>              predict(lm.object, newdata=newdata,
>                      se.fit=TRUE, conf.level=conf.level,
>                      interval = "confidence")
>            tmp <- new.p
>            tmp$ci.fit <- new.c$fit[,c("lwr","upr"),
drop=FALSE]
>            dimnames(tmp$ci.fit)[[2]] <-
c("lower","upper")
>            attr(tmp$ci.fit,"conf.level") <- conf.level
>            tmp$pi.fit <- new.p$fit[,c("lwr","upr"),
drop=FALSE]
>            dimnames(tmp$pi.fit)[[2]] <-
c("lower","upper")
>            attr(tmp$pi.fit,"conf.level") <- conf.level
>            tmp$fit <- tmp$fit[,"fit", drop=FALSE]
>            tmp
>          })
>   tpgsl <- trellis.par.get("superpose.line")
>   tpgsl <- Rows(tpgsl, 1:4)
>   tpgsl$col[1] <- 0
>   xyplot(formula.lm, data=data, newdata=newdata, newfit=newfit,
>          newdata.x=newdata.x,
>          xlim=xlim, ylim=ylim, pch=pch,
>          panel=function(..., newdata.x) {
>            panel.ci.plot(...)
>            if (length(newdata.x) > 0)
>              panel.rug(x=newdata.x)
>          },
>          main=main,
>          key=list(border=TRUE,
>            space="right",
>            text=list(c("observed","fit","conf
int","pred int")),
>            points=list(
>              pch=c(pch,32,32,32),
>              col=c(trellis.par.get("plot.symbol")$col,
tpgsl$col[2:4])
>              ),
>            lines=tpgsl),
>          ...)
> }
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

Thank you for the error report.  I have made the correction
which will be in HH_2.1-11 which will be posted in a few days.

Rich