POL90: Statistics

# POL90: Statistics
## Missing Data
### Prof. Wasow, Politics</br>Pomona College
### 2022-04-11

---

.regression12 table {
  font-size: 12px;     
}

.regression14 table {
  font-size: 14px;     
}

</style>

# Announcements

- Assignments
  - PS09
  - Report 3

- Missing Data
  - Handout

]

---
# Schedule

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Week </th>
   <th style="text-align:left;"> Date </th>
   <th style="text-align:left;"> Day </th>
   <th style="text-align:left;"> Title </th>
   <th style="text-align:right;"> Chapter </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 11 </td>
   <td style="text-align:left;"> Mar 28 </td>
   <td style="text-align:left;"> Mon </td>
   <td style="text-align:left;"> Interaction terms </td>
   <td style="text-align:right;"> 9 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 11 </td>
   <td style="text-align:left;"> Mar 30 </td>
   <td style="text-align:left;"> Wed </td>
   <td style="text-align:left;"> Interaction terms </td>
   <td style="text-align:right;"> 9 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 12 </td>
   <td style="text-align:left;"> Apr 4 </td>
   <td style="text-align:left;"> Mon </td>
   <td style="text-align:left;"> Logistic regression </td>
   <td style="text-align:right;"> 20 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 12 </td>
   <td style="text-align:left;"> Apr 6 </td>
   <td style="text-align:left;"> Wed </td>
   <td style="text-align:left;"> Logistic regression </td>
   <td style="text-align:right;"> 20 </td>
  </tr>
  <tr>
   <td style="text-align:right;color: black !important;background-color: yellow !important;"> 13 </td>
   <td style="text-align:left;color: black !important;background-color: yellow !important;"> Apr 11 </td>
   <td style="text-align:left;color: black !important;background-color: yellow !important;"> Mon </td>
   <td style="text-align:left;color: black !important;background-color: yellow !important;"> Missing data </td>
   <td style="text-align:right;color: black !important;background-color: yellow !important;"> Handout </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 13 </td>
   <td style="text-align:left;"> Apr 13 </td>
   <td style="text-align:left;"> Wed </td>
   <td style="text-align:left;"> Missing data </td>
   <td style="text-align:right;"> Handout </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 14 </td>
   <td style="text-align:left;"> Apr 18 </td>
   <td style="text-align:left;"> Mon </td>
   <td style="text-align:left;"> Matching </td>
   <td style="text-align:right;"> Handout </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 14 </td>
   <td style="text-align:left;"> Apr 20 </td>
   <td style="text-align:left;"> Wed </td>
   <td style="text-align:left;"> Matching </td>
   <td style="text-align:right;"> Handout </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 15 </td>
   <td style="text-align:left;"> Apr 25 </td>
   <td style="text-align:left;"> Mon </td>
   <td style="text-align:left;"> Causal inference: Panel data </td>
   <td style="text-align:right;"> Handout </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 15 </td>
   <td style="text-align:left;"> Apr 27 </td>
   <td style="text-align:left;"> Wed </td>
   <td style="text-align:left;"> Causal inference: Natural Experiments </td>
   <td style="text-align:right;"> Dunning </td>
  </tr>
</tbody>
</table>

---
## Assignment schedule

---

# Missing Data: Introduction

---
## Why do we care about missing data?

.vertical-center[
<img src="images/missingness_NA_drop.png" width="100%" style="display: block; margin: auto;" />
]

---
## Why do we care about missing data?

.vertical-center[
<img src="images/missingness_too_many_NAs.png" width="100%" style="display: block; margin: auto;" />
]

---
## Why do we care about missing data?

.vertical-center[
<img src="images/edstem_filtering_data.png" width="100%" style="display: block; margin: auto;" />
]

---
## Why do we care about missing data?
.vertical-center[
.large[

- Missingness is common to all real world data
  
  - What risks are there to "listwise deletion"?
]
]

---
.vertical-center[
<img src="images/listwise_deletion_is_evil.png" width="80%" style="display: block; margin: auto;" />
]

.footnote["Listwise deletion is evil: what to do about missing data in political science," http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.199.5654&rep=rep1&type=pdf]

---
## Why do we care about missing data?

- Different types of missingness
  - Missing at random might not bias results
  - Missingness that's systematic can lead to selection bias
    - Example: study of voting behavior biased by phone use
  
- Loss of useful data  
  - "Listwise deletion discards one-third of cases on average, which deletes both the few nonresponses and the many responses in those cases." 
    - King, Honaker, Joseph, Scheve (2001)

- Loss of *N*, loss of statistical power

]

---
class: center, middle

# Example: `airquality`

---
## Case `airquality` data

- This data set has daily air quality measurements in New York from May to September 1973 over a period of 5 months.
 
- Daily readings of the following air quality values for May 1, 1973 (a Tuesday) to September 30, 1973.
  - **Ozone**: Mean ozone in parts per billion from 1300 to 1500 hours at Roosevelt Island
  - **Solar.R**: Solar radiation in Langleys in the frequency band 4000–7700 Angstroms from 0800 to 1200 hours at Central Park
  - **Wind**: Average wind speed in miles per hour at 0700 and 1000 hours at LaGuardia Airport
  - **Temp**: Maximum daily temperature in degrees Fahrenheit at La Guardia Airport.

.footnote[Source: The data were obtained from the New York State Department of Conservation (ozone data) and the National Weather Service (meteorological data). Chambers, J. M., Cleveland, W. S., Kleiner, B. and Tukey, P. A. (1983) *Graphical Methods for Data Analysis*. Belmont, CA: Wadsworth.]

---
## Load `airquality` data

```r
# https://github.com/ropenscilabs/skimr
# quick overview of data
data(airquality) # load data
airquality <- airquality %>% clean_names()
dim(airquality)
```

```
[1] 153   6
```

```r
# look for missingness
head(airquality, 10)
```

```
   ozone solar_r wind temp month day
1     41     190  7.4   67     5   1
2     36     118  8.0   72     5   2
3     12     149 12.6   74     5   3
4     18     313 11.5   62     5   4
*5     NA      NA 14.3   56     5   5
*6     28      NA 14.9   66     5   6
7     23     299  8.6   65     5   7
8     19      99 13.8   59     5   8
9      8      19 20.1   61     5   9
*10    NA     194  8.6   69     5  10
```

---

# Missing Data: 
# Exploratory Analysis

---

## Many packages now help with missingness

```r
# viewing
install.packages("skimr")
install.packages("naniar")
install.packages("dataMaid")

# recoding
install.packages("tidyr")
install.packages("simputation")

# advanced: multiple imputation / panel imputation
install.packages("Amelia")
install.packages("mi")
install.packages("mice")
install.packages("panelView")
```

---
## Use `skim_without_charts` for Latex / pdfs

```r
# another way to summarize data
*skimr::skim_without_charts(airquality)
```

<table style='width: auto;'
        class='table table-condensed'>
<caption>Data summary</caption>
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:left;">   </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Name </td>
   <td style="text-align:left;"> airquality </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Number of rows </td>
   <td style="text-align:left;"> 153 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Number of columns </td>
   <td style="text-align:left;"> 6 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> _______________________ </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Column type frequency: </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> numeric </td>
   <td style="text-align:left;"> 6 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ________________________ </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Group variables </td>
   <td style="text-align:left;"> None </td>
  </tr>
</tbody>
</table>

**Variable type: numeric**

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> skim_variable </th>
   <th style="text-align:right;"> n_missing </th>
   <th style="text-align:right;"> complete_rate </th>
   <th style="text-align:right;"> mean </th>
   <th style="text-align:right;"> sd </th>
   <th style="text-align:right;"> p0 </th>
   <th style="text-align:right;"> p25 </th>
   <th style="text-align:right;"> p50 </th>
   <th style="text-align:right;"> p75 </th>
   <th style="text-align:right;"> p100 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> ozone </td>
   <td style="text-align:right;"> 37 </td>
   <td style="text-align:right;"> 0.76 </td>
   <td style="text-align:right;"> 42.13 </td>
   <td style="text-align:right;"> 32.99 </td>
   <td style="text-align:right;"> 1.0 </td>
   <td style="text-align:right;"> 18.0 </td>
   <td style="text-align:right;"> 31.5 </td>
   <td style="text-align:right;"> 63.25 </td>
   <td style="text-align:right;"> 168.0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> solar_r </td>
   <td style="text-align:right;"> 7 </td>
   <td style="text-align:right;"> 0.95 </td>
   <td style="text-align:right;"> 185.93 </td>
   <td style="text-align:right;"> 90.06 </td>
   <td style="text-align:right;"> 7.0 </td>
   <td style="text-align:right;"> 115.8 </td>
   <td style="text-align:right;"> 205.0 </td>
   <td style="text-align:right;"> 258.75 </td>
   <td style="text-align:right;"> 334.0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> wind </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1.00 </td>
   <td style="text-align:right;"> 9.96 </td>
   <td style="text-align:right;"> 3.52 </td>
   <td style="text-align:right;"> 1.7 </td>
   <td style="text-align:right;"> 7.4 </td>
   <td style="text-align:right;"> 9.7 </td>
   <td style="text-align:right;"> 11.50 </td>
   <td style="text-align:right;"> 20.7 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> temp </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1.00 </td>
   <td style="text-align:right;"> 77.88 </td>
   <td style="text-align:right;"> 9.47 </td>
   <td style="text-align:right;"> 56.0 </td>
   <td style="text-align:right;"> 72.0 </td>
   <td style="text-align:right;"> 79.0 </td>
   <td style="text-align:right;"> 85.00 </td>
   <td style="text-align:right;"> 97.0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> month </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1.00 </td>
   <td style="text-align:right;"> 6.99 </td>
   <td style="text-align:right;"> 1.42 </td>
   <td style="text-align:right;"> 5.0 </td>
   <td style="text-align:right;"> 6.0 </td>
   <td style="text-align:right;"> 7.0 </td>
   <td style="text-align:right;"> 8.00 </td>
   <td style="text-align:right;"> 9.0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> day </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1.00 </td>
   <td style="text-align:right;"> 15.80 </td>
   <td style="text-align:right;"> 8.86 </td>
   <td style="text-align:right;"> 1.0 </td>
   <td style="text-align:right;"> 8.0 </td>
   <td style="text-align:right;"> 16.0 </td>
   <td style="text-align:right;"> 23.00 </td>
   <td style="text-align:right;"> 31.0 </td>
  </tr>
</tbody>
</table>
]

---
## `naniar` also good for visualizing missingness

```r
naniar::gg_miss_var(airquality) + 
  theme_bw(base_size = 20)
```

---
## Why `naniar`? Nick Tierney says “NA-in-R”

.footnote[Photo credit: https://www.njtierney.com/about/]
---
## `naniar` offers multiple visulizations

```r
naniar::vis_miss(airquality)
```

---
## Can also easily facet (here by month)

```r
# What does facetting by month suggest about missingness?
naniar::gg_miss_var(airquality, facet = month) + theme_bw(base_size = 20)
```

---

# Missing Data: Modeling
# and Bivariate Relationships

---

## What happens with missingness in regression?

```r
# what happens with missingness in regression?
lm(ozone ~ solar_r, data = airquality) %>% summary()
```

```

Call:
lm(formula = ozone ~ solar_r, data = airquality)

Residuals:
   Min     1Q Median     3Q    Max 
-48.29 -21.36  -8.86  16.37 119.14

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  18.5987     6.7479    2.76  0.00686 ** 
solar_r       0.1272     0.0328    3.88  0.00018 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31.3 on 109 degrees of freedom
*  (42 observations deleted due to missingness)
Multiple R-squared:  0.121,	Adjusted R-squared:  0.113 
F-statistic: 15.1 on 1 and 109 DF,  p-value: 0.000179
```
]

---
## What about scatter plots? (Note warning)

```r
# https://github.com/njtierney/naniar
airquality %>% ggplot() + aes(x = solar_r, y = ozone) + geom_point()
```

```
Warning: Removed 42 rows containing missing values (geom_point).
```

---

## Alternate approaches with `naniar`

```r
# geom_miss_point shifts 
# missing values to 
# 10% below minimum 
airquality %>% 
  ggplot() +
    aes(x = solar_r, 
        y = ozone) +
*   geom_miss_point()
```
]

.right-plot[
<img src="week12_01_files/figure-html/visualizing_missingness2_plot-1.png" width="100%" style="display: block; margin: auto;" />
]

---
## How can `naniar` plot missing values?

```r
# geom_miss_point shifts 
# missing values to 
# 10% below minimum 
airquality %>% 
  ggplot() +
    aes(x = solar_r, 
        y = ozone) +
*   geom_miss_point()
```
]

.right-plot[
<img src="week12_01_files/figure-html/unnamed-chunk-20-1.png" width="100%" style="display: block; margin: auto;" />
]

---
## How can `naniar` plot missing values?

```r
# geom_miss_point shifts 
# missing values to 
# 10% below minimum 
airquality %>% 
  ggplot() +
    aes(x = solar_r, 
        y = ozone) +
*   geom_miss_point()
```
]

.right-plot[
<img src="week12_01_files/figure-html/unnamed-chunk-22-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Missing Data: Imputation

---
## What can we do about missingness?

```r
# create new data frame
air2 <- airquality %>% 
  mutate(
    ozone_orig   = ozone,
*   ozone_NA     = ifelse(is.na(ozone), "NA", "!NA"),
*   ozone_impute = ifelse(is.na(ozone), mean(ozone, na.rm = TRUE), ozone)
    )

# view mean 
mean(air2$ozone, na.rm = TRUE)
```

```
[1] 42.13
```

```r
# view imputation
air2 %>% select(ozone_orig, ozone_NA, ozone_impute) %>% head()
```

```
  ozone_orig ozone_NA ozone_impute
1         41      !NA        41.00
2         36      !NA        36.00
3         12      !NA        12.00
4         18      !NA        18.00
*5         NA       NA        42.13
6         28      !NA        28.00
```

---
## Pros & cons of imputing the mean?

```r
air2 %>% 
  ggplot() +
  aes(x = solar_r, 
      y = ozone_impute, 
      color = ozone_NA) +
    geom_point()
```
]

.right-plot[
<img src="week12_01_files/figure-html/unnamed-chunk-24-1.png" width="100%" style="display: block; margin: auto;" />
]

---
## How about a conditional mean?

```r
# could we just use regression?
lm1 <- lm(ozone ~ solar_r, data = air2)
summary(lm1)
```

```

Call:
lm(formula = ozone ~ solar_r, data = air2)

Residuals:
   Min     1Q Median     3Q    Max 
-48.29 -21.36  -8.86  16.37 119.14

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  18.5987     6.7479    2.76  0.00686 ** 
*solar_r       0.1272     0.0328    3.88  0.00018 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31.3 on 109 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared:  0.121,	Adjusted R-squared:  0.113 
F-statistic: 15.1 on 1 and 109 DF,  p-value: 0.000179
```

---
## Model allows us to impute some values for `ozone`

```r
equatiomatic::extract_eq(lm1, use_coefs = FALSE)
```

$$
\operatorname{ozone} = \alpha + \beta_{1}(\operatorname{solar\_r}) + \epsilon
$$

```r
equatiomatic::extract_eq(lm1, use_coefs = TRUE)
```

$$
\operatorname{\widehat{ozone}} = 18.6 + 0.13(\operatorname{solar\_r})
$$
]

---
## We can use a model to predict values of `ozone`

```r
air2 <- air2 %>% 
  mutate(
    ozone_bivariate = predict(lm1, data.frame(solar_r = air2$solar_r)),
    ozone_impute    = ifelse(is.na(ozone), ozone_bivariate, ozone)
)
air2 %>% select(solar_r, ozone, ozone_bivariate, ozone_impute) %>% head(10)
```

```
   solar_r ozone ozone_bivariate ozone_impute
1      190    41           42.76        41.00
2      118    36           33.60        36.00
3      149    12           37.55        12.00
4      313    18           58.40        18.00
*5       NA    NA              NA           NA
*6       NA    28              NA        28.00
7      299    23           56.62        23.00
8       99    19           31.19        19.00
9       19     8           21.01         8.00
*10     194    NA           43.27        43.27
```

---
## Visualizing `ozone` ~ `solar_r`

```r
# could we use regression?

air2 %>% 
  ggplot() +
    aes(x = solar_r, 
        y = ozone_impute,
        color = ozone_NA) +
    geom_smooth(method = "lm", se = FALSE, col = "gray80") +
    geom_point() 
```

]

.right-plot[
<img src="week12_01_files/figure-html/unnamed-chunk-29-1.png" width="100%" style="display: block; margin: auto;" />
]

---
## Bivariate linear model not great, why?

```r
# could we use regression?

air2 %>% 
  ggplot() +
    aes(x = solar_r, 
        y = ozone_impute,
        color = ozone_NA) +
    geom_smooth(method = "lm", se = FALSE, col = "gray80") +
    geom_point() 
```
]

.right-plot[
<img src="week12_01_files/figure-html/unnamed-chunk-31-1.png" width="100%" style="display: block; margin: auto;" />
]

---
## Bivariate linear model not great, why?

```r
# could we use regression?

air2 %>% 
  ggplot() +
    aes(x = solar_r, 
        y = ozone_impute,
        color = ozone_NA) +
    geom_smooth(method = "lm", se = FALSE, col = "gray80") +
    geom_point() 
```
]

.right-plot[
<img src="week12_01_files/figure-html/unnamed-chunk-32-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Missing Data: Imputation
# with Multiple Regression

---
## How else can we model `ozone`?

```r
lm2 <- 
  lm(ozone ~ solar_r + temp, 
     data = airquality)

stargazer(lm2,
  header = FALSE,
  omit.stat = c("f", "ser", "adj.rsq"),
  single.row = TRUE,
  digits = 2,
  type = 'html')
```
]
.right-table[

<table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
<tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
<tr><td style="text-align:left"></td><td>ozone</td></tr>
<tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">solar_r</td><td>0.06<sup>**</sup> (0.03)</td></tr>
<tr><td style="text-align:left">temp</td><td>2.28<sup>***</sup> (0.25)</td></tr>
<tr><td style="text-align:left">Constant</td><td>-145.70<sup>***</sup> (18.45)</td></tr>
<tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>111</td></tr>
<tr><td style="text-align:left">R<sup>2</sup></td><td>0.51</td></tr>
<tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
</table>
]

---
## Multiple regression model for `ozone`

```r
equatiomatic::extract_eq(lm2, use_coefs = FALSE)
```

`$$\text{ozone} = \alpha + \beta_{1}(\text{solar_r}) + \beta_{2}(\text{temp}) + \epsilon$$`

```r
equatiomatic::extract_eq(lm2, use_coefs = TRUE)
```

$$
\operatorname{\widehat{ozone}} = -145.7 + 0.06(\operatorname{solar\_r}) + 2.28(\operatorname{temp})
$$

---
## Try `simputation` with multiple regression

```r
# https://cran.r-project.org/web/packages/simputation/vignettes/intro.html

# add temp to model
air_imp <- air2 %>%
*   simputation::impute_lm(ozone ~ solar_r + temp)

# output is a data frame
dim(air_imp)
```

```
[1] 153  13
```

---
## Now `ozone` ~ `solar_r` + `temp`

```r
air_imp %>% 
* filter(ozone_NA == "NA") %>%
  select(solar_r, temp, ozone_orig, ozone) %>% head(15)
```

```
   solar_r temp ozone_orig   ozone
1       NA   56         NA      NA
2      194   69         NA  22.590
3       66   57         NA -12.061
4      266   58         NA   1.639
5       NA   57         NA      NA
6      286   78         NA  48.351
7      287   74         NA  39.294
8      242   67         NA  20.775
9      186   84         NA  56.310
10     220   85         NA  60.531
11     264   79         NA  49.373
12     273   87         NA  68.114
13     259   93         NA  80.986
14     250   92         NA  78.193
15     332   80         NA  55.535
```

---
## Try `simputation` with multiple regression

```r
air_imp %>%
    ggplot() + aes(x = solar_r, y = ozone, colour = ozone_NA) + 
    geom_miss_point()
```

---
## Could add additional terms like `wind`

```r
# impute data
air_imp <- air2 %>%
*   simputation::impute_lm(ozone ~ solar_r + temp + wind)

dim(air_imp)
```

```
[1] 153  13
```

```r
air_imp %>% select(solar_r, temp, wind, ozone, ozone_bivariate) %>% head(10)
```

```
   solar_r temp wind ozone ozone_bivariate
1      190   67  7.4 41.00           42.76
2      118   72  8.0 36.00           33.60
3      149   74 12.6 12.00           37.55
4      313   62 11.5 18.00           58.40
5       NA   56 14.3    NA              NA
6       NA   66 14.9 28.00              NA
7      299   65  8.6 23.00           56.62
8       99   59 13.8 19.00           31.19
9       19   61 20.1  8.00           21.01
10     194   69  8.6 32.59           43.27
```

---
## Scatterplot with imputed ozone