Regression standardization with the R package stdReg

Abstract

When studying the association between an exposure and an outcome, it is common to use regression models to adjust for measured confounders. The most common models in epidemiologic research are logistic regression and Cox regression, which estimate conditional (on the confounders) odds ratios and hazard ratios. When the model has been fitted, one can use regression standardization to estimate marginal measures of association. If the measured confounders are sufficient for confounding control, then the marginal association measures can be interpreted as poulation causal effects. In this paper we describe a new R package, stdReg, that carries out regression standardization with generalized linear models (e.g. logistic regression) and Cox regression models. We illustrate the package with several examples, using real data that are publicly available.

Appendix 1: Asymptotic distribution for standardized measures

For generalized linear models, let $x_0$ and $x_1$ be fixed constants. Let $\psi =g\{\theta (x_0),\theta (x_1)\}$ be a function of $\theta (x_0)$ and $\theta (x_1)$, e.g. $\theta (x_1)-\theta (x_0)$. Define $\nu =\{\beta ,\theta (x_0),\theta (x_1),\psi \}$. The estimator $\hat{\nu }=[\hat{\beta },\hat{\theta }(x_0),\hat{\theta }(x_1),g\{\hat{\theta } (x_0),\hat{\theta } (x_1)\}]$ is an M-estimator [5] that solves the estimating equation

$$\begin{aligned} \sum _{i=1}^nU_{\nu ,i}(\nu )=\sum _{i=1}^n \left[ \begin{array}{l} U_{\beta ,i}(\beta )\\ U_{\theta (x_0),i}\{\beta ,\theta (x_0)\}\\ U_{\theta (x_1),i}\{\beta ,\theta (x_1)\}\\ U_{\psi ,i}\{\theta (x_0),\theta (x_1),\psi \} \end{array}\right] =0, \end{aligned}$$

where $U_{\beta ,i}(\beta )$ is the contribution to the maximum likelihood score function from subject i, $U_{\theta (x),i}\{\beta ,\theta (x)\}=\eta ^{-1}\{h(X=x,Z_i;\beta )\}-\theta (x)$ for $x=x_1$ and $x=x_0$, and $U_{\psi ,i}\{\theta (x_0),\theta (x_1),\psi \}=g\{\theta (x_0),\theta (x_1)\}-\psi$.

For Cox regression models, let $x_0$, $x_1$ and t be fixed constants. Let $\psi =g\{\theta (t,x_0),\theta (t,x_1)\}$ be a function of $\theta (t,x_0)$ and $\theta (t,x_1)$, e.g. $\theta (t,x_1)-\theta (t,x_0)$. Define $\nu =\{\beta ,{\varLambda }_0(t),\theta (t,x_0),\theta (t,x_1),\psi \}$. The estimator $\hat{\nu }=[\hat{\beta },\hat{{\varLambda }}_0(t),\hat{\theta }(t,x_0),\hat{\theta }(t,x_1),g\{\hat{\theta } (t,x_0),\hat{\theta } (t,x_1)\}]$ is an M-estimator [5] that solves the estimating equation

$$\begin{aligned} \sum _{i=1}^nU_{\nu ,i}(\nu )=\sum _{i=1}^n\left[ \begin{array}{l}U_{\beta ,i}(\beta )\\ U_{{\varLambda }_0(t),i}\{\beta ,{\varLambda }_0(t)\}\\ U_{\theta (t,x_0),i}\{\beta ,{\varLambda }_0(t),\theta (t,x_0)\}\\ U_{\theta (t,x_1),i}\{\beta ,{\varLambda }_0(t),\theta (t,x_1)\}\\ U_{\psi ,i}\{\theta (t,x_0),\theta (t,x_1),\psi \} \end{array}\right] =0, \end{aligned}$$

where $U_{\beta ,i}(\beta )$ is the contribution to the Cox partial likelihood score function from subject i, $U_{{\varLambda }_0(t),i}\{\beta ,{\varLambda }_0(t)\}$ is the contribution to the estimating function for Breslow’s estimator of the cumulative baseline hazard from subject i, $U_{\theta (t,x),i}\{\beta ,{\varLambda }_0(t),\theta (t,x)\}=\text {exp}[-{\varLambda }_0(t)\text {exp}\{h(X=x,Z_i;\beta )\}]-\theta (t,x)$ for $x=x_1$ and $x=x_0$, and $U_{\psi ,i}\{\theta (t,x_0),\theta (t,x_1),\psi \}=g\{\theta (t,x_0),\theta (t,x_1)\}-\psi$.

For both generalized linear models and Cox regression models it now follows from standard theory for M-estimators [5] that $n^{1/2}(\hat{\nu }-\nu )$ is asymptotically normal with mean 0 and variance given by the ‘sandwich formula’

$$\begin{aligned} {\varSigma }=E^{\prime}\left\{ \frac{\partial U_{\nu ,i}(\nu )}{\partial \nu }\right\} ^{-1}\text {var}\{U_{\nu ,i}(\nu )\}E\left\{ \frac{\partial U_{\nu ,i}(\nu )}{\partial \nu }\right\} ^{-1}. \end{aligned}$$

(5)

A consistent estimate of the variance of $\hat{\nu }$ is obtained by replacing $\nu$ in (5) with $\hat{\nu }$, and the population moments in (5) by their sample counterparts.

The sandwich formula assumes that $U_{\nu ,i}(\nu )$ and $U_{\nu ,i^{\prime}}(\nu )$ are independent, for $i\ne i^{\prime}$. When data are clustered, as in the example in ‘Standardization with generalized linear models’ section, we may define $U_{\nu ,i}(\nu )=\sum _{j=1}^{n_i}U_{\nu ,ij}(\nu )$, where $U_{\nu ,ij}(\nu )$ is the contribution to the estimating equation from subject j within cluster i, and $n_i$ is the total number of subjects in cluster i. Provided that the clusters are independent we thus have that $U_{\nu ,i}(\nu )$ and $U_{\nu ,i^{\prime}}(\nu )$ are independent as well, for $i\ne i^{\prime}$, so that the sandwich formula still applies.