---
title: Hamermesh and Parker (2005)
author:
  - name: Andrew Pua
date: October 26, 2022
format:
  html:
    code-fold: true
    self-contained: true
---

## Preamble

Here I load the `foreign` package and the look at the names of the variables. I also calculate some summary statistics, but leave out the interpretation.

```{r}
rm(list=ls())
library(foreign)
ratings <- read.dta("TeachingRatings.dta")
names(ratings)
apply(ratings, 2, mean)
apply(ratings[, c("age", "beauty", "course_eval")], 2, sd)
```

### Remarks about summary statistics

In these solutions, I will focus on the interpretation of the regression coefficients and report the summary statistics. Note that standard deviations of dummy variables are harder to interpret. But the means of dummy variables are really proportions. It is also possible to show some histograms. 

### General remarks about the solutions

You were asked to provide the best interpretations possible. This does not necessarily mean that you have to talk about statistical significance. For those who have not encountered this phrase, skip this remark. For those who have encountered it, statistical significance requires assumptions which may not be satisfied by this application. Furthermore, interpreting coefficients and deciding whether the coefficients are significantly different from zero (which is not the most correct terminology either) are two different tasks!

## Regression 1

I recentered age at its mean to make the regression coefficients more interpretable. 

```{r}
ratings$c.age <- ratings$age - mean(ratings$age)
lm(course_eval ~ beauty + c.age, data = ratings)
```

1.    Intercept: The average course evaluation for instructors roughly `r round(mean(ratings$age),0)` years old having an average beauty rating is `r round(coef(lm(course_eval ~ beauty + c.age, data = ratings))[[1]],2)` points. 

2.    Coefficient of `beauty`: When we compare instructors with the same age but whose beauty ratings differ by a standard deviation, the instructors who have a higher beauty rating have average course evaluations that are higher by `r round(sd(ratings$beauty)*coef(lm(course_eval ~ beauty + c.age, data = ratings))[[2]],2)` points.

3.    Coefficient of `c.age`: When we compare instructors with the same beauty rating but whose ages differ by 10 years, the average course evaluation of older instructors are higher by `r round(10*coef(lm(course_eval ~ beauty + c.age, data = ratings))[[3]],3)` points.



### Remarks

1.    `beauty` is centered at 0 from the description. So `beauty = 0` means that the beauty rating is at the average. 
2.    The instructors you are comparing need not have to be `c.age = 0`. In fact, the resulting comparison will still be the same even if you fix  `c.age = 10` for example. Try it for yourself. One standard deviation for beauty ratings is about `r round(sd(ratings$beauty),2)` and you multiply this by the coefficient of `beauty` to get the result. 
3.    The instructors you are comparing need not have to be `beauty = 0`. I chose 10 years here because the coefficient of `age` is very small in magnitude and 10 years is roughly a standard deviation in the age distribution. 

## Regression 2

```{r}
lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings)
```

1.    Intercept: The average course evaluation for instructors roughly `r round(mean(ratings$age),0)` years old having an average beauty rating is `r round(coef(lm(course_eval ~ beauty + c.age, data = ratings))[[1]],2)` points. 
2.    Coefficient of `beauty`: When we compare instructors roughly `r round(mean(ratings$age),0)` years old but whose beauty ratings differ by a standard deviation, the instructors who have a higher beauty rating have average course evaluations that are higher by `r round(sd(ratings$beauty)*coef(lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings))[[2]],2)` points.
3.    Coefficient of `c.age`: When we compare instructors with average beauty rating but whose ages differ by 10 years, the average course evaluation of older instructors are higher by `r round(10*coef(lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings))[[3]],2)` points.
4.    Coefficient of `beauty:c.age`: When we compare instructors roughly `r round(mean(ratings$age),0)+10` years old but whose beauty ratings differ by a standard deviation, the instructors who have a higher beauty rating have average course evaluations that are higher by `r round(sd(ratings$beauty)*(coef(lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings))[[2]]+coef(lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings))[[4]]*10),2)` points.

### Remarks on the interaction term

1.    Based on the solutions to the related Technical Exercises, it is extremely difficult to interpret the coefficient of the interaction term alone. 
2.    I chose 58 years old because this is one standard deviation above the average age.
3.    To obtain the difference in average course evaluations, apply the result in the Technical Exercise if you wish. It is $(0.15+0.01\times(58-48))\times 0.79 \approx 0.2$. The value `r round(sd(ratings$beauty),2)` is the standard deviation for beauty ratings.
4.    To express in words the interpretation of the coefficient of the interaction term alone is extremely difficult because you will be comparing two comparisons (see the related Technical Exercise for this exercise set). 
5.    Some students used a causal interpretation to get around this difficulty. They state something like "The effect on course evaluation of a one point increase in beauty increases with age". This is ok if a causal interpretation could be justified, but in our data analysis so far, we are only interpreting what the data is telling us without any additional assumptions. 
6.    Another way to interpret the coefficient of `beauty:c.age` is as follows: When we compare instructors with beauty rating is a standard deviation above the average but whose ages differ by 10 years, older instructors have average course evaluations that are higher by `r round(10*(coef(lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings))[[3]]+coef(lm(course_eval ~ beauty + c.age + beauty:c.age, data = ratings))[[4]]*sd(ratings$beauty)),2)` points.