-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathptII_quan_Bayes_case_presidential-heights_helpfuncs.r
98 lines (87 loc) · 3.58 KB
/
ptII_quan_Bayes_case_presidential-heights_helpfuncs.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
### (C) 2005-2023 by Leo Guertler
### R-code supplement
### to the book
###
### "Subjektive Ansichten und objektive Betrachtungen"
###
### written by Gürtler & Huber (2023)
###
### All R-code is published under the GPL v3 license:
###
### https://www.gnu.org/licenses/gpl-3.0.en.html
###
### except for 'borrowed' code - see links and references.
### For this R-code the original license of the respective
### authors is valid.
###
### R-code published on
###
### https://osdn.net/projects/mixedmethod-rcode
### https://github.com/abcnorio/mixedmethod-rcode
# file:
# ptII_quan_Bayes_case_presidential-heights_helpfuncs.r
# location:
# chap. 6 [6.15.1]
# Präsidiale Höhenflüge
# HELPER FUNCTIONS
###### alternative function to bayes.binom.test from BayesianFirstAid
bayes.binom.alt <- function(x=NA, n=NA, n.chains=3, n.iter=15000, cred.mass=0.87, comp.theta=0.5,
x_name="", n_name="",
model_string=NA,
printout=TRUE, plotout=TRUE, prior.a=1, prior.b=1)
{
require(BayesianFirstAid)
if(is.na(model_string)) bbinom.model.string <- paste("model {\n x ~ dbinom(theta, n)\n theta ~ dbeta(",prior.a,", ",prior.b,")\n x_pred ~ dbinom(theta, n)\n}",sep="")
DNAME <- x_name
mcmc_samples <- BayesianFirstAid:::run_jags(
model_string=bbinom.model.string,
data=list(x=x, n=n), inits=list(theta = (n + 1) / (n + 2)),
params=c("theta", "x_pred"), n.chains=n.chains, n.adapt=0,
n.update=0, n.iter=n.iter, thin=1, progress.bar="none"
)
mcmc_samples
stats <- BayesianFirstAid:::mcmc_stats(mcmc_samples, cred_mass = cred.mass, comp_val = comp.theta)
bfa_object <- list(x = x, n = n, comp_theta = comp.theta, cred_mass = cred.mass,
x_name = x_name, n_name = n_name, data_name = DNAME,
mcmc_samples = mcmc_samples, stats = stats,
model=bbinom.model.string)
class(bfa_object) <- c("bayes_binom_test", "bayesian_first_aid")
if(printout) print(summary(bfa_object))
if(plotout) plot(bfa_object)
return(bfa_object)
}
# call:
# h.bbt.res <- bayes.binom.alt(x=x, n=n, cred.mass=0.87, comp.theta=0.5, x_name=x_name, n_name=n_name, model_string=bbinom.model.string)
# h.bbt.res <- bayes.binom.alt(x=x, n=n, cred.mass=0.87, comp.theta=0.5, x_name=x_name, n_name=n_name)
# h.bbt.res <- bayes.binom.alt(x=x, n=n, cred.mass=0.87, comp.theta=0.5, x_name=x_name, n_name=n_name, prior.a=2, prior.b=45)
########################## END OF FUNCTION
###### function to calculate posterior (Odds) Ratios
post.comp <- function(post1, post2, digits=3, RETURN=FALSE)
{
RR <- p1/p2
RR.inv <- 1-RR
OR <- (post1 * (1-post2)) / (post2 * (1-post1))
OR.inv <- 1-OR
res <- structure(list("post[1]"=post1, "post[2]"=post2,
"Ratio (p1 vs. p2)"=RR,
"Ratio (p2 vs. p1)"=1/RR,
"Odds Ratio (p1 vs. p2)"=OR,
"Odds Ratio (p2 vs. p1)"=1/OR,
WHATIS="Comparison of two posterior values", note=NULL), class="post.comp")
# taken from stats:::print.power.htest()
cat("\n ", res$WHATIS, "\n\n")
note <- res$note
res[c("WHATIS", "note")] <- NULL
cat(paste(format(names(res), width = 15L, justify = "right"),
format(res, digits = digits), sep = " = "), sep = "\n")
if(!is.null(note))
{
cat("\n", "NOTE: ", note, "\n\n", sep = "")
} else cat("\n")
if(RETURN) return(res) else invisible()
}
# call:
# p1 <- mean(meanDiff.diff > 0)
# p2 <- mean(meanDiff.diff < 0)
# p1vsp2 <- post.comp(post1=p1 , post2=p2, RETURN=TRUE)
########################## END OF FUNCTION