presid: Probability-scale residual

Description Usage Arguments Details Value References Examples

View source: R/presid.R

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

1
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. Probability-scale residuals for continuous, discrete, and censored data. Submitted.

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

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
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)

PResiduals documentation built on Oct. 6, 2017, 5:07 p.m.