Modelli multiparametrici

class: center, middle, inverse, title-slide

.title[
# Modelli multiparametrici
]
.subtitle[
## Normale e Multinomiale
]
.author[
### Matilde Trevisani
]
.date[
### 2022-12-01
]

---

## Modelli multiparametrici

Due esempi:

1. Il modello Normale univariato con media e  varianza non note

2. Il modello Multinomiale

---

## Tempi di sopravvivenza con il cancro

Congdon: esempio 2.4, pag. 21

 
*Program 2.4 Carcinoma Survival*

- Aitchison e Dunsmore (1975) presentano dati su
tempi di sopravvivenza `$y$` in settimane dopo che una combinazione di radioterapia e chirurgia sia applicata per un particolare carcinoma.

- A causa della forma asimmetrica dei dati originali, applichiamo una trasformazione logaritmica, i.e. supponiamo che `$y$` sia log-normale.

- Una questione di interesse è il tempo di sopravvivenza atteso per un nuovo paziente con questo carcinoma e sottoposto a questo tipo di trattamento.

- Un'altra domanda più specifica riguarda la probabilità che un nuovo paziente abbia un tempo di sopravvivenza superiore a `$150$` ( `$\log(150) = 5.01$` )

- Infine, può essere interessante esaminare gli effetti dell'età. Per esempio, possiamo dicotomizzare le età dei pazienti in età inferiori e superiori a `$50$` anni.

---

(esempio su tempi di spravvivenza continua)

- Supponiamo di assumere, in primo luogo, un'a priori non informativa.

- Secondariamenete, siamo disposti a specificare una convinzione a priori che la durata media di sopravvivenza prevista sia di `$30$` giorni
( `$\log(30)=3.4$` ) con un *numero di osservazioni pari a `$10$`* e che la varianza sia `$2$` con un *forza di grado `$10$`*.

---

### Inferenza su parametri normali: media e varianza non note

1. A priori non informative **indipendenti** 
`$p(\mu)\,p(\sigma^2)$`

*con R ...*

- A priori non informative: `$\,p(\mu,\sigma^2) \propto 1/\sigma^2$`

- calcolo della a posteriori attraverso la simulazione di marginale e condizionale: `$\,p(\mu,\sigma^2|y)=p(\mu|\sigma^2,y)p(\sigma^2|y)$`,
  ove

- distribuzione a posteriori marginale della varianza
  `$$p(\sigma^2|y)\sim \text{Inv-}{\chi^2}(n-1,s^2)$$`

- distribuzione a posteriori condizionale della media
  `$$p(\mu|\sigma^2,y)=\text{N}(\bar{y},\sigma^2/n)$$`
  
---

Inferenza su parametri normali: media e varianza non note (continua)

- distribuzione predittiva a posteriori per `$y_{n+1}$`:
  `$$p(y_{n+1}|{\bf y})=\int\,p(y_{n+1}|\mu,\sigma^2,{\bf y})\,p(\mu,\sigma^2|{\bf y})\,d\mu d\sigma^2$$`
    
    (i) simulazione della condizionale data l'estrazione dalla a posteriori di `$(\mu,\sigma^2)$`:
    `$$p(y_{n+1}|\mu,\sigma^2,{\bf y})=\text{N}(y_{n+1}|\mu,\sigma^2)$$`
    
    (ii) oppure, simulazione della marginale
    `$$p(y_{n+1}|{\bf y})= t_{n-1}\left(\bar{y},\left(1+\frac{1}{n}\right)s^2\right)$$`
    
  -  Probabilità di `$y_{n+1}>5.01 (=\log150)$`
  
    (i) simulazione della funzione `$f(y_{n+1})$`:
`$\,I_{[y_{n+1}>5.01]}(y_{n+1})$`

(ii) oppure, calcolo della probabilità: 
`$\,1-F_{t_{n-1}}((5.01-\bar{y})/s\sqrt{1+\frac{1}{n}})$`

---

#### Codice R (si veda il laboratorio 3)

.tiny[

```r
# 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 exceedes 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))
}
```
]

---

#### BUGS: a priori non informative indipendenti

BUGS richiede che sia definito un modello di probabilità (congiunto) completo, e quindi costringe **tutte le a priori** ad essere **proprie**
(in particolare le a priori sui "nodi fondatori" del grafico DAG).
 
Quindi, solitamente, vengono specificate a priori *appena proprie*. E.g., nel nostro esempio:

- a priori non informativa per la media:

`$N(0,\tau=0.0001)$` localmente uniforme su `$\pm 200$` ( `$\pm\,2$`sd )

- a priori non informativa per la precisione `$\tau$`:

`$G(0.0001,0.0001)$`  `$\, \sim 1/\tau$`

- calcolo delle a posteriori mediante **campionamento di Gibbs**.

---

### Inferenza su parametri normali: media e varianza non note

2. A priori **interdipendenti** 
`$p(\mu|\sigma^2)\,p(\sigma^2)$`

*con R ...*

- A priori coniugata: `$\text{N-Inv}_{\chi^2}(\mu_0,\sigma^2_0/m_0\,;\,\nu_0,\sigma^2_0)$`
  
`$$p(\mu|\sigma^2)p({\sigma^2})= \text{N}(\mu_0,\sigma^2/m_0)\;
\text{Inv}_{\chi^2}(\nu_0,\sigma^2_0)$$`

- Calcolo dalla a posteriori
`$\text{N-Inv}_{\chi^2}(\mu_n,\sigma^2_n/m_n\,;\,\nu_n,\sigma^2_n)$`
attraverso la simulazione marginale e condizionale:
`$p(\mu,\sigma|y)=p(\mu|\sigma,y)p(\sigma|y)$`

- Distribuzione marginale a posteriori della varianza:
`$$p(\sigma^2|y)\sim \text{Inv}_{\chi^2}(\nu_n,\sigma^2_n)$$`
- Distribuzione conditionale a posteriori della  media:
`$$p(\mu|\sigma,y)= \text{N}(\mu_n,\sigma^2/m_n)$$`
---

#### BUGS: a priori informative interdipendenti

- A priori per la media:  `$\,N(\mu_0,\,m_0\cdot\tau)$`

Nell'es. sui tempi di sopravvivenza: `$\mu_0=3.4$` ( `$=\log 30$` giorni), numerosità campionaria a priori `$m_0=10$`

- A priori per la precisione: `$\,G(\nu_0/2,\nu_0\,\sigma^2_0/2)$`

Nell'es.: `$\sigma^2_0=2$`, `$\nu_0=10$` a priori

- calcolo delle a posteriori con *Gibbs sampling*

---

## Tumori femminili

Congdon, Esempio 2.16, pag. 39
 
*Program 2.16 Female Cancers*

- Nel 1995, il cancro al seno nel Regno Unito rappresentava il `$18$`% delle
morti femminili per cancro.

- Supponiamo di essere interessati all'intervallo di credibilità per
questo tasso, in una situazione in cui abbiamo `$k=11$` tipi di cancro.

- Assumiamo una `$\text{Dirichlet}({\bf c})$` non informativa, con `$c_i=1, \,i=1,\ldots,k$`.

.small[Quindi, i dati del 1995 sul cancro femminile e altri tumori (il totale dei decessi per cancro è `$X=76680$`) tenderanno a “sommergere” l'a priori.]

---

### Dati multinomiali in Bugs

Possiamo usare
 - la distribuzione categoriale (discreta univariata)
`$$Y[i]\sim \text{dcat}(\pi[1:K])$$`
che genera una categoria `$j$` tra `$1$` e `$K$`,

- la distribuzione multinomiale (discreta multivariata)
`$$Z[i]\sim \text{dmulti}(\pi[1:K,1])$$`
che genera un vettore binario  `$K$`-dimensionale di  indicatori di categoria: 
    - `$Z_{ij}=1$` se `$j$` è la categoria scelta,   
    - `$Z_{ij}=0$` altrimenti

(i.e., un solo `$1$` e il resto `$0$`).

---

(esempio multinomiale continua)

- A parte la lo schema multinomiale, l'alternativa è
stimare i parametri multinomiali considerando ciascuno dei
`$k=11$` esiti come Poissons indipendenti.

- Possiamo assumere a priori Gamma non informative su ogni media di Poisson, e poi calcolare, per ogni estrazione
dalle medie a posteriori di queste `$11$` distribuzioni, le
stime `$\theta_i=\mu_i/\sum_i\mu_i$`.

---

#### Distribuzione multinomiale come prodotto di Poissons indipendenti

$$ p(y_1,y_2,\ldots, y_k)\propto \prod_j\pi_j^{y_j}
$$
come
`$$y_1\sim\,Po(\mu_1)\,,\;y_2\sim\,Po(\mu_2)\,,\ldots\ldots,\;y_k\sim\,Po(\mu_k)\,,$$`
sotto il vincolo che `$\sum_j\,y_j=n$`,
 
e con le probabilità multinomiali che si ottengono come `$\pi_j=\mu_j/\sum_j\mu_j$`

Dimostrazione
 
Se `$n$` è la somma di `$k$` Poissons, allora è una Poisson con media `$\sum_j\mu_j$`. Quindi la distribuzione di `$y_1,\ldots, y_k$` condizionale a `$n$` è

`\begin{align}
p(y_1,y_2,\ldots, y_k)/p(n)&=\\
&=\ldots
\end{align}`