Syntax

The primary fitting function, MRStdCRT_fit(), returns an object with estimates and inference for both cATE and iATE. The returned object also contains:

• The ICS test statistic and p-value in fit$ics,
• The number of clusters in fit$m,
• The observed cluster sizes in fit$N,
• The outcome-model method in fit$model and fitting option in fit$model_weights.

summary(fit) prints the estimates, inference, model details, number of clusters, and ICS test. Use fit$N to inspect the individual cluster sizes. Call MRStdCRT_fit() as follows:

MRStdCRT_fit(
  formula, data, cluster, trt, trtprob = rep(0.5, nrow(data)),
  method, family = gaussian(link = "identity"), corstr, scale,
  alpha = 0.05,
  model_weights = c("none", "inverse_cluster_size", "estimand_specific"),
  inference_scale = c("original", "log")
)

with the following arguments:

• formula: The outcome regression formula, with a single untransformed response column and an intercept. Supported terms are untransformed covariate main effects and treatment-by-covariate interactions, as described below.
• data: A data frame with categorical covariates converted to dummy variables. Cluster, treatment, outcome, and covariates must occupy distinct columns.
• cluster: The column name of the cluster identifier. Cluster IDs must not be missing.
• trt: The column name of the treatment variable. Values must be numeric or logical 0/1, without missing values, and constant within each cluster.
• trtprob: A numeric row-level or cluster-level vector of probabilities of assignment to treatment, \(P(A_i=1)\). Values must be finite, strictly between 0 and 1, and constant within each cluster. If NULL, a common probability is estimated as the proportion of clusters assigned to treatment, giving each cluster one vote. This probability is re-estimated in each jackknife sample. Use known probabilities from the randomization design when available.
• method: The method used for outcome regression model fitting, i.e. “GLM”, “LMM”, “GEE”, “GLMM”.
• family: A family object (default gaussian(link = "identity")). LMM uses a Gaussian identity working model. GLMM supports Gaussian identity, binomial logit, and Poisson log. Unweighted GLM and GEE pass the family to their fitting routines. The Gaussian-logit family is also an unweighted cluster-summary GLM working-model option for binary outcomes. Supported weighted combinations are described below.
• corstr: The working correlation structure, required for GEE and optional for other methods. Weighted GEE supports "independence" and "exchangeable".
• scale: The scale of estimand including risk difference (RD), risk ratio (RR), and odds ratio (OR).
• alpha: The significance level (default is 0.05).
• model_weights: The outcome-model fitting option. "none" (default) uses one unweighted model for both cATE and iATE. "inverse_cluster_size" uses one inverse-cluster-size weighted model for both estimands. "estimand_specific" uses a weighted model for cATE and an unweighted model for iATE, with both models refitted after each cluster deletion for paired inference and the ICS test. Both nondefault options support GEE, LMM, and GLMM; GLM supports only "none". See below for supported families, weight definitions, and executable examples.
• inference_scale: "original" (default) retains inference on the effect scale. For scale = "RR" or scale = "OR", "log" instead uses the log-transformed jackknife estimates for confidence intervals and p-values. It does not change the fitted model, point estimates, original-scale standard errors, or paired ICS test. It is not available for scale = "RD".

The treatment main effect is added automatically. Each interaction must involve treatment and a single covariate, and that covariate must also appear as a main effect. For example, with trt = "A", both y ~ x + A:x and y ~ A * x are supported. Transformations inside the formula, interactions between covariates, higher-order interactions, dot notation (.), no-intercept formulas, and term subtraction are not supported. To include a transformed covariate, first create a separate column and then refer to that column in the formula:

data$x_sq <- data$x^2
formula <- y ~ x + x_sq

summary() displays Estimate, SE, CI, and p-value for both estimands. The displayed confidence level is \(100(1-\alpha)\%\), using the alpha specified in MRStdCRT_fit(). The ics component is a named numeric vector: statistic is the paired ICS t statistic, and p.value is its two-sided p-value. These can be extracted as fit$ics[["statistic"]] and fit$ics[["p.value"]], respectively.

Illustrative Example

In this example, we demonstrate how to use the MRStdCRT_fit function to estimate treatment effects in a CRT using the ppact dataset. The goal is to estimate the cluster-averaged treatment effect (c-ATE) and the individual-averaged treatment effect (i-ATE) using the marginal model fitted by generalized estimating equation (GEE).

library(MRStdCRT)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

data(ppact)

prob <- rep(0.5, nrow(ppact))

cluster_assignments <- ppact %>%
  distinct(CLUST, INTERVENTION)

estimated_prob <- rep(
  mean(cluster_assignments$INTERVENTION == 1),
  nrow(ppact)
)

clusters <- sort(unique(ppact$CLUST))
n_clusters <- length(clusters) 
block_info <- data.frame(
  CLUST = clusters
) %>%
  mutate(
    # Create a 'block' identifier for each cluster
    block = case_when(
      row_number() <= 25 ~ 1,
      row_number() <= 66 ~ 2,
      TRUE ~ 3
    )
  ) %>%
  mutate(
    # Assign the desired hypothetical treatment probability to each block
    hypothetical_prob = case_when(
      block == 1 ~ 0.3,
      block == 2 ~ 0.5,
      block == 3 ~ 0.6
    )
  )
prob_lookup <- setNames(block_info$hypothetical_prob, block_info$CLUST)
probs_vector <- prob_lookup[as.character(ppact$CLUST)]

example <- MRStdCRT_fit(
  formula = PEGS ~ AGE + FEMALE + comorbid + Dep_OR_Anx + pain_count + PEGS_bl +
    BL_benzo_flag + BL_avg_daily + satisfied_primary + n,
  data = ppact,
  cluster = "CLUST",
  trt = "INTERVENTION",
  trtprob = prob,
  method = "GEE",
  corstr = "independence",
  scale = "RR"
)
## To view the summary, use the following command
summary(example)

## 
## Model-robust Standardization
## =========================================
##   Method   : GEE
##   Family   : gaussian (link = identity)
##   Clusters : 106
##   Scale    : Risk ratio
## 
## Estimates:
##       Estimate Std. Error         95% CI p-value
## c-ATE    0.907      0.028 (0.852, 0.962) 0.001**
## i-ATE    0.926      0.024 (0.879, 0.973) 0.002**
## 
## Test for no informative cluster size:
##   Statistic: -1.7204
##   p-value  : 0.0883

## One may also extract specific statistical feature, such as table for point and interval estimates, and standard errors
example$estimate

##       Estimate Std. Error  CI lower  CI upper     p-value
## cATE 0.9070587 0.02760612 0.8523209 0.9617965 0.001063919
## iATE 0.9262683 0.02369230 0.8792909 0.9732458 0.002393551

example$ics

##   statistic     p.value 
## -1.72035885  0.08831306

Fitting weights and supported models

For a cluster with \(N_i\) observed individuals, weighted GEE passes an individual fitting weight of \(1/N_i\) to geepack::geeglm(). Supported family/link combinations are Gaussian identity, binomial logit, and Poisson log, with independence or exchangeable working correlation.

Weighted LMM and GLMM use WeMix::mix() with a single cluster random intercept. The conditional weight for each individual is 1, and the cluster weight is \(1/N_i\). The weighted criterion is

\[ \ell_w(\theta)=\sum_i\frac{1}{N_i}\log L_i(\theta),\qquad L_i(\theta)=\int\left\{\prod_j f(Y_{ij}\mid b_i;\theta)\right\} \phi(b_i;0,\tau^2)\,db_i. \]

Thus the cluster’s likelihood is integrated over its random intercept before weighting its log contribution. This is different from raising each conditional individual likelihood contribution to \(1/N_i\) before integration. MRStdCRT constructs the two weight columns internally and uses the conditional-weight interface of WeMix.

Weighted LMM uses a Gaussian identity working model. This working model can also be applied to numeric binary or count outcomes. Weighted GLMM supports Gaussian identity, binomial logit, and Poisson log; the Gaussian case uses the same weighted linear-model fit as LMM. Weighted binomial models require outcomes coded 0/1, and Poisson models require nonnegative integer outcomes. The cluster-summary GLM supports only model_weights = "none".

Cluster sizes are calculated from the current analysis data. A column named N, Freq, or W1 is not required. No cluster-size main effect or treatment-by-cluster-size interaction is added to the outcome formula. The package retains its existing separation of individual covariates into within-cluster deviations and between-cluster means, and its handling of supported treatment interactions. Additional centering in WeMix is disabled. The weighted mixed-model path uses max_iteration = 100 and nQuad = 13 consistently for the full-data and jackknife fits. Gaussian models use the linear-model interface (family = NULL); numerical quadrature is used for nonlinear models. These fitting calculations are separate from the marginal-prediction approximations below. The weighted mixed-model criterion uses maximum likelihood. The default unweighted LMM retains its original REML fit. The "estimand_specific" option uses these same fitting procedures, so its LMM analysis combines a weighted maximum-likelihood fit for cATE with an unweighted REML fit for iATE.

For nonlinear mixed models, a random-intercept variance near zero can cause a numerical failure in WeMix. Fits that cannot be verified stop the analysis; the package does not remove the random intercept or impose a positive variance floor. A valid Gaussian variance boundary is allowed.

For the weighted path, necessary outcome and covariate data must be complete and finite. Backend failures, invalid parameter estimates, or invalid predictions stop the analysis. An error identifies the full-data fit or the deleted cluster in a jackknife refit. The analysis does not replace a failed weighted fit with an unweighted model or omit failed jackknife samples. With "estimand_specific", each component retains its corresponding fitting procedure. The unweighted component retains the existing warnings and nonparametric fallback after a model-fitting error. The selected weighted cATE, unweighted iATE, and their difference must be finite in the full-data analysis and every deletion sample. An error in either component identifies the weighted or unweighted fit and its full-data or deleted-cluster context.

Marginal predictions and inference

For LMM and Gaussian identity GLMM, the counterfactual prediction is \(X_{ij}(a)^\mathsf{T}\widehat\beta\), without fitted random intercepts or BLUPs. For binomial-logit GLMM, the package retains the variance-matching approximation used in the original MRS implementation:

\[ p_{ij}(a)=\operatorname{expit}\left\{ \frac{X_{ij}(a)^\mathsf{T}\widehat\beta} {\sqrt{1+3\widehat\tau^2/\pi^2}}\right\}. \]

Here \(\widehat\tau^2\) is the fitted random-intercept variance, and \(\pi\) is the mathematical constant. For Poisson-log GLMM, \(p_{ij}(a)=\exp\{X_{ij}(a)^\mathsf{T}\widehat\beta+\widehat\tau^2/2\}\). GEE uses response-scale predictions from the fitted marginal model.

The predictions are averaged within each cluster to obtain \(\widehat\eta_i(a)\). With "none" or "inverse_cluster_size", both estimands use the same augmented outcomes:

\[ R_i(a)=\widehat\eta_i(a)+\frac{I(A_i=a)}{\pi_i(a)} \{\overline Y_i-\widehat\eta_i(a)\},\qquad \widehat\mu_C(a)=\frac{1}{m}\sum_i R_i(a),\quad \widehat\mu_I(a)=\frac{\sum_i N_iR_i(a)}{\sum_i N_i}. \]

With "estimand_specific", let \(R_i^w(a)\) and \(R_i^u(a)\) denote the augmented outcomes from the weighted and unweighted models, respectively. The selected arm means are

\[ \widehat\mu_C^w(a)=\frac{1}{m}\sum_i R_i^w(a),\qquad \widehat\mu_I^u(a)=\frac{\sum_i N_iR_i^u(a)}{\sum_i N_i}. \]

The residual correction is always retained. Risk differences, risk ratios, and odds ratios are constructed from these arm-specific marginal means. The package does not average individual or cluster-specific odds ratios.

Every delete-one-cluster sample refits the selected model or models, including their correlation or variance parameters. For "estimand_specific", both fits omit the same cluster. The jackknife variance is centered at the mean of the deleted-cluster estimates:

\[ \widehat V_{\mathrm{jack}}(\widehat\theta) =\frac{m-1}{m}\sum_i \{\widehat\theta_{-i}-\overline\theta_{(-)}\}^2. \]

The ICS standard error uses the paired cATE-minus-iATE differences from the same deletions. For "estimand_specific", these differences are \(D_{-i}=\widehat\theta_{C,-i}^{w}-\widehat\theta_{I,-i}^{u}\). Applying the same jackknife formula to \(D_{-i}\) retains the covariance between the two estimates. The test statistic is the full-data difference \(\widehat\theta_C^w-\widehat\theta_I^u\) divided by this paired standard error. Inference uses \(m-1\) degrees of freedom. By default, confidence intervals and p-values are computed on the original RD/RR/OR scale. The optional log-scale procedure for RR and OR is described below; it leaves the paired ICS test unchanged. Weighted fitting is an additional working-model choice; it does not guarantee greater efficiency than the default.

Ordinary weighted-MRS examples

The following trial has unequal cluster sizes, a cluster-level covariate, and an individual-level covariate. The formulas do not include cluster size or its interaction with treatment. No weight columns are needed.

set.seed(68213)
sizes <- rep(c(18, 25, 33), 4)
m <- length(sizes)
trial <- data.frame(
  id = rep(seq_len(m), sizes),
  A = rep(rep(0:1, length.out = m), sizes)
)
random_intercept <- rep(rep(c(-1.4, -0.7, 0, 0.2, 0.8, 1.5), each = 2),
                        sizes)
trial$H <- rep(rnorm(m, sd = 0.5), sizes)
trial$x <- rep(rnorm(m, sd = 0.4), sizes) + rnorm(nrow(trial), sd = 0.7)
trial$continuous <- 3 + 0.5 * trial$A + 0.3 * trial$x + 0.2 * trial$H +
  random_intercept + rnorm(nrow(trial), sd = 0.8)
trial$binary <- rbinom(
  nrow(trial), 1,
  plogis(-0.2 + 0.5 * trial$A + 0.3 * trial$x + 0.2 * trial$H +
           random_intercept)
)
trial$count <- rpois(
  nrow(trial), exp(0.2 + 0.3 * trial$A + 0.15 * trial$x + 0.1 * trial$H +
                    0.7 * random_intercept)
)

A Gaussian GEE working model with exchangeable correlation uses an individual fitting weight of \(1/N_i\):

weighted_gee <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "inverse_cluster_size"
)
weighted_gee$estimate

##       Estimate Std. Error  CI lower CI upper   p-value
## cATE 0.2913647  0.6599991 -1.161283 1.744013 0.6674305
## iATE 0.3541473  0.7473106 -1.290672 1.998967 0.6448411

To use the weighted cATE together with the unweighted iATE, select "estimand_specific". The following comparison uses the same trial and Gaussian GEE specification as above. The selected rows agree with the corresponding rows from the two separate analyses; the ICS test uses the paired differences from refitting both models.

unweighted_gee <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "none"
)
estimand_specific_gee <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "estimand_specific"
)
stopifnot(
  isTRUE(all.equal(estimand_specific_gee$estimate["cATE", ],
                   weighted_gee$estimate["cATE", ])),
  isTRUE(all.equal(estimand_specific_gee$estimate["iATE", ],
                   unweighted_gee$estimate["iATE", ]))
)
estimand_specific_gee$estimate

##       Estimate Std. Error  CI lower CI upper   p-value
## cATE 0.2913647  0.6599991 -1.161283 1.744013 0.6674305
## iATE 0.3873097  0.7609663 -1.287566 2.062185 0.6208265

estimand_specific_gee$ics

##  statistic    p.value 
## -0.4829920  0.6385709

For \(m\) clusters, this option requires \(2(m+1)\) outcome-model fits: two full-data fits and two fits for each deletion. It provides the combined estimates and their paired ICS test in one call; an efficiency improvement over the other fitting options is not guaranteed.

A weighted LMM uses the same formula and a single cluster random intercept. Its predictions omit the fitted random effects:

weighted_lmm <- MRStdCRT_fit(
  continuous ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "LMM", family = gaussian(), scale = "RD",
  model_weights = "inverse_cluster_size"
)
weighted_lmm$estimate

##       Estimate Std. Error  CI lower CI upper   p-value
## cATE 0.2914823  0.6599850 -1.161135 1.744100 0.6672988
## iATE 0.3542419  0.7472853 -1.290522 1.999006 0.6447424

For a binary outcome, weighted GLMM combines the weighted fit with the marginal logit approximation described above:

weighted_glmm <- MRStdCRT_fit(
  binary ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GLMM", family = binomial(), scale = "OR",
  model_weights = "inverse_cluster_size"
)
summary(weighted_glmm)

## 
## Model-robust Standardization
## =========================================
##   Method   : GLMM
##   Family   : binomial (link = logit)
##   Clusters : 12
##   Scale    : Odds ratio
##   Fitting  : inverse-cluster-size weighted outcome model
## 
## Estimates:
##       Estimate Std. Error          95% CI p-value
## c-ATE    0.993      0.478 (-0.059, 2.044)   0.988
## i-ATE    1.010      0.541 (-0.181, 2.201)   0.986
## 
## Test for no informative cluster size:
##   Statistic: -0.1140
##   p-value  : 0.9113

For counts, the corresponding call uses count ~ A + x + H, family = poisson(), and, for example, scale = "RR". Setting method = "GLMM", family = gaussian() uses the weighted linear-model backend. For weighted GEE, all three supported families can be combined with either supported working correlation. A user-specified treatment interaction, such as continuous ~ A * x + H, remains available.

Optional log-scale inference for RR and OR

A symmetric confidence interval on the ratio scale can have a negative lower endpoint. With scale = "RR" or scale = "OR", inference_scale = "log" constructs the interval on the log scale and transforms its endpoints back. This option applies to all supported working-model methods and all three outcome-model weighting options. With "estimand_specific", the log procedure uses the selected weighted cATE and unweighted iATE estimates from each paired deletion. The default, inference_scale = "original", retains the original-scale inference.

Let \(\widehat\theta\) be the full-data RR or OR estimate for either cATE or iATE, and write \(z_{-i}=\log(\widehat\theta_{-i})\) for its estimate after deleting cluster \(i\) and refitting the model. The log-scale jackknife standard error is

\[ \widehat{\operatorname{SE}}_{\log} =\left\{\frac{m-1}{m}\sum_i(z_{-i}-\overline z_{(-)})^2\right\}^{1/2}. \]

The reported confidence interval is

\[ \exp\left\{\log(\widehat\theta)\mathbin{\pm} t_{m-1,1-\alpha/2}\widehat{\operatorname{SE}}_{\log}\right\}. \]

The two-sided p-value tests a ratio of 1, using \(\log(\widehat\theta)/\widehat{\operatorname{SE}}_{\log}\) and a \(t_{m-1}\) distribution. The original ratio estimate is retained; it is not replaced by an average of the deleted-cluster estimates. The log standard error is calculated directly from those estimates, rather than by dividing the original-scale standard error by the full-data estimate.

The following example compares the two inference options for a weighted binomial GEE on the trial data above:

ratio_original <- MRStdCRT_fit(
  binary ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = binomial(), corstr = "independence",
  scale = "OR", model_weights = "inverse_cluster_size"
)
ratio_log <- MRStdCRT_fit(
  binary ~ A + x + H, data = trial, cluster = "id", trt = "A",
  method = "GEE", family = binomial(), corstr = "independence",
  scale = "OR", model_weights = "inverse_cluster_size",
  inference_scale = "log"
)
stopifnot(
  identical(ratio_original$estimate[, c("Estimate", "Std. Error")],
            ratio_log$estimate[, c("Estimate", "Std. Error")]),
  identical(ratio_original$ics, ratio_log$ics)
)
ratio_original$estimate

##      Estimate Std. Error    CI lower CI upper   p-value
## cATE 1.015401  0.4875374 -0.05766194 2.088463 0.9753659
## iATE 1.029895  0.5435848 -0.16652675 2.226317 0.9571274

summary(ratio_log)

## 
## Model-robust Standardization
## =========================================
##   Method   : GEE
##   Family   : binomial (link = logit)
##   Clusters : 12
##   Scale    : Odds ratio
##   Fitting  : inverse-cluster-size weighted outcome model
##   Inference: log scale for CI and p-value; SE shown on original scale
## 
## Estimates:
##       Estimate Std. Error         95% CI p-value
## c-ATE    1.015      0.488 (0.364, 2.833)   0.974
## i-ATE    1.030      0.544 (0.327, 3.248)   0.956
## 
## Test for no informative cluster size:
##   Statistic: -0.0946
##   p-value  : 0.9263

ratio_log$log_se

##      cATE      iATE 
## 0.4662255 0.5217847

The Std. Error column of estimate and the SE printed by summary() remain on the original ratio scale. The separate log_se vector contains the log-scale standard errors for cATE and iATE when log inference is selected. summary() identifies this distinction: its confidence intervals and p-values use the log-scale procedure, while its estimates and displayed standard errors remain on the ratio scale. The paired ICS test continues to use the original-scale cATE-minus-iATE difference and its paired jackknife standard error.

Log inference requires valid marginal means for each reported estimand in the full-data fit and every deleted-cluster fit. With "estimand_specific", these are the weighted cATE arm means and unweighted iATE arm means. For OR, both arm means must be strictly between 0 and 1. For RR, both must be positive; for an observed binary outcome, they must also be at most 1. The outcome is treated as binary when all full-data values are 0/1; this classification is retained in every deletion sample. Ratios must be positive and finite. The analysis stops if these conditions fail, if the log-scale standard error is not positive and finite, or if the exponentiated interval cannot be represented with positive, finite endpoints. It does not truncate arm means or ratios, add small constants, or discard jackknife samples. The existing model-fitting procedures are unchanged. In particular, log inference cannot repair invalid marginal means caused by a working model or by the MRS residual correction. Positive confidence limits do not guarantee improved small-sample coverage.

A specific GEE comparison with cluster-size saturated g-computation

Lee, Harhay, and Li (2026) discuss models with a continuous cluster-size main effect and treatment-by-cluster-size interaction. Under suitable GEE conditions, using the same weighted fit and a common randomization probability, the cATE MRS residual correction vanishes and cATE agrees with g-computation. This is an additional property of these particular models; it is not required to use weighted MRS.

Here is a Gaussian identity example with exchangeable working correlation. Cluster size is explicitly included in the formula. The direct geeglm() call is only used to demonstrate the comparison.

set.seed(20261013)
cs_sizes <- rep(c(12, 20, 32, 48), 6)
cs_m <- length(cs_sizes)
assignment <- sample(rep(0:1, each = cs_m / 2))
cs_data <- data.frame(
  id = rep(seq_len(cs_m), cs_sizes), A = rep(assignment, cs_sizes),
  size = rep(cs_sizes, cs_sizes)
)
cs_data$y <- 2 + 0.02 * cs_data$size +
  cs_data$A * (0.4 + 0.01 * cs_data$size) +
  rep(rnorm(cs_m, sd = 0.6), cs_sizes) + rnorm(nrow(cs_data))

cs_mrs <- MRStdCRT_fit(
  y ~ A * size, data = cs_data, cluster = "id", trt = "A",
  method = "GEE", family = gaussian(), corstr = "exchangeable",
  scale = "RD", model_weights = "inverse_cluster_size"
)
cs_data$weight <- 1 / cs_data$size
cs_model <- geepack::geeglm(
  y ~ A * size, data = cs_data, id = id, weights = weight,
  family = gaussian(), corstr = "exchangeable"
)
stopifnot(cs_model$geese$error == 0)
ybar <- tapply(cs_data$y, cs_data$id, mean)
cluster_treatment <- tapply(cs_data$A, cs_data$id, mean)
g_mean <- correction <- numeric(2)
for (a in 0:1) {
  newdata <- cs_data
  newdata$A <- a
  eta <- tapply(predict(cs_model, newdata, type = "response"),
                cs_data$id, mean)
  g_mean[a + 1] <- mean(eta)
  correction[a + 1] <- mean((cluster_treatment == a) / 0.5 * (ybar - eta))
}
stopifnot(max(abs(correction)) < 1e-9,
          abs(cs_mrs$estimate[1, "Estimate"] - diff(g_mean)) < 1e-9)
rbind(g_computation = g_mean, residual_correction = correction)

##                             [,1]          [,2]
## g_computation       2.866507e+00  3.408728e+00
## residual_correction 1.048544e-16 -1.526557e-16

The analogous binomial-logit comparison is restricted here to cluster-level predictors and independence or exchangeable working correlation. These equivalence statements do not extend automatically to nonlinear exchangeable GEE with arbitrary individual covariates, cluster-varying treatment probabilities, or mixed models. The weighted iATE remains the MRS estimate calculated from the same weighted fit; no equality with an unweighted g-computation iATE is asserted.