[Rd] Bug in ci.plot(HH Package) (PR#11163)
ligges at statistik.tu-dortmund.de
ligges at statistik.tu-dortmund.de
Mon Apr 14 10:10:15 CEST 2008
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
More information about the R-devel
mailing list