### R code from vignette source '/home/zeileis/svn/projects/AER/Course/Ch-LinearRegression.Rnw' ################################################### ### code chunk number 1: setup ################################################### options(prompt = "R> ", continue = "+ ", width = 64, digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) library("AER") set.seed(1071) ################################################### ### code chunk number 2: journals-data ################################################### data("Journals", package = "AER") journals <- Journals[, c("subs", "price")] journals$citeprice <- Journals$price/Journals$citations summary(journals) ################################################### ### code chunk number 3: journals-lm (eval = FALSE) ################################################### ## plot(log(subs) ~ log(citeprice), data = journals) ## jour_lm <- lm(log(subs) ~ log(citeprice), data = journals) ## abline(jour_lm) ################################################### ### code chunk number 4: journals-lm1 ################################################### plot(log(subs) ~ log(citeprice), data = journals) jour_lm <- lm(log(subs) ~ log(citeprice), data = journals) abline(jour_lm) ################################################### ### code chunk number 5: class_names_lm ################################################### class(jour_lm) names(jour_lm) jour_lm$rank ################################################### ### code chunk number 6: journals-summary ################################################### summary(jour_lm) ################################################### ### code chunk number 7: journals-summary2 ################################################### jour_slm <- summary(jour_lm) class(jour_slm) names(jour_slm) jour_slm$coefficients ################################################### ### code chunk number 8: journals-anova ################################################### anova(jour_lm) ################################################### ### code chunk number 9: journals-coef ################################################### coef(jour_lm) ################################################### ### code chunk number 10: journals-confint ################################################### confint(jour_lm, level = 0.95) ################################################### ### code chunk number 11: journals-predict ################################################### predict(jour_lm, newdata = data.frame(citeprice = 2.11), interval = "confidence") predict(jour_lm, newdata = data.frame(citeprice = 2.11), interval = "prediction") ################################################### ### code chunk number 12: journals-predict-plot (eval = FALSE) ################################################### ## lciteprice <- seq(from = -6, to = 4, by = 0.25) ## jour_pred <- predict(jour_lm, interval = "prediction", ## newdata = data.frame(citeprice = exp(lciteprice))) ## plot(log(subs) ~ log(citeprice), data = journals) ## lines(jour_pred[, 1] ~ lciteprice, col = 1) ## lines(jour_pred[, 2] ~ lciteprice, col = 1, lty = 2) ## lines(jour_pred[, 3] ~ lciteprice, col = 1, lty = 2) ################################################### ### code chunk number 13: journals-predict-plot1 ################################################### lciteprice <- seq(from = -6, to = 4, by = 0.25) jour_pred <- predict(jour_lm, interval = "prediction", newdata = data.frame(citeprice = exp(lciteprice))) plot(log(subs) ~ log(citeprice), data = journals) lines(jour_pred[, 1] ~ lciteprice, col = 1) lines(jour_pred[, 2] ~ lciteprice, col = 1, lty = 2) lines(jour_pred[, 3] ~ lciteprice, col = 1, lty = 2) ################################################### ### code chunk number 14: journals-plot_lm (eval = FALSE) ################################################### ## plot(jour_lm) ################################################### ### code chunk number 15: journals-plot-individual (eval = FALSE) ################################################### ## plot(jour_lm, which = 2) ################################################### ### code chunk number 16: journals-plot_lm1 ################################################### par(mfrow = c(2, 2)) plot(jour_lm) par(mfrow = c(1, 1)) ################################################### ### code chunk number 17: linear_hypothesis1 ################################################### linearHypothesis(jour_lm, "log(citeprice) = -0.5") ################################################### ### code chunk number 18: linear_hypothesis2 (eval = FALSE) ################################################### ## linearHypothesis(jour_lm, hypothesis.matrix = c(0, 1), rhs = -0.5) ################################################### ### code chunk number 19: CPS-data ################################################### data("CPS1988", package = "AER") dim(CPS1988) summary(CPS1988) ################################################### ### code chunk number 20: CPS-lm ################################################### cps_lm <- lm(log(wage) ~ experience + I(experience^2) + education + ethnicity, data = CPS1988) ################################################### ### code chunk number 21: CPS-lm-summary ################################################### summary(cps_lm) ################################################### ### code chunk number 22: CPS-noeth ################################################### cps_noeth <- lm(log(wage) ~ experience + I(experience^2) + education, data = CPS1988) anova(cps_noeth, cps_lm) ################################################### ### code chunk number 23: CPS-anova ################################################### anova(cps_lm) ################################################### ### code chunk number 24: CPS-noeth2 (eval = FALSE) ################################################### ## cps_noeth <- update(cps_lm, formula = . ~ . - ethnicity) ################################################### ### code chunk number 25: CPS-waldtest1 (eval = FALSE) ################################################### ## waldtest(cps_lm, cps_noeth) ################################################### ### code chunk number 26: CPS-waldtest2 ################################################### waldtest(cps_lm, . ~ . - ethnicity) ################################################### ### code chunk number 27: CPS-spline ################################################### library("splines") cps_plm <- lm(log(wage) ~ bs(experience, df = 5) + education + ethnicity, data = CPS1988) summary(cps_plm) ################################################### ### code chunk number 28: CPS-spline-summary ################################################### summary(cps_plm) ################################################### ### code chunk number 29: CPS-spline-BIC ################################################### cps_bs <- lapply(3:10, function(i) lm(log(wage) ~ bs(experience, df = i) + education + ethnicity, data = CPS1988)) structure(sapply(cps_bs, AIC, k = log(nrow(CPS1988))), .Names = 3:10) ################################################### ### code chunk number 30: CPS1988-plm-plot (eval = FALSE) ################################################### ## cps <- data.frame(experience = -2:60, education = ## with(CPS1988, mean(education[ethnicity == "cauc"])), ## ethnicity = "cauc") ## cps$yhat1 <- predict(cps_lm, newdata = cps) ## cps$yhat2 <- predict(cps_plm, newdata = cps) ## ## plot(log(wage) ~ jitter(experience, factor = 3), pch = 19, ## cex = 1.5, col = rgb(0.5, 0.5, 0.5, alpha = 0.02), ## data = CPS1988) ## lines(yhat1 ~ experience, data = cps, lty = 2) ## lines(yhat2 ~ experience, data = cps) ## legend("topleft", c("quadratic", "spline"), ## lty = c(2, 1), bty = "n") ################################################### ### code chunk number 31: CPS1988-plm-plot1 ################################################### png(file = "Ch-LinearRegression-CPS1988-plm-plot.png", height = 500, width = 650) par(mar = c(5, 5, 2, 4)) cps <- data.frame(experience = -2:60, education = with(CPS1988, mean(education[ethnicity == "cauc"])), ethnicity = "cauc") cps$yhat1 <- predict(cps_lm, newdata = cps) cps$yhat2 <- predict(cps_plm, newdata = cps) plot(log(wage) ~ jitter(experience, factor = 3), pch = 19, cex = 1.5, col = rgb(0.5, 0.5, 0.5, alpha = 0.02), data = CPS1988) lines(yhat1 ~ experience, data = cps, lty = 2) lines(yhat2 ~ experience, data = cps) legend("topleft", c("quadratic", "spline"), lty = c(2, 1), bty = "n") dev.off() ################################################### ### code chunk number 32: CPS-interaction ################################################### cps_int <- lm(log(wage) ~ experience + I(experience^2) + education * ethnicity, data = CPS1988) coeftest(cps_int) ################################################### ### code chunk number 33: CPS-interaction2 (eval = FALSE) ################################################### ## cps_int <- lm(log(wage) ~ experience + I(experience^2) + ## education + ethnicity + education:ethnicity, ## data = CPS1988) ################################################### ### code chunk number 34: CPS-separate ################################################### cps_sep <- lm(log(wage) ~ ethnicity / (experience + I(experience^2) + education) - 1, data = CPS1988) ################################################### ### code chunk number 35: CPS-separate-coef ################################################### cps_sep_cf <- matrix(coef(cps_sep), nrow = 2) rownames(cps_sep_cf) <- levels(CPS1988$ethnicity) colnames(cps_sep_cf) <- names(coef(cps_lm))[1:4] cps_sep_cf anova(cps_sep, cps_lm) ################################################### ### code chunk number 36: CPS-relevel ################################################### CPS1988$region <- relevel(CPS1988$region, ref = "south") cps_region <- lm(log(wage) ~ ethnicity + education + experience + I(experience^2) + region, data = CPS1988) coef(cps_region) ################################################### ### code chunk number 37: journals-WLS ################################################### jour_wls1 <- lm(log(subs) ~ log(citeprice), data = journals, weights = 1/citeprice^2) jour_wls2 <- lm(log(subs) ~ log(citeprice), data = journals, weights = 1/citeprice) ################################################### ### code chunk number 38: journals-WLS-plot ################################################### plot(log(subs) ~ log(citeprice), data = journals) abline(jour_lm) abline(jour_wls1, lwd = 2, lty = 2) abline(jour_wls2, lwd = 2, lty = 3) legend("bottomleft", c("OLS", "WLS1", "WLS2"), lty = 1:3, lwd = 2, bty = "n") ################################################### ### code chunk number 39: journals-FGLS ################################################### auxreg <- lm(log(residuals(jour_lm)^2) ~ log(citeprice), data = journals) jour_fgls1 <- lm(log(subs) ~ log(citeprice), weights = 1/exp(fitted(auxreg)), data = journals) ################################################### ### code chunk number 40: journals-FGLS-iterated ################################################### gamma2i <- coef(auxreg)[2] gamma2 <- 0 while(abs((gamma2i - gamma2)/gamma2) > 1e-7) { gamma2 <- gamma2i fglsi <- lm(log(subs) ~ log(citeprice), data = journals, weights = 1/citeprice^gamma2) gamma2i <- coef(lm(log(residuals(fglsi)^2) ~ log(citeprice), data = journals))[2] } gamma2 jour_fgls2 <- lm(log(subs) ~ log(citeprice), data = journals, weights = 1/citeprice^gamma2) ################################################### ### code chunk number 41: journals-FGLS-plot ################################################### plot(log(subs) ~ log(citeprice), data = journals) abline(jour_lm) abline(jour_fgls2, lty = 2, lwd = 2) legend("bottomleft", c("OLS", "FGLS"), lty = 1:2, lwd = 2, bty = "n") ################################################### ### code chunk number 42: usmacro-plot (eval = FALSE) ################################################### ## data("USMacroG", package = "AER") ## plot(USMacroG[, c("dpi", "consumption")], lty = c(3, 1), ## lwd = 2, plot.type = "single", ylab = "") ## legend("topleft", legend = c("income", "consumption"), ## lwd = 2, lty = c(3, 1), bty = "n") ################################################### ### code chunk number 43: usmacro-plot1 ################################################### data("USMacroG", package = "AER") plot(USMacroG[, c("dpi", "consumption")], lty = c(3, 1), lwd = 2, plot.type = "single", ylab = "") legend("topleft", legend = c("income", "consumption"), lwd = 2, lty = c(3, 1), bty = "n") ################################################### ### code chunk number 44: usmacro-fit ################################################### library("dynlm") cons_lm1 <- dynlm(consumption ~ dpi + L(dpi), data = USMacroG) cons_lm2 <- dynlm(consumption ~ dpi + L(consumption), data = USMacroG) ################################################### ### code chunk number 45: usmacro-dl-summary ################################################### summary(cons_lm1) ################################################### ### code chunk number 46: usmacro-adl-summary ################################################### summary(cons_lm2) ################################################### ### code chunk number 47: usmacro-deviance ################################################### deviance(cons_lm1) deviance(cons_lm2) ################################################### ### code chunk number 48: usmacro-dynlm-plot (eval = FALSE) ################################################### ## plot(merge(as.zoo(USMacroG[,"consumption"]), fitted(cons_lm1), ## fitted(cons_lm2), 0, residuals(cons_lm1), ## residuals(cons_lm2)), screens = rep(1:2, c(3, 3)), ## col = rep(c(1, 2, 4), 2), xlab = "Time", ## ylab = c("Fitted values", "Residuals"), main = "") ## legend(0.05, 0.95, c("observed", "cons_lm1", "cons_lm2"), ## col = c(1, 2, 4), lty = 1, bty = "n") ################################################### ### code chunk number 49: usmacro-dynlm-plot1 ################################################### plot(merge(as.zoo(USMacroG[,"consumption"]), fitted(cons_lm1), fitted(cons_lm2), 0, residuals(cons_lm1), residuals(cons_lm2)), screens = rep(1:2, c(3, 3)), col = rep(c(1, 2, 4), 2), xlab = "Time", ylab = c("Fitted values", "Residuals"), main = "") legend(0.05, 0.95, c("observed", "cons_lm1", "cons_lm2"), col = c(1, 2, 4), lty = 1, bty = "n") ################################################### ### code chunk number 50: usmacro-encomp ################################################### cons_lmE <- dynlm(consumption ~ dpi + L(dpi) + L(consumption), data = USMacroG) ################################################### ### code chunk number 51: usmacro-encomp-anova ################################################### anova(cons_lm1, cons_lmE, cons_lm2) ################################################### ### code chunk number 52: usmacro-encomptest ################################################### encomptest(cons_lm1, cons_lm2) ################################################### ### code chunk number 53: Grunfeld-data ################################################### data("Grunfeld", package = "AER") library("plm") gr <- subset(Grunfeld, firm %in% c("General Electric", "General Motors", "IBM")) pgr <- plm.data(gr, index = c("firm", "year")) ################################################### ### code chunk number 54: Grunfeld-pool ################################################### gr_pool <- plm(invest ~ value + capital, data = pgr, model = "pooling") ################################################### ### code chunk number 55: Grunfeld-FE ################################################### gr_fe <- plm(invest ~ value + capital, data = pgr, model = "within") ################################################### ### code chunk number 56: Grunfeld-FE-summary (eval = FALSE) ################################################### ## summary(gr_fe) ################################################### ### code chunk number 57: Grunfeld-FE-summary-1 ################################################### out <- capture.output(summary(gr_fe)) writeLines(out[1:19]) ################################################### ### code chunk number 58: Grunfeld-FE-summary-1 ################################################### writeLines(out[-(1:19)]) ################################################### ### code chunk number 59: Grunfeld-pFtest ################################################### pFtest(gr_fe, gr_pool) ################################################### ### code chunk number 60: Grunfeld-RE ################################################### gr_re <- plm(invest ~ value + capital, data = pgr, model = "random", random.method = "walhus") ################################################### ### code chunk number 61: Grunfeld-RE-summary ################################################### summary(gr_re) ################################################### ### code chunk number 62: Grunfeld-RE-summary-pt1 ################################################### gr_re_out <- capture.output(summary(gr_re)) gr_re_out[1] <- gsub("(Wall", "\n(Wall", gr_re_out[1], fixed = TRUE) writeLines(gr_re_out[1:18]) ################################################### ### code chunk number 63: Grunfeld-RE-summary-pt2 ################################################### writeLines(gr_re_out[-(1:18)]) ################################################### ### code chunk number 64: Grunfeld-plmtest ################################################### plmtest(gr_pool) ################################################### ### code chunk number 65: Grunfeld-phtest ################################################### phtest(gr_re, gr_fe) ################################################### ### code chunk number 66: EmplUK-data ################################################### data("EmplUK", package = "plm") form <- log(emp) ~ log(wage) + log(capital) + log(output) ################################################### ### code chunk number 67: EmplUK-pgmm ################################################### empl_ab <- pgmm(dynformula(form, list(2, 1, 0, 1)), data = EmplUK, index = c("firm", "year"), effect = "twoways", model = "twosteps", gmm.inst = ~ log(emp), lag.gmm = list(c(2, 99))) ################################################### ### code chunk number 68: EmplUK-pgmm-summary ################################################### summary(empl_ab) ################################################### ### code chunk number 69: EmplUK-pgmm-summary-pt1 ################################################### empl_ab_out <- capture.output(summary(empl_ab)) writeLines(empl_ab_out[1:17]) ################################################### ### code chunk number 70: EmplUK-pgmm-summary-pt2 ################################################### writeLines(empl_ab_out[-(1:17)]) ################################################### ### code chunk number 71: Grunfeld-data2 ################################################### gr2 <- subset(Grunfeld, firm %in% c("Chrysler", "IBM")) pgr2 <- plm.data(gr2, c("firm", "year")) ################################################### ### code chunk number 72: Grunfeld-SUR ################################################### library("systemfit") gr_sur <- systemfit(invest ~ value + capital, method = "SUR", data = pgr2) ################################################### ### code chunk number 73: Grunfeld-SUR-summary ################################################### summary(gr_sur, residCov = FALSE, equations = FALSE)