The R Project SVN R

###-- Linear Models, basic functionality -- weights included.

## From John Maindonald :
roller <- data.frame(
 weight     = c(1.9,  3.1,  3.3,  4.8,  5.3,  6.1,  6.4,  7.6,  9.8, 12.4),
 depression = c(  2,  1,    5,    5,   20,   20,   23,   10,   30,   25))

roller.lmu  <-  lm(weight~depression, data=roller)
roller.lsfu <- lsfit(roller$depression, roller$weight)

roller.lsf  <- lsfit(roller$depression, roller$weight, wt = 1:10)
roller.lsf0 <- lsfit(roller$depression, roller$weight, wt = 0:9)
roller.lm  <-  lm(weight~depression, data=roller, weights= 1:10)
roller.lm0 <-  lm(weight~depression, data=roller, weights= 0:9)
roller.lm9 <-  lm(weight~depression, data=roller[-1,],weights= 1:9)
roller.glm <- glm(weight~depression, data=roller, weights= 1:10)
roller.glm0<- glm(weight~depression, data=roller, weights= 0:9)

predict(roller.glm0, type="terms")# failed till 2003-03-31

stopifnot(exprs = {
    all.equal(residuals(roller.glm0, type = "partial"),
              residuals(roller.lm0,  type = "partial"), tol = 1e-14)  # 4.0e-16
    all.equal(deviance(roller.lm),
              deviance(roller.glm), tol = 1e-14) # 2.4e-16
    all.equal(weighted.residuals(roller.lm),
              residuals         (roller.glm), tol = 2e-14) # 9.17e-16
    all.equal(deviance(roller.lm0),
              deviance(roller.glm0), tol = 1e-14) # 2.78e-16
    all.equal(weighted.residuals(roller.lm0, drop=FALSE),
              residuals         (roller.glm0),  tol = 2e-14) # 6.378e-16
})

(im.lm0 <- influence.measures(roller.lm0))

stopifnot(exprs = {
    all.equal(unname(im.lm0 $ infmat),
              unname(cbind(  dfbetas      (roller.lm0)
                         , dffits       (roller.lm0)
                         , covratio     (roller.lm0)
                          ,cooks.distance(roller.lm0)
                          ,lm.influence  (roller.lm0)$hat)
                     ))
    all.equal(rstandard(roller.lm9),
              rstandard(roller.lm0),tolerance = 1e-14)
    all.equal(rstudent(roller.lm9),
              rstudent(roller.lm0),tolerance = 1e-14)
    all.equal(rstudent(roller.lm),
              rstudent(roller.glm))
    all.equal(cooks.distance(roller.lm),
              cooks.distance(roller.glm))

    all.equal(summary(roller.lm0)$coefficients,
              summary(roller.lm9)$coefficients, tolerance = 1e-14)
    all.equal(print(anova(roller.lm0), signif.st=FALSE),
                    anova(roller.lm9), tolerance = 1e-14)
})


###  more regression tests for lm(), glm(), etc :

## moved from ?influence.measures:
lm.SR <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
(IM <- influence.measures(lm.SR))
summary(IM)
## colnames will differ in the next line
stopifnot(
    all.equal(dfbetas(lm.SR), IM$infmat[, 1:5], check.attributes = FALSE,
              tolerance = 1e-12)
)
signif(dfbeta(lm.SR), 3)
covratio (lm.SR)

## Multivariate lm ("mlm") --- Example from  ?SSD
reacttime <- matrix(c(
420, 420, 480, 480, 600, 780,
420, 480, 480, 360, 480, 600,
480, 480, 540, 660, 780, 780,
420, 540, 540, 480, 780, 900,
540, 660, 540, 480, 660, 720,
360, 420, 360, 360, 480, 540,
480, 480, 600, 540, 720, 840,
480, 600, 660, 540, 720, 900,
540, 600, 540, 480, 720, 780,
480, 420, 540, 540, 660, 780),
ncol = 6, byrow = TRUE,
dimnames = list(subj = 1:10,
              cond = c("deg0NA", "deg4NA", "deg8NA",
                       "deg0NP", "deg4NP", "deg8NP")))
mlmfit <- lm(reacttime ~ 1)
ImMLM <- influence.measures(mlmfit)## fails in R <= 3.5.1
## and the print() and summary() methods had failed additionally:
oo <- capture.output(ImMLM) # now ok
summary(ImMLM) # "ok"



## predict.lm(.)

all.equal(predict(roller.lm,                 se.fit=TRUE)$se.fit,
          predict(roller.lm, newdata=roller, se.fit=TRUE)$se.fit, tolerance =  1e-14)