From p.dalgaard@kubism.ku.dk  Mon Sep 15 22:57:02 1997
From: p.dalgaard@kubism.ku.dk (Peter Dalgaard BSA)
Date: 15 Sep 1997 23:57:02 +0200
Subject: R-alpha: predict.lm -- who ..?
In-Reply-To: Peter Dalgaard BSA's message of 15 Sep 1997 22:50:58 +0200
References: <199709151016.WAA19837@stat1.stat.auckland.ac.nz>
 <x24t7mfgi5.fsf@bush.kubism.ku.dk>
Message-ID: <x2yb4yckb5.fsf@bush.kubism.ku.dk>

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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-


R-alpha: predict.lm -- who ..?
Peter Dalgaard BSA