10  Step 7: Association

10.0.1 Set outcome formula

out.formula <- as.formula(paste0("outcome", "~", "exposure"))
out.formula
#> outcome ~ exposure

10.0.2 Obtain OR from unadjusted model

fit <- glm(out.formula,
            data = hdps.data,
            weights = W.out$weights,
            family= binomial(link = "logit"))
fit.summary <- summary(fit)$coef["exposure",
                                 c("Estimate", 
                                   "Std. Error", 
                                   "Pr(>|z|)")]
fit.summary[2] <- sqrt(sandwich::sandwich(fit)[2,2])
require(lmtest)
conf.int <- confint(fit, "exposure", level = 0.95, method = "hc1")
fit.summary_with_ci <- c(fit.summary, conf.int)
knitr::kable(t(round(fit.summary_with_ci,2)))
Estimate Std. Error Pr(>|z|) 2.5 % 97.5 %
0.42 0.08 0 0.35 0.49
  • We are using a crude outcome model here.
  • Somewhat controversial to adjust for all (investigator-specified and all 100 proxies) covariates.

10.0.3 Obtain RD from unadjusted model

fit <- glm(out.formula,
            data = hdps.data,
            weights = W.out$weights,
            family= gaussian(link = "identity"))
fit.summary <- summary(fit)$coef["exposure",
                                 c("Estimate", 
                                   "Std. Error", 
                                   "Pr(>|t|)")]
fit.summary[2] <- sqrt(sandwich::sandwich(fit)[2,2])
require(lmtest)
conf.int <- confint(fit, "exposure", level = 0.95, method = "hc1")
fit.summary <- c(fit.summary, conf.int)
knitr::kable(t(round(fit.summary, 2)))
Estimate Std. Error Pr(>|t|) 2.5 % 97.5 %
0.08 0.01 0 0.06 0.1

(Naimi and Whitcomb 2020)