Lab 1

Analisi Bayesiana di Dati Binomiali

library(dplyr)  # optional
library(kableExtra)

Esempio

  1. Stima della probabilità in un modello binomiale e prima analisi di sensitività
  2. idem ma con probabilità prossima a 0 e confrontando un campione numeroso con uno piccolo

Riferimento

  1. Gelman: sec. 2.4, 37-39
  2. Congdon: ex. 3.4, 77, Program 2.11

Per installare BUGS:
BUGS resources on-line (url: https://www.mrc-bsu.cam.ac.uk/software/bugs-project)
From there, go to links to WinBUGS or OpenBUGS.

Linguaggio: R

Soggetto

  1. La proporzione di nascite femminili nella popolazione di nati con placenta previa è inferiore a 0.485, la proporzione di nascite femminili nella popolazione generale?
  2. Qual è la probabilità che un adulto campionato casualmente risponda che le azioni del Presidente erano immorali?

Placenta previa

A priori uniforme

n <- 980
y <- 437
m <- 0.485
alpha <- beta <- 1
a <- alpha + y
b <- beta + n - y

theta <- seq(0.375, 0.525, 0.001)

# plot set
par(mar = c(4, 4, 1, 1))
plot(theta, dbeta(theta, a, b), type = "n", bty = "l", xlab = expression(theta), ylab = "density")

# likelihood.  It has been scaled (per 1000) to be comparable with posterior at
# alpha+beta=2.

lines(theta, 1000 * choose(n, y) * theta^y * (1 - theta)^(n - y), lwd = 2, col = "orange")
text(0.475, 1000 * choose(n, y) * 0.45^y * (1 - 0.45)^(n - y), label = "likelihood", col = "darkorange")

# prior

lines(theta, dbeta(theta, alpha, beta), lwd = 2, col = "aquamarine")
text(0.5, 2.25, label = "prior", col = "aquamarine3")

# posterior

lines(theta, dbeta(theta, a, b), lwd = 2, col = "dodgerblue")
text(0.48, dbeta(0.46, a, b), label = "posterior", col = "dodgerblue")

# posterior inference

ci <- qbeta(c(0.025, 0.975), a, b)

ind_ci <- theta >= ci[1] & theta <= ci[2]
polygon(c(ci[1], ci[1], theta[ind_ci], ci[2], ci[2]), c(0, dbeta(c(ci[1], theta[ind_ci],
    ci[2]), a, b), 0), col = "lightblue", border = F)

abline(v = m, lty = 2)

A priori Beta

Diverse a priori Beta sempre più concentrate intorno alla media fissata a 0.485, che è la media nota della popolazione generale (nascita femminile).

Dopo aver fissato diversi valori di \(\alpha + \beta\) (da cui i valori di \(\alpha\) e \(\beta\) risultano di conseguenza data la media assegnata a priori), il codice disegna, per ogni caso, l’a posteriori, la verosimiglianza e l’a priori.

Per valutare la sensibilità dell’inferenza a posteriori rispetto alla a priori assunta, gli intervalli a posteriori del 95% di \(\theta\) e la media della popolazione generale vengono tracciati.

n <- 980
y <- 437
m <- 0.485
theta <- seq(0.35, 0.55, 0.001)
posum <- c()

placenta.r <- function(pr_inf = c(0, 0, 10, 100, 1000, 10000)) {
    alpha_beta <- pr_inf + 2
    alpha <- c(1, m * alpha_beta[-1])
    beta <- c(1, alpha_beta[-1] - alpha[-1])
    a <- alpha + y
    b <- beta + n - y

    layout(matrix(1:6, 3, 2, byrow = T))
    par(mar = c(4, 4, 1, 1))

    for (i in 1:length(alpha)) {
        plot(theta, dbeta(theta, a[i], b[i]), type = "n", bty = "l", xlab = expression(theta),
            ylab = "density")

        text(x = 0.4, y = dbeta(a[i]/(a[i] + b[i]), a[i], b[i]) - 2.5, bquote(alpha +
            beta - 2 == .(pr_inf[i])), cex = 1.5)

        # likelihood It has been scaled (per 1000) to be comparable with posterior
        # at alpha+beta=2.

        lines(theta, 1000 * choose(n, y) * theta^y * (1 - theta)^(n - y), lwd = 2, col = "orange")

        # prior

        lines(theta, dbeta(theta, alpha[i], beta[i]), lwd = 2, col = "aquamarine")

        # posterior

        lines(theta, dbeta(theta, a[i], b[i]), lwd = 2, col = "dodgerblue")

        # posterior inference

        ci <- qbeta(c(0.025, 0.975), a[i], b[i])

        ind_ci <- theta >= ci[1] & theta <= ci[2]
        polygon(c(ci[1], ci[1], theta[ind_ci], ci[2], ci[2]), c(0, dbeta(c(ci[1], theta[ind_ci],
            ci[2]), a[i], b[i]), 0), col = "lightblue", border = F)

        abline(v = m, lty = 2)

        posum <- rbind(posum, c(pr_inf[i], alpha[i]/(alpha[i] + beta[i]), round(a[i]/(a[i] +
            b[i]), 3), round(qbeta(0.5, a[i], b[i]), 3), paste("[", round(ci[1], 3), ",",
            round(ci[2], 3), "]"), round(ci[2] - ci[1], 3)))

    }

    posum <- as.data.frame(posum)
    names(posum) <- c(expression(alpha + beta - 2), "mean", "mean", "median", "95% interval",
        "length")
    return(posum)
}

out <- placenta.r()

Analizziamo la sensitività o robustezza dell’inferenza rispetto alla scelta della a priori.
A former sensitivity analysis
Prior information
Posterior information
alpha + beta - 2 mean mean median 95% interval length
0 0.5 0.446 0.446 [ 0.415 , 0.477 ] 0.062
0 0.485 0.446 0.446 [ 0.415 , 0.477 ] 0.062
10 0.485 0.446 0.446 [ 0.416 , 0.477 ] 0.062
100 0.485 0.45 0.45 [ 0.42 , 0.479 ] 0.059
1000 0.485 0.466 0.466 [ 0.444 , 0.488 ] 0.044
10000 0.485 0.482 0.482 [ 0.472 , 0.491 ] 0.019

Si noti che nello studio \(n\approx 1000\) e \(\hat p=.446\).

Azioni immorali

# Internal calls in betabin.r are to the following functions:

logit.r <- function(p) {
    y <- log(p/(1 - p))
    return(y)
}


invlogit.r <- function(x) {
    y <- exp(x)/(1 + exp(x))
    return(y)
}

ci.r <- function(champion) {
    ci.mean <- mean(champion)
    ci.half <- 1.96 * sd(champion)
    return(c(ci.mean - ci.half, ci.mean + ci.half))
}
betabin.r <- function(y, n, a, b) {
    m <- length(a)  # m=1 or >1&even
    theta <- seq(0, 1, 0.001)

    upost <- rbeta(10000, y + 1, n - y + 1)
    uci <- round(quantile(upost, prob = c(0.025, 0.975), names = F), 3)

    maxupo <- max(hist(upost, nclass = 40, freq = F, plot = F)$density)
    maxli <- max(choose(n, y) * theta^y * (1 - theta)^(n - y))

    pr2po <- matrix(, nrow = m, ncol = 7)

    if (m > 1) {
        layout(matrix(1:m, nrow = m/2, byrow = T))
        par(mar = c(2, 2, 2, 2))
    }

    for (i in 1:m) {
        if (a[i] == 1 & b[i] == 1)
            post <- upost else post <- rbeta(10000, y + a[i], n - y + b[i])

        ci <- round(quantile(post, prob = c(0.025, 0.975), names = F), 3)
        nci <- round(ci.r(post), 3)
        logits <- logit.r(post)
        nlci <- ci.r(logits)
        nilci <- round(c(invlogit.r(nlci[1]), invlogit.r(nlci[2])), 3)

        maxpo <- max(hist(post, freq = F, plot = F)$density)
        dmax <- max(maxpo, maxupo)
        hist(post, freq = F, nclass = 40, col = "dodgerblue1", fg = gray(0.2), col.axis = gray(0.2),
            ylim = c(0, dmax), main = "")
        lines(theta, maxupo/maxli * choose(n, y) * theta^y * (1 - theta)^(n - y), lwd = 2,
            col = "coral")
        mtext(side = 4, col = "brown3", text = bquote(alpha + beta - 2 == .(round(a[i] +
            b[i] - 2, 4))), cex = 1, line = -1)
        mtext(side = 4, col = "brown3", text = bquote(alpha == .(a[i])), cex = 1, line = 0)

        abline(v = ci, col = "dodgerblue2", lwd = 3)
        abline(v = nilci, col = "violet", lwd = 2)
        abline(v = nci, col = "darkblue", lwd = 1)
        abline(v = uci, col = "aquamarine", lwd = 1.5)

        p <- a[i]/(a[i] + b[i])
        prv <- p * (1 - p)/(a[i] + b[i] + 1)
        pp <- (a[i] + y)/(a[i] + b[i] + n)
        # pp <- mean(post)
        pov <- pp * (1 - pp)/(a[i] + b[i] + n + 1)
        # pov <- var(post)

        pr2po[i, ] <- round(c(p, sqrt(prv), pp, sqrt(pov), ci, ci[2] - ci[1]), 3)
    }
    pr2po <- as.data.frame(pr2po)
    names(pr2po) <- c("mean", "sd", "mean", "sd", "95% L", "95% U", "length")
    return(pr2po)
}
alpha <- c(1, 0.001, 1, 1.8, 4.5, 45)
beta <- c(1, 0.001, 0.11, 0.2, 0.5, 5)

Le a priori

colori <- c("aquamarine3", "lightblue3", "goldenrod3", " orange", "coral", "brown3")
par(mar = c(2, 2, 1, 1))
for (i in length(alpha):1) {
    curve(dbeta(x, alpha[i], beta[i]), add = ifelse(i == 6, F, TRUE), lwd = 1 + 0.5 *
        i, ylim = c(0, 10), col = colori[i])
}

distinte per grado di “informatività”

par(mfrow = c(1, 3), mar = c(2, 2, 1, 1))
for (i in 1:2) {
    curve(dbeta(x, alpha[i], beta[i]), add = ifelse(i == 1, F, TRUE), col = colori[i],
        ylim = c(0, 1))
}
for (i in 3:4) {
    curve(dbeta(x, alpha[i], beta[i]), add = ifelse(i == 3, F, TRUE), col = colori[i],
        ylim = c(0, 9))
}
for (i in 5:6) {
    curve(dbeta(x, alpha[i], beta[i]), add = ifelse(i == 5, F, TRUE), col = colori[i])
}

n <- 751
y <- 150
out_1 <- betabin.r(y = y, n = n, a = alpha, b = beta)

out_1
   mean    sd  mean    sd 95% L 95% U length
1 0.500 0.289 0.201 0.015 0.173 0.230  0.057
2 0.500 0.500 0.200 0.015 0.172 0.229  0.057
3 0.901 0.206 0.201 0.015 0.172 0.231  0.059
4 0.900 0.173 0.202 0.015 0.173 0.231  0.058
5 0.900 0.122 0.204 0.015 0.176 0.234  0.058
6 0.900 0.042 0.243 0.015 0.215 0.274  0.059
Large sample
Prior information
Posterior information
mean sd mean sd 95% L 95% U length
0.500 0.289 0.201 0.015 0.173 0.230 0.057
0.500 0.500 0.200 0.015 0.172 0.229 0.057
0.901 0.206 0.201 0.015 0.172 0.231 0.059
0.900 0.173 0.202 0.015 0.173 0.231 0.058
0.900 0.122 0.204 0.015 0.176 0.234 0.058
0.900 0.042 0.243 0.015 0.215 0.274 0.059 
n <- 5
y <- 1
out_2 <- betabin.r(y = y, n = n, a = alpha, b = beta)

Tiny sample
Prior information
Posterior information
mean sd mean sd 95% L 95% U length
0.500 0.289 0.286 0.160 0.044 0.647 0.603
0.500 0.500 0.200 0.163 0.006 0.607 0.601
0.901 0.206 0.327 0.176 0.052 0.705 0.653
0.900 0.173 0.400 0.173 0.101 0.749 0.648
0.900 0.122 0.550 0.150 0.254 0.825 0.571
0.900 0.042 0.836 0.049 0.730 0.922 0.192