-
Notifications
You must be signed in to change notification settings - Fork 0
/
FixEff_dev_MCT.txt
74 lines (62 loc) · 1.31 KB
/
FixEff_dev_MCT.txt
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
model{
# Likelihood:
for (i in 1:50) {
r[i] ~ dbin(p[i], n[i])
# Sampling model:
logit(p[i]) <- mu[s[i]] + d[t[i]] - d[b[i]]
}
# Absolute treatment
# effect sampling
# model:
for (k in 1:4) {
logit(T[k]) <- A + d[k]
}
# Set d[1] to 0:
d[1] <- 0
# Vague Priors on baseline:
for (j in 1:24) {
mu[j] ~ dnorm(0, .0001)
}
# Vague Priors on treatment
# effects:
for (k in 2:4) {
d[k] ~ dnorm(0, .0001)
}
# Prior Absolute treatment
# effects:
A ~ dnorm(-2.6, precA)
# Transformation of
# Absolute treatment
# var to prec:
precA <- pow(.38, -2)
# Log-Odds calculations for
# each comparison
for (c in 1:3) {
# All pair-wise comparison
# log odds ratios:
for (k in (c + 1):4) {
# and single study comparison
# odds ratios:
OR[c, k] <- d[k] - d[c]
log(LOR[c, k]) <- OR[c,k]
}
}
# Rank treatment effect
# (where 1 = best):
rk <- 5 - rank(T[])
# & record the best treatment:
best <- equals(rk, 1)
# Model deviance
# calculations:
for (i in 1:50) {
# Predicted model
# deviance:
rhat[i] <- p[i] * n[i]
# Deviance of each
# data point:
dev[i] <- 2 * (r[i] * (log(r[i]) - log(rhat[i])) +
(n[i] - r[i]) * (log(n[i] - r[i]) - log(n[i] - rhat[i])))
}
# Residual deviance:
resdev <- sum(dev[])
}