# Duncan occupation data example # Chris Adolph # Clear memory rm(list=ls()) # Load the occupation data library(car) data(Duncan) attach(Duncan) # A regression analysis, with inference model1 <- prestige~income+education # model specification lm.res1 <- lm(model1, data=Duncan) # run the regression beta <- lm.res1\$coefficients # retrive the of betas vc <- vcov(lm.res1) # retrieve the var-cov matrix se <- sqrt(diag(vc)) # calc a vector of ses tstat <- beta/se # calc vector of tstats pval <- 2*(1-pt(tstat,42)) # calc p-values # Confidence intervals of betas by hand lower.95 <- beta - qt(0.025,42)*se upper.95 <- beta + qt(0.025,42)*se # Confidence intervals of betas the easy way library(stats) confint(lm.res1, level=0.95) ## Make a plot... ## CIs for yhat given a set of hypothetical income & education xhyp <- seq(min(income), max(income),1) zhyp <- rep(mean(education), length(xhyp)) hypo <- data.frame(income=xhyp, education=zhyp) pred <- predict(lm.res1, newdata=hypo, interval="confidence", level=0.95) yhat <- pred[,1] yhat.lower <- pred[,2] yhat.upper <- pred[,3] # Start the plot # Uncomment the below for pdf output (step 1) #pdf("yhatexample.pdf",width=5,height=4.5) plot(y=prestige, x=income, type="n") # Make the x-coord of a confidence envelope polygon xpoly <- c(xhyp, rev(xhyp), xhyp[1]) # Make the y-coord of a confidence envelope polygon ypoly <- c(yhat.lower, rev(yhat.upper), yhat.lower[1]) # Choose the color of the polygon col <- "gray" # Plot the polygon first, before the points & lines polygon(x=xpoly, y=ypoly, col=col, border=FALSE ) # Plot the fitted line lines(x=xhyp, y=yhat) # Uncomment the below for pdf output (step 2) #dev.off()