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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
## QRES.R
qresiduals <- qresid <- function(glm.obj, dispersion=NULL)
# Wrapper function for quantile residuals
# Peter K Dunn
# 28 Sep 2004. Last modified 5 Oct 2004.
{
glm.family <- glm.obj$family$family
if(substr(glm.family,1,17)=="Negative Binomial") glm.family <- "nbinom"
switch(glm.family,
binomial = qres.binom( glm.obj),
poisson = qres.pois(glm.obj),
Gamma = qres.gamma(glm.obj, dispersion),
inverse.gaussian = qres.invgauss(glm.obj, dispersion),
Tweedie = qres.tweedie(glm.obj, dispersion),
nbinom = qres.nbinom(glm.obj),
qres.default(glm.obj, dispersion))
}
qres.binom <- function(glm.obj)
# Randomized quantile residuals for binomial glm
# Gordon Smyth
# 20 Oct 96. Last modified 25 Jan 02.
{
p <- fitted(glm.obj)
y <- glm.obj$y
if(!is.null(glm.obj$prior.weights))
n <- glm.obj$prior.weights
else
n <- rep(1,length(y))
y <- n * y
a <- pbinom(y - 1, n, p)
b <- pbinom(y, n, p)
u <- runif(n = length(y), min = a, max = b)
qnorm(u)
}
qres.pois <- function(glm.obj)
# Quantile residuals for Poisson glm
# Gordon Smyth
# 28 Dec 96
{
y <- glm.obj$y
mu <- fitted(glm.obj)
a <- ppois(y - 1, mu)
b <- ppois(y, mu)
u <- runif(n = length(y), min = a, max = b)
qnorm(u)
}
qres.gamma <- function(glm.obj, dispersion = NULL)
# Quantile residuals for gamma glm
# Gordon Smyth
# 28 Dec 96. Last modified 5 Augusts 2016
{
mu <- fitted(glm.obj)
y <- glm.obj$y
df <- glm.obj$df.residual
w <- glm.obj$prior.weights
if(is.null(w)) w <- 1
if(is.null(dispersion)) dispersion <- sum(w * ((y - mu)/mu)^2)/df
logp <- pgamma((w * y)/mu/dispersion, w/dispersion, log.p=TRUE)
qnorm(logp, log.p=TRUE)
}
qres.invgauss <- function(glm.obj, dispersion = NULL)
# Quantile residuals for inverse Gaussian glm
# Gordon Smyth
# Created 15 Jan 98. Last modified 5 August 2016.
{
mu <- fitted(glm.obj)
y <- glm.obj$y
df <- glm.obj$df.residual
w <- glm.obj$prior.weights
if(is.null(w)) w <- 1
if(is.null(dispersion)) dispersion <- sum(w * (y - mu)^2 / (mu^2*y)) / df
up <- y>mu
down <- y<mu
if(any(down)) y[down] <- qnorm(pinvgauss(y[down],mean=mu[down],dispersion=dispersion,log.p=TRUE),log.p=TRUE)
if(any(up)) y[up] <- qnorm(pinvgauss(y[up],mean=mu[up],dispersion=dispersion,lower.tail=FALSE,log.p=TRUE),lower.tail=FALSE,log.p=TRUE)
y
}
qres.nbinom <- function(glm.obj)
{
# Quantile residuals for Negative Binomial glm
# Gordon Smyth
# 22 Jun 97. Last modified 18 Nov 2008.
#
y <- glm.obj$y
if(is.null(glm.obj$theta)) {
size <- glm.obj$call$family[[2]]
} else {
size <- glm.obj$theta
}
mu <- fitted(glm.obj)
p <- size/(mu + size)
a <- ifelse(y > 0, pbeta(p, size, pmax(y, 1)), 0)
b <- pbeta(p, size, y + 1)
u <- runif(n = length(y), min = a, max = b)
qnorm(u)
}
qres.tweedie <- function(glm.obj, dispersion = NULL)
# Quantile residuals for Tweedie glms
# Gordon Smyth
# Created 29 April 1998. Last modified 30 March 2015.
{
requireNamespace("tweedie")
mu <- fitted(glm.obj)
y <- glm.obj$y
df <- glm.obj$df.residual
w <- glm.obj$prior.weights
if(is.null(w))
w <- 1
p <- get("p",envir=environment(glm.obj$family$variance))
if(is.null(dispersion))
dispersion <- sum((w * (y - mu)^2)/mu^p)/df
u <- tweedie::ptweedie(q=y, power=p, mu=fitted(glm.obj), phi=dispersion/w)
if(p>1&&p<2)
u[y == 0] <- runif(sum(y == 0), min = 0, max = u[y == 0])
qnorm(u)
}
qres.default <- function(glm.obj, dispersion=NULL)
# Quantile residuals for Gaussian and default glms
# Gordon Smyth
# 5 Oct 2004.
{
r <- residuals(glm.obj, type="deviance")
if(is.null(dispersion)) {
df.r <- glm.obj$df.residual
if(df.r > 0) {
if(any(glm.obj$weights==0)) warning("observations with zero weight ", "not used for calculating dispersion")
dispersion <- sum(glm.obj$weights*glm.obj$residuals^2)/df.r
} else
dispersion <- 1
}
r/sqrt(dispersion)
}
|
