## ----, warning = FALSE, message = FALSE---------------------------------- slid.dat <- read.csv("http://www.jkarreth.net/files/car_SLID.csv") summary(slid.dat) ## ----, warning = FALSE, message = FALSE---------------------------------- mod <- lm(wages ~ education + age + sex, data = slid.dat) summary(mod) ## ----, warning = FALSE, message = FALSE---------------------------------- summary(resid(mod)) ## ----, warning = FALSE, message = FALSE---------------------------------- plot(density(resid(mod))) ## ----, warning = FALSE, message = FALSE---------------------------------- summary(rstudent(mod)) plot(density(rstudent(mod))) ## ----, warning = FALSE, message = FALSE---------------------------------- qqnorm(rstudent(mod)) qqline(rstudent(mod)) ## ----, warning = FALSE, message = FALSE---------------------------------- hist(slid.dat$wages) ## ----, warning = FALSE, message = FALSE---------------------------------- summary(slid.dat$wages) slid.dat$wages.ln <- log(slid.dat$wages) ## ----, warning = FALSE, message = FALSE---------------------------------- mod.ln <- lm(wages.ln ~ education + age + sex, data = slid.dat) summary(mod.ln) qqnorm(rstudent(mod.ln)) qqline(rstudent(mod.ln)) ## ----, warning = FALSE, message = FALSE---------------------------------- plot(density(rstudent(mod.ln))) ## ----, warning = FALSE, message = FALSE---------------------------------- plot(x = fitted.values(mod.ln), y = rstudent(mod)) ## ----, warning = FALSE, message = FALSE---------------------------------- library(car) residualPlots(mod.ln) ## ----, warning = FALSE, message = FALSE---------------------------------- summary(mod) ## ----, warning = FALSE, message = FALSE---------------------------------- crPlots(mod.ln) ## ----, warning = FALSE, message = FALSE---------------------------------- crPlots(model = mod.ln, terms = ~ . - sex) ## ----, warning = FALSE, message = FALSE---------------------------------- vif(mod.ln) ## ----, warning = FALSE, message = FALSE---------------------------------- set.seed(123) n.sample <- 200 x1 <- rnorm(n.sample, 2, 5) x2 <- x1 + rnorm(n.sample, 0, 1) x3 <- runif(n.sample, min = -10, max = 10) a <- 0.2 b1 <- 1.1 b2 <- 0.2 b3 <- -0.7 e <- rnorm(n.sample, 0, 10) y <- a + b1 * x1 + b2 * x2 + b3 * x3 + e mod.sim <- lm(y ~ x1 + x2 + x3) summary(mod.sim) vif(mod.sim) sqrt(vif(mod.sim)) mod.sim2 <- lm(y ~ x1 + x3) summary(mod.sim2) vif(mod.sim2) ## ----, warning = FALSE, message = FALSE---------------------------------- robSE <- function(mod){ require(lmtest, quiet = TRUE) X.mat <- model.matrix(mod) # note that X includes 1 for the intercept y <- mod$mod[1] # Response variable e <- resid(mod) # Vector of residuals n <- dim(X.mat)[1] # Number of observations k <- dim(X.mat)[2] # Number of predictors (including intercept) vcov.mat <- vcov(mod) # Variance-covariance matrix SE <- sqrt(diag(vcov.mat)) # Standard errors E.mat <- matrix(0, ncol = n, nrow = n) # n-by-n Matrix of residuals diag(E.mat) <- e^2 # Squared residuals in the diagonal sum.mat <- t(X.mat) %*% E.mat %*% X.mat vcov.mat.rob <- abs(solve(t(X.mat) %*% X.mat) %*% sum.mat %*% solve(t(X.mat) %*% X.mat)) rob.SE <- sqrt(diag(vcov.mat.rob)) ## Add the degrees-of-freedom correction dfc <- sqrt(nrow(X.mat)) / sqrt(nrow(X.mat) - ncol(X.mat)) rob.SE.dfc <- dfc * sqrt(diag(vcov.mat.rob)) return(rob.SE.dfc) } ## ----, warning = FALSE, message = FALSE---------------------------------- summary(mod.ln) robSE(mod.ln) ## ----, warning=FALSE, message=FALSE-------------------------------------- library(texreg) screenreg(list(mod.ln, mod.ln), override.se = list(NA, robSE(mod.ln)), custom.model.names = c("Original SEs", "Robust SEs"), digits = 5 )