Description Usage Arguments Value See Also Examples
Using matrices of bootstrapped multinomial logistic regression coefficients, calculate and plot predicted probabilities and 95% confidence limits for each outcome at a given level of a continuous covariate. Also includes Wald p-values for each comparison. This function can handle covariates which are either linear or are modeled using restricted cubic splines with three knots. This function is commonly used with two model fits and corresponding coefficient matrices, representing the outcome with the lowest and highest level as reference, respectively. This allows for calculation of confidence limits for all levels of the outcome. Note that in this case, the coefficient matrices, variance-covariance matrices, and model objects should all be specified in a consistent order (eg, reference = lowest, followed by reference = highest).
1 2 3 4 5 |
xval |
String; name of variable to plot on X axis. |
xval.limits |
Numeric vector; percentiles of |
xval.knots |
Numeric vector of length 3; percentile values of knot placement if
|
data.set |
Data frame to use for quantiles and label of |
adjto.vals |
Named vector of values to adjust covariates to; length must equal the number of coefficients in the model, and names must match the column names of the elements of coef.list, ignoring the ":[0-9]" suffix. |
mod.objs |
List of |
coef.list |
List of matrices containing bootstrapped coefficient estimates. |
vcov.list |
List of variance-covariance matrices of coefficient estimates. Defaults to a list of var(coef.list[[x]]). |
plot.raw |
Add raw data to plots? Defaults to TRUE. |
xval.lab |
String to combine with Wald p-values for X axis label. If not specified, defaults to either label of data.set[,xval] if labelled, or xval if not labeled. |
pvals.oneline |
Whether to print Wald p-values on a single line or separate lines. Defaults to separate lines. |
yval.lab |
String to use for Y axis label. Defaults to 'Adjusted Probability of Outcome.' |
yval.limits |
Numeric vector for Y axis limits. Defaults to c(0, 1). |
alpha.lev |
Numeric value for transparency level of CI ribbons. Defaults to 1. |
List of 1) Data frame of all predicted probabilities and confidence
limits (pp.data
); 2) Data frame of Wald test results for each outcome
comparison (results.data
); 3) ggplot2 object which contains plots of
xval
vs. predicted probability of each outcome level.
vglm
, which this function assumes you are using;
rcs()
; Hmisc (label()
; aod (wald.test()
); ggplot2; dplyr; tidyr.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | ## Create data frame
df <- data.frame(id = sample(1:20, size = 100, replace = TRUE),
x1 = rnorm(n = 100),
x2 = rbinom(p = 0.75, n = 100, size = 1),
y = sample(LETTERS[1:3], size = 100, replace = TRUE))
df <- df[order(df$id),]
df$time <- unlist(lapply(1:length(unique(df$id)),
FUN = function(idnum){ 1:nrow(df[df$id == unique(df$id)[idnum],]) }))
## Using create.sampdata(), generate list of cluster bootstrapped data sets
bootdata.list <- create.sampdata(org.data = df,
id.var = 'id',
n.sets = 25)
## Fit model to original and bootstrapped data frame, saving errors and warnings to .txt file
boot.fits.a <- multi.bootstrap(org.data = df,
data.sets = bootdata.list,
ref.outcome = grep('A', levels(df$y)),
multi.form = as.formula('y ~ x1 + x2'))
## Create matrices of coefficients for all bootstrap fits
boot.matrix.a <- do.call(rbind,
lapply(boot.fits.a$boot.models,
FUN = function(x){ x@coefficients }))
## Calculate predicted probs and CIs for x1 at outcomes B, C
## Values to adjust to: intercept = 1, X2 = 1
vals.adjust <- c("(Intercept)" = 1, "x2" = 1)
## Calculate and plot predicted probabilities for outcomes at levels of x1
## To get probabilities and CIs for all outcome levels, run models twice,
## with lowest and highest outcome levels as reference
## Fit model to original and bootstrapped data frame, saving errors and warnings to .txt file
boot.fits.c <- multi.bootstrap(org.data = df,
data.sets = bootdata.list,
ref.outcome = grep('C', levels(df$y)),
multi.form = as.formula('y ~ x1 + x2'))
boot.matrix.c <- do.call(rbind,
lapply(boot.fits.c$boot.models,
FUN = function(x){ x@coefficients }))
x1.prob.results <- multi.plot.probs(xval = 'x1',
data.set = df,
adjto.vals = vals.adjust,
mod.objs = list(boot.fits.a$org.mod,
boot.fits.c$org.mod),
coef.list = list(boot.matrix.a,
boot.matrix.c),
vcov.list = list(var(boot.matrix.a),
var(boot.matrix.c)))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.