presid: Probability-scale residual

In PResiduals: Probability-Scale Residuals and Residual Correlations

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

object	The model object for which the probability-scale residual is calculated
...	Additional arguements passed to methods

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)

Example output

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.3491  0.2995  0.5296  0.5138  0.7393  1.0000 
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.983745 -0.497222 -0.042775 -0.004053  0.461956  0.999351 

Re-fitting to get Hessian

Call:
polr(formula = y ~ x, method = "logistic")

Coefficients:
   Value Std. Error t value
x -1.046    0.06776  -15.44

Intercepts:
    Value    Std. Error t value 
0|1  -0.9866   0.0752   -13.1210
1|2  -0.0838   0.0686    -1.2213
2|3   0.9422   0.0747    12.6071
3|4   2.0529   0.0977    21.0099

Residual Deviance: 2868.486 
AIC: 2878.486 
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.971660 -0.504768 -0.008136  0.000000  0.489763  0.979312 

PResiduals documentation built on June 24, 2021, 9:10 a.m.