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

Il modello Normale univariato con media e varianza non note e a priori non-informativa
Stima della velocità della luce (impianto teorico come sopra)
Il modello Multinomiale

Riferimento

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

Linguaggio

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

Soggetto

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?
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?
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

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