Lab on Binomial model

Bayesian Analysis of Binomial Data

library(dplyr) # optional
library(kableExtra)
library(tidyr) # gather
library(gridExtra)

Exercises on Bayesian inference with binomial data, taken from the Bayesian Data Analysis course site (associated to BDA textbook).

Visualization by ggplot2

Probabilità di una nascita femminile data la condizione placenta previa

(BDA3 p. 37).

library(ggplot2)
theme_set(theme_minimal())

L’a posteriori è una Beta(438,544)

# seq creates evenly spaced values
df1 <- data.frame(theta = seq(0.375, 0.525, 0.001)) 
a <- 438
b <- 544
# dbeta computes the posterior density
df1$p <- dbeta(df1$theta, a, b)

Calcoliamo anche l’intervallo centrale del 95%

# seq creates evenly spaced values from 2.5% quantile to 97.5% quantile (i.e., 95% central interval)
# qbeta computes the value for a given quantile given parameters a and b
df2 <- data.frame(theta = seq(qbeta(0.025, a, b), qbeta(0.975, a, b), length.out = 100))
# compute the posterior density
df2$p <- dbeta(df2$theta, a, b)

Si disegni l’a posteriori (Beta(438,544)) e la linea verticale per la media della popolazione (48.8%)

ggplot2

One of the key ideas behind ggplot2 is that it allows you to easily iterate, building up a complex plot a layer at a time. Each layer can come from a different dataset and have a different aesthetic mappings, making it possible to create sophisticated plots that display data from multiple sources.

ggplot(mapping = aes(theta, p)) +
  geom_line(data = df1) + 
  # Add a layer of colorized 95% posterior interval
  geom_area(data = df2, aes(fill='1')) +
  # Add the proportion of girl babies in general population
  geom_vline(xintercept = 0.488, linetype='dotted', linewidth=1.25) +
  # Decorate the plot a little
  labs(title='Uniform prior -> Posterior is Beta(438,544)',x=expression(theta), y = '') +
  scale_y_continuous(expand = c(0, 0.1), breaks = NULL) +
  scale_fill_manual(values = 'lightblue', labels = '95% posterior interval') +
  theme(legend.position = 'bottom', legend.title = element_blank())

change from ggplot2 3.4.0

Warning: Using size aesthetic for lines was deprecated in ggplot2 3.4.0. ℹ Please use linewidth instead.

Effetto di diverse a priori in un modello binomiale

Effetto dell’a priori e confronto delle a posteriori con valori dei parametri differenti per le a priori beta

Oltre a ggplot2 per la visualizzazione, tidyr è usato per lavorare sui data frame

Dati osservati: 437 femmine e 543 maschi

a <- 437
b <- 543

La densità viene valutata nei punti equamente distribuiti tra 0.375 e 0.525

df1 <- data.frame(theta = seq(0.375, 0.525, 0.001))

A posteriori con Beta(1,1), ie. a priori uniforme

df1$pu <- dbeta(df1$theta, a+1, b+1)

Tre scelte differenti per le a priori

Beta(\(0.488*2\),\((1-0.488)*2\))
Beta(\(0.488*20\),\((1-0.488)*20\))
Beta(\(0.488*200\),\((1-0.488)*200\))

n <- c(2, 20, 200) # conteggi a priori
apr <- 0.488 # prop di successo a priori

# helperf ritorna, dati un numero di osservazioni a priori, una prop di successo a priori, 
# un numero di successi e fallimenti osservati e un dataframe con i valori per theta, 
# un dataframe nuovo con valori a priori e a posteriori calcolati nei punti theta.
# 
helperf <- function(n, apr, a, b, df)
  cbind(df, pr = dbeta(df$theta, n*apr, n*(1-apr)), po = dbeta(df$theta, n*apr + a, n*(1-apr) + b), n = n)

# lapply the above function over prior counts n, then bind by row, and lastly gather results into key-value pairs.
df2 <- lapply(n, helperf, apr, a, b, df1) %>% do.call(rbind, args = .)  %>%
  gather(key = grp, value = p, -c(theta, n), factor_key = T)
df2 %>% head; df2 %>% tail
  theta n grp            p
1 0.375 2  pu 0.0008507679
2 0.376 2  pu 0.0011418738
3 0.377 2  pu 0.0015257278
4 0.378 2  pu 0.0020295190
5 0.379 2  pu 0.0026876416
6 0.380 2  pu 0.0035433725
     theta   n grp            p
1354 0.520 200  po 0.0006885263
1355 0.521 200  po 0.0005006071
1356 0.522 200  po 0.0003622430
1357 0.523 200  po 0.0002608721
1358 0.524 200  po 0.0001869727
1359 0.525 200  po 0.0001333673
# add correct labels for plotting
df2$title <- factor(paste0('alpha/(alpha+beta)=0.488, alpha+beta=',df2$n))
levels(df2$grp)
[1] "pu" "pr" "po"
levels(df2$grp) <- c('Posterior with unif prior', 'Prior', 'Posterior')
df2 %>% head
  theta n                       grp            p
1 0.375 2 Posterior with unif prior 0.0008507679
2 0.376 2 Posterior with unif prior 0.0011418738
3 0.377 2 Posterior with unif prior 0.0015257278
4 0.378 2 Posterior with unif prior 0.0020295190
5 0.379 2 Posterior with unif prior 0.0026876416
6 0.380 2 Posterior with unif prior 0.0035433725
                                   title
1 alpha/(alpha+beta)=0.488, alpha+beta=2
2 alpha/(alpha+beta)=0.488, alpha+beta=2
3 alpha/(alpha+beta)=0.488, alpha+beta=2
4 alpha/(alpha+beta)=0.488, alpha+beta=2
5 alpha/(alpha+beta)=0.488, alpha+beta=2
6 alpha/(alpha+beta)=0.488, alpha+beta=2

Grafico delle distribuzioni

ggplot(data = df2) +
  geom_line(aes(theta, p, color = grp)) +
  geom_vline(xintercept = apr, linetype = 'dotted') +
  facet_wrap(~ title, ncol = 1) +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Inferenza basata sulla simulazione

Simuliamo campioni da una Beta(438,544), disegniamo un istogramma (con quantili), e facciamo lo stesso per una variabile trasformata.

Campioniamo dalla a posteriori Beta(438,544). Otteniamo tutte le estrazioni in una volta e salviamole in un vettore ‘theta’

a <- 438
b <- 544
theta <- rbeta(10000, a, b)

Calcoliamo l’odds ratio per tutte le estrazioni

phi <- theta/(1 - theta)

Calcoliamo i quantili 2.5% e 97.5% aprossimati usando i campioni

quantiles <- c(0.025, 0.975)
thetaq <- quantile(theta, quantiles)
phiq <- quantile(phi, quantiles)

Istogrammi con 30 bins per theta e phi

# merge the data into one data frame for plotting
df1 <- data.frame(phi,theta) %>% gather()
#df1 %>% head
# merge quantiles into one data frame for plotting
df2 <- data.frame(phi = phiq, theta = thetaq) %>% gather()
#df2
ggplot() +
  geom_histogram(data = df1, aes(value), bins = 30) +
  geom_vline(data = df2, aes(xintercept = value), linetype = 'dotted', linewidth = 1.25) +
  facet_wrap(~key, ncol = 1, scales = 'free_x')  +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL)

Se vogliamo aggoiungere la densità

ggplot(data = df1, aes(value)) +
  geom_histogram(aes(y = after_stat(density)), fill="white", col="black", bins = 30) +
  geom_density(col = 2, fill = 2, alpha=.2) +
  geom_vline(data = df2, aes(xintercept = value), linetype = 'dotted', linewidth = 1.25) +
  facet_wrap(~key, ncol = 1, scales = 'free')  +
  labs(x = '', y = '') +
  scale_y_continuous(breaks = NULL)

Grid e inverse-cdf sampling

Campionamento con griglia e con l’inversa della funzione di ripartizione

  • Calcoliamo la distribuzione a posteriori su una griglia discreta di punti moltiplicando la verosimiglianza e un’a priori non coniugata in ogni punto e poi normalizzando sui punti.
  • Simuliamo campioni dalla distribuzione a posteriori non standard risultante utilizzando l’inversa della cdf e la griglia discreta.

A posteriori con osservazioni (437,543) e a priori uniforme (Beta(1,1))

a <- 438
b <- 544

df1 <- data.frame(theta = seq(0.1, 1, 0.001))
df1$con <- dbeta(df1$theta, a, b)

Calcoliamo la densità di un’a priori non-coniugata nei punti discreti di una griglia.

E.g., un’a priori lineare a pezzi con una certa % di probabilità al di fuori dell’intervallo (0.388, 0.588) (si veda la figura 2.4 in BDA3)

pp <- rep(1, nrow(df1))
pi <- sapply(c(0.388, 0.488, 0.588), function(pi) which(df1$theta == pi))
pi
[1] 289 389 489
pm <- 11
pp[pi[1]:pi[2]] <- seq(1, pm, length.out = length(pi[1]:pi[2])) #%>% print
pp[pi[3]:pi[2]] <- seq(1, pm, length.out = length(pi[3]:pi[2])) #%>% print

Normalizziamo la a priori

df1$nc_p <- pp / sum(pp)

Calcoliamo l’a posteriori non-normalizzata e non-coniugata sulla griglia, quindi la normalizziamo

po <- dbeta(df1$theta, a, b) * pp

df1$nc_po <- po / sum(po)

Disegniamo l’a posteriori con a priori uniforme, l’a priori non-coniugata e la corrispondente a posteriori non-coniugata

# gather the data frame into key-value pairs
# and change variable names for plotting
df2 <- gather(df1, grp, p, -theta, factor_key = T) #%>%
levels(df2$grp) <- c('Posterior with uniform prior', 'Non-conjugate prior', 'Non-conjugate posterior')
ggplot(data = df2) +
  geom_line(aes(theta, p)) +
  facet_wrap(~grp, ncol = 1, scales = 'free_y') +
  coord_cartesian(xlim = c(0.35,0.6)) +
  scale_y_continuous(breaks=NULL) +
  labs(x = '', y = '')

Inverse cdf sampling

calcoliamo la densità cumulata (funzione di ripartizione) sulla griglia

df1$cs_po <- cumsum(df1$nc_po)

Campioniamo dalla distr uniforme U(0,1)

# set.seed(seed) is used to set seed for the random number generator
set.seed(2601)
# runif(k) returns k uniform random numbers from interval [0,1]
r <- runif(10000)

Campionamento tramite l’inversa della cdf

# function to find the smallest value theta at which the cumulative sum of the posterior densities is greater than r.
invcdf <- function(r, df) df$theta[sum(df$cs_po < r) + 1]
# sapply function for each sample r. The returned values s are now random draws from the distribution.
s <- sapply(r, invcdf, df1)

Creiamo tre grafici: p1 è l’a posteriori, p2 è la cdf dell’a posteriori e p3 è l’istogramma delle estrazioni dall’a posteriori (estratte usando l’inversa-cdf)

p1 <- ggplot(data = df1) +
  geom_line(aes(theta, nc_po)) +
  coord_cartesian(xlim = c(0.35, 0.6)) +
  labs(title = 'Non-conjugate posterior', x = '', y = '') +
  scale_y_continuous(breaks = NULL)
p2 <- ggplot(data = df1) +
  geom_line(aes(theta, cs_po)) +
  coord_cartesian(xlim = c(0.35, 0.6)) +
  labs(title = 'Posterior-cdf', x = '', y = '') +
  scale_y_continuous(breaks = NULL)
p3 <- ggplot() +
  geom_histogram(aes(s, y = after_stat(density)), binwidth = 0.003) +
  geom_line(aes(s, y = after_stat(density)), stat = "bin", binwidth = 0.003, col='red') +
  coord_cartesian(xlim = c(0.35, 0.6)) +
  labs(title = 'Histogram of posterior samples', x = '', y = '') +
  scale_y_continuous(breaks = NULL)
# combine the plots
grid.arrange(p1, p2, p3)