Fit the core outcome-model and MSM components for two-mixture g-computation
Source:R/qgcompmulti_msm_fit.R
qgcompmulti_msm_fit.RdFits the outcome regression, computes predicted potential outcomes under joint interventions on two exposure mixtures, and estimates the marginal structural model (MSM) coefficients that summarize the resulting two-dimensional intervention-response surface.
Arguments
- f
A model formula for the outcome regression. The formula should include the outcome and any baseline covariates. Mixture variables listed in
mix1andmix2should also appear in the formula if they are to be included in the outcome model.- data
A data frame containing the outcome, exposure variables, and any covariates in the model.
- mix1
A character vector giving the names of the variables in the first exposure mixture.
- mix2
A character vector giving the names of the variables in the second exposure mixture.
- interaction
Logical; if
TRUE, the package includes an interaction term in both the outcome regression and the fitted MSM. In the current implementation, the outcome model is augmented with a cross-product between the sums of the components inmix1andmix2, and the MSM includes thepsi1 * psi2interaction term. IfFALSE, both models are fit without that interaction.- family
A GLM family object (e.g.,
gaussian(),binomial(),poisson()) specifying the outcome model.- estimand_scale
Optional character string naming the fitted MSM estimand scale. Supported values depend on
familyand its link. IfNULL, the Version0.5.0defaults are used.- q
Integer greater than or equal to 2 giving the number of quantiles used to discretize the exposure variables, or
NULLto skip quantization and fit the outcome model on the original exposure scale. Whenq = NULL, the fit-time intervention grid is defined by the pooled 25th, 50th, and 75th percentile values within each mixture, and under each intervention every component in a mixture is set to the same pooled mixture-specific value.- centering
Character string controlling how the marginal structural model intervention variables are coded when
q = NULL. Must be one of"none"or"median". Centering affects only the MSM predictors and does not change the outcome regression fit. This argument is ignored whenqis numeric.- id
Optional character string giving the name of a cluster identifier variable. If supplied, Monte Carlo subsampling is performed at the cluster level rather than the observation level.
- MCsize
Optional integer controlling the Monte Carlo sample size used to approximate the marginalization step in g-computation. If
MCsizeis greater than or equal to the current analysis sample size, all observations are used. Smaller values compute predicted outcomes over a random subsample drawn from the empirical distribution. Whenidis supplied andMCsize < nrow(data), the approximation is implemented by samplingMCsizeclusters with replacement.- seed
Optional integer random seed used to make Monte Carlo subsampling reproducible when
MCsize < nrow(data). IfNULL, the current RNG state is used and not modified byqgcompmulti_msm_fit().
Value
A list with components:
outcome_fitThe fitted outcome regression model object.
msm_fitThe fitted marginal structural model object.
coefficientsA named vector of MSM coefficients.
n_usedThe number of observations used in the g-computation prediction step.
intervention_gridThe fit-time intervention grid on the intervention-value scale.
msm_gridThe corresponding fit-time grid on the MSM coding scale.
counterfactual_surfaceThe exact fit-time counterfactual mean surface implied by the fitted outcome model.
msm_surfaceThe fitted MSM surface evaluated on the common fit-time grid.
surface_comparisonA direct exact-versus-MSM comparison object on the common fit-time grid.
counterfactual_surface_targetThe transformed fit-time MSM target surface used in the MSM fit.
msm_surface_targetThe fitted MSM target-scale surface evaluated on the common fit-time grid.
surface_comparison_targetA direct exact-versus-MSM comparison object on the target scale.
Details
This is a lower-level fitting helper used internally by
qgcomp.glm.multi(). It is exported because it can be useful for method
development, testing, and direct inspection of the fitted outcome model, MSM,
and stored fit-time surfaces. Most users will want qgcomp.glm.multi()
instead.
This function carries out the core g-computation step for the two-mixture
extension. It first fits a generalized linear model for the outcome using
glm(). It then constructs a grid of intervention levels over the two
mixtures, replaces the observed mixture values with the intervention values,
and computes predicted outcomes for each observation (or for a Monte Carlo
subsample if MCsize < nrow(data)). These predicted outcomes are stacked
into a pseudo-dataset and used to fit a marginal structural model that
summarizes the dose-response surface. Note that predicted potential
outcomes are computed on the response scale. The marginal structural model
is then fit either on that response-scale surface or on a transformed target
surface implied by estimand_scale.
When q is an integer, the intervention grid is 0, 1, ..., q - 1 for each
mixture, corresponding to simultaneous quantile increases in all components
of that mixture.
When q = NULL, each mixture is instead set to common pooled percentile
values on the original exposure scale. The resulting MSM coefficients are
therefore scale-dependent and should be interpreted in the units of the
underlying exposures. If centering = "median", the intercept corresponds
to the pooled-median intervention for both mixtures.
Because the number of intervention combinations grows as q^2, this step
can become computationally expensive for large datasets. The MCsize
argument reduces this burden by approximating the empirical covariate
distribution using a random subset while leaving the outcome model fit
unchanged.
Examples
dat <- sim_mixture_data(
n = 500,
pA = 3,
pB = 3,
rho_within_A = 0.3,
rho_within_B = 0.3,
rho_between = 0.2,
psi1 = 0.5,
psi2 = 0.3,
psi12 = 0.2,
return_quantized = FALSE,
seed = 123
)
qgcompmulti_msm_fit(
f = Y ~ X1 + X2 + X3 + W1 + W2 + W3 + C,
data = dat,
mix1 = c("X1", "X2", "X3"),
mix2 = c("W1", "W2", "W3"),
interaction = TRUE,
q = 4,
MCsize = nrow(dat),
seed = 13
)
#> $outcome_fit
#>
#> Call: glm(formula = outcome_formula, family = family, data = data)
#>
#> Coefficients:
#> (Intercept) X1
#> 1.63670 0.26724
#> X2 X3
#> 0.34837 0.28444
#> W1 W2
#> 0.22404 0.26326
#> W3 C
#> 0.11939 1.01466
#> I((X1 + X2 + X3) * (W1 + W2 + W3))
#> 0.01217
#>
#> Degrees of Freedom: 499 Total (i.e. Null); 491 Residual
#> Null Deviance: 1463
#> Residual Deviance: 513.6 AIC: 1452
#>
#> $msm_fit
#>
#> Call: glm(formula = msm_formula, data = msmdat)
#>
#> Coefficients:
#> (Intercept) psi1 psi2 psi1:psi2
#> 1.6138 0.9000 0.6067 0.1095
#>
#> Degrees of Freedom: 7999 Total (i.e. Null); 7996 Residual
#> Null Deviance: 17420
#> Residual Deviance: 1.569e-24 AIC: -487700
#>
#> $coefficients
#> (Intercept) psi1 psi2 psi1:psi2
#> 1.6137564 0.9000449 0.6066983 0.1095483
#>
#> $n_used
#> [1] 500
#>
#> $intervention_grid
#> grid_id psi1 psi2
#> 1 1 0 0
#> 2 2 1 0
#> 3 3 2 0
#> 4 4 3 0
#> 5 5 0 1
#> 6 6 1 1
#> 7 7 2 1
#> 8 8 3 1
#> 9 9 0 2
#> 10 10 1 2
#> 11 11 2 2
#> 12 12 3 2
#> 13 13 0 3
#> 14 14 1 3
#> 15 15 2 3
#> 16 16 3 3
#>
#> $msm_grid
#> grid_id psi1 psi2
#> 1 1 0 0
#> 2 2 1 0
#> 3 3 2 0
#> 4 4 3 0
#> 5 5 0 1
#> 6 6 1 1
#> 7 7 2 1
#> 8 8 3 1
#> 9 9 0 2
#> 10 10 1 2
#> 11 11 2 2
#> 12 12 3 2
#> 13 13 0 3
#> 14 14 1 3
#> 15 15 2 3
#> 16 16 3 3
#>
#> $counterfactual_surface
#> grid_id intervention_psi1 intervention_psi2 msm_psi1 msm_psi2 exact_mean
#> 1 1 0 0 0 0 1.613756
#> 2 2 1 0 1 0 2.513801
#> 3 3 2 0 2 0 3.413846
#> 4 4 3 0 3 0 4.313891
#> 5 5 0 1 0 1 2.220455
#> 6 6 1 1 1 1 3.230048
#> 7 7 2 1 2 1 4.239641
#> 8 8 3 1 3 1 5.249235
#> 9 9 0 2 0 2 2.827153
#> 10 10 1 2 1 2 3.946295
#> 11 11 2 2 2 2 5.065436
#> 12 12 3 2 3 2 6.184578
#> 13 13 0 3 0 3 3.433851
#> 14 14 1 3 1 3 4.662541
#> 15 15 2 3 2 3 5.891231
#> 16 16 3 3 3 3 7.119921
#>
#> $msm_surface
#> grid_id intervention_psi1 intervention_psi2 msm_psi1 msm_psi2 msm_mean
#> 1 1 0 0 0 0 1.613756
#> 2 2 1 0 1 0 2.513801
#> 3 3 2 0 2 0 3.413846
#> 4 4 3 0 3 0 4.313891
#> 5 5 0 1 0 1 2.220455
#> 6 6 1 1 1 1 3.230048
#> 7 7 2 1 2 1 4.239641
#> 8 8 3 1 3 1 5.249235
#> 9 9 0 2 0 2 2.827153
#> 10 10 1 2 1 2 3.946295
#> 11 11 2 2 2 2 5.065436
#> 12 12 3 2 3 2 6.184578
#> 13 13 0 3 0 3 3.433851
#> 14 14 1 3 1 3 4.662541
#> 15 15 2 3 2 3 5.891231
#> 16 16 3 3 3 3 7.119921
#>
#> $surface_comparison
#> grid_id intervention_psi1 intervention_psi2 msm_psi1 msm_psi2 exact_mean
#> 1 1 0 0 0 0 1.613756
#> 2 2 1 0 1 0 2.513801
#> 3 3 2 0 2 0 3.413846
#> 4 4 3 0 3 0 4.313891
#> 5 5 0 1 0 1 2.220455
#> 6 6 1 1 1 1 3.230048
#> 7 7 2 1 2 1 4.239641
#> 8 8 3 1 3 1 5.249235
#> 9 9 0 2 0 2 2.827153
#> 10 10 1 2 1 2 3.946295
#> 11 11 2 2 2 2 5.065436
#> 12 12 3 2 3 2 6.184578
#> 13 13 0 3 0 3 3.433851
#> 14 14 1 3 1 3 4.662541
#> 15 15 2 3 2 3 5.891231
#> 16 16 3 3 3 3 7.119921
#> msm_mean residual
#> 1 1.613756 8.437695e-15
#> 2 2.513801 1.021405e-14
#> 3 3.413846 1.199041e-14
#> 4 4.313891 1.509903e-14
#> 5 2.220455 9.325873e-15
#> 6 3.230048 1.110223e-14
#> 7 4.239641 1.332268e-14
#> 8 5.249235 1.687539e-14
#> 9 2.827153 1.021405e-14
#> 10 3.946295 1.287859e-14
#> 11 5.065436 1.509903e-14
#> 12 6.184578 1.953993e-14
#> 13 3.433851 1.154632e-14
#> 14 4.662541 1.421085e-14
#> 15 5.891231 1.687539e-14
#> 16 7.119921 2.042810e-14
#>
#> $counterfactual_surface_target
#> grid_id intervention_psi1 intervention_psi2 msm_psi1 msm_psi2 exact_target
#> 1 1 0 0 0 0 1.613756
#> 2 2 1 0 1 0 2.513801
#> 3 3 2 0 2 0 3.413846
#> 4 4 3 0 3 0 4.313891
#> 5 5 0 1 0 1 2.220455
#> 6 6 1 1 1 1 3.230048
#> 7 7 2 1 2 1 4.239641
#> 8 8 3 1 3 1 5.249235
#> 9 9 0 2 0 2 2.827153
#> 10 10 1 2 1 2 3.946295
#> 11 11 2 2 2 2 5.065436
#> 12 12 3 2 3 2 6.184578
#> 13 13 0 3 0 3 3.433851
#> 14 14 1 3 1 3 4.662541
#> 15 15 2 3 2 3 5.891231
#> 16 16 3 3 3 3 7.119921
#>
#> $msm_surface_target
#> grid_id intervention_psi1 intervention_psi2 msm_psi1 msm_psi2 msm_target
#> 1 1 0 0 0 0 1.613756
#> 2 2 1 0 1 0 2.513801
#> 3 3 2 0 2 0 3.413846
#> 4 4 3 0 3 0 4.313891
#> 5 5 0 1 0 1 2.220455
#> 6 6 1 1 1 1 3.230048
#> 7 7 2 1 2 1 4.239641
#> 8 8 3 1 3 1 5.249235
#> 9 9 0 2 0 2 2.827153
#> 10 10 1 2 1 2 3.946295
#> 11 11 2 2 2 2 5.065436
#> 12 12 3 2 3 2 6.184578
#> 13 13 0 3 0 3 3.433851
#> 14 14 1 3 1 3 4.662541
#> 15 15 2 3 2 3 5.891231
#> 16 16 3 3 3 3 7.119921
#>
#> $surface_comparison_target
#> grid_id intervention_psi1 intervention_psi2 msm_psi1 msm_psi2 exact_target
#> 1 1 0 0 0 0 1.613756
#> 2 2 1 0 1 0 2.513801
#> 3 3 2 0 2 0 3.413846
#> 4 4 3 0 3 0 4.313891
#> 5 5 0 1 0 1 2.220455
#> 6 6 1 1 1 1 3.230048
#> 7 7 2 1 2 1 4.239641
#> 8 8 3 1 3 1 5.249235
#> 9 9 0 2 0 2 2.827153
#> 10 10 1 2 1 2 3.946295
#> 11 11 2 2 2 2 5.065436
#> 12 12 3 2 3 2 6.184578
#> 13 13 0 3 0 3 3.433851
#> 14 14 1 3 1 3 4.662541
#> 15 15 2 3 2 3 5.891231
#> 16 16 3 3 3 3 7.119921
#> msm_target residual_target
#> 1 1.613756 8.437695e-15
#> 2 2.513801 1.021405e-14
#> 3 3.413846 1.199041e-14
#> 4 4.313891 1.509903e-14
#> 5 2.220455 9.325873e-15
#> 6 3.230048 1.110223e-14
#> 7 4.239641 1.332268e-14
#> 8 5.249235 1.687539e-14
#> 9 2.827153 1.021405e-14
#> 10 3.946295 1.287859e-14
#> 11 5.065436 1.509903e-14
#> 12 6.184578 1.953993e-14
#> 13 3.433851 1.154632e-14
#> 14 4.662541 1.421085e-14
#> 15 5.891231 1.687539e-14
#> 16 7.119921 2.042810e-14
#>