[Rd] Wrong SEs in predict.lm(..., type="terms") (PR#528)

[john.maindonald@anu.edu.au]

>From e980153 Tue Apr 25 14:42:27 2000
To: r-help@stat.math.ethz.ch
Subject: Wrong SEs in predict.lm(..., type="terms")

For what it is worth, I am using RW-1.0.0 under Windows 98.
I submitted this earlier to r-help.  There is one change
below to my proposed corrected code:

predict.lm(..., type="terms") gives wrong standard errors.

Below, I have provided what I believe are the necessary fixes.
However, there are subtleties, and the code needs careful
checking.  Some of the looping is surely not necessary, but
it is surely best to begin with the minimum necessary changes.

My tests, including checks against S-PLUS, have extended to fitting
spline curves.  I have not tested the code with any example where
the pivot sequence does not run sequentially from 1 to p1.

df<-data.frame(x=1:9,y=c(2,3,6,4,8,10,12,14,15))
> df.lm_lm(y~x,data=df)
> as.vector(predict(df.lm,se=T,type="terms")$se.fit)
[1] 0.9904885 0.7428664 0.4952442 0.2476221 0.0000000 0.2476221 0.4952442
[8] 0.7428664 0.9904885

Here is what one should get:
> abs((df$x-mean(df$x)))*summary(df.lm)$coef[2,2]
[1] 0.5769836 0.4327377 0.2884918 0.1442459 0.0000000 0.1442459 0.2884918
[8] 0.4327377 0.5769836

Here is the relevant section of corrected code:


    if (type=="terms"){
      asgn <- attrassign(object)
      beta<-coef(object)
      hasintercept<-attr(tt,"intercept")>0
      if (hasintercept)
        asgn$"(Intercept)"<-NULL
      nterms<-length(asgn)
      predictor<-matrix(ncol=nterms,nrow=NROW(X))
      dimnames(predictor)<-list(rownames(X),names(asgn))

      if (se.fit){
        ip<-matrix(ncol=nterms,nrow=NROW(X))
        dimnames(ip)<-list(rownames(X),names(asgn))
      }
      if (hasintercept){     # This was omitted in the code on r-help
        avX<-apply(X,2,mean)  
        X<-sweep(X,2,avX)    # We'd best sweep out column means
                             # before we start
      }
      for (i in seq(1,nterms,length=nterms)){
        ii<-piv[asgn[[i]]]
	        predictor[,i]<-X[,ii,drop=F]%*%(beta[ii])
        if (se.fit){
          ii2_asgn[[i]]      # X[,piv] matches rows & columns of R
          vci<-R[ii2,ii2]*res.var
          for(j in (1:NROW(X))){
            xi<-X[j,ii,drop=F]  # Do not multiply by beta[ii]
           ip[j,i]<-sum(xi%*% vci %*%t(xi))

          }
        }
      }

Here is the existing section of the (1.0.0) code:

    if (type=="terms"){
      asgn <- attrassign(object)
      beta<-coef(object)
      hasintercept<-attr(tt,"intercept")>0
      if (hasintercept)
        asgn$"(Intercept)"<-NULL
      nterms<-length(asgn)
      predictor<-matrix(ncol=nterms,nrow=NROW(X))
      dimnames(predictor)<-list(rownames(X),names(asgn))
      if (se.fit){
        ip<-matrix(ncol=nterms,nrow=NROW(X))
        dimnames(ip)<-list(rownames(X),names(asgn))
      }
      for (i in seq(1,nterms,length=nterms)){
        if (hasintercept)  # Redundant
          i0<-1            # Redundant
        else               # Redundant
          i0<-NULL         # Redundant
        ii<-piv[asgn[[i]]]
        predictor[,i]<-X[,ii,drop=F]%*%(beta[ii])  # Should be centred
        X[,ii]<-X[,ii]-mean(X[,ii]) # Columns must be centred individually
        if (se.fit){
          vci<-R[ii,ii]*res.var    # The pivot should not be applied to R
          for(j in (1:NROW(X))){
            xi<-X[j,ii,drop=F]*(beta[ii])  # Do not Xply by beta
            ip[j,i]<-sum(xi%*% vci %*%t(xi))
          }
        }
      }


John Maindonald               email : john.maindonald@anu.edu.au        
Statistical Consulting Unit,  phone : (6249)3998        
c/o CMA, SMS,                 fax   : (6249)5549  
John Dedman Mathematical Sciences Building
Australian National University
Canberra ACT 0200
Australia


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._