-
-
Notifications
You must be signed in to change notification settings - Fork 478
/
magnesium.stan
88 lines (76 loc) · 2.17 KB
/
magnesium.stan
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
// Sensitivity to prior distributions: application to Magnesium meta-analysis
// http://www.openbugs.net/Examples/Magnesium.html
// prior specification 1, 2, 3, 4, 5, and 6 combined
// in one posterior distribution.
// For individual model specification, see magnesium$i.stan ($i=1,...,6)
//
// Note: tau is really the standard deviation parameter in the normal distribution.
//
data {
int N_studies;
array[N_studies] int rt;
array[N_studies] int nt;
array[N_studies] int rc;
array[N_studies] int nc;
}
transformed data {
int N_priors;
real<lower=0> s0_sq;
real<lower=0> p0_sigma;
N_priors = 6;
s0_sq = 0.1272041;
p0_sigma = 1 / sqrt(Phi(0.75) / s0_sq);
}
parameters {
array[N_priors] real<lower=-10, upper=10> mu;
array[N_priors, N_studies] real theta;
array[N_priors, N_studies] real<lower=0, upper=1> pc;
real<lower=0> inv_tau_sq_1;
real<lower=0, upper=50> tau_sq_2;
real<lower=0, upper=50> tau_3;
real<lower=0, upper=1> B0;
real<lower=0, upper=1> D0;
real<lower=0> tau_sq_6;
}
transformed parameters {
array[N_priors] real<lower=0> tau;
tau[1] = 1 / sqrt(inv_tau_sq_1);
tau[2] = sqrt(tau_sq_2);
tau[3] = tau_3;
tau[4] = sqrt(s0_sq * (1 - B0) / B0);
tau[5] = sqrt(s0_sq) * (1 - D0) / D0;
tau[6] = sqrt(tau_sq_6);
}
model {
// prior 1: gamma(0.001, 0.001) on inv_tau_sq
inv_tau_sq_1 ~ gamma(0.001, 0.001);
// Prior 2: Uniform(0, 50) on tau.sqrd
tau_sq_2 ~ uniform(0, 50);
// Prior 3: Uniform(0, 50) on tau
tau_3 ~ uniform(0, 50);
// Prior 4: Uniform shrinkage on tau.sqrd
B0 ~ uniform(0, 1);
// Prior 5: Dumouchel on tau
D0 ~ uniform(0, 1);
// Prior 6: Half-Normal on tau.sqrd
tau_sq_6 ~ normal(0, p0_sigma) T[0, ];
mu ~ uniform(-10, 10);
for (prior in 1 : N_priors) {
pc[prior] ~ uniform(0, 1);
theta[prior] ~ normal(mu[prior], tau[prior]);
}
for (prior in 1 : N_priors) {
vector[N_studies] tmpm;
for (i in 1 : N_studies) {
tmpm[i] = theta[prior, i] + logit(pc[prior, i]);
}
rc ~ binomial(nc, pc[prior]);
rt ~ binomial_logit(nt, tmpm);
}
}
generated quantities {
array[N_priors] real OR;
for (prior in 1 : N_priors) {
OR[prior] = exp(mu[prior]);
}
}