###-- 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)