R allows interaction spline functions, wide variety of predictor parameterizations, wide variety of models, unifying model formula language, model validation by resampling.

R is comprehensive:

6.1 The R Modeling Language

R statistical modeling language:

Code 
response ~ terms

y ~ age + sex            # age + sex main effects
y ~ age + sex + age:sex  # add second-order interaction
y ~ age*sex              # second-order interaction +
                         # all main effects
y ~ (age + sex + pressure)^2
                         # age+sex+pressure+age:sex+age:pressure...
y ~ (age + sex + pressure)^2 - sex:pressure
                         # all main effects and all 2nd order
                         # interactions except sex:pressure
y ~ (age + race)*sex     # age+race+sex+age:sex+race:sex
y ~ treatment*(age*race + age*sex) # no interact. with race,sex
sqrt(y) ~ sex*sqrt(age) + race
# functions, with dummy variables generated if
# race is an R factor (classification) variable
y ~ sex + poly(age,2)    # poly generates orthogonal polynomials
race.sex <- interaction(race,sex)
y ~ age + race.sex       # for when you want dummy variables for
                         # all combinations of the factors

The formula for a regression model is given to a modeling function, e.g.

Code 
lrm(y ~ rcs(x,4))

is read “use a logistic regression model to model y as a function of x, representing x by a restricted cubic spline with 4 default knots”1.

  • 1 lrm andrcs are in the rms package.

    • update function: re-fit model with changes in terms or data:

    Code 
    f  <- lrm(y ~ rcs(x,4) + x2 + x3)
f2 <- update(f, subset=sex=="male")
f3 <- update(f, .~.-x2)         # remove x2 from model
f4 <- update(f, .~. + rcs(x5,5))# add rcs(x5,5) to model
f5 <- update(f, y2 ~ .)         # same terms, new response var.

    6.2 User-Contributed Functions

    • R is high-level object-oriented language.
    • R runs on all platforms, and can be run on chromebooks using RStudio Server
    • Multitude of user-contributed functions freely available
    • International community of users

    Some R functions:

    • See Venables and Ripley
    • Hierarchical clustering: hclust
    • Principal components: princomp, prcomp
    • Canonical correlation: cancor
    • Nonparametric transform-both-sides additive models:
      ace, avas
    • Parametric transform-both-sides additive models:
      areg, areg.boot (Hmisc package in R))
    • Rank correlation methods:
      rcorr, hoeffd, spearman2 (Hmisc)
    • Variable clustering: varclus (Hmisc)
    • Single imputation: transcan (Hmisc)
    • Multiple imputation: aregImpute (Hmisc)
    • Restricted cubic splines:
      rcspline.eval (Hmisc)
    • Re-state restricted spline in simpler form:
      rcspline.restate (Hmisc)

    6.3 The rms Package

    FAQ
    • datadist function to compute predictor distribution summaries
    Code 
    y ~ sex + lsp(age,c(20,30,40,50,60)) +
    sex %ia% lsp(age,c(20,30,40,50,60))

    E.g. restrict age \(\times\) cholesterol interaction to be of form \(AF(B) + BG(A)\):

    Code 
    y ~ lsp(age,30) + rcs(cholesterol,4) +
    lsp(age,30) %ia% rcs(cholesterol,4)

    Special fitting functions by Harrell to simplify procedures described in these notes:

    rms Fitting Functions
    Function Purpose Related R Functions
    ols Ordinary least squares linear model lm
    lrm Binary and ordinal logistic regression model with options for penalize MLE glm
    orm Ordinal semi-parametric regression model for continuous \(Y\) and several link functions polr,lrm
    psm Accelerated failure time parametric survival models survreg
    cph Cox proportional hazards regression coxph
    bj Buckley-James censored least squares model survreg,lm
    Glm rms version of glm glm
    Gls rms version of gls gls (nlme package)
    Rq rms version of rq rq (quantreg package)
    rms Transformation Functions
    Function Purpose
    asis No post-transformation (seldom used explicitly)
    rcs Restricted cubic splines
    pol Polynomial using standard notation
    lsp Linear spline
    catg Categorical predictor (seldom)
    scored Ordinal categorical variables
    matrx Keep variables as group for anova and fastbw
    strat Non-modeled stratification factors (used for cph only)

    The transformation functions work also with regular R functions, e.g. when predict() is called the predicted values are computed by looking up the knot locations for rcs.

    Below notice that there are three graphic models implemented for depicting the effects of predictors in the fitted model: lattice graphics, a ggplot method using the ggplot2 package (which has an option to convert the result to plotly), and a direct plotly method. plotly is used to create somewhat interactive graphics with drill-down capability, and the rms package takes advantage of this capability. plotly graphics are best used with RStudio Rmarkdown html output.

    rms After-fit Functions
    Function Purpose Related Functions
    print Print parameters and statistics of fit
    coef Fitted regression coefficients
    formula Formula used in the fit
    specs Detailed specifications of fit
    vcov Fetch covariance matrix
    logLik Fetch maximized log-likelihood
    AIC Fetch AIC with option to put on chi-square basis
    lrtest Likelihood ratio test for two nested models
    univarLR Compute all univariable LR \(\chi^{2}\)
    robcov Robust covariance matrix estimates
    bootcov Bootstrap covariance matrix estimates and bootstrap distribution of estimates
    pentrace Find optimum penalty factors by tracing effective AIC for a grid of penalties
    effective.df Print effective d.f. for each type of variable in model, for penalized fit or pentrace results
    summary Summary of effects of predictors
    plot.summary Plot continuously shaded confidence bars for results of summary
    anova Wald and LR tests of most meaningful hypotheses
    plot.anova Graphical depiction of anova
    contrast General contrasts, C.L., tests
    gendata Easily generate predictor combinations
    predict Obtain predicted values or design matrix
    Predict Obtain predicted values and confidence limits easily varying a subset of predictors and others set at default values
    plot.Predict Plot the result of Predict using lattice
    ggplot.Predict Plot the result of Predict using ggplot2
    plotp.Predict Plot the result of Predict using plotly
    fastbw Fast backward step-down variable selection step
    residuals (or resid) Residuals, influence stats from fit
    sensuc Sensitivity analysis for unmeasured confounder
    which.influence Which observations are overly influential residuals
    latex \(\LaTeX\) representation of fitted model Function
    Function R function analytic representation of \(X\hat{\beta}\) from a fitted regression model latex
    Hazard R function analytic representation of a fitted hazard function (for psm)
    Survival R function analytic representation of fitted survival function (for psm, cph)
    Quantile R function analytic representation of fitted function for quantiles of survival time (for psm, cph)
    Mean R function analytic representation of fitted function for mean survival time or for ordinal logistic
    nomogram Draws a nomogram for the fitted model latex, plot
    survest Estimate survival probabilities (psm, cph) survfit
    survplot Plot survival curves (psm, cph) plot.survfit
    survplotp Plot survival curves with plotly features survplot
    validate Validate indexes of model fit using resampling
    val.prob External validation of a probability model lrm
    val.surv External validation of a survival model calibrate
    calibrate Estimate calibration curve using resampling val.prob
    vif Variance inflation factors for fitted model
    naresid Bring elements corresponding to missing data back into predictions and residuals
    naprint Print summary of missing values
    impute Impute missing values aregImpute
    rmsb: Bayesian Regression Modeling Strategies Package, Focusing on Semiparametric Univariate and Longitudinal Models
    Function Purpose
    blrm Bayesian binary and ordinal logistic model
    stackMI Bayesian posterior stacking for multiple imputation
    stanDx Stan diagnostics on fit
    stanDxplot Trace plots to check posterior sampling convergence
    PostF Creates R function for computing posterior probabilities
    plot.rmsb Plot posterior densities, intervals, point summaries
    compareBmods Compare two models using LOO-cv
    HPDint Compute highest posterior density interval
    distSym Compute meaure of symmetry of posterior distribution

    An extensive overview of Bayesian capabilities of the rmsb package may be found at hbiostat.org/R/rmsb/blrm.html.

    Global options prType and grType control printed and some graphical output, respectively as shown in example code below. The default is plain output and static graphics. If using plotly interactive graphics through ggplot or plotp or with anova or summary functions it is best to do so with RStudio html output or html notebooks. If using html output you must be producing an html document or notebook. When setting grType to use \(\LaTeX\) or html it is highly recommended that you use the knitr package.

    Example:

    • treat: categorical variable with levels "a","b","c"
    • num.diseases: ordinal variable, 0-4
    • age: continuous
      Restricted cubic spline
    • cholesterol: continuous
      (3 missings; use median)
      log(cholesterol+10)
    • Allow treat \(\times\) cholesterol interaction
    • Program to fit logistic model, test all effects in design, estimate effects (e.g. inter-quartile range odds ratios), plot estimated transformations
    Code 
    require(rms)                     # make new functions available
options(prType='html')           # print, summary, anova, validate: html output
                                 # others: 'latex', 'plain'
options(grType='plotly')         # plotly graphics for ggplot, anova, summary
                                 # default is 'base' for static graphics
ddist <- datadist(cholesterol, treat, num.diseases, age)
# Could have used ddist <- datadist(data.frame.name)
options(datadist="ddist")        # defines data dist. to rms
cholesterol <- impute(cholesterol)
fit <- lrm(y ~ treat + scored(num.diseases) + rcs(age) +
                 log(cholesterol+10) + treat:log(cholesterol+10),
                     x=TRUE, y=TRUE)  # needed for robcov, anova test='LR'
fit    # outputs plain, LaTeX, or html markup
describe(y ~ treat + scored(num.diseases) + rcs(age))
# or use describe(formula(fit)) for all variables used in fit
# describe function (in Hmisc) gets simple statistics on variables
# fit <- robcov(fit)             # Would make all statistics that follow
                                 # use a robust covariance matrix
                                 # would need x=T, y=T in lrm()
specs(fit)                       # Describe the design characteristics
anova(fit)                       # Wald tests, plain, LaTex, or html
anova(fit, test='LR')            # Likelihood ratio tests
anova(fit, treat, cholesterol)   # Test these 2 by themselves
plot(anova(fit))                 # Summarize anova graphically
summary(fit)                     # Estimate effects using default ranges
                                 # prints plain, LaTeX, or html
plot(summary(fit))               # Graphical display of effects with C.I.
summary(fit, treat="b", age=60)  # Specify reference cell and adjustment val
summary(fit, age=c(50,70))       # Estimate effect of increasing age from
                                 # 50 to 70
summary(fit, age=c(50,60,70))    # Increase age from 50 to 70, adjust to
                                 # 60 when estimating effects of other
                                 # factors
# If had not defined datadist, would have to define ranges for all var.

# Estimate and test treatment (b-a) effect averaged over 3 cholesterols
contrast(fit, list(treat='b', cholesterol=c(150,200,250)),
              list(treat='a', cholesterol=c(150,200,250)),
         type='average')
# See the help file for contrast.rms for several examples of
# how to obtain joint tests of multiple contrasts and how to get
# double differences (interaction contrasts)

p <- Predict(fit, age=seq(20,80,length=100), treat, conf.int=FALSE)
plot(p)                          # Plot relationship between age and log
# or ggplot(p), plotp(p)         # odds, separate curve for each treat,
                                 # no C.I.
plot(p, ~ age | treat)           # Same but 2 panels
ggplot(p, groups=FALSE)
bplot(Predict(fit, age, cholesterol, np=50))
                                 # 3-dimensional perspective plot for age,
                                 # cholesterol, and log odds using default
                                 # ranges for both variables
plot(Predict(fit, num.diseases, fun=function(x) 1/(1+exp(-x)), conf.int=.9),
     ylab="Prob")                # Plot estimated probabilities instead of
                                 # log odds (or use ggplot())
                                 # can also use plotp() for plotly
# Again, if no datadist were defined, would have to tell plot all limits
logit <- predict(fit, expand.grid(treat="b",num.dis=1:3,age=c(20,40,60),
                 cholesterol=seq(100,300,length=10)))
# Could also obtain list of predictor settings interactively}
logit <- predict(fit, gendata(fit, nobs=12))

# Since age doesn't interact with anything, we can quickly and
# interactively try various transformations of age, taking the spline
# function of age as the gold standard. We are seeking a linearizing
# transformation.

ag <- 10:80
logit <- predict(fit, expand.grid(treat="a", num.dis=0, age=ag,
                 cholesterol=median(cholesterol)), type="terms")[,"age"]
# Note: if age interacted with anything, this would be the age
#        "main effect" ignoring interaction terms
# Could also use
#    logit <- Predict(f, age=ag, ...)$yhat,
# which allows evaluation of the shape for any level of interacting
#  factors.  When age does not interact with anything, the result from
# predict(f, ..., type="terms") would equal the result from
# Predict if all other terms were ignored

# Could also specify
#    logit <- predict(fit, gendata(fit, age=ag, cholesterol=...))
# Un-mentioned variables set to reference values

plot(ag^.5, logit)               # try square root vs. spline transform.
plot(ag^1.5, logit)              # try 1.5 power

latex(fit)                       # fit in math notation
# Draw a nomogram for the model fit
plot(nomogram(fit))

# Compose R function to evaluate linear predictors analytically
g <- Function(fit)
g(treat='b', cholesterol=260, age=50)
# Letting num.diseases default to reference value

    To examine interactions in a simpler way, you may want to group age into tertiles:

    Code 
    age.tertile <- cut2(age, g=3)
# For automatic ranges later, add age.tertile to datadist input
fit <- lrm(y ~ age.tertile * rcs(cholesterol))

    6.4 Other Functions

    • processMI: works with Hmisc::fit.mult.impute to process resampling validation statistics and likelihood ratio \(\chi^2\) statistics accounting for multiple imputation
    • supsmu: Friedman’s “super smoother”
    • lowess: Cleveland’s scatterplot smoother
    • glm: generalized linear models (see Glm)
    • gam: Generalized additive models
    • rpart: Like original CART with surrogate splits for missings, censored data extension (Atkinson & Therneau)
    • validate.rpart: in rms; validates recursive partitioning with respect to certain accuracy indexes
    • loess: multi-dimensional scatterplot smoother
    Code 
    f <- loess(y ~ age * pressure)
plot(f)                           # cross-sectional plots
ages <- seq(20,70,length=40)
pressures <- seq(80,200,length=40)
pred <- predict(f, expand.grid(age=ages, pressure=pressures))
persp(ages, pressures, pred)      # 3-d plot