class: center, middle, inverse, title-slide # Fitting the outcome model ### Lucy D’Agostino McGowan ### Wake Forest University ### 2020-07-29 (updated: 2020-07-29) --- ## Outcome Model ```r library(broom) lm(outcome ~ exposure, data = df, weights = wts) %>% tidy() ``` -- ✅ This will get us the point estimate -- ❌ This will get NOT us the correct confidence intervals -- 📦 {rsample} --- <span class = "num">1</span> <h3> Create a function to run your analysis once on a sample of your data</h3> .small[ ```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 = 1 / ifelse(exposure == 0, 1 - .fitted, .fitted)) # fit correctly bootsrapped ipw model lm(outcome ~ exposure, data = .df, weights = wts) %>% tidy() } ``` ] --- <span class = "num">1</span> <h3> Create a function to run your analysis once on a sample of your data</h3> .small[ ```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 = 1 / ifelse(exposure == 0, 1 - .fitted, .fitted)) # fit correctly bootsrapped ipw model lm(outcome ~ exposure, data = .df, weights = wts) %>% tidy() } ``` ] --- <span class = "num">1</span> <h3> Create a function to run your analysis once on a sample of your data</h3> .small[ ```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 = 1 / ifelse(exposure == 0, 1 - .fitted, .fitted)) # fit correctly bootsrapped ipw model lm(outcome ~ exposure, data = .df, weights = wts) %>% tidy() } ``` ] --- <span class = "num">1</span> <h3> Create a function to run your analysis once on a sample of your data</h3> .small[ ```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 = 1 / ifelse(exposure == 0, 1 - .fitted, .fitted)) # fit correctly bootsrapped ipw model lm(outcome ~ exposure, data = .df, weights = wts) %>% tidy() } ``` ] --- <span class = "num">1</span> <h3> Create a function to run your analysis once on a sample of your data</h3> .small[ ```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 = 1 / ifelse(exposure == 0, 1 - .fitted, .fitted)) * # fit correctly bootsrapped ipw model * lm(outcome ~ exposure, data = .df, weights = wts) %>% * tidy() } ``` ] --- <span class = "num">2</span> <h3> Use {rsample} to bootstrap our causal effect</h3> ```r library(rsample) # fit ipw model to bootstrapped samples ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>% mutate(results = map(splits, fit_ipw)) ``` --- <span class = "num">2</span> <h3> Use {rsample} to bootstrap our causal effect</h3> ```r library(rsample) # fit ipw model to bootstrapped samples *ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>% mutate(results = map(splits, fit_ipw)) ``` --- <span class = "num">2</span> <h3> Use {rsample} to bootstrap our causal effect</h3> ```r library(rsample) # fit ipw model to bootstrapped samples ipw_results <- bootstraps(df, 1000, apparent = TRUE) %>% * mutate(results = map(splits, fit_ipw)) ``` --- <span class = "num">3</span> <h3> Pull out the causal effect</h3> ```r # get t-statistic-based CIs *boot_estimate <- int_t(ipw_results, results) %>% filter(term == "exposure") ``` --- ## Your Turn <div class="countdown" id="timer_5f2190f8" style="right:0;bottom:0;" data-warnwhen="0"> <code class="countdown-time"><span class="countdown-digits minutes">07</span><span class="countdown-digits colon">:</span><span class="countdown-digits seconds">00</span></code> </div> 1. Create a function called `ipw_fit` that fits the propensity score model and the weighted outcome model for the effect between `qsmk` and `wt82_71` 2. Using the `bootstraps()` and `int_t()` functions to estimate the final effect.