# Basic regression commands and # plots of estimated regression models using # predict() and the base graphics package # # Code and example from Iversen & Soskice 2003 # # Chris Adolph # 25 June 2015 # Clear memory of all objects rm(list=ls()) # Load libraries library(RColorBrewer) # For nice colors # Load data file <- "iverRevised.csv" iversen <- read.csv(file,header=TRUE) # Create dummy variables for each party system iversen$majoritarian <- as.numeric(iversen$partySystem=="Majoritarian") iversen$proportional <- as.numeric(iversen$partySystem=="Proportional") iversen$unanimity <- as.numeric(iversen$partySystem=="Unanimity") # A bivariate model, using a formula to log transform a variable model1 <- povertyReduction ~ log(effectiveParties) lm.res1 <- lm(model1, data=iversen) summary(lm.res1) # A new model with multiple regressors model2 <- povertyReduction ~ log(effectiveParties) + majoritarian + proportional lm.res2 <- lm(model2, data=iversen) summary(lm.res2) # A new model with multiple regressors and no constant model3 <- povertyReduction ~ log(effectiveParties) + majoritarian + proportional + unanimity - 1 lm.res3 <- lm(model3, data=iversen) summary(lm.res3) # A new model with multiple regressors and an interaction model4 <- povertyReduction ~ log(effectiveParties) + majoritarian + proportional + log(effectiveParties):majoritarian lm.res4 <- lm(model4, data=iversen) summary(lm.res4) # A more efficient way to specify an interaction model5 <- povertyReduction ~ log(effectiveParties)*majoritarian + proportional lm.res5 <- lm(model5, data=iversen) summary(lm.res5) # Make a plot of the data (automatic axes, etc) plot(x=iversen$effectiveParties, y=iversen$povertyReduction, log="x", xlab="Effective Number of Parties", ylab="Poverty Reduction") # One way to add a regression line to the plot abline(lm.res1$coefficients[1], # Intercept lm.res1$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 from model 1 for # povertyReduction at each level of effectiveParties # Counterfactual levels of effectiveParties effectiveParties.hyp <- seq(from=min(iversen$effectiveParties), to=max(iversen$effectiveParties), by=0.1) # The list below must contain a vector of hypothetical values # for every covariate in the model formula, with names to match xnew <- list(effectiveParties=effectiveParties.hyp ) # Expected level of poverty reduction yhyp <- predict(lm.res1, 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=7, # width=7) # Create a new plot plot.new() # Get nice colors col <- brewer.pal(3, "Set1") # Set the plotting region limits par(usr=c(min(iversen$effectiveParties) - 0.25, # x minimum of graphic max(iversen$effectiveParties) + 0.25, # x maximum of graphic 0, # y minimum 100) # y maximum ) # Create the x-axis x.ticks <- c(2, 3, 4, 5, 6, 7) axis(1, # Which axis to make (1 indicates x) at=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(effectiveParties.hyp, rev(effectiveParties.hyp), effectiveParties.hyp[1]) # Make the y-coord of a confidence envelope polygon ypoly <- c(yhyp[,2], rev(yhyp[,3]), yhyp[1,2]) # Plot the polygon first, before the points & lines polygon(x=xpoly, y=ypoly, col="gray80", border=FALSE ) # Plot the expected values for the regression model lines(x=effectiveParties.hyp, y=yhyp[,1], col="black") # Plot the data for the regression model # Uncomment this to do all points in same color and symbol #points(x=iversen$effectiveParties[iversen$majoritarian==1], # y=iversen$povertyReduction[iversen$majoritarian==1], # col="black", # pch=1) # #text(x=iversen$effectiveParties, # y=iversen$povertyReduction - 2.5, # labels=iversen$country, # col="black", # cex=0.7 # ) # Plot majoritarian cases with labels points(x=iversen$effectiveParties[iversen$majoritarian==1], y=iversen$povertyReduction[iversen$majoritarian==1], col=col[1], # from RColorBrewer above; # see example(brewer.pal) for ideas pch=17) # see example(points) for symbols text(x=iversen$effectiveParties[iversen$majoritarian==1], y=iversen$povertyReduction[iversen$majoritarian==1] - 2.5, labels=iversen$country[iversen$majoritarian==1], col=col[1], cex=0.7 ) # Plot proportional cases with labels points(x=iversen$effectiveParties[iversen$proportional==1], y=iversen$povertyReduction[iversen$proportional==1], col=col[2], pch=15) text(x=iversen$effectiveParties[iversen$proportional==1], y=iversen$povertyReduction[iversen$proportional==1] - 2.5, labels=iversen$country[iversen$proportional==1], col=col[2], cex=0.7 ) # Plot unanimous cases with labels points(x=iversen$effectiveParties[iversen$unanimity==1], y=iversen$povertyReduction[iversen$unanimity==1], col=col[3], pch=16) text(x=iversen$effectiveParties[iversen$unanimity==1], y=iversen$povertyReduction[iversen$unanimity==1] - 2.5, labels=iversen$country[iversen$unanimity==1], col=col[3], cex=0.7 ) # Finish drawing the box around the plot area box() # Close the device (that is, save the graph) # Uncomment this for pdf output #dev.off() #################################################################### # For an alternative approach using the tile package, # see: http://faculty.washington.edu/cadolph/vis/inequalityScatter.R