Plot Effects of Variables Estimated by a Regression Model Fit Using plotly
Uses plotly
graphics (without using ggplot2) to plot the effect
of one or two predictors
on the linear predictor or X beta scale, or on some transformation of
that scale. The first argument specifies the result of the
Predict
function. The predictor is always plotted in its
original coding. Hover text shows point estimates, confidence
intervals, and on the leftmost x-point, adjustment variable settings.
If Predict
was run with no variable settings, so that each
predictor is varied one at a time, the result of plotp.Predict
is a list with two elements. The first, named Continuous
, is a
plotly
object containing a single graphic with all the
continuous predictors varying. The second, named Categorical
,
is a plotly
object containing a single graphic with all the
categorical predictors varying. If there are no categorical
predictors, the value returned by by plotp.Predict
is a single
plotly
object and not a list of objects.
If rdata
is given, a spike histogram is drawn showing
the location/density of data values for the x-axis variable. If
there is a superposition variable that generated separate
curves, the data density specific to each class of points is shown.
The histograms are drawn by histSpikeg
.
To plot effects instead of estimates (e.g., treatment differences as a
function of interacting factors) see contrast.rms
and
summary.rms
.
Unlike ggplot.Predict
, plotp.Predict
does not handle
groups
, anova
, or perim
arguments.
## S3 method for class 'Predict' plotp(data, subset, xlim, ylim, xlab, ylab, rdata=NULL, nlevels=3, vnames=c('labels','names'), histSpike.opts=list(frac=function(f) 0.01 + 0.02 * sqrt(f - 1)/sqrt(max(f, 2) - 1), side=1, nint=100), ncols=3, width=800, ...)
data |
a data frame created by |
subset |
a subsetting expression for restricting the rows of
|
xlim |
ignored unless predictors were specified to |
ylim |
Range for plotting on response variable axis. Computed by default and includes the confidence limits. |
xlab |
Label for |
ylab |
Label for |
rdata |
a data frame containing the original raw data on which the
regression model were based, or at least containing the x-axis
and grouping variable. If |
nlevels |
A non-numeric x-axis variable with |
vnames |
applies to the case where multiple plots are produced
separately by predictor. Set to |
histSpike.opts |
a list containing named elements that specifies
parameters to |
ncols |
number of columns of plots to use when plotting multiple continuous predictors |
width |
width in pixels for |
... |
ignored |
a plotly
object or a list containing two elements, each
one a plotly
object
Frank Harrell
Department of Biostatistics, Vanderbilt University
fh@fharrell.com
Fox J, Hong J (2009): Effect displays in R for multinomial and proportional-odds logit models: Extensions to the effects package. J Stat Software 32 No. 1.
## Not run: n <- 350 # define sample size set.seed(17) # so can reproduce the results age <- rnorm(n, 50, 10) blood.pressure <- rnorm(n, 120, 15) cholesterol <- rnorm(n, 200, 25) sex <- factor(sample(c('female','male'), n,TRUE)) label(age) <- 'Age' # label is in Hmisc label(cholesterol) <- 'Total Cholesterol' label(blood.pressure) <- 'Systolic Blood Pressure' label(sex) <- 'Sex' units(cholesterol) <- 'mg/dl' # uses units.default in Hmisc units(blood.pressure) <- 'mmHg' # Specify population model for log odds that Y=1 L <- .4*(sex=='male') + .045*(age-50) + (log(cholesterol - 10)-5.2)*(-2*(sex=='female') + 2*(sex=='male')) + .01 * (blood.pressure - 120) # Simulate binary y to have Prob(y=1) = 1/[1+exp(-L)] y <- ifelse(runif(n) < plogis(L), 1, 0) ddist <- datadist(age, blood.pressure, cholesterol, sex) options(datadist='ddist') fit <- lrm(y ~ blood.pressure + sex * (age + rcs(cholesterol,4)), x=TRUE, y=TRUE) p <- plotp(Predict(fit)) p$Continuous p$Categorical # When using Rmarkdown html notebook, best to use # prList(p) to render the two objects plotp(Predict(fit), rdata=llist(blood.pressure, age))$Continuous # spike histogram plot for two of the predictors p <- Predict(fit, name=c('age','cholesterol')) # Make 2 plots plotp(p) p <- Predict(fit, age, sex) plotp(p, rdata=llist(age,sex)) # rdata= allows rug plots (1-dimensional scatterplots) # on each sex's curve, with sex- # specific density of age # If data were in data frame could have used that p <- Predict(fit, age=seq(20,80,length=100), sex='male', fun=plogis) # works if datadist not used plotp(p, ylab='P') # plot predicted probability in place of log odds # Compute predictions for three predictors, with superpositioning or # conditioning on sex, combined into one graph p1 <- Predict(fit, age, sex) p2 <- Predict(fit, cholesterol, sex) p3 <- Predict(fit, blood.pressure, sex) p <- rbind(age=p1, cholesterol=p2, blood.pressure=p3) plotp(p, ncols=2, rdata=llist(age, cholesterol, sex)) ## End(Not run)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.