Plot Effects of Variables Estimated by a Regression Model Fit Using plotly
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.
subset: a subsetting expression for restricting the rows of data that are used in plotting. For example, predictions may have been requested for males and females but one wants to plot only females.
xlim: ignored unless predictors were specified to Predict. Specifies the x-axis limits of the single plot produced.
ylim: Range for plotting on response variable axis. Computed by default and includes the confidence limits.
xlab: Label for x-axis when a single plot is made, i.e., when a predictor is specified to Predict. Default is one given to asis, rcs, etc., which may have been the "label" attribute of the variable.
ylab: Label for y-axis. If fun is not given, default is "log Odds" for lrm, "log Relative Hazard" for cph, name of the response variable for ols, TRUE or log(TRUE) for psm, or "X * Beta" otherwise. Specify ylab=NULL to omit y-axis labels.
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 rdata is present and contains the needed variables, the original data are added to the graph in the form of a spike histogram using histSpikeg in the Hmisc package.
nlevels: A non-numeric x-axis variable with nlevels or fewer unique values will cause a horizontal dot plot to be drawn instead of an x-y plot.
vnames: applies to the case where multiple plots are produced separately by predictor. Set to 'names' to use variable names instead of labels for these small plots.
histSpike.opts: a list containing named elements that specifies parameters to histSpikeg when rdata is given. The col parameter is usually derived from other plotting information and not specified by the user.
ncols: number of columns of plots to use when plotting multiple continuous predictors
width: width in pixels for plotly graphics
...: ignored
Returns
a plotly object or a list containing two elements, each one a plotly object
Author(s)
Frank Harrell
Department of Biostatistics, Vanderbilt University
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 sizeset.seed(17)# so can reproduce the resultsage <- 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 Hmisclabel(cholesterol)<-'Total Cholesterol'label(blood.pressure)<-'Systolic Blood Pressure'label(sex)<-'Sex'units(cholesterol)<-'mg/dl'# uses units.default in Hmiscunits(blood.pressure)<-'mmHg'# Specify population model for log odds that Y=1L <-.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 objectsplotp(Predict(fit), rdata=llist(blood.pressure, age))$Continuous
# spike histogram plot for two of the predictorsp <- Predict(fit, name=c('age','cholesterol'))# Make 2 plotsplotp(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 thatp <- Predict(fit, age=seq(20,80,length=100), sex='male', fun=plogis)# works if datadist not usedplotp(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 graphp1 <- 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)