Effect modification is specified inside dose() with the
modifier argument. Formula modifiers can be binary
numeric/logical variables or factors. Character variables are not
accepted as modifiers; convert them to factors first so the level order
is explicit.
There are two useful ways to code effect modification:
modifier = ~ M uses reference-plus-contrast
coding.modifier = ~ 0 + M or modifier = ~ M - 1
uses subgroup-specific coding.These are reparameterizations of the same model. They give the same fitted likelihood when the same model and data are used, but the coefficients answer different questions directly.
Reference-plus-contrast coding is the default formula-modifier
coding. For a binary modifier M1, the coefficient named
dose is the dose effect in the reference group, and
dose:M1 is the difference between the non-reference group
and the reference group.
fit_contrast <- ameras(
Y.gaussian ~ dose(V1:V10, modifier = ~ M1) + X1 + X2,
data = data,
family = "gaussian",
methods = "RC"
)
coef(fit_contrast)
#> RC
#> (Intercept) -1.3602199
#> X1 0.4788215
#> X2 -0.5192733
#> dose 1.0836911
#> dose:M1 0.1566173
#> sigma 1.1021197This coding is convenient when the effect-modification contrast is the main quantity of interest.
Subgroup-specific coding reports one dose-effect parameter per
subgroup. For the same binary modifier M1, the coefficient
names are dose[M1=0] and dose[M1=1].
fit_subgroup <- ameras(
Y.gaussian ~ dose(V1:V10, modifier = ~ 0 + M1) + X1 + X2,
data = data,
family = "gaussian",
methods = "RC"
)
coef(fit_subgroup)
#> RC
#> (Intercept) -1.3602200
#> X1 0.4788212
#> X2 -0.5192734
#> dose[M1=0] 1.0836912
#> dose[M1=1] 1.2403085
#> sigma 1.1021195The two codings are equivalent. The non-reference subgroup effect is the reference effect plus the contrast:
contrast_coef <- fit_contrast$RC$coefficients
subgroup_coef <- fit_subgroup$RC$coefficients
data.frame(
term = c("reference group", "non-reference group"),
from_contrast = c(
contrast_coef["dose"],
contrast_coef["dose"] + contrast_coef["dose:M1"]
),
from_subgroup = c(
subgroup_coef["dose[M1=0]"],
subgroup_coef["dose[M1=1]"]
),
row.names = NULL
)
#> term from_contrast from_subgroup
#> 1 reference group 1.083691 1.083691
#> 2 non-reference group 1.240308 1.240309Subgroup-specific coding is usually the cleaner choice when subgroup-specific effect estimates and confidence intervals are the main goal. For linear ERR and linear-exponential models, it can also be more numerically robust because default bounds are applied directly to each subgroup dose-effect parameter. With reference-plus-contrast coding, contrast parameters are not bounded by the default transformation.
Factor modifiers use the existing level order of the factor. The first level is the reference level.
data$M_factor <- factor(
ifelse(data$M1 == 0, "unexposed modifier group", "modifier group"),
levels = c("unexposed modifier group", "modifier group")
)With modifier = ~ M_factor, ameras reports contrast
terms for the non-reference levels.
fit_factor_contrast <- ameras(
Y.gaussian ~ dose(V1:V10, modifier = ~ M_factor) + X1 + X2,
data = data,
family = "gaussian",
methods = "RC"
)
coef(fit_factor_contrast)
#> RC
#> (Intercept) -1.3602199
#> X1 0.4788215
#> X2 -0.5192733
#> dose 1.0836911
#> dose:M_factor=modifier group 0.1566173
#> sigma 1.1021197With modifier = ~ 0 + M_factor, ameras reports
subgroup-specific effects for all levels.
fit_factor_subgroup <- ameras(
Y.gaussian ~ dose(V1:V10, modifier = ~ 0 + M_factor) + X1 + X2,
data = data,
family = "gaussian",
methods = "RC"
)
coef(fit_factor_subgroup)
#> RC
#> (Intercept) -1.3602200
#> X1 0.4788212
#> X2 -0.5192734
#> dose[M_factor=unexposed modifier group] 1.0836912
#> dose[M_factor=modifier group] 1.2403085
#> sigma 1.1021195Multi-level factors are supported. With reference-plus-contrast coding, ameras reports one contrast for each non-reference level. With subgroup-specific coding, ameras reports one dose-effect parameter for each factor level.
data$M3 <- factor(
rep(c("low", "middle", "high"), length.out = nrow(data)),
levels = c("low", "middle", "high")
)
fit_three_level <- ameras(
Y.gaussian ~ dose(V1:V10, modifier = ~ 0 + M3) + X1 + X2,
data = data,
family = "gaussian",
methods = "RC"
)
coef(fit_three_level)
#> RC
#> (Intercept) -1.3619937
#> X1 0.4813357
#> X2 -0.5196477
#> dose[M3=low] 1.1433969
#> dose[M3=middle] 1.1715892
#> dose[M3=high] 1.1807593
#> sigma 1.1073548For subgroup-specific confidence intervals, fit the model using
subgroup-specific coding and then call confint() as usual.
For RC, ERC, and MCML, profile
likelihood intervals treat subgroup effects as ordinary coefficient
rows.
For FMA and BMA, sample-based confidence or
credible intervals are computed from the sampled coefficients. When
subgroup-specific coding is used, BMA priors and MCMC sampling remain on
the internal reference-plus-contrast scale, but stored samples,
summaries, diagnostics, trace plots, and sample-based intervals are
reported on the subgroup-specific scale.
Use reference-plus-contrast coding when the contrast itself is the main estimand. Use subgroup-specific coding when subgroup effects are the main estimands, when subgroup-specific profile likelihood intervals are desired, or when ERR/LINEXP optimization is sensitive to the unbounded contrast parameterization.
Testing whether effect modification improves the model is a
model-comparison question. Fit a model without the modifier and a model
with the modifier, then compare the fitted likelihoods using the
appropriate degrees of freedom. dose_lrt() tests
dose-related parameters, but its global dose test is not the same as a
test for effect modification.