Plot Effects of Variables

Description

Plots the effect of one or two predictors on the linear predictor or X beta scale, or on some transformation of that scale. The predictor is always plotted in its original coding on the x or y-axis. perimeter is a function used to generate the boundary of data to plot when a 3-d plot is made. It finds the area where there are sufficient data to generate believable interaction fits.

Legend is a generic function for adding legends to an existing graph according to the specific plot made by plot.Design. The specific Legend method for plot.Design is Legend.plot.Design. It handles legends for image plots. For other plots with one or more curves, make legends using the label.curves parameter.

datadensity is a function for showing the data density (raw data) on each curve generated for curve-type plots. This is a rug plot showing the location/density of data values for the x-axis variable. If there was a second variable specified to plot that generated separate curves, the data density specific to each class of points is shown. This assumes that the second variable was categorical. The rug plots are drawn by scat1d.

To plot effects instead of estimates (e.g., treatment differences as a function of interacting factors) see contrast.Design and summary.Design.

Usage

perimeter(x, y, xinc=diff(range(x))/10, 
          n=10, lowess.=TRUE)

## S3 method for class 'Design':
plot(x, ..., xlim, ylim, fun, xlab, ylab, 
     conf.int=.95, conf.type=c('mean','individual'),
     add=FALSE, label.curves=TRUE,
     eye, theta=0, phi=15, perspArgs=NULL,
     lty, col=1, lwd=par('lwd'), lwd.conf=1, pch=1, 
     adj.zero=FALSE, ref.zero=FALSE, adj.subtitle, cex.adj,
     non.slopes, time=NULL, loglog=FALSE, val.lev=FALSE,
     digits=4, log="", perim,
     method=c("persp","contour","image","dotchart","default"),
     sortdot=c('neither','ascending','descending'),
     nlevels=10, name, zlim=range(zmat,na.rm=TRUE),
     vnames=c('labels','names'), abbrev=FALSE)
# Just say plot(fit, ...)

## S3 method for class 'plot.Design':
print(x, ...)

## S3 method for class 'perimeter':
lines(x, ...)

Legend(object, ...)
## S3 method for class 'plot.Design':
Legend(object, x, y, size=c(1,1), horizontal=TRUE,
       nint=50, fun, at, zlab, ...)

## S3 method for class 'plot.Design':
datadensity(object, x1, x2, ...)

Arguments

Details

When there are no intercepts in the fitted model, plot subtracts adjustment values from each factor while computing variances for confidence limits.

perimeter is a kind of generalization of datadist for 2 continuous variables. First, the n smallest and largest x values are determined. These form the lowest and highest possible xs to display. Then x is grouped into intervals bounded by these two numbers, with the interval widths defined by xinc. Within each interval, y is sorted and the nth smallest and largest y are taken as the interval containing sufficient data density to plot interaction surfaces. The interval is ignored when there are insufficient y values. When plot.Design readies the data for persp, it uses the approx function to do linear interpolation of the y-boundaries as a function of the x values actually used in forming the grid (the values of the first variable specified to plot). To make the perimeter smooth, specify lowess.=TRUE to perimeter.

Specifying time will not work for Cox models with time-dependent covariables. Use survest or survfit for that purpose.

Use ps.slide, win.slide, gs.slide to set up nice defaults for plotting. These also set a system option mgp.axis.labels to allow x and y-axes to have differing mgp graphical parameters (see par). This is important when labels for y-axis tick marks are to be written horizontally (par(las=1)), as a larger gap between the labels and the tick marks are needed. You can set the axis-specific 2nd components of mgp using mgp.axis.labels(c(xvalue,yvalue)).

Note that because the generic plot method has the variable x as its first argument, you cannot explicitly specify that you want to plot the effect of a predictor named x.

Value

perimeter returns a matrix of class perimeter. This outline can be conveniently plotted by lines.perimeter. Legend.plot.Design invisibly returns the position of the legend. plot.Design invisibly returns an invisible object of class "plot.Design" with the following components (use print.plot.Design to get a nice printout of the object):

Author(s)

Frank Harrell
Department of Biostatistics, Vanderbilt University
f.harrell@vanderbilt.edu

See Also

datadist, predict.Design, contrast.Design, summary.Design, persp, Design, Design.trans, survest, survplot, Design.Misc, contour, image, labcurve, scat1d, dotchart2, mgp.axis.labels Overview, par, ps.slide, xYplot, smearingEst

Examples

n <- 1000    # 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'))
# 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)

par(mfrow=c(2,2))
plot(fit)                # Plot effects of all 4 predictors
par(mfrow=c(1,2))
plot(fit, name=c('age','cholesterol'))   # Make 2 plots
par(mfrow=c(1,1))
plot(fit, age=seq(20,80,length=100), sex=NA, conf.int=FALSE)
                         # Plot relationship between age and log
                         # odds, separate curve for each sex,
                         # no C.I.
z <- plot(fit, age=NA, sex=NA, label.curves=FALSE)
                         # use label.curves=list(keys=c('a','b'))'
                         # to use 1-letter abbreviations
datadensity(z, age, sex) # rug plots (1-dimensional scatterplots)
                         # on each treatment curve, with treatment-
                         # specific density of age
plot(fit, age=seq(20,80,length=100), sex='male')  # works if datadist not used
plot(fit, age=NA, cholesterol=NA)# 3-dimensional perspective plot for age,
                         # cholesterol, and log odds using default
                         # ranges for both variables
boundaries <- perimeter(age, cholesterol, lowess=TRUE)
plot(age, cholesterol)   # show bivariate data density
lines(boundaries)        # and perimeter that will be used for 3-D plot
z <- plot(fit, age=NA, cholesterol=NA, perim=boundaries, method='image')
                         # draws image() plot
                         # don't show estimates where data are sparse
                         # doesn't make sense here since vars don't interact
if(!.R.)Legend(z, fun=plogis, at=qlogis(c(.01,.05,.1,.2,.3,.4,.5)),
               zlab='Probability')   # gray scale or color legend for prob.
plot(fit, age=NA, fun=function(x) 1/(1+exp(-x)) , # or fun=plogis
     ylab="Prob", conf.int=.9)    # Plot estimated probabilities instead of
                                  # log odds

# Plot the age effect as an odds ratio
# comparing the age shown on the x-axis to age=30 years

ddist$limits$age[2] <- 30    # make 30 the reference value for age
# Could also do: ddist$limits["Adjust to","age"] <- 30
fit <- update(fit)   # make new reference value take effect
plot(fit, age=NA, ref.zero=TRUE, fun=exp, ylab='Age=x:Age=30 Odds Ratio')
abline(h=1, lty=2, col=2); abline(v=30, lty=2, col=2)

# Make two curves, and plot the predicted curves as two trellis panels
w <- plot(fit, age=NA, sex=NA)   # Would be nice if a pl=FALSE option was avail.
z <- data.frame(w$x.xbeta)     # Makes variable names legal
if(.R.) library(lattice)
xyplot(log.odds ~ age | sex, data=z, type='l')
# To add confidence bands we need to use the Hmisc xYplot function in
# place of xyplot
xYplot(Cbind(log.odds,lower,upper) ~ age | sex, data=z, 
       method='bands', type='l')
# If non-displayed variables were in the model, add a subtitle to show
# their settings using title(sub=paste('Adjusted to',w$adjust),adj=0)
# See predict.Design for an example using predict and xYplot without plot()



# Plots for a parametric survival model
n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "Age"
sex <- factor(sample(c('Male','Female'), n, 
              rep=TRUE, prob=c(.6, .4)))
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
t <- -log(runif(n))/h
label(t) <- 'Follow-up Time'
e <- ifelse(t<=cens,1,0)
t <- pmin(t, cens)
units(t) <- "Year"
ddist <- datadist(age, sex)
Srv <- Surv(t,e)

# Fit log-normal survival model and plot median survival time vs. age
f <- psm(Surv(t, e) ~ rcs(age), dist=if(.R.)'lognormal' else 'gaussian')
med <- Quantile(f)       # Creates function to compute quantiles
                         # (median by default)
plot(f, age=NA, fun=function(x)med(lp=x), ylab="Median Survival Time")
# Note: This works because med() expects the linear predictor (X*beta)
#       as an argument.  Would not work if use 
#       plot(..., ref.zero=TRUE or adj.zero=TRUE)
# Also, confidence intervals from this method are approximate since
# they don't take into account estimation of scale parameter

# Fit an ols model to log(y) and plot the relationship between x1
# and the predicted mean(y) on the original scale without assuming
# normality of residuals; use the smearing estimator
set.seed(1)
x1 <- runif(300)
x2 <- runif(300)
ddist <- datadist(x1,x2)
y  <- exp(x1+x2-1+rnorm(300))
f <- ols(log(y) ~ pol(x1,2)+x2)
r <- resid(f)
smean <- function(yhat)smearingEst(yhat, exp, res, statistic='mean')
formals(smean) <- list(yhat=numeric(0), res=r[!is.na(r)])
#smean$res <- r[!is.na(r)]   # define default res argument to function
plot(f, x1=NA, fun=smean, ylab='Predicted Mean on y-scale')

options(datadist=NULL)

## Not run: 
# Example in which separate curves are shown for 4 income values
# For each curve the estimated percentage of voters voting for
# the democratic party is plotted against the percent of voters
# who graduated from college.  scat1d is used to indicate
# the income-interval-specific data density for college.  For
# this purpose show the distribution of percent in college for
# those having an income level within +/- the half-width of
# the income interval.  scat1d shows the rug plot superimposed
# on the estimated curve.  Data are county-level percents.
# This can't be done automatically using datadensity on the object
# returned by plot.Design, as the variable representing different
# curves (income) is a continuous variable.

incomes <- seq(22900, 32800, length=4)  
# equally spaced to outer quintiles
pl <- plot(f, college=NA, income=incomes, 
           conf.int=FALSE, xlim=c(0,35), ylim=c(30,55),
           lty=1, lwd=c(.25,1.5,3.5,6), col=c(1,1,2,2))
graph.points <- pl$x.xbeta
for(i in 1:4) {
  college.in.income.group <- college[abs(income-incomes[i]) < 1650]
  this.income <- graph.points[,'income']==incomes[i]
  scat1d(college.in.income.group,
         curve=list(x=graph.points[this.income,'college'],
           y=graph.points[this.income,'democrat']))
}

# Instead of showing a rug plot on each curve, erase end portions
# of each curve where there are fewer than 10 counties having
# % college graduates to the left of the x-coordinate being plotted,
# for the subset of counties having median family income with 1650
# of the target income for the curve

show.pts <- function(college.pts, income.pt) {
  s <- abs(income - income.pt) < 1650  #assumes income known to top frame
  x <- college[s]
  x <- sort(x[!is.na(x)])
  n <- length(x)
  low <- x[10]; high <- x[n-9]
  college.pts >= low & college.pts <= high
}

plot(f, college=NA, income=incomes,
     conf.int=FALSE, xlim=c(0,35), ylim=c(30,55),
     lty=1, lwd=c(.25,1.5,3.5,6), col=c(1,1,2,2),
     perim=show.pts)
## End(Not run)