multi.plot.probs: Calculate and Plot Predicted Probabilities and 95% Confidence...

In jenniferthompson/ClusterBootMultinom: Multinomial Logistic Regression Using Clustered Bootstrapping

Description

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).

Usage

multi.plot.probs(xval, xval.limits = c(0.1, 0.9), xval.knots = c(0.05, 0.5,
  0.95), data.set, adjto.vals, mod.objs, coef.list, vcov.list = NULL,
  plot.raw = TRUE, xval.lab = NULL, pvals.oneline = FALSE,
  yval.lab = "Adjusted Probability of Outcome", yval.limits = c(0, 1),
  alpha.lev = 1)

Arguments

xval	String; name of variable to plot on X axis.
xval.limits	Numeric vector; percentiles of xval to include. Defaults to excluding lower and upper 10% of X values.
xval.knots	Numeric vector of length 3; percentile values of knot placement if xval is modeled with restricted cubic splines. Defaults to c(0.05, 0.50, 0.95) - rcspline.eval() defaults.
data.set	Data frame to use for quantiles and label of xval.
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 vglm() model objects - original model fit(s).
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.

Value

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.

See Also

vglm, which this function assumes you are using; rcs(); Hmisc (label(); aod (wald.test()); ggplot2; dplyr; tidyr.

Examples

## 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)))

jenniferthompson/ClusterBootMultinom documentation built on May 19, 2019, 4:03 a.m.