Outcome Models

# Outcome Models

**Session 16**

]

---

# Weighted Estimator

.box-1[
`$$\hat{\tau} = \frac{\sum_{i=1}^n Y_i X_i w_i}{\sum_{i=1}^nX_iw_i}-\frac{\sum_{i=1}^n Y_i (1-X_i) w_i}{\sum_{i=1}^n(1-X_i)w_i}$$`
]

# Example

```
## # A tibble: 1,100 × 5
##    coffee_town cups_of_coffee   age job   smoke
##          <dbl>          <dbl> <dbl> <chr> <chr>
##  1           1              0  12.1 none  never
##  2           0              0  14.6 none  never
##  3           1              0  11.5 none  never
##  4           0              1  18.8 easy  never
##  5           0              0  15.0 none  never
##  6           1              0  11.5 none  never
##  7           1              0  15.5 none  never
##  8           1              0  16.2 easy  never
##  9           1              0  15.7 none  never
## 10           0              0  13.1 none  never
## # … with 1,090 more rows
```
]

---

# Example

.box-1.medium[Causal Question]

.box-inv-1[What is the average causal effect of living in coffee town on the number of cups of coffee consumed across the whole population?]

---

# <svg aria-hidden="true" role="img" viewBox="0 0 640 512" style="height:1em;width:1.25em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M624 416H381.54c-.74 19.81-14.71 32-32.74 32H288c-18.69 0-33.02-17.47-32.77-32H16c-8.8 0-16 7.2-16 16v16c0 35.2 28.8 64 64 64h512c35.2 0 64-28.8 64-64v-16c0-8.8-7.2-16-16-16zM576 48c0-26.4-21.6-48-48-48H112C85.6 0 64 21.6 64 48v336h512V48zm-64 272H128V64h384v256z"/></svg> Application Exercise

1. Calculate a propensity score model estimating for the exposure `coffee_town`
2. Calculate `ate` weights using the propensity score
3. Knit commit & push to GitHub

---

# Example

```r
library(broom)

coffee_town_df <- 
  glm(coffee_town ~ age + job + smoke, 
      data = coffee_town_df,
      family = binomial()) %>%
  augment(data = coffee_town_df, type.predict = "response") %>%
  mutate(
    ate = coffee_town / .fitted + (1 - coffee_town)  / (1 - .fitted)
  )
```
]

---

# Estimating the ATE

```r
coffee_town_df %>%
  summarise(
    ate = sum(cups_of_coffee * coffee_town * ate) / 
      sum(coffee_town * ate) - 
      sum(cups_of_coffee * (1 - coffee_town) * ate) / 
      sum((1 - coffee_town) * ate))
```

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

---

# Estimating the ATE

```r
coffee_town_df %>%
  summarise(
*   ate = sum(cups_of_coffee * coffee_town * ate) /
      sum(coffee_town * ate) - 
      sum(cups_of_coffee * (1 - coffee_town) * ate) / 
      sum((1 - coffee_town) * ate))
```

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

---

# Estimating the ATE

```r
coffee_town_df %>%
  summarise(
    ate = sum(cups_of_coffee * coffee_town * ate) / 
*     sum(coffee_town * ate) -
      sum(cups_of_coffee * (1 - coffee_town) * ate) / 
      sum((1 - coffee_town) * ate))
```

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

---

# Estimating the ATE

```r
coffee_town_df %>%
  summarise(
    ate = sum(cups_of_coffee * coffee_town * ate) / 
      sum(coffee_town * ate) - 
*     sum(cups_of_coffee * (1 - coffee_town) * ate) /
      sum((1 - coffee_town) * ate))
```

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

---

# Estimating the ATE

```r
coffee_town_df %>%
  summarise(
    ate = sum(cups_of_coffee * coffee_town * ate) / 
      sum(coffee_town * ate) - 
      sum(cups_of_coffee * (1 - coffee_town) * ate) / 
*     sum((1 - coffee_town) * ate))
```

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

---

# Estimating the ATE

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

```r
lm(cups_of_coffee ~ coffee_town, 
   weights = ate, 
   data = coffee_town_df)
```

```
## 
## Call:
## lm(formula = cups_of_coffee ~ coffee_town, data = coffee_town_df, 
##     weights = ate)
## 
## Coefficients:
## (Intercept)  coffee_town  
##      2.2064       0.5656
```
]
]

---
class: subtitle subtitle-1 middl

# Warning: This may not be the right estimate! We haven't checked / iterated, this is just for example purposes

---
class: title-1 title

# Estimating the ATE

```
## # A tibble: 1 × 1
##     ate
##   <dbl>
## 1 0.566
```
]
]

```r
lm(cups_of_coffee ~ coffee_town, 
*  weights = ate,
   data = coffee_town_df)
```

```
## 
## Call:
## lm(formula = cups_of_coffee ~ coffee_town, data = coffee_town_df, 
##     weights = ate)
## 
## Coefficients:
## (Intercept)  coffee_town  
##      2.2064       0.5656
```
]
]
---

# Can we just use the CIs from `lm`?

```r
lm(cups_of_coffee ~ coffee_town, 
*  weights = ate,
   data = coffee_town_df) %>%
  tidy(conf.int = TRUE)
```

```
## # A tibble: 2 × 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)    2.21     0.0593     37.2  1.04e-196    2.09      2.32 
## 2 coffee_town    0.566    0.0869      6.51 1.15e- 10    0.395     0.736
```
]

---

# Option 1: Bootstrap

```r
fit_ipw <- function(split, ...) {
  .df <- analysis(split)
  # fit propensity score model
  propensity_model <- glm(
    exposure ~ confounder_1 + confounder_2 + ...
    family = binomial(), 
    data = .df
  )
  # calculate inverse probability weights
  .df <- propensity_model %>% 
    augment(type.predict = "response", data = .df) %>% 
    mutate(wts = exposure / .fitted + (1 - exposure) / (1  - .fitted)) 
  # fit correctly bootsrapped ipw model
  lm(outcome ~ exposure, data = .df, weights = wts) %>% 
    tidy()
}
```
]
---

# Option 1: Bootstrap

```r
*fit_ipw <- function(split, ...) {
  .df <- analysis(split)
  # fit propensity score model
  propensity_model <- glm(
    exposure ~ confounder_1 + confounder_2 + ...
    family = binomial(), 
    data = .df
  )
  # calculate inverse probability weights
  .df <- propensity_model %>% 
    augment(type.predict = "response", data = .df) %>% 
    mutate(wts = exposure / .fitted + (1 - exposure) / (1  - .fitted)) 
  # fit correctly bootsrapped ipw model
  lm(outcome ~ exposure, data = .df, weights = wts) %>% 
    tidy()
*} 
```
]

---

# Option 1: Bootstrap

```r
fit_ipw <- function(split, ...) {
* .df <- analysis(split)
  # fit propensity score model
  propensity_model <- glm(
    exposure ~ confounder_1 + confounder_2 + ...
    family = binomial(), 
    data = .df
  )
  # calculate inverse probability weights
  .df <- propensity_model %>% 
    augment(type.predict = "response", data = .df) %>% 
    mutate(wts = exposure / .fitted + (1 - exposure) / (1  - .fitted)) 
  # fit correctly bootsrapped ipw model
  lm(outcome ~ exposure, data = .df, weights = wts) %>% 
    tidy()
}
```
]

---

# Option 1: Bootstrap

```r
fit_ipw <- function(split, ...) {
  .df <- analysis(split)
* # fit propensity score model
* propensity_model <- glm(
*   exposure ~ confounder_1 + confounder_2 + ...
*   family = binomial(),
*   data = .df
* )
  # calculate inverse probability weights
  .df <- propensity_model %>% 
    augment(type.predict = "response", data = .df) %>% 
    mutate(wts = exposure / .fitted + (1 - exposure) / (1  - .fitted)) 
  # fit correctly bootsrapped ipw model
  lm(outcome ~ exposure, data = .df, weights = wts) %>% 
    tidy()
}
```
]

---

# Option 1: Bootstrap

```r
fit_ipw <- function(split, ...) {
  .df <- analysis(split)
  # fit propensity score model
  propensity_model <- glm(
    exposure ~ confounder_1 + confounder_2 + ...
    family = binomial(), 
    data = .df
  )
* # calculate inverse probability weights
* .df <- propensity_model %>%
*   augment(type.predict = "response", data = .df) %>%
*   mutate(wts = exposure / .fitted + (1 - exposure) / (1  - .fitted))
  # fit correctly bootsrapped ipw model
  lm(outcome ~ exposure, data = .df, weights = wts) %>% 
    tidy()
}
```
]

---

# Option 1: Bootstrap

```r
fit_ipw <- function(split, ...) {
  .df <- analysis(split)
  # fit propensity score model
  propensity_model <- glm(
    exposure ~ confounder_1 + confounder_2 + ...
    family = binomial(), 
    data = .df
  )
  # calculate inverse probability weights
  .df <- propensity_model %>% 
    augment(type.predict = "response", data = .df) %>% 
    mutate(wts = exposure / .fitted + (1 - exposure) / (1  - .fitted)) 
* # fit correctly bootsrapped ipw model
* lm(outcome ~ exposure, data = .df, weights = wts) %>%
*   tidy()
}
```
]

---

# Option 1: Bootstrap

```r
library(rsample)
# fit ipw model to bootstrapped samples
ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>% 
  mutate(results = map(splits, fit_ipw))
```
]
---

# Option 1: Bootstrap

```r
*library(rsample)
# fit ipw model to bootstrapped samples
ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>% 
  mutate(results = map(splits, fit_ipw))
```
]

---

# Option 1: Bootstrap

```r
library(rsample)
# fit ipw model to bootstrapped samples
*ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>%
  mutate(results = map(splits, fit_ipw))
```
]

---

# Option 1: Bootstrap

```r
library(rsample)
# fit ipw model to bootstrapped samples
ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>% 
* mutate(results = map(splits, fit_ipw))
```
]

---

# Option 1: Bootstrap

```r
# get t-statistic-based CIs
boot_estimate <- int_t(ipw_results, results) %>%
  filter(term == "exposure")
```

---

# Option 1: Bootstrap

```r
# get t-statistic-based CIs
*boot_estimate <- int_t(ipw_results, results) %>%
  filter(term == "exposure")
```

---

# Option 1: Bootstrap

```r
# get t-statistic-based CIs
boot_estimate <- int_t(ipw_results, results) %>%
* filter(term == "exposure")
```

---

1. `install.packages("rsample")`
2. Create a function called that fits the propensity score model and the weighted outcome model for the effect between `coffee_town` and `cups_of_coffee`
3. Using the `bootstraps()` and `int_t()` functions to estimate the final effect.
4. Knit commit & push to GitHub

---

# Bootstrap estimate

```
## # A tibble: 1 × 6
##   term        .lower .estimate .upper .alpha .method  
##   <chr>        <dbl>     <dbl>  <dbl>  <dbl> <chr>    
## 1 coffee_town 0.0402     0.538   1.17   0.05 student-t
```

.box-1.medium[Naive estimate]

```
## # A tibble: 2 × 4
##   term        estimate conf.low conf.high
##   <chr>          <dbl>    <dbl>     <dbl>
## 1 (Intercept)    2.21     2.09      2.32 
## 2 coffee_town    0.566    0.395     0.736
```

---

# Option 2: Robust Standard Errors

`$$u(\theta)=\begin{pmatrix}(Y-\mu_1)Xe^{-1}\\(Y-\mu_0)(1-X)(1-e)^{-1}\end{pmatrix}\\\theta = (\mu_1, \mu_0)$$`

]

--
<br>
.box-1[
`$e$` = the propensity score, this method assumes this is **fixed**
]

---

# M-Estimation

.box-1[
`$$\sqrt{n}(\hat\theta-\theta)\rightarrow^d N(0, \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-T})$$`
]

---

# Sandwich estimator

`$$\textrm{var}(\theta) = \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-T}$$`
]

<br>

`$$\mathbf{A} = -E\left[\frac{\partial \mathbf{u}(\theta)}{\partial\theta}\right]$$`
]
---

# Sandwich estimator

`$$\textrm{var}(\theta) = \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-T}$$`
]

<br>

`$$\mathbf{B} = E[\mathbf{u}(\theta)\mathbf{u}(\theta)^T]$$`
]

---

# Estimating **A**

--
<br>

# Estimating **B**

`$$\mathbf{B} = E[\mathbf{u}(\theta)\mathbf{u}(\theta)^T]$$`
]

--
<br>

.box-1[
This is your "robust" estimate of the information (the variance of the *score function* if you remember from statistical inference!)
]

---

# Option 2: Robust Standard Errors

```r
library(survey)
des <- svydesign(
  ids = ~1, 
  data = coffee_town_df,
  weights = ~ate
)
svyglm(cups_of_coffee ~ coffee_town, des) %>%
  tidy(conf.int = TRUE)
```
---

# Option 2: Robust Standard Errors

```r
*library(survey)
des <- svydesign(
  ids = ~1, 
  data = coffee_town_df,
  weights = ~ate
)
svyglm(cups_of_coffee ~ coffee_town, des) %>%
  tidy(conf.int = TRUE)
```

---
class: title title-1

# Option 2: Robust Standard Errors

```r
library(survey)
*des <- svydesign(
* ids = ~1,
* data = coffee_town_df,
* weights = ~ate
*) 
svyglm(cups_of_coffee ~ coffee_town, des) %>%
  tidy(conf.int = TRUE)
```

---

# Option 2: Robust Standard Errors

```r
library(survey)
des <- svydesign(
  ids = ~1, 
  data = coffee_town_df,
  weights = ~ate
)
*svyglm(cups_of_coffee ~ coffee_town, des) %>%
  tidy(conf.int = TRUE)
```

---

1. Use the `svyglm` function from the survey package to model the average treatment effect with the robust standard errors.
2. Knit commit & push to GitHub

---

# Option 2: Robust Standard Errors

```
## # A tibble: 2 × 4
##   term        estimate conf.low conf.high
##   <chr>          <dbl>    <dbl>     <dbl>
## 1 (Intercept)    2.21    2.15        2.26
## 2 coffee_town    0.566   0.0857      1.05
```

---
class: title title-1

# Option 3: Sandwich estimator

.box-1[
`$$u(\theta)=\begin{pmatrix}(Y-\mu_1)Xe(\mathbf{Z}, \beta)^{-1}\\(Y-\mu_0)(1-X)(1-e(\mathbf{Z}, \beta))^{-1}\\\mathbf{Z}(X-e(\mathbf{Z}, \beta))\end{pmatrix}\\\theta = (\mu_1, \mu_0, \beta)$$`
]

--
<br>
.box-1[
`$e(\mathbf{Z}, \beta)$` = the propensity score
]

---

# Option 3: Sandwich estimator

.box-1[
`$$u(\theta)=\begin{pmatrix}(Y-\mu_1)\color{orange}{Xe(\mathbf{Z}, \beta)^{-1}}\\(Y-\mu_0)\color{orange}{(1-X)(1-e(\mathbf{Z}, \beta))^{-1}}\\\mathbf{Z}(X-e(\mathbf{Z}, \beta))\end{pmatrix}\\\theta = (\mu_1, \mu_0, \beta)$$`
]

<br>

---
class: title title-1

# Option 3: Sandwich estimator

.box-1[
`$$u(\theta)=\begin{pmatrix}(Y-\mu_1)Xe(\mathbf{Z}, \beta)^{-1}\\(Y-\mu_0)(1-X)(1-e(\mathbf{Z}, \beta))^{-1}\\\color{orange}{\mathbf{Z}(X-e(\mathbf{Z}, \beta))}\end{pmatrix}\\\theta = (\mu_1, \mu_0, \beta)$$`
]

--
<br>
.box-1[
Estimating equation for the propensity score model
]

---

# M-Estimation

.box-1[
`$$\sqrt{n}(\hat\theta-\theta)\rightarrow^d N(0, \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-T})$$`
]

---

# Sandwich estimator

`$$\textrm{var}(\theta) = \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-T}$$`
]

<br>

`$$\mathbf{A} = -E\left[\frac{\partial \mathbf{u}(\theta)}{\partial\theta}\right]$$`
]
---

# Sandwich estimator

`$$\textrm{var}(\theta) = \mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-T}$$`
]

<br>

`$$\mathbf{B} = E[\mathbf{u}(\theta)\mathbf{u}(\theta)^T]$$`
]

---

# Estimating **A**

--
<br>

`$$\mathbf{A} = \begin{pmatrix}a_{11}& 0& a_{13}\\ 0 & a_{22} & a_{23}\\0&0&a_{33}\end{pmatrix}$$`
]

---

# Estimating **A**

--
<br>
.box-1[Why?]

.box-1[
`$$-E\left[\frac{\partial(Y-\mu_1)Xe(\mathbf{Z}, \beta)^{-1}}{\partial \mu_1}\right]\\=-E[-Xe(\mathbf{Z}, \beta)^{-1}]$$`
]
---

# Estimating **A**

--
<br>
.box-1[
`$a_{22} = E[(1-X) (1 - e(\mathbf{Z}, \beta))^{-1}]$`
]

--
<br>
.box-1[
`$a_{33} = E[\mathbf{ZZ}^Te(\mathbf{Z}, \beta)(1-e(\mathbf{Z}, \beta))]$`

]

---

# Estimating **A**

<br>
.box-1[
`$a_{23} = E[\mathbf{ZZ}^T(X-e(\mathbf{Z}, \beta))^2]$`
]

---

# Estimating **A**

---

# Estimating **A**

---

# Estimating **A**

---

# Estimating **A**

<br>
.box-1[
`$\hat{a}_{33} = \mathbf{ZZ}^T\hat{e}(1-\hat{e})$`

]

---

# Estimating **B**

`$$\mathbf{B} = E[\mathbf{u}(\theta)\mathbf{u}(\theta)^T]$$`
]

--
<br>
.box-1[
`$$\mathbf{B}=\begin{pmatrix}b_{11}&0&b_{13}\\0&b_{22}&b_{23}\\b_{13}^T&b_{23}^T& b_{33}\end{pmatrix}$$`
]

---
class: title-1 title
# Estimating **B**

<br>

---

<br>

<br>
.box-1[
`$b_{33}=E[\mathbf{ZZ}^Te(1-e)]$`
]

---

# Putting the variance together

<br>

<br>

---

# Doing it in R!

```r
library(PSW) 
coffee_town_df <- coffee_town_df %>%
  mutate(job_hard = ifelse(job == "hard", 1, 0),
         job_easy = ifelse(job == "easy", 1, 0),
         smoke_former = ifelse(smoke == "former", 1, 0),
         smoke_never = ifelse(smoke == "never", 1, 0)
  ) %>%
  as.data.frame()
```

---

# Doing it in R!

```r
*library(PSW)
coffee_town_df <- coffee_town_df %>%
  mutate(job_hard = ifelse(job == "hard", 1, 0),
         job_easy = ifelse(job == "easy", 1, 0),
         smoke_former = ifelse(smoke == "former", 1, 0),
         smoke_never = ifelse(smoke == "never", 1, 0)
  ) %>%
  as.data.frame()
```

---

# Doing it in R!

```r
library(PSW) 
coffee_town_df <- coffee_town_df %>%
* mutate(job_hard = ifelse(job == "hard", 1, 0),
         job_easy = ifelse(job == "easy", 1, 0),
         smoke_former = ifelse(smoke == "former", 1, 0),
         smoke_never = ifelse(smoke == "never", 1, 0)
  ) %>%
  as.data.frame()
```

---
class: title title-1

# Doing it in R!

```r
library(PSW) 
coffee_town_df <- coffee_town_df %>%
  mutate(job_hard = ifelse(job == "hard", 1, 0),
*        job_easy = ifelse(job == "easy", 1, 0),
         smoke_former = ifelse(smoke == "former", 1, 0),
         smoke_never = ifelse(smoke == "never", 1, 0)
  ) %>%
  as.data.frame()
```

---class: title title-1

# Doing it in R!

```r
library(PSW) 
coffee_town_df <- coffee_town_df %>%
  mutate(job_hard = ifelse(job == "hard", 1, 0),
         job_easy = ifelse(job == "easy", 1, 0),
*        smoke_former = ifelse(smoke == "former", 1, 0),
         smoke_never = ifelse(smoke == "never", 1, 0)
  ) %>%
  as.data.frame()
```

---
class: title title-1

# Doing it in R!

```r
library(PSW) 
coffee_town_df <- coffee_town_df %>%
  mutate(job_hard = ifelse(job == "hard", 1, 0),
         job_easy = ifelse(job == "easy", 1, 0),
         smoke_former = ifelse(smoke == "former", 1, 0),
*        smoke_never = ifelse(smoke == "never", 1, 0)
  ) %>%
  as.data.frame()
```

---
class: title title-1

# Doing it in R!

---

# Doing it in R!
.small[

```r
mod <- psw(data = coffee_town_df, 
    form.ps = 
      "coffee_town ~ age + job_hard + job_easy +
       smoke_former + smoke_never",
    weight = "ATE",
    wt = TRUE,
    out.var = "cups_of_coffee", 
    family = "gaussian")
```

]

---

# Doing it in R!
.small[

```r
*mod <- psw(data = coffee_town_df,
    form.ps = 
      "coffee_town ~ age + job_hard + job_easy +
       smoke_former + smoke_never",
    weight = "ATE",
    wt = TRUE,
    out.var = "cups_of_coffee", 
    family = "gaussian")
```

]
---

# Doing it in R!
.small[

```r
mod <- psw(data = coffee_town_df, 
*   form.ps =
*     "coffee_town ~ age + job_hard + job_easy +
*      smoke_former + smoke_never",
    weight = "ATE",
    wt = TRUE,
    out.var = "cups_of_coffee", 
    family = "gaussian")
```

]

---

# Doing it in R!
.small[

```r
mod <- psw(data = coffee_town_df, 
    form.ps = 
      "coffee_town ~ age + job_hard + job_easy +
       smoke_former + smoke_never",
*   weight = "ATE",
    wt = TRUE,
    out.var = "cups_of_coffee", 
    family = "gaussian")
```

]
---

# Doing it in R!
.small[

```r
mod <- psw(data = coffee_town_df, 
    form.ps = 
      "coffee_town ~ age + job_hard + job_easy +
       smoke_former + smoke_never",
    weight = "ATE",
*   wt = TRUE,
    out.var = "cups_of_coffee", 
    family = "gaussian")
```

]

---

# Doing it in R!
.small[

```r
mod <- psw(data = coffee_town_df, 
    form.ps = 
      "coffee_town ~ age + job_hard + job_easy +
       smoke_former + smoke_never",
    weight = "ATE",
    wt = TRUE,
*   out.var = "cups_of_coffee",
    family = "gaussian")
```

]

---

# Doing it in R!
.small[

```r
mod <- psw(data = coffee_town_df, 
    form.ps = 
      "coffee_town ~ age + job_hard + job_easy +
       smoke_former + smoke_never",
    weight = "ATE",
    wt = TRUE,
    out.var = "cups_of_coffee", 
*   family = "gaussian")
```

]

---

# Doing it in R!
.small[

```r
mod$est.wt
```

```
## [1] 0.5655656
```

```r
# lb
mod$est.wt - 1.96 * mod$std.wt
```

```
## [1] 0.03946954
```

```r
# ub
mod$est.wt + 1.96 * mod$std.wt
```

```
## [1] 1.091662
```
]

---

1. `install.packages("PSW")`
2. Calculate the large sample variance sandwich estimator using the `psw` function
3. Knit commit & push to GitHub