---
title: "Demonstrating some of the math behind `lm()`"
subtitle: Solutions
format: html
self-contained: true
editor: visual
---

## Preamble

```{r}
rm(list=ls())
library(foreign)
ratings <- read.dta("TeachingRatings.dta")
```

## Item 1

```{r}
lm(course_eval ~ beauty - 1, data = ratings)
mean(residuals(lm(course_eval ~ beauty - 1, data = ratings)))
```

The residuals do not have mean zero. 

## Item 2

```{r}
adj.course_eval <- residuals(lm(course_eval ~ beauty, data = ratings))
adj.female <- residuals(lm(female ~ beauty, data = ratings))
lm(adj.course_eval ~ adj.female - 1, data = ratings)
lm(course_eval ~ beauty + female, data = ratings)
```

Thus, the slope of `adj.female` matches the slope obtained from the slope of `female` in the linear regression of course evaluations on beauty and female. 

When we compare courses which have instructors with the same level of beauty rating, the average course evaluations of female instructors is `r -round(coef(lm(course_eval ~ beauty + female, data = ratings))[[3]], 2)` points lower than the average course evaluations of male instructors

## Item 3

```{r}
results <- lm(course_eval ~ beauty + female, data = ratings)
mean(fitted(results))
mean(ratings$course_eval)
mean(residuals(results))
cor(fitted(results), residuals(results))
cor(ratings$beauty, residuals(results))
cor(ratings$female, residuals(results))
```

### Remarks for Item 3

These do not have to be part of the solution. You may choose to assign names to the objects generated by `fitted()` and `residuals()`. You can also make these fitted values and residuals part of the original data set. You can also use the `with()` to help you reduce typing `ratings` so many times. Code related to these remarks, which I left unevaluated, is presented below. 


```{r}
#| eval: false
results <- lm(course_eval ~ beauty + female, data = ratings)
results.f <- fitted(results)
# or something like ratings$fitted <- fitted(results)
mean(results.f)
results.r <- residuals(results)
# or something like ratings$residuals <- residuals(results)
with(ratings, mean(course_eval))
mean(results.r)
cor(results.f, results.r)
# or something like with(ratings, cor(fitted, residuals))
cor(ratings$beauty, results.r)
cor(ratings$female, results.r)
```
