#--------------------------------- # Three simple models #--------------------------------- # The main points with respect to programming R are: # (i) lm() and other modeling functions return # "lm" objects, aka "fit" objects, # (ii) the language for specifying formulas is a # little different than the language for doing # math, and the two can be used simultaneously # (iii) so-called "extractor" functions, aka "accessor" # functions work on "lm" objects, and # (iv) many generic functions have "lm" methods. #--------------------------------- # Example 1 - One sample t-test # and its related linear model #--------------------------------- y <- rnorm(25) + 5 ### y <- rnorm(25, mean=5) hist(y) t.test(y, mu=5) t.test(y) # the mu=0 model fit1 <- lm(y ~ 1) # In the related linear model, we estimate mu fit1 # from print.lm() coef(fit1) confint(fit1) summary(fit1) # from summary.lm() str(fit1) plot(fit1) # from plot.lm() qqnorm(residuals(fit1)) qqline(residuals(fit1)) #--------------------------------- # Example 2 - Two group t-test # paired t-test, and related # linear models #--------------------------------- sleep ?sleep with(sleep, tapply(extra, group, mean)) with(sleep, tapply(extra, group, sd)) t.test(extra ~ group, data=sleep, var.equal=TRUE) summary(fit2 <- lm(extra ~ group, data=sleep)) summary(fit3 <- update(fit2, . ~ . -group)) #plot(fit2) # These data are actually paired/repeated measures plot(extra ~ ID, data=sleep) plot(extra ~ as.integer(ID), data=sleep) points(as.integer(sleep\$ID)[1:10], sleep\$extra[1:10], col="red") t.test(extra ~ group, data=sleep, paired=TRUE) # Doing math within a model formula summary(fit2a <- lm(extra[group==1]-extra[group==2] ~ 1, data=sleep)) # A similar model summary(fit2b <- lm(extra ~ ID+group, data=sleep)) # How do the residuals from fit2a relate to the residuals in # fit2b? abs(abs(residuals(fit2a)) - 2*abs(residuals(fit2b))) abs(abs(residuals(fit2a)) - 2*abs(residuals(fit2b))) < .Machine\$double.eps #--------------------------------- # Example 3 - Classic Bivariate # Regression #--------------------------------- library(faraway) ?nels88 str(nels88) cor(nels88\$math, nels88\$ses) plot(math ~ ses, data=nels88) summary(fit3 <- lm(math ~ ses, data=nels88)) plot(math ~ ses, data=nels88) abline(fit3) segments(x0=nels88\$ses, y0=predict(fit3), y1=nels88\$math) #--------------------------------- # Example 4 - A Factor and a # Covariate #--------------------------------- x <- rep(c(0, 1), 20) z <- runif(40, -5, 5) y <- 3 + 6*x + 2.5*z + rnorm(40, sd=2.5) plot(y ~ z) points(z[x==1], y[x==1], col="blue") plot(y ~ factor(x)) fit4 <- lm(y ~ 1 + x +z) fit4 summary(fit4) anova(fit4) plot(y ~ z) points(predict(fit4) ~ z, pch="+") segments(z, predict(fit4), y1=y) fit5 <- lm(y ~ x * z) # crossed summary(fit5) plot(y ~ z) points(predict(fit5) ~ z, pch="+") fit6 <- lm(math ~ (sex + race)*ses, data=nels88) summary(fit6) anova(fit6) fit7 <- lm(math ~ (sex + race)*ses - sex:ses - sex, data=nels88) anova(fit7) fit7a <- update(fit6, .~.-sex-sex:ses) anova(fit7a) summary(fit7a) anova(fit7a, fit6) plot(math ~ ses, data=nels88) points(predict(fit7) ~ ses, data=nels88, pch="+") levels(nels88\$race) raceb <- relevel(nels88\$race, "Black") summary(update(fit7, .~.-race+raceb))