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.ci <- confint(fit, level = 0.95)["exposure", ]
fit.summary_with_ci <- c(fit.summary, fit.ci)
round(fit.summary_with_ci,2) 
#>   Estimate Std. Error   Pr(>|z|)      2.5 %     97.5 % 
#>       0.42       0.04       0.00       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)
round(fit.summary, 2) 
#>   Estimate Std. Error   Pr(>|t|)      2.5 %     97.5 % 
#>       0.08       0.01       0.00       0.06       0.10

(Naimi and Whitcomb 2020)