# Linear Models #-------------- # Simulated Example #------------------ y <- rep(1:4, each=25) # note this is numeric, not a factor z <- 5 + 2*y + rnorm(100) # first response var a <- rep(1:5, times=20) za <- 4 + 3*y +2*a +rnorm(100) # a second response var table(y,a) # balanced simulation ds <- data.frame(y,a,z,za) # Some model specification #------------------------- (reg1 <- lm(z ~ 1 + y + a, data=ds)) # Explicit intercept lm(z ~ y + a, data=ds) # Implicit intercept lm(z ~ a + 1 + y, data=ds) # Somewhat arbitrary order (reg2 <- lm(z ~ 0 + y + a, data=ds)) # Zero intercept lm(z ~ 5 + y + a, data=ds) # ERROR, 5 is not a valid term ds\$int <- 5 # recycling! lm(z ~ offset(int) + y + a, data=ds) # Offset the model lm(z ~ offset(int) + 0 + y + a, data=ds) # Offset the zero dose model lm(z ~ 1 + y + a - 1, data=ds) # Remove a term lm(z ~ (1 + y) + a, data=ds) # Group terms, visually lm(z ~ I(1 + y) + a, data=ds) # A calculated term, "protected" lm(z ~ I(y^2), data=ds) # another lm(log(z) ~ y + a, data=ds) # another lm(z ~ 1 + (y + a), data=ds) # Another grouping lm(z ~ 1 + I(y + a), data=ds) # Another calculated term lm(z ~ 1 + 7*(y + a), data=ds) # ERROR, confuses crossing and multiplication lm(z ~ 1 + I(7*(y + a)), data=ds) # Calculated term lm(z ~ y*a, data=ds) # crossed term lm(z ~ y + a + y:a, data=ds) # explicit main and interaction terms m <- formula(z~y) # formulas can be objects str(m) lm(m) plot(m) # Extracting information from a model fit #---------------------------------------- reg1 <- lm(z ~ y, data=ds) str(reg1) summary(reg1) anova(reg1) model1 <- lm(za ~ y, data=ds) model2 <- lm(za ~ a + y, data=ds) anova(model1, model2) coef(model1) plot(za ~ y) abline(coef(model1)) confint(model2) #################### # Example Regression #################### library(lattice) xyplot(dist ~ speed, data=cars, type=c("p", "r")) regcar <- lm(dist ~ speed, data=cars) summary(regcar) # Simulated ANOVA #---------------- yf <- as.factor(y) bwplot(z~yf) aov1 <- lm(z ~ yf) summary(aov1) anova(aov1) anova(reg1) xyplot(residuals(reg1)~predict(reg1)) xyplot(residuals(aov1)~predict(aov1)) ############# # Example AOV ############# aov2 <- lm(weight ~ feed, chickwts) summary(aov2) plot(chickwts) xyplot(weight ~ feed, data=chickwts) xyplot(residuals(aov2) ~ predict(aov2), group=chickwts\$feed, auto.key=TRUE, xlab="Predicted values", ylab="Residuals")