forked from gertvv/gemtc
-
Notifications
You must be signed in to change notification settings - Fork 0
/
deviance-example.R
138 lines (117 loc) · 3.73 KB
/
deviance-example.R
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
data <- read.table(textConnection('
study treatment responders sampleSize
1 1 3 39
1 2 3 38
2 1 14 116
2 2 7 114
3 1 11 93
3 2 5 69
4 1 127 1520
4 2 102 1533
5 1 27 365
5 2 28 355
6 1 6 52
6 2 4 59
7 1 152 939
7 2 98 945
8 1 48 471
8 2 60 632
9 1 37 282
9 2 25 278
10 1 188 1921
10 2 138 1916
11 1 52 583
11 2 64 873
12 1 47 266
12 2 45 263
13 1 16 293
13 2 9 291
14 1 45 883
14 2 57 858
15 1 31 147
15 2 25 154
16 1 38 213
16 2 33 207
17 1 12 122
17 2 28 251
18 1 6 154
18 2 8 151
19 1 3 134
19 2 6 174
20 1 40 218
20 2 32 209
21 1 43 364
21 2 27 391
22 1 39 674
22 2 22 680'), header=T)
# Binom/logit (blockers)
#network <- mtc.network(data.ab = data)
#model <- mtc.model(network, linearModel='fixed')
# Binom/logit (smoking)
#network <- read.mtc.network(system.file('extdata/luades-smoking.gemtc', package='gemtc'))
#model <- mtc.model(network, linearModel='fixed')
# Relative effect data
#network <- mtc.network(data.re=read.table('gemtc/tests/data/parkinson-diff.data.txt', header=TRUE))
#model <- mtc.model(network, linearModel='fixed', likelihood='normal', link='identity')
#model <- mtc.model(network, linearModel='random', likelihood='normal', link='identity')
# Mixed data
#model <- mtc.model(parkinson_shared, linearModel='fixed', likelihood='normal', link='identity')
#model <- mtc.model(parkinson_shared, linearModel='random', likelihood='normal', link='identity')
# Poisson/log (dietary fat, 2b)
#network <- mtc.network(dget('gemtc/tests/data/fat-survival.data.txt'))
# ... edit
#model <- mtc.model(network, linearModel='fixed', likelihood='poisson', link='log')
model <- mtc.model(certolizumab, linearModel='fixed', likelihood='binom', link='logit')
result <- mtc.run(model)
summary(result)
print(result$deviance)
## Below not generalized for multi-arm data
# w <- sign(r - rfit) * sqrt(devbar)
# plot(w, leverage, xlim=c(-3,3), ylim=c(0, 4.5))
# x <- seq(from=-3, to=3, by=0.05)
# for (c in 1:4) { lines(x, c - x^2) }
## residual deviance plot
if (length(model$data$studies.a) == length(model$data$studies)) {
tpl <- gemtc:::arm.index.matrix(model[['network']])
study <- matrix(rep(1:nrow(tpl), times=ncol(tpl)), nrow=nrow(tpl), ncol=ncol(tpl))
study <- t(study)[t(!is.na(tpl))]
devbar <- t(result$deviance$dev.ab)[t(!is.na(tpl))]
title <- "Per-arm residual deviance"
xlab <- "Arm"
} else {
nd <- model$data$na
nd[-model$data$studies.a] <- nd[-model$data$studies.a] - 1
devbar <- c(apply(result$deviance$dev.ab, 1, sum, na.rm=TRUE), result$deviance$dev.re) / nd
study <- 1:length(devbar)
title <- "Per-study mean per-datapoint residual deviance"
xlab <- "Study"
}
plot(devbar, ylim=c(0,max(devbar, na.rm=TRUE)),
ylab="Residual deviance", xlab=xlab,
main=title, pch=c(1, 22)[(study%%2)+1])
for (i in 1:length(devbar)) {
lines(c(i, i), c(0, devbar[i]))
}
# w <- sign(r - rfit) * sqrt(devbar)
# plot(w, leverage, xlim=c(-3,3), ylim=c(0, 4.5))
# x <- seq(from=-3, to=3, by=0.05)
# for (c in 1:4) { lines(x, c - x^2) }
fit.ab <- apply(result$deviance$fit.ab, 1, sum, na.rm=TRUE)
dev.ab <- apply(result$deviance$dev.ab, 1, sum, na.rm=TRUE)
lev.ab <- dev.ab - fit.ab
fit.re <- result$deviance$fit.re
dev.re <- result$deviance$dev.re
lev.re <- dev.re - fit.re
nd <- model$data$na
nd[-model$data$studies.a] <- nd[-model$data$studies.a] - 1
w <- sqrt(c(dev.ab, dev.re) / nd)
lev <- c(lev.ab, lev.re) / nd
plot(w, lev, xlim=c(0, max(c(w, 2.5))), ylim=c(0, max(c(lev, 4))),
xlab="Square root of residual deviance", ylab="Leverage",
main="Leverage versus residual deviance")
mtext("Per-study mean per-datapoint contribution")
x <- seq(from=0, to=3, by=0.05)
for (c in 1:4) { lines(x, c - x^2) }
rhat <- gemtc:::arrayize(result$deviance$fitted)
plot(model$data$r/model$data$n, rhat$rhat/model$data$n, xlab="Observed", ylab="Fitted")
abline(a=0,b=1,lty=3)