## ----opts, echo=F, include=F--------------------------------------------- options(useFancyQuotes=FALSE) ## ----mods, echo=F, results='asis'---------------------------------------- library(car) library(apsrtable) data(Prestige) Prestige$income <- Prestige$income/1000 mod1 <- lm(prestige ~ income, data=Prestige) mod2 <- lm(prestige ~ education, data=Prestige) mod3 <- lm(prestige ~ women, data=Prestige) mod4 <- lm(prestige ~ income + education + women, data=Prestige) apsrtable(mod1, mod2, mod3, mod4, digits=3, caption="OLS Regressions of Occupational Prestige", notes=list("Main entries are OLS coefficients", se.note, "$^{*}$ indicates singificance at $p < 0.05$ (two-tailed)"), coef.names = c("Intercept", "Mean income of incumbents (in \\$1000)", "Mean incumbent years of education", "Percentage of women incumbents")) ## ----partial, echo=T----------------------------------------------------- EP <- lm(prestige ~ education + women, data=Prestige)$residuals EI <- lm(income ~ education + women, data=Prestige)$residuals partial.mod <- lm(EP ~ EI) coef(mod4) coef(partial.mod) ## ----vecs, echo=T-------------------------------------------------------- u <- c(3,3,3,3) v <- c(1,2,3,4) u+v u-v ## ----inner, echo=T, include=T-------------------------------------------- u %*% v ## ----outer, echo=T------------------------------------------------------- u %o% v outer(u, v, "*") ## ----mats, echo=T-------------------------------------------------------- x <- matrix(c(1,3,2,4), ncol=2, nrow=2) y <- matrix(c(5,7,6,8), ncol=2, nrow=2) x + y ## ----matmult, echo=T----------------------------------------------------- x %*% y ## ----matinv, echo=T------------------------------------------------------ x solve(x) x %*% solve(x) ## ----mod, echo=T--------------------------------------------------------- library(car) data(Duncan) mod <- lm(prestige ~ income + education, data=Duncan) summary(mod) ## ----betavar, echo=T----------------------------------------------------- X <- with(Duncan, cbind(1,income, education)) y <- matrix(Duncan[["prestige"]], ncol=1) b <- solve(t(X)%*%X)%*%t(X)%*%y b coef(mod) e <- matrix(y - X%*%b, ncol=1) Vb <- c((t(e)%*%e)/(nrow(X) - 2 - 1))* solve(t(X)%*%X) Vb vcov(mod) ## ----ftest, echo=T------------------------------------------------------- restricted.mod <- lm(prestige ~ 1, data=Duncan) anova(restricted.mod, mod, test="F") xtx <- solve(t(X)%*%X) s2e <- (t(e) %*% e)/(nrow(X)-3) F0 <- (t(b[2:3]) %*% solve(xtx[2:3,2:3]) %*% b[2:3])/(2*s2e) F0 pf(F0, 2, nrow(X)-3, lower.tail=F) ## ----glht, echo=T-------------------------------------------------------- linearHypothesis(mod, c("education = .5", "income = .5")) ## ----vapred_hide, echo=F, results="hide"--------------------------------- library(car) data(Duncan) X <- with(Duncan, cbind(1,income, education)) y <- matrix(Duncan[["prestige"]], ncol=1) b <- solve(t(X)%*%X)%*%t(X)%*%y e <- matrix(y - X%*%b, ncol=1) s2e <- c((t(e)%*%e)/(nrow(X) - 2 - 1)) Vb <- s2e* solve(t(X)%*%X) ## ----varpred, echo=T----------------------------------------------------- x0 <- cbind(c(1, 21, 26), c(1,64,84), c(1,41.87, 52.56)) pred0 <- t(x0) %*% b pred0 v0 <- as.numeric(s2e) * (1 + t(x0) %*% solve(t(X)%*% X) %*% x0) sqrt(diag(v0)) vhat0 <- as.numeric(s2e) * (t(x0) %*% solve(t(X)%*% X) %*% x0) sqrt(diag(vhat0)) ## ----varpred2, echo=T---------------------------------------------------- mod <- lm(prestige ~ income + education, data=Duncan) newdat <- t(x0) colnames(newdat) <- c("int", "income", "education") newdat <- as.data.frame(newdat) predict(mod, newdat, interval="confidence") predict(mod, newdat, interval="prediction")