Robust variance estimation for covariate-adjusted unconditional treatment effect in randomized clinical trials with binary outcomes

Abstract

To improve the precision of estimation and power of testing hypothesis for an unconditional treatment effect in randomized clinical trials with binary outcomes, researchers and regulatory agencies recommend using g-computation as a reliable method of covariate adjustment. However, the practical application of g-computation is hindered by the lack of an explicit robust variance formula that can be used for different unconditional treatment effects of interest. To fill this gap, we provide explicit and robust variance estimators for g-computation estimators and demonstrate through simulations that the variance estimators can be reliably applied in practice.

Keywords:

• G-computation
• model-assisted
• nonlinear covariate adjustment
• risk difference
• logistic regression
• standardization

1. Introduction

In randomized clinical trials, adjusting for baseline covariates has been advocated as a way to improve the precision of estimating and power of testing treatment effects (Freedman, 2008; Lin, 2013; Tsiatis et al., 2008; Yang & Tsiatis, 2001; Ye et al., 2023, 2022). We focus on binary outcomes in this article. When a logistic model is used as a working model for baseline covariate adjustment, the g-computation (Freedman, 2008; Moore & van der Laan, 2009) provides asymptotically normal estimators of unconditional treatment effects such as the risk difference, relative risk and odds ratio, regardless of whether the logistic model is correct or not. In May 2021, the US Food and Drug Administration released a draft guidance (FDA, 2021) for the use of covariates in the analysis of randomized clinical trials, and recommended the g-computation as a ‘statistically reliable method of covariate adjustment for an unconditional treatment effect with binary outcomes '.

However, to the best of our knowledge, no explicit robust variance estimation formula for g-computation is currently available that can be used for inference on different unconditional treatment effects of interest. Moreover, some existing variance estimation formulas in the literature, such as the formula in Ge et al. (2011) for risk difference and two treatment arms, are model-based and do not fit the model-robust inference paradigm. Additionally, the formula in Ge et al. (2011) does not take into account a source of variability due to covariates and nonlinearity of logistic model, which can lead to confidence intervals with insufficient coverage probabilities.

The purpose of this article is to fill this gap by providing explicit and robust variance estimators for g-computation estimators. Our simulations demonstrate that the provided variance estimators can be reliably applied in practice.

2. Robust variance estimation

Consider a k-arm trial with n subjects. For each subject i, let ${\mathit{A}}_i$ be the k-dimensional treatment indicator vector that equals ${\mathit{a}}_t$ if patient i receives treatment t for $t=1,\dots ,k$, where ${\mathit{a}}_t$ denotes the k-dimensional vector whose tth component is 1 and other components are 0, $Y_i^{\left(t\right)}$ be the binary potential outcome under treatment t, and ${\mathit{X}}_i$ be the baseline covariate vector for adjustment. The observed outcome is $Y_i=Y_i^{\left(t\right)}$ if and only if ${\mathit{A}}_i={\mathit{a}}_t$. We consider simple randomization where ${\mathit{A}}_i$ is completely random with known ${\pi }_t=P({\mathit{A}}_i={\mathit{a}}_t)$, ${\pi }_t>0$ and $\sum _{t=1}^k{\pi }_t=1$. We assume that $(Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)},{\mathit{A}}_i,{\mathit{X}}_i),i=1,\dots ,n$, are independent and identically distributed with finite second order moments. To simplify the notation, we drop the subscript i when referring to a generic subject from the population. Write the unconditional response means as ${\theta }_t=\mathrm{E}\left(Y^{\left(t\right)}\right)$ and $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_k{)}^{\mathrm{\top }}$, where the superscript ⊤ denotes the transpose of a vector throughout. The target parameter is a given contrast of the unconditional response mean vector $\mathit{\theta }$ denoted as $f\left(\mathit{\theta }\right)$, such as the risk difference ${\theta }_t-{\theta }_s$, risk ratio ${\theta }_t/{\theta }_s$, and odds ratio $\frac{{\theta }_t/(1-{\theta }_t)}{{\theta }_s/(1-{\theta }_s)}$ between two treatment arms t and s.

Throughout this article, we consider the g-computation procedure that fits a working logistic model $E(Y\mid \mathit{A},\mathit{X})=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}({\text{β}}_A^{\mathrm{\top }}\mathit{A}+{\text{β}}_X^{\mathrm{\top }}\mathit{X})$, where $\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}\left(x\right)=\exp \left(x\right)/\{1+\exp (x\left)\right\}$, and ${\text{β}}_A$ and ${\text{β}}_X$ are unknown parameter vectors (FDA, 2021). The logistic model does not need to be correct and is only used as an intermediate step to obtain g-computation estimators. Let ${\widehat{\text{β}}}_A$ and ${\widehat{\text{β}}}_X$ be the maximum likelihood estimators of ${\text{β}}_A$ and ${\text{β}}_X$, respectively, under the working logistic model. Then, ${\hat{\mu }}_t\left({\mathit{X}}_i\right)=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}({\widehat{\text{β}}}_A^{\mathrm{\top }}{\mathit{a}}_t+{\widehat{\text{β}}}_X^{\mathrm{\top }}{\mathit{X}}_i)$ is the predicted probability of response under treatment t. The g-computation estimator of $\mathit{\theta }$ is $\hat{\mathit{\theta }}=({\hat{\theta }}_1,\dots ,{\hat{\theta }}_k{)}^{\mathrm{\top }}$ with ${\hat{\theta }}_t=n^{-1}\sum _{i=1}^n{\hat{\mu }}_t\left({\mathit{X}}_i\right)$, and of a given contrast $f\left(\mathit{\theta }\right)$ is $f\left(\hat{\mathit{\theta }}\right)$. Hence, the g-computation takes a summary-then-contrast approach ((2019), R1).

Next, we derive the asymptotic distribution of the g-computation estimator $\hat{\mathit{\theta }}$ and apply the delta method to obtain the asymptotic distribution of the g-computation estimator $f\left(\hat{\mathit{\theta }}\right)$. As the logistic regression uses a canonical link, the first-order conditions of the maximum likelihood estimation ensure that, for $t=1,\dots ,k$, $\sum _{i=1}^{n}I({\mathit{A}}_i={\mathit{a}}_t)\{Y_i^{\left(t\right)}-{\hat{\mu }}_t({\mathit{X}}_i\left)\right\}=0,$ where $I({\mathit{A}}_i={\mathit{a}}_t)$ is the indicator of ${\mathit{A}}_i={\mathit{a}}_t$. Hence, the g-computation estimator is equal to ${\hat{\theta }}_t=\frac{1}{n}\sum _{i=1}^n{\hat{\mu }}_t\left({\mathit{X}}_i\right)=\frac{1}{n}\sum _{i=1}^n[\frac{I({\mathit{A}}_i={\mathit{a}}_t)}{{\hat{\pi }}_t}\{Y_i^{\left(t\right)}-{\hat{\mu }}_t({\mathit{X}}_i\left)\right\}+{\hat{\mu }}_t({\mathit{X}}_i\left)\right],$ where ${\hat{\pi }}_t=n_t/n$ and $n_t$ is the number of subjects assigned to treatment t. Since ${\mathit{A}}_i$'s are assigned completely at random, ${\hat{\pi }}_t$ and ${\hat{\mu }}_t\left(\mathit{x}\right)$ can converge to ${\pi }_t$ and ${\mu }_t\left(\mathit{x}\right)$ with $n^{-1/2}$ rate, respectively, where $\mathit{x}$ is a fixed point and ${\mu }_t\left(\mathit{x}\right)$ is a function not necessarily equal to $E(Y^{\left(t\right)}\mid \mathit{X}=\mathit{x})$ under model misspecification but satisfies $E\{Y_i^{\left(t\right)}-{\mu }_t({\mathit{X}}_i\left)\right\}=0$ due to the above first-order conditions, $t=1,\dots ,k$. Then, by Kennedy (2016) and Chernozhukov et al. (2017), ${\hat{\theta }}_t=\frac{1}{n}\sum _{i=1}^n[\frac{I({\mathit{A}}_i={\mathit{a}}_t)}{{\pi }_t}\{Y_i^{\left(t\right)}-{\mu }_t({\mathit{X}}_i\left)\right\}+{\mu }_t({\mathit{X}}_i\left)\right]+o_p\left(n^{-1/2}\right),$ where $o_p\left(n^{-1/2}\right)$ denotes the remaining term multiplied by $n^{1/2}$ converges to 0 in probability. Therefore, an application of the central limit theorem shows that, regardless of whether the working model is correct or not, $\begin{array}{r}\sqrt{n}(\hat{\mathit{\theta }}-\mathit{\theta })⟹dN(\mathbf{0},\mathit{V}),\mathit{V}=\left(\begin{array}{cccc}{v}_{11}& v_{12}& \dots & v_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ v_{1k}& v_{2k}& \dots & v_{kk}\end{array}\right),\end{array}$ where $⟹d$ denotes convergence in distribution, $\mathbf{0}$ is the k-dimensional vector of zeros, and $\begin{array}{rl}{v}_{tt}& ={\pi }_t^{-1}\mathrm{v}\mathrm{a}\mathrm{r}\{Y^{\left(t\right)}-{\mu }_t(\mathit{X}\left)\right\}+2\mathrm{c}\mathrm{o}\mathrm{v}\{Y^{\left(t\right)},{\mu }_t(\mathit{X}\left)\right\}-\mathrm{v}\mathrm{a}\mathrm{r}\left\{{\mu }_t\right(\mathit{X}\left)\right\},t=1,\dots ,k,\\ v_{ts}& =\mathrm{c}\mathrm{o}\mathrm{v}\{Y^{\left(t\right)},{\mu }_s(\mathit{X}\left)\right\}+\mathrm{c}\mathrm{o}\mathrm{v}\{Y^{\left(s\right)},{\mu }_t(\mathit{X}\left)\right\}-\mathrm{c}\mathrm{o}\mathrm{v}\left\{{\mu }_t\right(\mathit{X}),{\mu }_s(\mathit{X}\left)\right\},1\le t<s\le k.\end{array}$ By the delta method, when $f\left(\mathit{\theta }\right)$ is differentiable at $\mathit{\theta }$ with partial derivative vector $\mathrm{\nabla }f\left(\mathit{\theta }\right)$, we have $\sqrt{n}\left\{f\right(\hat{\mathit{\theta }})-f(\mathit{\theta }\left)\right\}⟹dN(0,\{\mathrm{\nabla }f\left(\mathit{\theta }\right){\}}^{\mathrm{\top }}\mathit{V}\left\{\mathrm{\nabla }f\right(\mathit{\theta }\left)\right\}).$ Some examples are: $\begin{array}{r}\left\{\mathrm{\nabla }f\right(\mathit{\theta }\left){\}}^{\mathrm{\top }}\mathit{V}\right\{\mathrm{\nabla }f\left(\mathit{\theta }\right)\}=\{\begin{array}{ll}{v}_{tt}-2v_{ts}+v_{ss},& f\left(\mathit{\theta }\right)={\theta }_t-{\theta }_s,\\ \frac{v_{tt}}{{\theta }_t^2}-\frac{2v_{ts}}{{\theta }_t{\theta }_s}+\frac{v_{ss}}{{\theta }_s^2},& f\left(\mathit{\theta }\right)=\log \frac{{\theta }_t}{{\theta }_s},\\ \frac{v_{tt}}{{\theta }_t^2(1-{\theta }_t{)}^2}-\frac{2v_{ts}}{{\theta }_t(1-{\theta }_t){\theta }_s(1-{\theta }_s)}+\frac{v_{ss}}{{\theta }_s^2(1-{\theta }_s{)}^2}& f\left(\mathit{\theta }\right)=\log \frac{{\theta }_t/(1-{\theta }_t)}{{\theta }_s/(1-{\theta }_s)}.\end{array}\end{array}$ Note that we apply normal approximation for the log transformed risk ratio and odds ratio because the log transformation typically can improve the performance of normal approximation (Haldane, 1956; Woolf, 1955).

For robust inference, we propose the following variance estimator for $f\left(\hat{\mathit{\theta }}\right)$ that is always consistent regardless of model misspecification: (1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) where $\begin{array}{rl}{\hat{v}}_{tt}& ={\pi }_t^{-1}{S}_{rt}^2+2Q_{ytt}-S_{\mu t}^2,t=1,\dots ,k,\\ {\hat{v}}_{ts}& =Q_{yts}+Q_{yst}-Q_{\mu ts},1\le t<s\le k,\end{array}$ $S_{rt}^2$ is the sample variance of $Y_i-{\hat{\mu }}_t\left({\mathit{X}}_i\right)$ for subjects with $A_i={\mathit{a}}_t$, $Q_{ytt}$ is the sample covariance of $Y_i$ and ${\hat{\mu }}_t\left({\mathit{X}}_i\right)$ for subjects with $A_i={\mathit{a}}_t$, $S_{\mu t}^2$ is the sample variance of ${\hat{\mu }}_t\left({\mathit{X}}_i\right)$ for all subjects, $Q_{yts}$ is the sample covariance of $Y_i$ and ${\hat{\mu }}_s\left({\mathit{X}}_i\right)$ for subjects with $A_i={\mathit{a}}_t$, and $Q_{\mu ts}$ is the sample covariance of ${\hat{\mu }}_t\left({\mathit{X}}_i\right)$ and ${\hat{\mu }}_s\left({\mathit{X}}_i\right)$ for all subjects. These robust variance estimators can be directly calculated using our R package RobinCar that is publicly available at https://github.com/tye27/RobinCar.

To end this section, we describe the variance estimator in Ge et al. (2011) for the g-computation estimator of risk difference ${\hat{\theta }}_2-{\hat{\theta }}_1$ in a two-arm trial and discuss why it can be inconsistent and underestimate the true variance. In our notation, Ge et al. (2011) wrote the g-computation estimator ${\hat{\theta }}_2-{\hat{\theta }}_1$ as $g_n\left(\widehat{\text{β}}\right)$, where $g_n\left(\widehat{\text{β}}\right)=\frac{1}{n}\sum _{i=1}^n\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}({\widehat{\text{β}}}_A^{\mathrm{\top }}{\mathit{a}}_2+{\widehat{\text{β}}}_X^{\mathrm{\top }}{\mathit{X}}_i)-\frac{1}{n}\sum _{i=1}^n\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}({\widehat{\text{β}}}_A^{\mathrm{\top }}{\mathit{a}}_1+{\widehat{\text{β}}}_X^{\mathrm{\top }}{\mathit{X}}_i)$ and $\widehat{\text{β}}=({\widehat{\text{β}}}_A^{\mathrm{\top }},{\widehat{\text{β}}}_X^{\mathrm{\top }}{)}^{\mathrm{\top }}$. Then they applied the Taylor expansion $g_n\left(\widehat{\text{β}}\right)-g_n\left(\text{β}\right)=\left\{\mathrm{\nabla }{g}_n\right(\text{β}\left){\}}^{\mathrm{\top }}\right(\widehat{\text{β}}-\text{β})+o_p(n^{-1/2}),$ where $\text{β}$ is the probability limit of $\widehat{\text{β}}$, and proposed $n^{-1}\left\{\mathrm{\nabla }{g}_n\right(\widehat{\text{β}}\left){\}}^{\mathrm{\top }}{\hat{\mathit{V}}}_{\mathrm{M}}\right\{\mathrm{\nabla }{g}_n\left(\widehat{\text{β}}\right)\}$ as a variance estimator for ${\hat{\theta }}_2-{\hat{\theta }}_1$, where ${\hat{\mathit{V}}}_{\mathrm{M}}$ is the model-based variance estimator for $\sqrt{n}(\widehat{\text{β}}-\text{β})$ from the standard maximum likelihood approach. This approach has two problems. First, it uses the model-based variance estimator ${\hat{\mathit{V}}}_{\mathrm{M}}$, which may be inconsistent to the true variance of $\widehat{\text{β}}$ under model misspecification. Second, from $({\hat{\theta }}_2-{\hat{\theta }}_1)-({\theta }_2-{\theta }_1)=\left\{g_n\right(\widehat{\text{β}})-g_n(\text{β}\left)\right\}+\left\{g_n\right(\text{β})-({\theta }_2-{\theta }_1\left)\right\},$ the variance estimator proposed by Ge et al. (2011) only accounts for the variance of the first term $g_n\left(\widehat{\text{β}}\right)-g_n\left(\text{β}\right)$ but misses the variability from $g_n\left(\text{β}\right)-({\theta }_2-{\theta }_1)$ that is not 0 as the function $\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}(\cdot )$ is nonlinear. This second problem can lead to a confidence interval with too low coverage probability, which can be seen from the simulation results in the following section.

3. Simulations

We conduct simulations to evaluate the finite-sample performance of our robust variance estimator in (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) ). We consider two arms or three arms, simple randomization for treatment assignments with equal allocation (i.e., ${\pi }_1={\pi }_2=1/2$ for two arms and ${\pi }_1={\pi }_2={\pi }_3=1/3$ for three arms), a one-dimensional covariate $X\sim N(0,3^2)$, and n = 200 or 500.

We consider the following three outcome data generating processes.

• Case I: $P(Y=1\mid \mathit{A},X)=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}\{-2+5I(\mathit{A}={\mathit{a}}_2)+X\}$.

• Case II: $P(Y=1\mid \mathit{A}={\mathit{a}}_1,X)=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}(-2+X)$ and $P(Y=1\mid \mathit{A}={\mathit{a}}_2,X)=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}(3+1.5X-0.01X^2)$.

• Case III: $P(Y=1\mid \mathit{A},X)=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}(-2+2I(\mathit{A}={\mathit{a}}_2)+4I(\mathit{A}={\mathit{a}}_3)+X)$.

In order to determine the true values of the unconditional response means, we simulate a large dataset of sample size ${10}^7$ for each case and obtain that $({\theta }_1,{\theta }_2)=(0.2830,0.8057)$ for Case I, $({\theta }_1,{\theta }_2)=(0.2830,0.7297)$ for Case II, and $({\theta }_1,{\theta }_2,{\theta }_3)=(0.2827,0.5004,0.7172)$ for Case III. In each case, the g-computation estimator is based on fitting a working logistic model $P(Y=1\mid \mathit{A},X)=\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{i}\mathrm{t}({\text{β}}_A^{\mathrm{\top }}\mathit{A}+{\text{β}}_{X}X)$, which is correctly specified under Case I and Case III, but is misspecified under Case II.

For Cases I–II, which have two arms, we focus on estimating ${\theta }_2-{\theta }_1$ and also include the variance estimator in Ge et al. (2011). For Case III, which has three arms, we evaluate our robust variance estimators for three common unconditional treatment effects for binary outcomes. The results for Cases I–II are in Table 1 and for Case III are in Table 2, which include (i) the true parameter value, (ii) Monte Carlo mean and standard deviation (SD) of g-computation point estimators, (iii) average of standard error (SE), and (iv) coverage probability (CP) of 95% confidence intervals. We use sample size $n=$ 200 or 500, and 10000 simulation runs.

Table 1. Simulation mean and standard deviation (SD) of ${\hat{\theta }}_2-{\hat{\theta }}_1$, average standard error (SE), and coverage probability (CP) of 95% asymptotic confidence interval for ${\theta }_2-{\theta }_1$ under Cases I–II and simple randomization.

			${\hat{\theta }}_2-{\hat{\theta }}_1$		Robust SE in (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) )		SE in Ge et al. (2011)
Case	${\theta }_2-{\theta }_1$	n	Mean	SD	SE	CP (%)	SE	CP (%)
I	0.5227	200	0.5228	0.0464	0.0464	94.44	0.0415	91.27
		500	0.5227	0.0295	0.0294	94.70	0.0264	91.94
II	0.4467	200	0.4469	0.0457	0.0458	94.56	0.0404	91.08
		500	0.4463	0.0289	0.0290	94.90	0.0257	91.77

Table 2. Simulation mean and standard deviation (SD) of g-computation estimators, average standard error (SE), and coverage probability (CP) of 95% asymptotic confidence interval based on robust SE (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) ) under Case III and simple randomization.

		n = 200				n = 500
Parameter	Truth	Mean	SD	SE	CP	Mean	SD	SE	CP
${\theta }_2-{\theta }_1$	0.2177	0.2176	0.0578	0.0573	94.34	0.2170	0.0366	0.0363	94.82
$\log \left({\theta }_2/{\theta }_1\right)$	0.5711	0.5798	0.1701	0.1664	94.50	0.5726	0.1053	0.1042	94.59
$\log \frac{{\theta }_2/(1-{\theta }_2)}{{\theta }_1/(1-{\theta }_1)}$	0.9328	0.9440	0.2620	0.2586	94.63	0.9341	0.1637	0.1624	94.79
${\theta }_3-{\theta }_1$	0.4346	0.4348	0.0581	0.0568	94.15	0.4347	0.0360	0.0360	94.84
$\log \left({\theta }_3/{\theta }_1\right)$	0.9311	0.9432	0.1653	0.1611	94.43	0.9353	0.1018	0.1009	94.92
$\log \frac{{\theta }_3/(1-{\theta }_3)}{{\theta }_1/(1-{\theta }_1)}$	1.8621	1.8852	0.2920	0.2851	94.57	1.8712	0.1791	0.1788	95.01

From Tables 1–2, we see that the g-computation estimators have negligible biases compared to the standard deviations. Our robust standard error, which is the squared root of variance estimator in (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) ), is always very close to the actual standard deviation, and the related confidence interval has nominal coverage across all settings. In contrast, the standard error in Ge et al. (2011) underestimates the actual standard deviation under Case I when there is no model misspecification, as well as under Case II when there is model misspecification, and the related confidence intervals have too low coverage probabilities in both cases.

4. Summary and discussion

In this article, we provide an explicit robust variance estimator formula for g-computation estimators, which can be used for different unconditional treatment effects of interest and clinical trials with two or more arms. Our simulations demonstrate that the variance estimator can be reliably used in practice.

In this article, for the purpose of being specific, we focus on the logistic model that regresses the outcome on the treatment indicators and covariates, which is arguably the most widely used model for binary outcomes. However, our robust variance estimation formula in (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) ) is not limited to this model and can be used with different specifications of the working model (e.g., fitting a separate logistic model for each treatment arm) or with other generalized linear models using a canonical link for non-binary outcomes (e.g., Poisson regression for count outcomes). Additionally, although our article considers simple randomization, our robust variance formula in (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) ) can also be used for a complete randomization scheme where the sample size in every group t is fixed to be $n{\pi }_t$, because this randomization scheme leads to the same asymptotic distribution as the simple randomization (Ye et al., 2023). Simulation results under this randomization scheme are similar to those under simple randomization; see Tables 3-4 in the Appendix.

We implement an R package called RobinCar to conveniently compute the g-computation estimator and our robust variance estimators, which is publicly available at https://github.com/tye27/RobinCar.

Appendix

In Tables 3–4, we include simulation results under a complete randomization scheme where the sample size in every group t is fixed to be $n{\pi }_t$.

Table 3. Simulation mean and standard deviation (SD) of ${\hat{\theta }}_2-{\hat{\theta }}_1$, average standard error (SE) and coverage probability (CP) of 95% asymptotic confidence interval for ${\theta }_2-{\theta }_1$ under Cases I–II and complete randomization that fixes $n_t=n{\pi }_t$.

			${\hat{\theta }}_2-{\hat{\theta }}_1$		Robust SE in (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) )		SE in Ge et al. (2011)
Case	${\theta }_2-{\theta }_1$	n	Mean	SD	SE	CP (%)	SE	CP (%)
I	0.5227	200	0.5231	0.0464	0.0462	94.62	0.0414	91.54
	500	0.5230	0.0298	0.0294	94.67	0.0264	91.70
II	0.4467	200	0.4469	0.0457	0.0456	94.68	0.0403	91.29
	500	0.4471	0.0290	0.0290	94.74	0.0257	91.45

Table 4. Simulation mean and standard deviation (SD) of g-computation estimators, average standard error (SE), and coverage probability (CP) of 95% asymptotic confidence interval based on robust SE (1(1) $\begin{array}{r}{n}^{-1}\left\{\mathrm{\nabla }f\right(\hat{\mathit{\theta }}\left){\}}^{\mathrm{\top }}\hat{\mathit{V}}\right\{\mathrm{\nabla }f\left(\hat{\mathit{\theta }}\right)\},\hat{\mathit{V}}=\left(\begin{array}{cccc}{\hat{v}}_{11}& {\hat{v}}_{12}& \dots & {\hat{v}}_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{v}}_{1k}& {\hat{v}}_{2k}& \dots & {\hat{v}}_{kk}\end{array}\right),\end{array}$(1) ) under Case III and complete randomization that fixes $n_t=n{\pi }_t$.

		n = 200				n = 500
Parameter	Truth	Mean	SD	SE	CP	Mean	SD	SE	CP
${\theta }_2-{\theta }_1$	0.2177	0.2177	0.0580	0.0570	94.25	0.2182	0.0364	0.0362	94.61
$\log \left({\theta }_2/{\theta }_1\right)$	0.5711	0.5790	0.1687	0.1652	94.60	0.5755	0.1046	0.1040	95.02
$\log \frac{{\theta }_2/(1-{\theta }_2)}{{\theta }_1/(1-{\theta }_1)}$	0.9328	0.9437	0.2616	0.2567	94.78	0.9391	0.1627	0.1621	94.73
${\theta }_3-{\theta }_1$	0.4346	0.4349	0.0579	0.0567	93.96	0.4354	0.0364	0.0360	94.23
$\log \left({\theta }_3/{\theta }_1\right)$	0.9311	0.9424	0.1640	0.1602	94.02	0.9371	0.1017	0.1008	94.99
$\log \frac{{\theta }_3/(1-{\theta }_3)}{{\theta }_1/(1-{\theta }_1)}$	1.8621	1.8853	0.2910	0.2844	94.31	1.8751	0.1811	0.1788	94.48

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Institute of Allergy and Infectious Diseases [NIAID 5 UM1 AI068617].