R code for the beauty and teaching evaluations example

# Load in the "arm" package. (Always do this at the start of your R session.)
library (arm)

# The data came as an Excel file (ProfEvaltnsBeautyPublic.xls)
# I went into Excel and saved it as a .csv file (comma-separated version)

# Read the data into R, including the variable names (headers)
beauty.data <- read.table ("http://www.stat.columbia.edu/~gelman/arm/examples/beauty/ProfEvaltnsBeautyPublic.csv", header=TRUE, sep=",")
#beauty.data <- read.table ("ProfEvaltnsBeautyPublic.csv", header=TRUE,sep=",")

# Attach the data so they are accessible using their variable names
attach.all (beauty.data)

###############################################################################
# Do more beautiful profs get higher evaluations?
###############################################################################

# Rename the two variables for convenience
beauty <- btystdave
eval <- courseevaluation

# Make a scatterplot
plot (beauty, eval)

# Make a scatterplot, labeling the axes
plot (beauty, eval, xlab="beauty", ylab="average teaching evaluation")

# Fit a linear regression
lm.1 <- lm (eval ~ beauty)
display (lm.1)

# Extract the information: fitted line is y = a + b*x + error
# Errors have mean 0, standard deviation sigma
a <- coef (lm.1)["(Intercept)"]
b <- coef (lm.1)["beauty"]
sigma <- sigma.hat (lm.1)

# Plot the regression line on top of the scatterplot
curve (a + b*x, add=TRUE)

# Add dotted lines to show +/- 1 standard deviation
curve (a + b*x + sigma, lty=2, add=TRUE)
curve (a + b*x - sigma, lty=2, add=TRUE)

###############################################################################
# Do things differ for male and female profs? Parallel regression lines
###############################################################################

lm.2 <- lm (eval ~ beauty + female)
display (lm.2)

# Set up a 2x2 grid of plots
par (mfrow=c(2,2))

# Make separate plots for men and women
b0 <- coef (lm.2)["(Intercept)"]
b1 <- coef (lm.2)["beauty"]
b2 <- coef (lm.2)["female"]

plot (beauty[female==0], eval[female==0], xlim=range(beauty), ylim=range(eval),
xlab="beauty", ylab="average teaching evaluation", main="Men")
curve (b0 + b1*x + b2*0, add=TRUE)

plot (beauty[female==1], eval[female==1], xlim=range(beauty), ylim=range(eval),
xlab="beauty", ylab="average teaching evaluation", main="Women")
curve (b0 + b1*x + b2*1, add=TRUE)

# Display both sexes on the same plot
# First make the plot with type="n" (which displays axes but does not plot
# the points), then plot the points and lines separately for each sex
plot (beauty, eval, xlab="beauty", ylab="average teaching evaluation",
main="Both sexes", type="n")
points (beauty[female==0], eval[female==0], col="blue")
curve (b0 + b1*x + b2*0, add=TRUE, col="blue")
points (beauty[female==1], eval[female==1], col="red")
curve (b0 + b1*x + b2*1, add=TRUE, col="red")

###############################################################################
# Do things differ for male and female profs? Non-parallel regression lines
###############################################################################

lm.3 <- lm (eval ~ beauty + female + beauty*female)
display (lm.3)

# Set up a new 2x2 grid of plots
par (mfrow=c(2,2))

# Display the parallel regression lines in gray and the non-parallel lines
# in heavy black
b0.new <- coef (lm.3)["(Intercept)"]
b1.new <- coef (lm.3)["beauty"]
b2.new <- coef (lm.3)["female"]
b3.new <- coef (lm.3)["beauty:female"]

plot (beauty[female==0], eval[female==0], xlim=range(beauty), ylim=range(eval),
xlab="beauty", ylab="average teaching evaluation", main="Men")
curve (b0 + b1*x + b2*0, lwd=.5, col="gray", add=TRUE)
curve (b0.new + b1.new*x + b2.new*0 + b3.new*x*0, lwd=2, col="black", add=TRUE)

plot (beauty[female==1], eval[female==1], xlim=range(beauty), ylim=range(eval),
xlab="beauty", ylab="average teaching evaluation", main="Women")
curve (b0 + b1*x + b2*1, lwd=.5, col="gray", add=TRUE)
curve (b0.new + b1.new*x + b2.new*1 + b3.new*x*1, lwd=2, col="black", add=TRUE)

###############################################################################
# More models
###############################################################################

lm.4 <- lm (eval ~ beauty + female + beauty*female + age)
display (lm.4)

lm.5 <- lm (eval ~ beauty + female + beauty*female + minority)
display (lm.5)

lm.6 <- lm (eval ~ beauty + female + beauty*female + nonenglish)
display (lm.6)

lm.7 <- lm (eval ~ beauty + female + beauty*female + nonenglish + lower)
display (lm.7)