The null hypothesis is that the infection is not endemic (prevalence <= 50%).
dat <- read.table(head = TRUE, text = "
theta hypothesis
0.0 H0
0.1 H0
0.2 H0
0.3 H0
0.4 H0
0.5 H0
0.6 H1
0.7 H1
0.8 H1
0.9 H1
1.0 H1
")
Prior probabilities
The investigator is 80% sure it is endemic, thus 0.80 is distributed uniformly over the prevalence parameter values consistent with \( H_1 \). Similarly, 0.20 is distributed over the prevalence parameter values consistent with \( H_0 \).
dat$prior <- NA
dat$prior[1:6] <- 0.20 / 6
dat$prior[7:11] <- 0.80 / 5
dat
theta hypothesis prior
1 0.0 H0 0.03333
2 0.1 H0 0.03333
3 0.2 H0 0.03333
4 0.3 H0 0.03333
5 0.4 H0 0.03333
6 0.5 H0 0.03333
7 0.6 H1 0.16000
8 0.7 H1 0.16000
9 0.8 H1 0.16000
10 0.9 H1 0.16000
11 1.0 H1 0.16000
barplot(dat$prior, name = dat$theta, ylim = c(0,1), main = "Prior")
Likelihood
The likelihood function is binomial distribution. The data is fixed, thus the number of positive results is 7 out of 10. The value of theta (parameter) is changed over values between 0.0 to 1.0.
dat$likelihood <- dbinom(x = 7, size = 10, prob = dat$theta)
dat
theta hypothesis prior likelihood
1 0.0 H0 0.03333 0.000000000
2 0.1 H0 0.03333 0.000008748
3 0.2 H0 0.03333 0.000786432
4 0.3 H0 0.03333 0.009001692
5 0.4 H0 0.03333 0.042467328
6 0.5 H0 0.03333 0.117187500
7 0.6 H1 0.16000 0.214990848
8 0.7 H1 0.16000 0.266827932
9 0.8 H1 0.16000 0.201326592
10 0.9 H1 0.16000 0.057395628
11 1.0 H1 0.16000 0.000000000
barplot(dat$likelihood, name = dat$theta, ylim = c(0,1), main = "Likelihood")
Posterior probabilities
dat <- within(dat, {
likelihood.x.prior <- likelihood * prior
})
dat$posterior <- prop.table(dat$likelihood.x.prior)
dat
theta hypothesis prior likelihood likelihood.x.prior posterior
1 0.0 H0 0.03333 0.000000000 0.0000000000 0.000000000
2 0.1 H0 0.03333 0.000008748 0.0000002916 0.000002349
3 0.2 H0 0.03333 0.000786432 0.0000262144 0.000211177
4 0.3 H0 0.03333 0.009001692 0.0003000564 0.002417179
5 0.4 H0 0.03333 0.042467328 0.0014155776 0.011403538
6 0.5 H0 0.03333 0.117187500 0.0039062500 0.031467770
7 0.6 H1 0.16000 0.214990848 0.0343985357 0.277105970
8 0.7 H1 0.16000 0.266827932 0.0426924691 0.343919816
9 0.8 H1 0.16000 0.201326592 0.0322122547 0.259493839
10 0.9 H1 0.16000 0.057395628 0.0091833005 0.073978364
11 1.0 H1 0.16000 0.000000000 0.0000000000 0.000000000
barplot(dat$posterior, name = dat$theta, ylim = c(0,1), main = "Posterior")
Uninformative prior
The posterior distribution has the same shape as the likelihood function, i.e., the same result as the frequentist approach.
dat$prior <- 1 / length(dat$prior)
dat <- within(dat, {
likelihood.x.prior <- likelihood * prior
})
dat$posterior <- prop.table(dat$likelihood.x.prior)
dat
theta hypothesis prior likelihood likelihood.x.prior posterior
1 0.0 H0 0.09091 0.000000000 0.0000000000 0.000000000
2 0.1 H0 0.09091 0.000008748 0.0000007953 0.000009613
3 0.2 H0 0.09091 0.000786432 0.0000714938 0.000864218
4 0.3 H0 0.09091 0.009001692 0.0008183356 0.009892049
5 0.4 H0 0.09091 0.042467328 0.0038606662 0.046667768
6 0.5 H0 0.09091 0.117187500 0.0106534091 0.128778506
7 0.6 H1 0.09091 0.214990848 0.0195446225 0.236255574
8 0.7 H1 0.09091 0.266827932 0.0242570847 0.293219860
9 0.8 H1 0.09091 0.201326592 0.0183024175 0.221239788
10 0.9 H1 0.09091 0.057395628 0.0052177844 0.063072625
11 1.0 H1 0.09091 0.000000000 0.0000000000 0.000000000
barplot(dat$prior, name = dat$theta, ylim = c(0,1), main = "Prior")
barplot(dat$likelihood, name = dat$theta, ylim = c(0,1), main = "Likelihood")
barplot(dat$posterior, name = dat$theta, ylim = c(0,1), main = "Posterior")
