 
---
title: "Lab 12: Inference for Linear Regression"
author: "Insert Name"
output:
  pdf_document
 
---

## Instructions

**Optional: If for credit due by Friday May 7th**

Read through the following document (it may be helpful to have both the .rmd and .pdf files open side-by-side). Interspersed throughout the text are 23 **Exercises**. Answer these questions in the .rmd file by inserting a code chunk or typing text between the stars directly below the exercise statement.


```{r, include=FALSE}
# Do not edit this code block/chunk
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning = FALSE, fig.width = 5, fig.height = 3)
```

## Load packages

For this lab we will be using the `tidyverse` for data wrangling and visualization, `moderndive` for linear modeling, `infer` for statistical inference, and a data set will come from the `oilabs` package.

```{r, include=FALSE, message = FALSE, warning = FALSE}
library(tidyverse)
library(knitr)
library(moderndive)
library(infer)
library(oilabs)
```


 

## Twins and IQ

In the mid 20th century, a study was conducted that tracked down
identical twins that were separated at birth: one child was raised in
the home of their biological parents and the other in a foster home. In
an attempt to answer the question of whether intelligence is the result
of nature or nurture, both children were given IQ tests.

The researcher that collected this data was named Cyril
Burt, and his work was the subject of some controversy, including the claim that he falsified some of his research data.  One of the suspect datasets is the identical twins data that we will look at.  This has come to be called "The Burt Affair".  Take a moment and read about it before beginning this lab.

<https://en.wikipedia.org/wiki/Cyril_Burt#%22The_Burt_Affair%22>

This data is in a `.csv` file, which is a very general format for data.
You can load our usual packages as well as the data (into the `twins`
data frame) with the following commands:

```{r}
 
twins <- read.csv("http://andrewpbray.github.io/data/twins.csv")
options(scipen = 99)

set.seed(0909)
```

To reduce scientific notation appearing in the output, we can also
specify the option `scipen` to be a large positive value. More info is
available via `?options`.

Since this data was read in from a `.csv` file, there is no help file
associated with it, but the IQ of the foster
child is recorded in `Foster` and the IQ of the biologically-raised
child in `Biological`. The socieoeconomic status of the biological
family was also recorded in `Social`.

### Exercise 1

Create a plot of the relationship between the IQs of the two
    children and describe the relationship in words (there is no clear
    independent variable here, so just put `Biological` on the x-axis).
    
**********************************************************************

**********************************************************************

### Exercise 2

We will use `infer`to generate a plot with 12 facets, each one a
    scatterplot illustrating what the data would look like if the IQ of
    the foster child is independent of the IQ of the child raised by the
    biological parents. 
    
Run the following code to generate the plots. Be sure you understand what each line of code is doing.

```{r eval = T}
twins %>%
  specify(Foster ~ Biological) %>%
  hypothesize(null = "independence" ) %>%
  generate(reps = 12, type = "permute" ) %>%
  ggplot(aes(x = Biological , y = Foster)) +
  geom_point() +
  facet_wrap(~ replicate)
```

Does your plot from question 1 look like these other plots or is it 
  systematically different? How so?


**********************************************************************

**********************************************************************

### Exercise 3

Fit a linear model to this data and have a look at the regression table
for the model. Which one of the following is *false* and why?
    
  A. For each 10 point increase in the biological twin’s IQ, we
        expect the foster twin’s IQ to increase on average by roughly 9
        points.
        
  B. The linear model (the estimated mean function) is approximately
        $\hat{y} =9.2+0.9$x.
        
  C. Foster children with higher than average IQs (average of the
        foster children) will *always* have a biological twins with higher than
        average IQs as well (average of the biological children).


**********************************************************************

**********************************************************************

### Exercise 4

Is it reasonable to draw conclusions from this model? Assess the
    conditions for a valid model based on the context of this problem
    using appropriate plots. Recall that it is only appropriate to 
    make predictions based on a regression model if:
    
  1. The data should show a linear trend.
  
  2. The distribution of residuals is approximately Normal.
  
  3. The residuals have constant variability.

**********************************************************************

**********************************************************************

### Exercise 5

This study was used to weigh in on the question of whether IQ is a
    result of nature (your genes) or nurture (your environment). If IQ
    was purely a result of nurture, what slope would you expect to see
    in your linear model?
    
Run the following code to generate a null distribution for slopes, assuming
    Foster IQ and Biological IQ are independent. Then use `infer` to visualize
    the distribution and compute the p_value. Compare the P-value you get here
    to the one given for slope in the regression table in Exercise 2.

```{r cache = T}
twins_null<-twins %>% 
  specify(Foster ~ Biological) %>% 
  hypothesize(null = "independence") %>% 
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "slope")
```

**********************************************************************

**********************************************************************

### Exercise 6

At the other end of the argument, the distribution of IQ might be
    entirely determined by genes. What slope would you expect to see
    if this were true?

Run the following code to generate a bootstrap distribution for slopes. 
    Then use `infer` and the percentile method to create a 95% confidence interval
    on your slope. Compare the confidence interval here to the one given in the
    regression table in Exercise 2. Finally, assess whether this claim is reasonable
    given your data set.
    
```{r}
twins_boot<-twins %>% 
  specify(Foster ~ Biological) %>% 
  generate(reps = 1000, type = "bootstrap") %>% 
  calculate(stat = "slope")

```


**********************************************************************

**********************************************************************

### Exercise 7

Now consider how the association between IQs would change if social
    class were taken into account. Create a plot to illustrate the
    relationship between all three variables (recall the use of the
    `color` argument in `aes()`). Just by looking at the plot and
    without fitting any models, does it appear that the relationship
    between the IQs is the same for high social class children as it is
    for low social class children?

**********************************************************************

**********************************************************************

### Exercise 8

Considering the results that you found and what you know about the controversy surrounding the data, would you report these results?  Why, or why not?

**********************************************************************

**********************************************************************



### Exercise 9

Begin by reading the following blog post:

<https://marginalrevolution.com/marginalrevolution/2020/03/what-worked-in-1918-1919.html>

Then access the original paper that seeded this post, which was published in the Journal of the American Medical Association (Don't be discuouraged if some of the statistical terminology is unfamiliar; this is research-level paper using the full toolbox of statistical analysis)

<https://biotech.law.lsu.edu/blog/joc70085_644_654.pdf>

A. What were the objectives of the study? Write this in a jargon-free way another person not familiar with the terminology or content of our course could understand.

*********************************************************************************

*********************************************************************************

B. For the main data set, featured in Table 1, what is the unit of observation?

*********************************************************************************

*********************************************************************************


C. Classify each type of variable in Table 1. Think carefully about this one: this contains one data type that we haven’t explicitly discussed in class.

*********************************************************************************

*********************************************************************************

D. Read the descriptions of the various response and explanatory variables that are featured in figure 1. Which set seem to you to best capture what we’re interested in here: a sense of the safety measures taken by the city and the effect it had on public health.

*********************************************************************************

*********************************************************************************

E. Of the four relationships shown in figure 1, which one is strongest? How can you tell?

*********************************************************************************

*********************************************************************************




### Exercise 10

Use the article refernced in Problem 1 to answer these questions.

A. While we can make predictions off of any regression line, in order to perform inference accurately, we need our mathematical assumptions to be reasonable. Do you have reservations about the accuracy of inference in any of these four settings?

*********************************************************************************

*********************************************************************************

B. Write out the null hypothesis corresponding the p-values listed in these plots.

*********************************************************************************

*********************************************************************************

C. Two cities, highlighted in blue, are what the authors flag as “outliers”. If they were removed, which of the correlations and p-values do you expect would change the most? Which would be most unaffected?


*********************************************************************************

*********************************************************************************

D. Summarize in 1-2 paragrpahs why the authors were justified in deducing a **causal** link between nonpharmaceutical intervention and excessive death rate, despite the fact that the investigation was an observational study and not an experiment.

*********************************************************************************

*********************************************************************************

## Grading the professor

Many college courses conclude by giving students the opportunity to
evaluate the course and the instructor anonymously. However, the use of
these student evaluations as an indicator of course quality and teaching
effectiveness is often criticized because these measures may reflect the
influence of non-teaching related characteristics, such as the physical
appearance of the instructor. The article titled, “Beauty in the
classroom: instructors’ pulchritude and putative pedagogical
productivity” by Hamermesh and Parker found that instructors who are
viewed to be better looking receive higher instructional ratings.

In this lab we will analyze the data from this study in order to learn
what goes into a positive professor evaluation.

### The data

The data were gathered from end of semester student evaluations for a
large sample of professors from the University of Texas at Austin. In
addition, six students rated the professors’ physical appearance. The
result is a data frame where each row contains a different course and
columns represent variables about the courses and professors.

Let’s load the data:

```{r}
data(evals)
```

We have observations on 21 different variables, some categorical and
some numerical. The meaning of each variable can be found by bringing up
the help file (select the help file in the oilabs package):

```{r eval = F}
?evals
```

## Exploring the data

### Exercise 11

Is this an observational study or an experiment? The original
research question posed in the paper is whether beauty causes
a difference in course evaluations. Given the investigation design,
is it possible to answer this question as it is phrased? If not,
rephrase the question.
    
************************************************************************

************************************************************************

### Exercise 12

Describe the shape, center and spread of the distribution of `score`.
What does that tell you about how students rate courses? Is this
what you expected to see? Explain.
    
    
************************************************************************

************************************************************************


## Multiple linear regression

In order to see if beauty is a significant predictor of professor
score after we’ve accounted for the gender of the professor, we can build a model that includes both beauty and gender terms.

Run the following code to fit a multiple regression model to the data predicting `score` as a function of `bty_avg` and `gender`. Investigate the results in the regression table and regression summary.

```{r}
m_bty_gen <- lm(score ~ bty_avg + gender, data = evals)

get_regression_table(m_bty_gen)
get_regression_summaries(m_bty_gen)
```




### Exercise 13

Is `bty_avg` a significant predictor of `score` in a model that controls for `gender`? 

************************************************************************

************************************************************************

Note that the estimate for `gender` is now called `gendermale`. You’ll
see this name change whenever you introduce a categorical variable. The
reason is that R recodes `gender` from having the values of `female` and
`male` to being an indicator variable called `gendermale` that takes a
value of 0 for females and a value of 1 for males. (Such
variables are often referred to as “dummy” variables.)

As a result, for females, the parameter estimate is multiplied by zero,
leaving the intercept and slope form familiar from simple regression.

$$\begin{aligned} \widehat{score} &= \hat{\beta}_0 + \hat{\beta}_1 \times bty\_avg + \hat{\beta}_2 \times (0) \\ &= \hat{\beta}_0 + \hat{\beta}_1 \times bty\_avg\end{aligned}$$

We can plot this line and the line corresponding to males with the
following tweak by using `geom_parallel_slopes()` instead of `geom_smooth()`.

```{r}
ggplot(evals, aes(x = bty_avg, y = score, color = gender)) +
  geom_point() +
  geom_parallel_slopes(se = FALSE)
```

### Exercise 14

What is the equation of the line corresponding to males? For two
professors who received the same beauty rating, which gender tends
to have the higher course evaluation score? Explain why the two lines
should have the same slope, based on the equation of the multiple 
regression model.
    
************************************************************************

************************************************************************

The decision to call the indicator variable `gendermale` instead
of`genderfemale` has no deeper meaning. R simply codes the category that
comes first alphabetically as a \(0\). (You can change the reference
level of a categorical variable, which is the level that is coded as a
0, using the`relevel` function. Use `?relevel` to learn more.)

### Exercise 15

Create a new model called `m_bty_rank` with `gender` removed and
    `rank` added in. How does R appear to handle categorical variables
    that have more than two levels? Note that the rank variable has
    three levels: `teaching`, `tenure track`, `tenured`.
    
************************************************************************

************************************************************************

The interpretation of the coefficients in multiple regression is
slightly different from that of simple regression. The estimate for
`bty_avg` reflects how much higher a group of professors is expected to
score if they have a beauty rating that is one point higher *while
holding all other variables constant*. In this case, that translates
into considering only professors of the same rank with `bty_avg` scores
that are one point apart.

## The search for the best model

We will start with a full model that predicts professor score based on
rank, ethnicity, gender, language of the university where they got their
degree, age, proportion of students that filled out evaluations, class
size, course level, number of professors, number of credits, average
beauty rating, outfit, and picture color.

### Exercise 16

Which variable would you expect to have the highest p-value in this
model? Why? *Hint:* Think about which variable would you expect to
not have any association with the professor score.
    
    
************************************************************************

************************************************************************

Let’s run the model...

```{r}
m_full <- lm(score ~ rank + ethnicity + gender + language + 
               age + cls_perc_eval + cls_students + cls_level +
               cls_profs + cls_credits + bty_avg + pic_outfit +
               pic_color, data = evals)

get_regression_table(m_full)
get_regression_summaries(m_full)
```

### Exercise 17

Check your suspicions from the previous exercise. And support your conclusions 
with reference to the regression table.
    
************************************************************************

************************************************************************

### Exercise 18

Interpret the coefficient associated with the ethnicity variable.

************************************************************************

************************************************************************

### Exercise 19

Drop the variable with the highest p-value and re-fit the model. Did
the coefficients and significance of the other explanatory variables
change? (One of the things that makes multiple regression
interesting is that coefficient estimates depend on the other
variables that are included in the model.)
    
************************************************************************

************************************************************************

### Exercise 20


Backward-selection using p-value is a model selection process in which
we start with the full model and consecutively remove the variable 
from the model with highest p-value and refit the model, until all 
variables have p-value at most 0.05 (or another discipline specific
significance level),

Using backward-selection and p-value as the selection criterion,
determine the best model. You do not need to show all steps in your
answer, just the output for the final model. Also, write out the
linear model for predicting score based on the final model you
settle on.
    
************************************************************************

************************************************************************



### Exercise 21

The original paper describes how these data were gathered by taking
a sample of professors from the University of Texas at Austin and
including all courses that they have taught. Considering that each
row represents a course, could this new information have an impact
on any of the conditions of linear regression?


************************************************************************

************************************************************************

### Exercise 22

Based on your final model, describe the characteristics of a
professor and course at University of Texas at Austin that would be
associated with a high evaluation score.
    
************************************************************************

************************************************************************

### Exercise 23

Would you be comfortable generalizing your conclusions to apply to
professors generally (at any university)? Why or why not?
    
************************************************************************

************************************************************************
