modav | R Documentation |
Fits separate (marginal) models for each candidate subgroup, i.e. including the subgroup as a main effect and interaction with treatment for each model. These models are used predict the treatment effect for the subgroup of interest by predicting the effect for all patients in the subgroup and then averaging. These subgroup effects are then calculated for all models and then averaged according to posterior model weights. Details of the procedure are explained in Bornkamp et al. (2017) and Thomas and Bornkamp (2017).
modav(resp, trt, subgr, covars = NULL, data, fitfunc = c("lm", "glm", "glm.nb", "survreg", "coxph", "rlm"), event, exposure, level=0.1, prior = 1, nullprior = 0, ...)
resp |
Character giving the name of the response variable. The variable can be
either defined in the global environment or in the data-set |
trt |
Character giving the name of the treatment variable. The variable can
be either defined in the global environment or in the data-set
|
subgr |
Character vector giving the variable names in |
covars |
Formula, specifying additional (prognostic) covariates to be included
in the models (need to be available in |
data |
Data frame containing the variables referenced in |
fitfunc |
Model fitting functions. Currrently one of |
event |
Character giving the name of the event variable. Has to be specified
when using fit functions |
exposure |
Character giving the name of the exposure variable, needed for
negative binomial regression, when using fit functions
|
level |
Significance level for confidence intervals will be calculated for treatment effect estimates. |
prior |
Numeric vector of prior model/subgroup probabilities of the same
length as |
nullprior |
Numeric giving the prior model probability of the model without any
subgroup effect. This needs to be specified on the same scale as
the prior argument. E.g. if there are 2 subgroups, |
... |
other arguments passed to the model fitting function. |
In the simple linear case (e.g when using fitfunc lm
) for each of the P candidate subgroups the fitted model is of the form
M_p : y_i ~ N(μ_i^(p), σ_p^2), i=1,...,n
where
μ_i^{(p)} = α_p + β_p z_i + (γ_p + δ_p z_i) s_{pi} + ∑_{k = 1}^{K} τ_k x_{ik}
where s_i denotes the subgroup indicators (the column vectors of subgr
), z_i is the treatment indicator (from trt
) and x_{i1}, ..., x_{iK} are additional covariates as specified in covars
.
For other fitting functions the models are of similar form, including prognostic and predictive effects of subgroups.
A treatment effect (on the scale determined by fitfunc
) of the candidate subgroups can be derived naively as \hat{β_p} + \hat{δ_p} and a treatment effect estimate for the complement is given by \hat{β_p}.
Note that choosing subgroups based on these unadjusted treatment effect estimates may lead to overoptimistic conclusions in regards to the treatment effect in that subgroup. Naive estimates do not consider model
selection uncertainty and will often suffer from selection bias.
For each subgroup a treatment effect is obtained by estimating the treatment effect for that subgroup under all models (by averaging the individual predictions in that subgroup) and approximating the resulting estimate within each model by a normal distribution for details see Bornkamp et al, 2017. Posterior model weights are obtained using BIC model weights (Raftery, 1995), so that overall a normal mixture is used to approximate the posterior distribution for every subgroup effect.
The returned treatment effect estimates are based on the median of this distribution, credible bounds are based on posterior quantiles.
Estimates of the interaction (difference in treatment effect between subgroup and complement) are also derived as the median and quantiles of the corresponding mixture distribution.
A list (object of class subtee
). The most important entries are (i) fitmods
containing all
fitted subgroup models and the overall model (ii) trtEff
containing the treatment effect estimates and CI for subgroup and subgroup
complements. (iii) trtEffDiff
containing the differences in
treatment effect estimates (subgroup vs complement) and CI.
Ballarini, N. Thomas, M., Rosenkranz, K. and Bornkamp, B. (2021) "subtee: An R Package for Subgroup Treatment Effect Estimation in Clinical Trials" Journal of Statistical Software, 99, 14, 1-17, doi: 10.18637/jss.v099.i14
Thomas, M., and Bornkamp, B. (2017) "Comparing Approaches to Treatment Effect Estimation for Subgroups in Early Phase Clinical Trials." Statistics in Biopharmaceutical Research, 9, 160-171, doi: 10.1080/19466315.2016.1251490
Bornkamp, B., Ohlssen, D., Magnusson, B. P., and Schmidli, H. (2017) "Model averaging for treatment effect estimation in subgroups." Pharmaceutical Statistics, 16, 133-142, doi: 10.1002/pst.1796
Raftery, A. E. (1995) "Bayesian model selection in social research." Sociological Methodology, 25, 111-163.
summary.subtee
, plot.subtee
,
lm
, glm
, glm.nb
,
survreg
, coxph
## toy example calls using the simulated datnorm data-set without ## treatment and subgroup effect, see ?datnorm for details data(datnorm) head(datnorm) ## first need to create candidate subgroups (if not already defined in data-set) ## here generate candidate subgroups manually (need to be numeric 0-1 variables) groups <- data.frame(labvalL.5=as.numeric(datnorm$labvalue < 0.5), regUS=as.numeric(datnorm$region == "US"), hgtL175=as.numeric(datnorm$height < 175)) fitdat <- cbind(datnorm, groups) # bind subgroup variables to main data ## subgroups of interest subgr <- c("labvalL.5", "regUS", "hgtL175") res <- modav(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "lm") summary(res) plot(res, show.compl=TRUE) ## compare to unadjusted analysis res <- unadj(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "lm") summary(res) plot(res) ## generate candidate subgroups using the subbuild function ## semi-automatically i.e. some groups specified directly (height and ## smoker), for region and labvalue subbuild generates subgroups (see ## ?subbuild). cand.groups <- subbuild(datnorm, height < 175, smoker == 1, region, labvalue) head(cand.groups) fitdat <- cbind(datnorm, cand.groups) subgr <- colnames(cand.groups) resMA <- modav(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "lm") summary(resMA) plot(resMA, show.compl = TRUE) ## toy example call for binary data on simulated datbin data-set data(datbin) cand.groups <- subbuild(datbin, height < 175, smoker == 1, region, labvalue) fitdat <- cbind(datbin, cand.groups) subgr <- colnames(cand.groups) res <- modav(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "glm", family = binomial(link = "logit")) ## scale of the treatment effect estimate: difference on log-odds scale summary(res) plot(res) ## toy example call for parametric and semi-parametric survival data on ## datsurv data-set data(datsurv) cand.groups <- subbuild(datsurv, height < 175, smoker == 1, region, labvalue) fitdat <- cbind(datsurv, cand.groups) subgr <- colnames(cand.groups) res.survreg <- modav(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "survreg", event = "event", dist = "exponential") ## parametric survival model (here exponential distribution) ## scale of treatment effect estimate: log scale (see ?survreg for details) summary(res.survreg) plot(res.survreg) res.cox <- modav(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "coxph", event = "event") ## scale of treatment effect estimate: difference in log-hazard rate summary(res.cox) plot(res.cox) ## toy example call overdispersed count data on datcount data-set data(datcount) cand.groups <- subbuild(datcount, height < 175, smoker == 1, region, labvalue) fitdat <- cbind(datcount, cand.groups) subgr <- colnames(cand.groups) res <- modav(resp = "y", trt = "treat", subgr = subgr, data = fitdat, covars = ~ x1 + x2, fitfunc = "glm.nb", exposure = "exposure") ## scale of treatment effect estimate: difference on log scale summary(res) plot(res)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.