Package 'sure'

Title: Surrogate Residuals for Ordinal and General Regression Models
Description: An implementation of the surrogate approach to residuals and diagnostics for ordinal and general regression models; for details, see Liu and Zhang (2017, <doi:https://doi.org/10.1080/01621459.2017.1292915>) and Greenwell et al. (2017, <https://journal.r-project.org/archive/2018/RJ-2018-004/index.html>). These residuals can be used to construct standard residual plots for model diagnostics (e.g., residual-vs-fitted value plots, residual-vs-covariate plots, Q-Q plots, etc.). The package also provides 'plot' methods for producing standard diagnostic plots using lightweight base R graphics (via 'tinyplot'). The package currently supports cumulative link models from packages 'MASS', 'ordinal', 'rms', and 'VGAM'. Support for binary regression models using the standard 'glm' function is also available.
Authors: Brandon Greenwell [aut, cre], Andrew McCarthy [aut], Boehmke Bradley [aut], Liz Glaser [ctb]
Maintainer: Brandon Greenwell <[email protected]>
License: GPL (>= 2)
Version: 0.3.0
Built: 2026-07-07 04:53:42 UTC
Source: https://github.com/bgreenwell/sure

Help Index

Calculate surrogate residuals for multinomial/discrete choice models

Description

Calculate surrogate residuals for multinomial logit and general discrete choice models using a fast, exact conditional Gumbel simulation method.

Usage

calc_resid(y, proba)

Arguments

y

Vector of observed classes (either a factor or a 1-based integer vector).
proba

Matrix of predicted probabilities, where rows correspond to observations and columns correspond to classes.

Value

A matrix of surrogate residuals of the same dimensions as proba.

Goodness-of-Fit Simulation

Description

Simulate p-values from a goodness-of-fit test.

Usage

gof(object, nsim = 10, test = c("ks", "ad", "cvm"), ...)

## Default S3 method:
gof(object, nsim = 10, test = c("ks", "ad", "cvm"), ...)

## S3 method for class 'gof'
plot(x, ...)

Arguments

object

An object of class ordinal::clm(), stats::glm(), rms::lrm(), rms::orm(), MASS::polr(), or VGAM::vglm().
nsim

Integer specifying the number of bootstrap replicates to use.
test

Character string specifying which goodness-of-fit test to use. Current options include: "ks" for the Kolmogorov-Smirnov test, "ad" for the Anderson-Darling test, and "cvm" for the Cramer-Von Mises test. Default is "ks".
...

Additional optional arguments. (Currently ignored.)
x

An object of class "gof".

Details

Under the null hypothesis, the distribution of the p-values should appear uniformly distributed on the interval ⁠[0, 1]⁠. This can be visually investigated using the plot method. A 45 degree line is indicative of a "good" fit.

Value

A numeric vector of class ⁠"gof", "numeric"⁠ containing the simulated p-values.

Examples

# See ?resids for an example
?resids

Plot Residual Diagnostics

Description

Residual-based diagnostic plots for cumulative link and general regression models using tinyplot::tinyplot() graphics.

Usage

## S3 method for class 'resid'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

## S3 method for class 'clm'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

## S3 method for class 'glm'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

## S3 method for class 'lrm'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

## S3 method for class 'orm'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

## S3 method for class 'polr'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

## S3 method for class 'vglm'
plot(
  x,
  y = NULL,
  what = c("qq", "fitted", "covariate"),
  covariate = NULL,
  fit = NULL,
  distribution = qnorm,
  ncol = NULL,
  alpha = 1,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  smooth = TRUE,
  ...
)

Arguments

x

An object of class "resid" (returned by resids()), or a fitted model object of class ordinal::clm(), stats::glm(), rms::lrm(), rms::orm(), MASS::polr(), or VGAM::vglm().
y

Ignored.
what

Character string specifying what to plot. Default is "qq" which produces a quantile-quantile plot of the residuals.
covariate

A vector giving the covariate values to use for residual-by- covariate plots (i.e., when what = "covariate").
fit

The fitted model from which the residuals were extracted. (Only required if what = "fitted" and x inherits from class "resid".)
distribution

Function that computes the quantiles for the reference distribution to use in the quantile-quantile plot. Default is qnorm which is only appropriate for models using a probit link function. When jitter_scale = "probability", the reference distribution is always U(-0.5, 0.5).
ncol

Integer specifying the number of columns to use for the plot layout (if requesting multiple plots). Default is NULL.
alpha

A single value in the interval ⁠[0, 1]⁠ controlling the opacity alpha of the plotted points. Only used when nsim > 1.
xlab

Character string giving the text to use for the x-axis label. Default is NULL.
ylab

Character string giving the text to use for the y-axis label. Default is NULL.
main

Character string giving the plot title. Default is NULL.
smooth

Logical indicating whether or not to add a nonparametric smooth to certain plots. Default is TRUE.
...

Additional optional arguments to be passed onto tinyplot::tinyplot() (e.g., col, pch, cex, etc.) or resids() (if calling the plot method directly on a fitted model object).

Value

No return value, called for plotting.

Note

For vglm (and vgam) objects, the VGAM package defines its own S4 plot method which will override this S3 method if VGAM is loaded. In such cases, we recommend extracting the residuals first using resids() and then plotting the resulting resid object (e.g. plot(resids(fit))), which will correctly dispatch to plot.resid().

Surrogate residuals

Description

Simulate surrogate response values for cumulative link regression models using the latent method described in Liu and Zhang (2017).

Usage

resids(object, ...)

## Default S3 method:
resids(
  object,
  nsim = 1L,
  method = c("latent", "jitter"),
  jitter_scale = c("response", "probability"),
  ...
)

## S3 method for class 'multinom'
resids(object, ...)

## S3 method for class 'matrix'
resids(object, y, ...)

Arguments

object

An object of class ordinal::clm(), stats::glm(), rms::lrm(), rms::orm(), MASS::polr(), or VGAM::vglm(). Or a matrix of fitted probabilities (e.g., from a multinomial model).
...

Additional optional arguments. (Currently ignored.)
nsim

Integer specifying the number of bootstrap replicates to use. Default is 1L meaning no bootstrap samples.
method

Character string specifying which method to use to generate the surrogate response values. Current options are "latent" and "jitter". Default is "latent".
jitter_scale

Character string specifyint the scale on which to perform the jittering whenever method = "jitter". Current options are "response" and "probability". Default is "response".
y

A vector of the response variable (only used when object is a matrix of probabilities).

Value

A numeric vector of class c("numeric", "surrogate") containing the simulated surrogate response values. Additionally, if nsim > 1, then the result will contain the attributes:

  • boot_reps: A matrix with nsim columns, one for each bootstrap replicate of the surrogate values. Note, these are random and do not correspond to the original ordering of the data.

  • boot_id: A matrix with nsim columns. Each column contains the observation number each surrogate value corresponds to in boot_reps. (This is used for plotting purposes.)

Note

Surrogate response values require sampling from a continuous distribution; consequently, the result will be different with every call to surrogate. The internal functions used for sampling from truncated distributions are based on modified versions of truncdist::rtrunc() and truncdist::qtrunc().

For "glm" objects, only the binomial() family is supported.

References

Liu, Dungang and Zhang, Heping. Residuals and Diagnostics for Ordinal Regression Models: A Surrogate Approach. Journal of the American Statistical Association (accepted). URL http://www.tandfonline.com/doi/abs/10.1080/01621459.2017.1292915?journalCode=uasa20

Nadarajah, Saralees and Kotz, Samuel. R Programs for Truncated Distributions. Journal of Statistical Software, Code Snippet, 16(2), 1-8, 2006. URL https://www.jstatsoft.org/v016/c02.

Examples

# Generate data from a quadratic probit model
set.seed(101)
n <- 2000
x <- runif(n, min = -3, max = 6)
z <- 10 + 3 * x - 1 * x^2 + rnorm(n)
y <- ifelse(z <= 0, yes = 0, no = 1)

# Scatterplot matrix
pairs(~ x + y + z)

# Setup for side-by-side plots
par(mfrow = c(1, 2))

# Misspecified mean structure
fm1 <- glm(y ~ x, family = binomial(link = "probit"))
scatter.smooth(x, y = resids(fm1),
               main = "Misspecified model",
               ylab = "Surrogate residual",
               lpars = list(lwd = 3, col = "red2"))
abline(h = 0, lty = 2, col = "blue2")

# Correctly specified mean structure
fm2 <- glm(y ~ x + I(x ^ 2), family = binomial(link = "probit"))
scatter.smooth(x, y = resids(fm2),
               main = "Correctly specified model",
               ylab = "Surrogate residual",
               lpars = list(lwd = 3, col = "red2"))
abline(h = 0, lty = 2, col = "blue2")

Simulate Data for Diagnostic Plots

Description

Generate simulated data sets designed to illustrate various model misspecifications and diagnostics for ordinal regression models, as described in Liu and Zhang (2017).

Usage

sim_data(
  n = 2000,
  type = c("quadratic", "heteroscedastic", "gumbel", "proportionality", "interaction"),
  ...
)

Arguments

n

Integer specifying the number of observations to simulate. Default is 2000.
type

Character string specifying the type of data/model to simulate. Default is "quadratic" (linear relation with a quadratic trend). Other options include: "heteroscedastic" (non-constant variance), "gumbel" (Gumbel error distribution), "proportionality" (non-proportional hazards), and "interaction" (interaction effect).
...

Additional optional arguments (currently ignored).

Value

A data frame containing the simulated predictor(s) and the ordered factor response variable y.

References

Liu, Dungang and Zhang, Heping. Residuals and Diagnostics for Ordinal Regression Models: A Surrogate Approach. Journal of the American Statistical Association (accepted).

Examples

# Simulate quadratic data
set.seed(101)
df <- sim_data(n = 500, type = "quadratic")
head(df)

Surrogate response

Description

Simulate surrogate response values for cumulative link regression models using the latent method described in Liu and Zhang (2017).

Usage

surrogate(
  object,
  nsim = 1L,
  method = c("latent", "jitter"),
  jitter_scale = c("response", "probability"),
  ...
)

Arguments

object

An object of class ordinal::clm(), stats::glm() rms::lrm(), rms::orm(), MASS::polr(), or VGAM::vglm().
nsim

Integer specifying the number of bootstrap replicates to use. Default is 1L meaning no bootstrap samples.
method

Character string specifying which method to use to generate the surrogate response values. Current options are "latent" and "jitter". Default is "latent".
jitter_scale

Character string specifyint the scale on which to perform the jittering whenever method = "jitter". Current options are "response" and "probability". Default is "response".
...

Additional optional arguments. (Currently ignored.)

Value

A numeric vector of class c("numeric", "surrogate") containing the simulated surrogate response values. Additionally, if nsim > 1, then the result will contain the attributes:

  • boot_reps: A matrix with nsim columns, one for each bootstrap replicate of the surrogate values. Note, these are random and do not correspond to the original ordering of the data.

  • boot_id: A matrix with nsim columns. Each column contains the observation number each surrogate value corresponds to in boot_reps. (This is used for plotting purposes.)

Note

Surrogate response values require sampling from a continuous distribution; consequently, the result will be different with every call to surrogate. The internal functions used for sampling from truncated distributions are based on modified versions of truncdist::rtrunc() and truncdist::qtrunc().

For "glm" objects, only the binomial() family is supported.

References

Liu, Dungang and Zhang, Heping. Residuals and Diagnostics for Ordinal Regression Models: A Surrogate Approach. Journal of the American Statistical Association (accepted). URL http://www.tandfonline.com/doi/abs/10.1080/01621459.2017.1292915?journalCode=uasa20

Nadarajah, Saralees and Kotz, Samuel. R Programs for Truncated Distributions. Journal of Statistical Software, Code Snippet, 16(2), 1-8, 2006. URL https://www.jstatsoft.org/v016/c02.

Examples

# Generate data from a quadratic probit model
set.seed(101)
n <- 2000
x <- runif(n, min = -3, max = 6)
z <- 10 + 3*x - 1*x^2 + rnorm(n)
y <- ifelse(z <= 0, yes = 0, no = 1)

# Scatterplot matrix
pairs(~ x + y + z)

# Setup for side-by-side plots
par(mfrow = c(1, 2))

# Misspecified mean structure
fm1 <- glm(y ~ x, family = binomial(link = "probit"))
s1 <- surrogate(fm1)
scatter.smooth(x, s1 - fm1$linear.predictors,
               main = "Misspecified model",
               ylab = "Surrogate residual",
               lpars = list(lwd = 3, col = "red2"))
abline(h = 0, lty = 2, col = "blue2")

# Correctly specified mean structure
fm2 <- glm(y ~ x + I(x ^ 2), family = binomial(link = "probit"))
s2 <- surrogate(fm2)
scatter.smooth(x, s2 - fm2$linear.predictors,
               main = "Correctly specified model",
               ylab = "Surrogate residual",
               lpars = list(lwd = 3, col = "red2"))
abline(h = 0, lty = 2, col = "blue2")