# Basic regression commands, and two ways to make # plots of estimated regression models using # predict() and the base graphics package # # Code and example from Iversen & Soskice 2003 # # Chris Adolph # April 12, 2011 # Clear memory of all objects rm(list=ls()) # Load data file <- "iver.csv" data <- read.csv(file,header=TRUE) attach(data) # A bivariate model lm.result <- lm(povred~lnenp) print(summary(lm.result)) # A new model with multiple regressors lm.result2 <- lm(povred~lnenp+maj+pr) print(summary(lm.result2)) # A new model with multiple regressors and no constant lm.result3 <- lm(povred~lnenp+maj+pr+unam-1) print(summary(lm.result3)) # A new model with multiple regressors and an interaction lm.result4 <- lm(povred~lnenp+maj+pr+lnenp:maj) print(summary(lm.result4)) # A different way to specify an interaction lm.result5 <- lm(povred~pr+lnenp*maj) print(summary(lm.result5)) # Make a plot of the data (automatic axes, etc) plot(x=lnenp, y=povred, xlab="Log Effective Number of Parties", ylab="Poverty Reduction") # One way to add a regression line to the plot abline(lm.result\$coefficients[1], # Intercept lm.result\$coefficients[2], # Slope col="black") # The above is easy for bivariate models # For multivariate models, you need to calculate # an appropriate intercept to take account # of all the other covariates # Generate expected values & CIs for povred at each lnenp lnenp.hyp <- seq(min(lnenp),max(lnenp),0.1) xnew <- list(lnenp=lnenp.hyp ) povred.pred <- predict(lm.result, newdata=xnew, interval="confidence", level=0.95 ) # Open a pdf file for plotting # Uncomment this for pdf output instead of screen #pdf("redist.pdf", # height=5, # width=5) # Create a new plot plot.new() # Set the plotting region limits par(usr=c(0.5,2,0,100)) # Create the x-axis x.ticks <- c(2,3,4,5,6,7) axis(1, # Which axis to make (1 indicates x) at=log(x.ticks), # Where to put the ticks labels=x.ticks # How to label the ticks ) # Create the y-axis axis(2,at=seq(0,100,10)) # Add plot titles title(xlab="Effective Number of Parties", ylab="Poverty Reduction" ) # Plot ci for the regression line # Make the x-coord of a confidence envelope polygon xpoly <- c(lnenp.hyp, rev(lnenp.hyp), lnenp.hyp[1]) # Make the y-coord of a confidence envelope polygon ypoly <- c(povred.pred[,2], rev(povred.pred[,3]), povred.pred[1,2]) # Choose the color of the polygon col <- "gray70" # Plot the polygon first, before the points & lines polygon(x=xpoly, y=ypoly, col=col, border=FALSE ) # Plot the expected values for the regression model lines(x=lnenp.hyp, y=povred.pred[,1], col="black") # Plot the data for the regression model # Uncomment this to do all points in same color and symbol #points(x=lnenp, # y=povred, # col="black", # see colors() for color names # pch=1) # see example(points) for symbols points(x=lnenp[maj==1], y=povred[maj==1], col="blue", # see colors() for color names pch=17) # see example(points) for symbols points(x=lnenp[pr==1], y=povred[pr==1], col="green", # see colors() for color names pch=15) # see example(points) for symbols points(x=lnenp[unam==1], y=povred[unam==1], col="red", # see colors() for color names pch=16) # see example(points) for symbols text(x=lnenp, y=povred-3, labels=cty, col="black", cex=0.5 ) # Finish drawing the box around the plot area box() # Close the device (ie, save the graph) # Uncomment this for pdf output #dev.off()