Estimate, compare, and visualize exposure histories
Source:R/compareHistories.R
compareHistories.Rd
Takes fitted model output to created predicted values for user-specified histories (pooling for imputed data), before conducting contrast comparisons (pooling for imputed data), correcting for multiple comparisons, and then plotting results.
Usage
compareHistories(
fit,
hi_lo_cut,
dose_level = "h",
reference = NULL,
comparison = NULL,
mc_comp_method = "BH",
verbose = FALSE,
save.out = FALSE
)
# S3 method for class 'devMSM_comparisons'
print(x, save.out = FALSE, ...)
# S3 method for class 'devMSM_comparisons'
plot(
x,
colors = "Dark2",
exp_lab = NULL,
out_lab = NULL,
save.out = FALSE,
...
)
# S3 method for class 'devMSM_comparisons'
summary(object, type = "comps", ...)
Arguments
- fit
list of model outputs from
fitModel()
- hi_lo_cut
list of two numbers indicating quantile values that reflect high and low values, respectively, for continuous exposure
- dose_level
(optional) "l" or "h" indicating whether low or high doses should be tallied in tables and plots (default is high "h")
- reference
lists of one or more strings of "-"-separated "l" and "h" values indicative of a reference exposure history to which to compare comparison, required if comparison is supplied
- comparison
(optional) list of one or more strings of "-"-separated "l" and "h" values indicative of comparison history/histories to compare to reference, required if reference is supplied
- mc_comp_method
(optional) character abbreviation for multiple comparison correction method for stats::p.adjust, default is Benjamini-Hochburg ("BH")
- verbose
(optional) TRUE or FALSE indicator for printing output to console. default is FALSE.
- save.out
(optional) Either logical or a character string. If
TRUE
, it will output the result to a default file name withinhome_dir
set ininitMSM()
. You can load the data withx <- readRDS(file)
. To use a non-default file name, specify a character string with the file name. It will save relative tohome_dir
. There might be naming conflicts where two objects get saved to the same file. In these cases, users should specify a custom name. default is FALSE.- x
devMSM_histories object from
compareHistories()
- ...
ignored
- colors
(optional) character specifying Brewer palette or list of colors (n(epochs)+1) for plotting (default is "Dark2" palette)
- exp_lab
(optional) character label for exposure variable in plots (default is variable name)
- out_lab
(optional) character label for outcome variable in plots (default is variable name)
- object
devMSM_histories object from
compareHistories()
- type
Either "preds" or "comps" corresponding to the results of
marginaleffects::avg_predictions()
at low and high dosages ormarginaleffects::avg_comparisons()
respectively
Value
list containing two dataframes: preds
with predictions from
marginaleffects::avg_predictions()
containing average expected outcome
for different exposure histories and comps
with contrasts from
marginaleffects::comparisons()
comparing different exposure history
See also
marginaleffects::avg_predictions()
,
https://cran.r-project.org/web/packages/marginaleffects/marginaleffects.pdf;
marginaleffects::hypotheses()
,
https://cran.r-project.org/web/packages/marginaleffects/marginaleffects.pdf;
stats::p.adjust()
,
https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/p.adjust;
Examples
library(devMSMs)
set.seed(123)
data <- data.frame(
ID = 1:50,
A.1 = rnorm(n = 50),
A.2 = rnorm(n = 50),
A.3 = rnorm(n = 50),
B.1 = rnorm(n = 50),
B.2 = rnorm(n = 50),
B.3 = rnorm(n = 50),
C = rnorm(n = 50),
D.3 = rnorm(n = 50)
)
obj <- initMSM(
data,
exposure = c("A.1", "A.2", "A.3"),
ti_conf = c("C"),
tv_conf = c("B.1", "B.2", "B.3", "D.3")
)
f <- createFormulas(obj, type = "short")
w <- createWeights(data = data, formulas = f)
fit <- fitModel(
data = data, weights = w,
outcome = "D.3", model = "m0"
)
comp <- compareHistories(
fit = fit,
hi_lo_cut = c(0.3, 0.6)
)
#> Error: Assertion failed. One of the following must apply:
#> * checkmate::check_character(hypothesis): Must comply to pattern
#> * '=|<=|>='
#> * checkmate::check_numeric(hypothesis): Must be of type 'numeric', not
#> * 'character'
#> * checkmate::check_formula(hypothesis): Must be a formula, not
#> * character
#> * checkmate::check_matrix(hypothesis): Must be of type 'matrix', not
#> * 'character'
#> * checkmate::check_function(hypothesis): Must be a function, not
#> * 'character'
#> * checkmate::check_null(hypothesis): Must be NULL
print(comp)
#> Error: object 'comp' not found
plot(comp)
#> Error: object 'comp' not found
summary(comp, "preds")
#> Error: object 'comp' not found
summary(comp, "comps")
#> Error: object 'comp' not found
comp2 <- compareHistories(
fit = fit,
hi_lo_cut = c(0.3, 0.6),
reference = "l-l-l",
comparison = c("h-h-h", "h-h-l")
)
print(comp2)
#> Summary of Exposure Main Effects:
#> Warning: There are no participants in your sample in the following histories: l-l-l.
#> Please revise your reference/comparison histories and/or the high/low cutoffs, if applicable.
#> USER ALERT: Out of the total of 50 individuals in the sample, below is the distribution of the 6 (12%) individuals that fall into 2 user-selected exposure histories (out of the 24 total) created from 30th and 60th percentile values for low and high levels of exposure-epoch A.1, A.2, A.3.
#> USER ALERT: Please inspect the distribution of the sample across the following exposure histories and ensure there is sufficient spread to avoid extrapolation and low precision:
#>
#> +---------------+---+
#> | epoch_history | n |
#> +===============+===+
#> | h-h-h | 3 |
#> +---------------+---+
#> | h-h-l | 3 |
#> +---------------+---+
#>
#> Table: Summary of user-selected exposure histories based on exposure main effects A.1, A.2, A.3:
#>
#> Below are the pooled average predictions by user-specified history:
#> +-------+-------+------+--------+----------+-----------+----------+-----------+
#> | term | A.1 | A.2 | A.3 | estimate | std.error | conf.low | conf.high |
#> +=======+=======+======+========+==========+===========+==========+===========+
#> | l-l-l | -0.47 | -0.3 | -0.804 | -0.15 | 0.2 | -0.54 | 0.24 |
#> +-------+-------+------+--------+----------+-----------+----------+-----------+
#> | h-h-l | 0.24 | 0.31 | -0.804 | -0.179 | 0.13 | -0.43 | 0.07 |
#> +-------+-------+------+--------+----------+-----------+----------+-----------+
#> | h-h-h | 0.24 | 0.31 | -0.051 | -0.093 | 0.14 | -0.36 | 0.18 |
#> +-------+-------+------+--------+----------+-----------+----------+-----------+
#>
#> Conducting multiple comparison correction for all pairings between comparison histories and each reference history using the BH method.
#>
#>
#> +-------------------+-------------+-----------+-----------+------------+-----------+--------------+
#> | term | estimate | std.error | p.value | conf.low | conf.high | p.value_corr |
#> +===================+=============+===========+===========+============+===========+==============+
#> | (h-h-h) - (l-l-l) | 0.05723992 | 0.1571878 | 0.7157462 | -0.2508426 | 0.3653224 | 0.8193787 |
#> +-------------------+-------------+-----------+-----------+------------+-----------+--------------+
#> | (h-h-l) - (l-l-l) | -0.02878467 | 0.1260583 | 0.8193787 | -0.2758544 | 0.2182850 | 0.8193787 |
#> +-------------------+-------------+-----------+-----------+------------+-----------+--------------+
plot(comp2)
summary(comp2, "preds")
#> term A.1 A.2 A.3 estimate std.error statistic
#> 1 l-l-l -0.4684962 -0.2971206 -0.80434443 -0.14979959 0.1983536 -0.7552150
#> 7 h-h-l 0.2359494 0.3148300 -0.80434443 -0.17858426 0.1267806 -1.4086089
#> 8 h-h-h 0.2359494 0.3148300 -0.05134827 -0.09255967 0.1371104 -0.6750739
#> p.value s.value conf.low conf.high df dose
#> 1 0.4501200 1.151619 -0.5385654 0.23896625 Inf 0
#> 7 0.1589509 2.653347 -0.4270696 0.06990112 Inf 2
#> 8 0.4996288 1.001071 -0.3612912 0.17617185 Inf 3
summary(comp2, "comps")
#> term estimate std.error statistic p.value s.value
#> 1 (h-h-h) - (l-l-l) 0.05723992 0.1571878 0.3641498 0.7157462 0.4824800
#> 2 (h-h-l) - (l-l-l) -0.02878467 0.1260583 -0.2283442 0.8193787 0.2873977
#> conf.low conf.high dose p.value_corr
#> 1 -0.2508426 0.3653224 3 - 0 0.8193787
#> 2 -0.2758544 0.2182850 2 - 0 0.8193787