R-alpha: predict.lm -- who ..?

[Peter Dalgaard BSA]

Peter Dalgaard BSA <p.dalgaard@kubism.ku.dk> writes:

> Hum, maybe my version was closer after all. Thing was that I tried to
...
> However I did get the SE's done and soforth. Here's my effort:

OK. A little thought and some cleanup later, I have this one. I'm a
little reluctant to go completely Splus compatible (return a list with
components $fit $se.fit $residual.scale $df if se.fit=T, otherwise
just the predicted values) and losing the confidence/prediction
intervals. Opinions?

"predict.lm" <-
  function (object, newdata = model.frame(object)) 
{
  form<-formula(object)
  form<-as.call(list(form[[1]],form[[3]])) # x~a+b becomes ~a+b
  class(form)<-"formula"
  X <- model.matrix(form,newdata)
  n <- NROW(object$qr$qr)
  p <- object$rank
  p1 <- 1:p
  piv <- object$qr$pivot[p1]
  r <- resid(object)
  f <- fitted(object)
  w <- weights(object)
  if (is.null(w)) 
    rss <- sum(r^2)
  else 
    rss <- sum(r^2*w)
  R <- chol2inv(object$qr$qr[p1, p1, drop = FALSE])
  est <- object$coefficients[piv]
  predictor <- c(X[,piv] %*% est)
  ip <- real(NROW(X))
  resvar <- rss/(n - p)
  vcov <- resvar * R
  for (i in (1:NROW(X))) {
    xi <- X[i,piv]
    ip[i] <- xi %*% vcov %*% xi
  }
  stderr1 <- sqrt(ip)
  stderr2 <- sqrt(resvar + ip)
  tt <- qt(0.975, n - p)
  conf.l <- predictor - tt * stderr1
  conf.u <- predictor + tt * stderr1
  pred.l <- predictor - tt * stderr2
  pred.u <- predictor + tt * stderr2
  data.frame(predictor=predictor, conf.l=conf.l, conf.u=conf.u,
	pred.l=pred.l,pred.u=pred.u,row.names=rownames(newdata))
}
-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-