From john.maindonald at anu.edu.au  Mon May 19 07:16:09 2003
From: john.maindonald at anu.edu.au (john.maindonald@anu.edu.au)
Date: Mon May 19 06:16:22 2003
Subject: [Rd] predict.lm(), se=TRUE, with weights; erroneous SEs (PR#3043)
Message-ID: <200305190416.h4J4G92b006930@pubhealth.ku.dk>

<<insert bug report here>>
"xy" <-
     structure(list(x = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
                    y = c(1.036, 1.883, 3.11, 3.148, 4.226,
                    6.291, 7.312, 7.796, 9.563, 9.563),
                    w = c(0.067, 0.06, 1.588, 0.239, 0.015,
                    0.938, 0.041, 0.473, 0.483, 0.044)),
               .Names = c("x", "y", "w"),
               row.names = c("1", "2", "3", "4", "5",
               "6", "7", "8", "9", "10"), class = "data.frame")
 > xy.lm <- lm(y~x, weights=w, data=xy)  # Result is wrong
 > predict(xy.lm, se=T)$se
  [1] 0.05895647 0.04606261 0.19247671 0.06140987 0.01380818 0.11510861
  [7] 0.02865223 0.11991962 0.14786285 0.05331879
 > predict(xy.lm, newdata=xy, se=T)$se  # This gives the correct result
         1         2         3         4         5         6         7   
       8
0.2277687 0.1880498 0.1527401 0.1256144 0.1127433 0.1188520 0.1415033 
0.1743652
         9        10
0.2127578 0.2541874

Find the lines
             if (type != "terms") {
             if (missing(newdata))

The fix is to replace the next line (around line 47):
                 XRinv <- qr.Q(object$qr)[,  p1,  drop = FALSE]
     by
                 XRinv <- if(is.null(w)) qr.Q(object$qr)[,  p1,  drop = 
FALSE]
             else qr.Q(object$qr)[,  p1,  drop = FALSE]/sqrt(w)}

As far as I can discover, this bug is left over from extensive fixes
that I contributed to predict.lm() just after version 1.0.  If so, this
has never worked correctly.

I note that further down the code, following     if (interval != 
"none") {
the symbol w is used for another purpose.  To avoid  confusing
casual perusers of the code, this second use of w might be
changed to, e.g., halfwid

--please do not edit the information below--

Version:
  platform = powerpc-apple-darwin6.5
  arch = powerpc
  os = darwin6.5
  system = powerpc, darwin6.5
  status =
  major = 1
  minor = 7.0
  year = 2003
  month = 04
  day = 16
  language = R

Search Path:
  .GlobalEnv, package:methods, package:ctest, package:mva, 
package:modreg, package:nls, package:ts, Autoloads, package:base

John Maindonald             email: john.maindonald@anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Bioinformation Science, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.


[Rd] predict.lm(), se=TRUE, with weights; erroneous SEs (PR#3043)
john.maindonald@anu.edu.au