Computes probability-scale residuals and residual correlations.

Description

This package outputs probability-scale residuals from multiple models and computes residual correlation. Probability-scale residual can be computed for continuous, ordinal, binary, count, and time-to-event data (although the current implementation is only for ordinal variables). Plots of probability-scale residuals can be useful for model diagnostics. Residual correlation can be used to test for conditional independence between multiple types of variables.

Author(s)

Bryan Shepherd bryan.shepherd@vanderbilt.edu

Chun Li cxl791@case.edu

Qi Liu qi.liu4@merck.com

Charles Dupont charles.dupont@vanderbilt.edu

Jeffrey Horner jeffrey.horner@vanderbilt.edu

Goodman-Kruskal's \gamma

Description

Computes Goodman-Kruskal's \gamma

Usage

GKGamma(M)

Arguments

Value

References

Goodman LA, Kruskal WH (1954) Measures of association for cross classifications, Journal of the American Statistical Association, 49, 732-764.

Example Dataset for PResiduals Package

Description

This is a dataset used in Examples Section of PResiduals package help files.

Usage

PResidData

Format

A data frame with 200 rows and 5 variables:

x

an ordered categorical variable with 5 levels

y

an ordered categorical variable with 4 levels

z

a continuous variable

w

a continuous variable

c

a count variable

Source

Simulated

Conditional ordinal by ordinal tests for association.

Description

cobot tests for independence between two ordered categorical variables, X and Y conditional on other variables, Z. The basic approach involves fitting models of X on Z and Y on Z and determining whether there is any remaining information between X and Y. This is done by computing one of 3 test statistics. T1 compares empirical distribution of X and Y with the joint fitted distribution of X and Y under independence conditional on Z. T2 computes the correlation between ordinal (probability-scale) residuals from both models and tests the null of no residual correlation. T3 evaluates the concordance–disconcordance of data drawn from the joint fitted distribution of X and Y under conditional independence with the empirical distribution. Details are given in Li C and Shepherd BE, Test of association between two ordinal variables while adjusting for covariates. Journal of the American Statistical Association 2010, 105:612-620.

Usage

cobot(
  formula,
  link = c("logit", "probit", "cloglog", "loglog", "cauchit"),
  link.x = link,
  link.y = link,
  data,
  subset,
  na.action = na.fail,
  fisher = TRUE,
  conf.int = 0.95
)

Arguments

Details

formula is specified as X | Y ~ Z. This indicates that models of X ~ Z and Y ~ Z will be fit. The null hypothsis to be tested is H_0 : X independant of Y conditional on Z.

Note that T2 can be thought of as an adjusted rank correlation.(Li C and Shepherd BE, A new residual for ordinal outcomes. Biometrika 2012; 99:473-480)

Value

object of ‘⁠cobot⁠’ class.

References

Li C and Shepherd BE, Test of association between two ordinal variables while adjusting for covariates. Journal of the American Statistical Association 2010, 105:612-620.

Li C and Shepherd BE, A new residual for ordinal outcomes. Biometrika 2012; 99:473-480

See Also

Formula, as.data.frame

Examples

## The code is commented out because it will give a Fortran 90 runtime
## error due to the outdated version of Fortran 90 on CRAN's Debian test
## environment (11/2025).
#data(PResidData)
#cobot(x|y~z, data=PResidData)

Conditional continuous by ordinal tests for association.

Description

cocobot tests for independence between an ordered categorical variable, X, and a continuous variable, Y, conditional on other variables, Z. The basic approach involves fitting an ordinal model of X on Z, a linear model of Y on Z, and then determining whether there is any residual information between X and Y. This is done by computing residuals for both models, calculating their correlation, and testing the null of no residual correlation. This procedure is analogous to test statistic T2 in cobot. Two test statistics (correlations) are currently output. The first is the correlation between probability-scale residuals. The second is the correlation between the observed-minus-expected residual for the continuous outcome model and a latent variable residual for the ordinal model (Li C and Shepherd BE, 2012).

Usage

cocobot(
  formula,
  data,
  link = c("logit", "probit", "cloglog", "loglog", "cauchit"),
  subset,
  na.action = getOption("na.action"),
  emp = TRUE,
  fisher = TRUE,
  conf.int = 0.95
)

Arguments

Details

Formula is specified as X | Y ~ Z. This indicates that models of X ~ Z and Y ~ Z will be fit. The null hypothsis to be tested is H_0 : X independant of Y conditional on Z. The ordinal variable, X, must precede the | and be a factor variable, and Y must be continuous.

Value

object of ‘⁠cocobot⁠’ class.

References

Li C and Shepherd BE (2012) A new residual for ordinal outcomes. Biometrika. 99: 473–480.

Shepherd BE, Li C, Liu Q (2016) Probability-scale residuals for continuous, discrete, and censored data. The Canadian Journal of Statistics. 44: 463–479.

Examples

data(PResidData)
cocobot(y|w ~ z, data=PResidData)

Conditional Partial Spearman's Rank Correlation

Description

conditional_Spearman computes the partial Spearman's rank correlation between variable X and variable Y adjusting for variable Z conditional on Zc. X and Y can be any orderable variables, including continuous and discrete variables. Covariate Z can be multidimensional. X, Y, and Z are specified by the argument ‘⁠formula⁠’. Zc is a one-dimensional covariate, specified by the argument ‘⁠conditional.by⁠’. The basic approach involves fitting a specified model of X on Z, a specified model of Y on Z, obtaining the probability-scale residuals, Xres and Yres, from both models, and then modeling their Pearson's correlation conditional on Zc. Different methods are provided to model the Pearson's correlation between the two sets of probability-scale residuals. See details in ‘⁠conditional.method⁠’. As in ‘⁠partial.Spearman⁠’, by default conditional_Spearman uses cumulative link models for both continous and discrete ordinal variables X and Y to preserve the rank-based nature of Spearman's correlation. For some specific types of variables, options of fitting parametric models are also available. See details in ‘⁠fit.x⁠’ and ‘⁠fit.y⁠’.

Usage

conditional_Spearman(
  formula,
  conditional.by,
  data,
  conditional.method = c("lm", "kernel", "stratification"),
  conditional.formula = paste("~", conditional.by, sep = ""),
  kernel.function = c("normal", "gaussian", "triweight", "quartic", "biweight",
    "epanechnikov", "uniform", "triangle"),
  kernel.bandwidth = "silverman",
  fit.x = "orm",
  fit.y = "orm",
  link.x = c("logit", "probit", "cloglog", "loglog", "cauchit", "logistic"),
  link.y = c("logit", "probit", "cloglog", "loglog", "cauchit", "logistic"),
  subset,
  na.action = getOption("na.action"),
  fisher = TRUE,
  conf.int = 0.95
)

Arguments

Details

To compute the partial Spearman's rank correlation between X and Y adjusting for Z conditional on Zc, ‘⁠formula⁠’ is specified as X | Y ~ Z and ‘⁠conditional.by⁠’ is specified as Zc. This indicates that models of X ~ Z and Y ~ Z will be fit, and the correlation between the probability-scale residuals from these two models will be modeled conditional on Zc.

Value

object of ‘⁠conditional_Spearman⁠’ class.

References

Li C and Shepherd BE (2012) A new residual for ordinal outcomes. Biometrika. 99: 473–480.

Shepherd BE, Li C, Liu Q (2016) Probability-scale residuals for continuous, discrete, and censored data. The Canadian Jouranl of Statistics. 44:463–476.

Liu Q, Shepherd BE, Wanga V, Li C (2018) Covariate-Adjusted Spearman's Rank Correlation with Probability-Scale Residuals. Biometrics. 74:595–605.

See Also

print.conditional_Spearman,print.conditional_Spearman

Examples

data(PResidData)
library(rms)
#### fitting cumulative link models for both Y and W
result <- conditional_Spearman(c|y~ x + w, conditional.by="w",
                                       conditional.method="lm", conditional.formula="~rcs(w)",
                                       fit.x="poisson",fit.y="orm",
                                       data=PResidData, fisher=TRUE)
plot(result)

Calculates the weighted correlation given a data set and a set of weights.

Description

This is a copy of corr function from the boot package. It calculates the correlation coefficient in weighted form.

Usage

corr(d, w = rep(1, nrow(d))/nrow(d))

Arguments

Value

the correlation coefficient between d[,1] and d[,2].

Conditional count by ordinal tests for association.

Description

countbot tests for independence between an ordered categorical variable, X, and a count variable, Y, conditional on other variables, Z. The basic approach involves fitting an ordinal model of X on Z, a Poisson or Negative Binomial model of Y on Z, and then determining whether there is any residual information between X and Y. This is done by computing residuals for both models, calculating their correlation, and testing the null of no residual correlation. This procedure is analogous to test statistic T2 in cobot. Two test statistics (correlations) are currently output. The first is the correlation between probability-scale residuals. The second is the correlation between the Pearson residual for the count outcome model and a latent variable residual for the ordinal model (Li C and Shepherd BE, 2012).

Usage

countbot(
  formula,
  data,
  link.x = c("logit", "probit", "loglog", "cloglog", "cauchit"),
  fit.y = c("poisson", "negative binomial"),
  subset,
  na.action = getOption("na.action"),
  fisher = TRUE,
  conf.int = 0.95
)

Arguments

Details

Formula is specified as X | Y ~ Z. This indicates that models of X ~ Z and Y ~ Z will be fit. The null hypothesis to be tested is H_0 : X independent of Y conditional on Z. The ordinal variable, X, must precede the | and be a factor variable, and Y must be an integer.

Value

object of ‘⁠cocobot⁠’ class.

References

Li C and Shepherd BE (2012) A new residual for ordinal outcomes. Biometrika. 99: 473–480.

Shepherd BE, Li C, Liu Q (2016) Probability-scale residuals for continuous, discrete, and censored data. The Canadian Journal of Statistics. 44: 463–479.

Examples

data(PResidData)
countbot(x|c ~z, fit.y="poisson",data=PResidData)
countbot(x|c ~z, fit.y="negative binomial",data=PResidData)

Extract or construct a diagonal matrix.

Description

This works like diag except when x is a single integer value. If x is a single integer value then it assumes that you want a 1 by 1 matrix with the value set to x

Usage

diagn(x = 1, nrow = length(x), ncol = nrow)

Arguments

Value

matrix with diagonal elements set to x

See Also

diag

Examples

diag(5)
diagn(5)

kernel.function

Description

kernel.function calculates several kernel functions (uniform, triangle, epanechnikov, biweight, triweight, gaussian).

Usage

kernel.function(u, kernel = "normal", product = TRUE)

Arguments

Details

slightly modified version of the kernel.function from the gplm package. The kernel parameter is a text string specifying the univariate kernel function which is either the gaussian pdf or proportional to (1-|u|^p)^q. Possible text strings are "triangle" (p=q=1), "uniform" (p=1, q=0), "epanechnikov" (p=2, q=1), "biweight" or "quartic" (p=q=2), "triweight" (p=2, q=3), "gaussian" or "normal" (gaussian pdf). The multivariate kernels are obtained by a product of unvariate kernels K(u_1)...K(u_d) or by a spherical (radially symmetric) kernel proportional to K(||u||). (The resulting kernel is a density, i.e. integrates to 1.)

Value

matrix with diagonal elements set to x

Conditional tests for association.

Description

megabot tests for correlation between a variable, X, and another variable, Y, conditional on other variables, Z. The basic approach involves fitting an specified model of X on Z, a specified model of Y on Z, and then determining whether there is any remaining information between X and Y. This is done by computing residuals for both models, calculating their correlation, and testing the null of no residual correlation. The test statistic output is the correlation between probability-scale residuals. X and Y can be continous or ordered discrete variables. megabot replicates the functionality of cobot, cocobot, and countbot

Usage

megabot(
  formula,
  data,
  fit.x,
  fit.y,
  link.x = c("logit", "probit", "cloglog", "loglog", "cauchit", "logistic"),
  link.y = c("logit", "probit", "cloglog", "loglog", "cauchit", "logistic"),
  subset,
  na.action = getOption("na.action"),
  fisher = TRUE,
  conf.int = 0.95
)

Arguments

Details

Formula is specified as X | Y ~ Z. This indicates that models of X ~ Z and Y ~ Z will be fit. The null hypothesis to be tested is H_0 : X independent of Y conditional on Z.

Value

object of ‘⁠cocobot⁠’ class.

References

Li C and Shepherd BE (2012) A new residual for ordinal outcomes. Biometrika. 99: 473–480.

Shepherd BE, Li C, Liu Q (2016) Probability-scale residuals for continuous, discrete, and censored data. The Canadian Journal of Statistics. 44: 463–479.

Examples

data(PResidData)
megabot(y|w ~ z, fit.x="ordinal", fit.y="lm.emp", data=PResidData)

slightly modified version of polr from MASS

Description

slightly modified version of polr from MASS

Usage

newpolr(
  formula,
  data,
  weights,
  start,
  ...,
  subset,
  na.action,
  contrasts = NULL,
  Hess = FALSE,
  model = TRUE,
  method = c("logit", "probit", "cloglog", "loglog", "cauchit")
)

Arguments

Value

A object of class "polr". This has components

References

polr from MASS

See Also

optim, glm, multinom

Partial Spearman's Rank Correlation

Description

partial_Spearman computes the partial Spearman's rank correlation between variable X and variable Y adjusting for other variables, Z. The basic approach involves fitting a specified model of X on Z, a specified model of Y on Z, obtaining the probability-scale residuals from both models, and then calculating their Pearson's correlation. X and Y can be any orderable variables, including continuous or discrete variables. By default, partial_Spearman uses cumulative probability models (also referred as cumulative link models in literature) for both X on Z and Y on Z to preserve the rank-based nature of Spearman's correlation, since the model fit of cumulative probability models only depends on the order information of variables. However, for some specific types of variables, options of fitting parametric models are also available. See details in fit.x and fit.y

Usage

partial_Spearman(
  formula,
  data,
  fit.x = "orm",
  fit.y = "orm",
  link.x = c("logit", "probit", "cloglog", "loglog", "cauchit", "logistic"),
  link.y = c("logit", "probit", "cloglog", "loglog", "cauchit", "logistic"),
  subset,
  na.action = getOption("na.action"),
  fisher = TRUE,
  conf.int = 0.95
)

Arguments

Details

To compute the partial Spearman's rank correlation between X and Y adjusting for Z, ‘⁠formula⁠’ is specified as X | Y ~ Z. This indicates that models of X ~ Z and Y ~ Z will be fit.

Value

object of ‘⁠partial_Spearman⁠’ class.

References

Li C and Shepherd BE (2012) A new residual for ordinal outcomes. Biometrika. 99: 473–480.

Shepherd BE, Li C, Liu Q (2016) Probability-scale residuals for continuous, discrete, and censored data. The Canadian Jouranl of Statistics. 44:463–476.

Liu Q, Shepherd BE, Wanga V, Li C (2018) Covariate-Adjusted Spearman's Rank Correlation with Probability-Scale Residuals. Biometrics. 74:595–605.

See Also

print.partial_Spearman

Examples

data(PResidData)
library(rms)
#### fitting cumulative probability models for both Y and W
partial_Spearman(c|w ~ z,data=PResidData)
#### fitting a cumulative probability model for W and a poisson model for c
partial_Spearman(c|w~z, fit.x="poisson",data=PResidData)
partial_Spearman(c|w~z, fit.x="poisson", fit.y="lm.emp", data=PResidData )

conditional_Spearman class plot method

Description

conditional_Spearman class plot method

Usage

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

Arguments

Value

No return value, called for side effects

Probability-scale residual

Description

presid calculates the probability-scale residual for various model function objects. Currently supported models include glm (Poisson, binomial, and gaussian families), lm in the stats library; survreg (Weibull, exponential, gaussian, logistic, and lognormal distributions) and coxph in the survival library; polr and glm.nb in the MASS library; and ols, cph, lrm, orm, psm, and Glm in the rms library.

Usage

presid(object, ...)

Arguments

Details

Probability-scale residual is P(Y < y) - P(Y > y) where y is the observed outcome and Y is a random variable from the fitted distribution.

Value

The probability-scale residual for the model

References

Shepherd BE, Li C, Liu Q (2016) Probability-scale residuals for continuous, discrete, and censored data. The Canadian Jouranl of Statistics. 44:463–476.

Li C and Shepherd BE (2012) A new residual for ordinal outcomes. Biometrika. 99: 473–480.

Examples

library(survival)
library(stats)

set.seed(100)
n <- 1000
x <- rnorm(n)
t <- rweibull(n, shape=1/3, scale=exp(x))
c <- rexp(n, 1/3)
y <- pmin(t, c)
d <- ifelse(t<=c, 1, 0)

mod.survreg <- survreg(Surv(y, d) ~ x, dist="weibull")
summary(presid(mod.survreg))
plot(x, presid(mod.survreg))

##### example for proprotional hazards model
n <- 1000
x <- rnorm(n)
beta0 <- 1
beta1 <- 0.5
t <- rexp(n, rate = exp(beta0 + beta1*x))
c <- rexp(n, rate=1)
y <- ifelse(t<c, t, c)
delta <- as.integer(t<c)

mod.coxph <- coxph(Surv(y, delta) ~ x)
presid <- presid(mod.coxph)
plot(x, presid, cex=0.4, col=delta+2)

#### example for Negative Binomial regression
library(MASS)

n <- 1000
beta0 <- 1
beta1 <- 0.5
x <- runif(n, min=-3, max=3)
y <- rnbinom(n, mu=exp(beta0 + beta1*x), size=3)

mod.glm.nb <- glm.nb(y~x)
presid <- presid(mod.glm.nb)
summary(presid)
plot(x, presid, cex=0.4)

##### example for proportional odds model
library(MASS)

n <- 1000
x  <- rnorm(n)
y  <- numeric(n)
alpha = c(-1, 0, 1, 2)
beta <- 1
py  <-  (1 + exp(- outer(alpha, beta*x, "+"))) ^ (-1)
aa = runif(n)
for(i in 1:n)
  y[i] = sum(aa[i] > py[,i])
y <-  as.factor(y)


mod.polr <- polr(y~x, method="logistic")
summary(mod.polr)
presid <- presid(mod.polr)
summary(presid)
plot(x, presid, cex=0.4)

cobot class print method

Description

cobot class print method

Usage

## S3 method for class 'cobot'
print(x, ...)

Arguments

Value

No return value, called for side effects

cocobot class print method

Description

cocobot class print method

Usage

## S3 method for class 'cocobot'
print(x, ...)

Arguments

Value

No return value, called for side effects

conditional_Spearman class print method

Description

conditional_Spearman class print method

Usage

## S3 method for class 'conditional_Spearman'
print(x, ...)

Arguments

Value

No return value, called for side effects

partial_Spearman class print method

Description

partial_Spearman class print method

Usage

## S3 method for class 'partial_Spearman'
print(x, ...)

Arguments

Value

No return value, called for side effects