-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathbayes.test.R
75 lines (75 loc) · 3.5 KB
/
bayes.test.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
function (rhos, marg, res.se = 1, n, singular, surv = 0.975)
{
prior_jeffrey = dget("prior_jeffrey.R")
Berger = dget("berger.R")
lubrano = dget("lubrano.R")
(beta = dbeta(rhos, shape1 = 10, shape2 = 2))
odds = BF1 = BF2 = BF3 = BF4 = postH0 = point.est = NULL
instat = which(rhos >= surv)
(prior = prior_jeffrey(rhos = rhos, TT = n, sigma = res.se))
(output1 = prior * marg)
(const = sintegral(x = rhos, output1)$value)
(output1 = output1/const)
(point1 = sintegral(x = rhos, fx = rhos * output1)$value)
(numJ = sintegral(rhos[instat], output1[instat])$value)
(prior = dnorm(rhos, mean = 1, sd = res.se))
(output2 = marg * prior)
(const = sintegral(x = rhos, output2)$value)
(output2 = output2/const)
(point2 = sintegral(x = rhos, fx = rhos * output2)$value)
(numN = sintegral(rhos[instat], output2[instat])$value)
(prior = Berger(rhos, singular))
(output3 = marg * prior)
(const = sintegral(x = rhos, output3)$value)
(output3 = output3/const)
(point3 = sintegral(x = rhos, fx = rhos * output3)$value)
(numB = sintegral(rhos[instat], output3[instat])$value)
(prior = lubrano(rho = rhos, nu = 0.333))
(output4 = marg * prior)
(const = sintegral(x = rhos, output4)$value)
(output4 = output4/const)
(point4 = sintegral(x = rhos, fx = rhos * output4)$value)
(numL = sintegral(rhos[instat], output4[instat])$value)
stat_H1 = which(rhos >= 0 & rhos < 1)
(denomJ = sintegral(rhos[stat_H1], output1[stat_H1])$value)
(denomN = sintegral(rhos[stat_H1], output2[stat_H1])$value)
(denomB = sintegral(rhos[stat_H1], output3[stat_H1])$value)
(denomL = sintegral(rhos[stat_H1], output4[stat_H1])$value)
(PO_J = numJ/denomJ)
(PO_L = numL/denomL)
(PO_N = numN/denomN)
(PO_B = numB/denomB)
odds = rbind(odds, c(Ratios.Jeff = PO_J, Ratios.Norm = PO_N,
Ratios.BY = PO_B, Ratios.LU = PO_L))
(const = sintegral(x = rhos, marg)$value)
mlike.norm = marg/const
(BFnum = mlike.norm[which(rhos == 1)])
pp = 0.5
(BFdenom = sintegral(rhos[stat_H1], output1[stat_H1])$value)
(BFJ = BFnum/BFdenom)
(pH0_J = 1/(1 + (1 - pp)/pp * (1/BFJ)))
BFJ.trans = round(-2 * log(BFJ), digits = 4)
(BF1 = rbind(BF1, c(BFJ, BFJ.trans, H0rejectJ = BFJ < 1)))
(BFdenom = sintegral(rhos[stat_H1], output2[stat_H1])$value)
(BFN = BFnum/BFdenom)
(pH0_N = 1/(1 + (1 - pp)/pp * (1/BFN)))
BFN.trans = round(-2 * log(BFN), digits = 4)
(BF3 = rbind(BF3, c(BFN, BFN.trans, H0rejectN = BFN < 1)))
(BFdenom = sintegral(rhos[stat_H1], output3[stat_H1])$value)
(BFB = BFnum/BFdenom)
(pH0_BY = 1/(1 + (1 - pp)/pp * (1/BFB)))
BFB.trans = round(-2 * log(BFB), digits = 4)
(BF4 = rbind(BF4, c(BFB, BFB.trans, H0rejectB = BFB < 1)))
(BFdenom = sintegral(rhos[stat_H1], output4[stat_H1])$value)
(BFL = BFnum/BFdenom)
(pH0_LU = 1/(1 + (1 - pp)/pp * (1/BFL)))
BFL.trans = round(-2 * log(BFL), digits = 4)
(BF2 = rbind(BF2, c(BFL, BFL.trans, H0rejectL = BFL < 1)))
(postH0 = rbind(postH0, c(pH0_N, pH0_BY, pH0_J, pH0_LU)))
(point.est = rbind(point.est, c(Jeff = point1, Normal = point2,
BY = point3, LU = point4)))
(probs = t(data.frame(postH0.N = postH0[, 1], postH0.BY = postH0[,
2], postH0.J = postH0[, 3], postH0.LU = postH0[, 4],
pt.est = point.est, probJ = numJ, probN = numN)))
return(list(test_results = probs, Jeff.mpost = output1, Normal.mpost = output2))
}