Lab 3

Modelli multiparametrici: Normale e Multinomiale

library(dplyr) # optional
library(kableExtra)
library(ggplot2) # the last three libr for ex 2 
library(gridExtra)
library(tidyr)

Esempi

  1. Il modello Normale univariato con media e varianza non note e a priori non-informativa

  2. Stima della velocità della luce (impianto teorico come sopra)

  3. Il modello Multinomiale

Riferimento

  1. Congdon: ex. 2.4, 21, Survival times from carcinoma
  2. Gelman: sec. 3.2, p. 66, Estimating the speed of light
  3. Congdon: ex. 2.16, 39, Cancers in women

Linguaggio

  1. BUGS, Program 2.4 Carcinoma survival
    R, norm2.r
  2. R
  3. BUGS, Program 2.16 Female Cancers
    R, mnomial.r

Soggetto

  1. Qual è la durata della sopravvivenza prevista per un nuovo paziente nelle stesse condizioni dei pazienti dell’indagine? E la probabilità di un tempo di sopravvivenza superiore a \(150\) settimane? E l’età incide su questi valori?
  2. Newcomb misurò (\(n=66\)) il tempo necessario alla luce per viaggiare dal suo laboratorio sul fiume Potomac a uno specchio alla base del Washington Monument e ritorno, una distanza totale di 7422 metri. Sulla base di questi dati, qual è una stima della velocità della luce?
  3. Nel 1995, il cancro al seno nel Regno Unito ha provocato il \(18\)% delle morti femminili per cancro. Qual è l’intervallo di credibilità per questo tasso in una situazione in cui abbiamo 11 tipi di cancro?

Esempio su normale con media e varianza non note e a priori non-informativa

Dati

dat <- read.table("norm2.dat", header=T)

par(mfrow=c(1,2))
for(i in 1:2) hist(dat[,i], main="", xlab=names(dat)[i])

Trasformiamo \(y\) in \(\log(y)\)

y <- dat[,1] %>% log()
y %>% hist(main="", xlab=expression(log(y)))

Analisi Bayesiana con a priori non-informativa

# Codice R per la generazione di campioni dalla densità a posteriori congiunta di (mu, sigma2) nel modello di dati Normale(mu, sigma2) con a priori non informativa p(mu,sigma2) proporzionale a 1/sigma2.
# 
norm2.r <- function(y, m=10000)
{
  n <- length(y)
  # s2: sample variance
  s2 <- 1/(n-1)*sum((y-mean(y))^2)
  ybar <- mean(y)
  
  # Simulation from the marginal posterior density for sigma2 
  chi <- rchisq(m, n-1)
  sigma2post <- (n-1)*s2/chi
  
  # Simulation from conditional posterior density for mu 
  # and from the (conditional) posterior predictive density for ynew. 
  # Simulation from a transformation of ynew: the probability that survival time exceeds 150 (log150=5.01) 
  
  mupostcond <- ynew <- g150 <- vector(length=m)

  for(i in 1:m){
    mupostcond[i] <- rnorm(1, ybar, sqrt(sigma2post[i]/n))
    ynew[i] <- rnorm(1, mupostcond[i], sqrt(sigma2post[i])) 
    g150[i] <- as.numeric(ynew[i]>5.01)
    }
  
  # Alternatively: simulation from posterior predictive distribution for a future observation ynew; 
  # probability computation of g.150 mean value

  ti <- rt(m,n-1)
  ynew_d <- sqrt(1+1/n)*sqrt(s2)*ti+ybar

  # Posterior probability of survival time exceeding 150

  g150_exact <- 1-pt((log(150)-ybar)/(sqrt(1+1/n)*sqrt(s2)),n-1)
  

    return(list(sigma=sigma2post,mu=mupostcond,ynew=ynew,g150=g150,ynew_d=ynew_d,g150_exact=g150_exact))
#  return(list(sigma=sigma2post,mu=mupostcond,ynew=ynew,g150=g150))
}
out <- norm2.r(y)
summary(out$mu)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.475   3.687   3.860   3.861   4.029   4.936 
sd(out$mu)
[1] 0.2580708
summary(out$ynew)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -3.579   3.062   3.846   3.846   4.619   9.679 
sd(out$ynew)
[1] 1.199971
summary(out$ynew_d)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -1.366   3.082   3.862   3.862   4.643  11.491 
summary(out$sigma)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4115  1.0137  1.2602  1.3558  1.5812  8.0492 
sd(out$sigma)
[1] 0.4957498
summary(out$g150)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.0000  0.1581  0.0000  1.0000 
out$g150_exact
[1] 0.1598639
bru=20
par(mfrow=c(2,2), mar=c(2,2,2,0.5))
hist(out$mu, prob=T, br=bru, main=expression(mu))
lines(density(out$mu))
hist(out$sigma, prob=T, br=bru, main=expression(sigma^2))
lines(density(out$sigma))
hist(out$ynew, prob=T, br=bru, main=expression(y[new]))
lines(density(out$ynew))
#hist(out$ynew_d)
lines(density(out$ynew_d), col=2)
hist(out$mu, prob=T, br=bru, main=bquote(mu~","~ y[new]), xlim=c(0,8))
lines(density(out$mu))
lines(density(out$ynew_d), col=2)

Stima della velocità della luce utilizzando il modello normale

Dati

Terza serie di misurazioni di Newcomb del tempo di passaggio della luce, effettuate dal 24 luglio 1882 al 5 settembre 1882. I valori indicati *0,001+24,8 sono le misurazioni di Newcomb, lette lungo le colonne, registrate in milionesimi di secondo. La tabella mostra tre set di dati raccolti in giorni diversi.

Newcomb data
Newcomb data

In Moodle i dati sono contenuti nel file light.txt

y <- read.table("light.txt")$V1
range(y)
[1] -44  40
n <- length(y) %>% print
[1] 66

Nota: il tempo effettivo misurato da Newcomb si ottiene a partire dai valori nei dati addizionando \(24800\) nanosecondi (\(10^{−9}\) secondi) (e.g., il primo valore osservato corrisponde a \(0.000024828\) secondi).

The data are recorded as deviations from 24800 nanoseconds.

sort(y)[1]
[1] -44
(24800+sort(y)[1])/10^9
[1] 2.4756e-05

Statistiche sufficienti

s2 <- var(y)
my <- mean(y)
cat(my, sqrt(s2))
26.21212 10.74532

Se eliminiamo i due valori anomali (che hanno valore negativo e riteniamo così solo valori positivi) le statistiche sufficienti diventano

y_pos <- y[y > 0]
range(y_pos)
[1] 16 40
n_pos <- length(y_pos) %>% print()
[1] 64

s2_pos <- var(y_pos)
my_pos <- mean(y_pos)
cat(my_pos, sqrt(s2_pos))
27.75 5.083431

Calcoliamo la densità di mu su questi punti

tl1 <- c(18, 34)
df1 <- data.frame(t1 = seq(tl1[1], tl1[2], length.out = 1000))

Calcoliamo la densità marginale esatta per mu (si usa la relazione delle densità di variabili trasformate: la densità della var su scala naturale è uguale alla densità della variabile standardizzata moltiplicata per lo Jacobiano) (Si faccia un confronto con l’es sopra)

# multiplication by 1./sqrt(s2/n) is due to the transformation of variable
# z=(x-mean(y))/sqrt(s2/n), see BDA3 p. 21
df1$pm_mu <- dt((df1$t1 - my)/sqrt(s2/n), n-1) / sqrt(s2/n)

Calcoliamo la densità marginale esatta per mu per i dati positivi

df1$pm_mu_pos = dt((df1$t1 - my_pos) / sqrt(s2_pos/n_pos), n_pos-1) / sqrt(s2_pos/n_pos)
theme_set(theme_minimal())

Creiamo un istogramma delle misure

p1 <- ggplot() +
  geom_histogram(aes(y), binwidth = 2, fill = 'steelblue', color = 'black') +
  coord_cartesian(xlim = c(-42, 42)) +
  scale_y_continuous(breaks = NULL) +
  labs(title = 'Newcomb\'s measurements', x = 'y')

Creiamo un grafico del modello normale

# gather the data points into key-value pairs
df1 %>% head
        t1        pm_mu    pm_mu_pos
1 18.00000 6.369248e-08 1.512572e-22
2 18.01602 6.681809e-08 1.643435e-22
3 18.03203 7.009541e-08 1.785760e-22
4 18.04805 7.353169e-08 1.940564e-22
5 18.06406 7.713453e-08 2.108955e-22
6 18.08008 8.091191e-08 2.292138e-22
df2 <- gather(df1, grp, p, -t1)
df2 %>% head
        t1   grp            p
1 18.00000 pm_mu 6.369248e-08
2 18.01602 pm_mu 6.681809e-08
3 18.03203 pm_mu 7.009541e-08
4 18.04805 pm_mu 7.353169e-08
5 18.06406 pm_mu 7.713453e-08
6 18.08008 pm_mu 8.091191e-08
# legend labels
labs2 <- c('Posterior of mu', 'Posterior of mu given y > 0', 'Modern estimate')
p2 <- ggplot(data = df2) +
  geom_line(aes(t1, p, color = grp)) +
  geom_vline(aes(xintercept = 33),
             linetype = 'dashed') +
  coord_cartesian(xlim = c(-42, 42)) +
  scale_y_continuous(breaks = NULL) +
  labs(title = 'Normal model', x = 'mu', y = '') +
  scale_color_manual(values = c('blue', 'darkgreen')) +
  guides(color="none") +
  annotate("text", x=24, y=0.26, label=labs2[1], hjust="right", size=5) +
  annotate("text", x=26, y=0.58, label=labs2[2], hjust="right", size=5) +
  annotate("text", x=32, y=0.7, label=labs2[3], hjust="right", size=5)

Combiniamo i grafici

grid.arrange(p1, p2)

Come risulta la stima della velocità della luce dai dati completi e da quelli senza i due valori anomali rispetto alla stima correntemente accettata della velocità della luce (nelle date condizioni dell’esperimento) corrispondente a \(33\) (\(0.00002483302\) secondi)?

Calcolare l’intervallo centrale per mu dai due datasets.

All 66 measurements:  23.57059 28.85365
Only positive data:   26.4802 29.0198

… This reinforces the fact that posterior inferences are only as good as the model and the experiment that produced the data. (Gelman)

Per quanto concerne i dati

sum(y>=33)
[1] 10
(1:n/n)[sort(y)==33]
[1] 0.8636364 0.8787879

il valore \(33\) corrisponde all’\(86\)esimo percentile.

Per quanto riguarda il modello, vedere nel capitolo 6 del Gelman le misure test per verificare l’adeguatezza del modello.

Dopo averne verificato l’inadeguatezza, possibili alternative sono un modello normale contaminato o un modello \(t\) di Student.

Esempio su multinomiale

Dati sui tumori femminili in UK nel 1995

# Data taken from Congdon ex. 2.16, Cancers in women
x <- c(14080,12990,6440,4350,3420,3190,2600,2420,1820,1760,23610)

Consideriamo un modello Dirichlet-Multinomial con a priori uniforme

mnomial.r <- function(x, m=10000)
{
  # Simulation from the posterior density for (theta_1,...,theta_k)
  pi_post <- rdiric.r(m,x+1)
  
  # Alternatively use: rdirichlet(n, alpha)
  
  # posterior Mean, Median and 95%CI for the category-1 probability
  mm <- summary(pi_post[,1])[3:4]
  ci <- quantile(pi_post[,1],c(.025,.975)) 
  
  # Alternatively, analytical computation of the posterior mean 
  # ((x[1]+1)/sum(x+1)) 
  
  return(list(mm,ci))
}


# R code for generating n random samples from a Dirichlet(a)
# It uses the fact that a Dirichlet can be obtained by sampling independent gamma densities (with equal scale parameter beta, e.g., beta=1, and shape parameter alpha equal to the corresponding element of the Dirichlet a vector), and then dividing such draws by their sum.

rdiric.r <- function(n,a) {
  p <- length(a)
  m <- matrix(nrow=n,ncol=p)
  for (i in 1:p) {
    m[,i] <- rgamma(n,a[i])
  }
  sumvec <- m %*% rep(1,p) # apply(m, 1 sum)
  m/as.vector(sumvec)
}
# Let try rdiric.r

op <- options()
options(digits=3)

rdiric.r(1,rep(1,10))
       [,1]   [,2]   [,3]   [,4]  [,5]   [,6]  [,7] [,8]    [,9] [,10]
[1,] 0.0504 0.0131 0.0505 0.0649 0.226 0.0512 0.302 0.11 0.00899 0.124

rdiric.r(4,rep(1,10))
        [,1]    [,2]     [,3]    [,4]   [,5]   [,6]   [,7]    [,8]   [,9]   [,10]
[1,] 0.09656 0.00629 0.265019 0.00599 0.1064 0.1702 0.0500 0.24636 0.0507 0.00257
[2,] 0.10848 0.33637 0.035592 0.04567 0.0981 0.0642 0.0453 0.00881 0.0422 0.21537
[3,] 0.00433 0.10030 0.000196 0.20476 0.1071 0.1610 0.1232 0.19716 0.0351 0.06688
[4,] 0.04727 0.12242 0.079743 0.02898 0.1103 0.1798 0.1749 0.03958 0.1037 0.11333

rdiric.r(4,rep(1,10)) %>% apply(1,sum)
[1] 1 1 1 1

options(op)     # reset (all) initial options
options("digits")
$digits
[1] 7
out <- mnomial.r(x)
out
[[1]]
   Median      Mean 
0.1836013 0.1836218 

[[2]]
     2.5%     97.5% 
0.1808959 0.1863826 
(x[1]+1)/(sum(x+1))
[1] 0.1836069