class: center, middle, inverse, title-slide .title[ # Modelli a parametro singolo ] .subtitle[ ## <br><br> Un’introduzione a BUGS e all’apprendimento bayesiano in due modelli di dati standard: normale e Poisson ] .author[ ### <br><br>Matilde Trevisani ] .date[ ### 2024-10-21 ] --- class: middle ## Introduzione a BUGS --- ## BUGS BUGS project at https://www.mrc-bsu.cam.ac.uk/software Il progetto BUGS (**Bayesian Inference Using Gibbs Sampling**) si occupa di software flessibili per l'analisi bayesiana di modelli statistici complessi utilizzando metodi Markov chain Monte Carlo (MCMC). Il progetto è iniziato nel 1989 nell'unità di biostatistica dell'MRC, Cambridge, e ha portato inizialmente al programma BUGS "Classic", quindi al software WinBUGS sviluppato in collaborazione con l'Imperial College School of Medicine di St Mary's di Londra. Gli sviluppi si sono successivamente concentrati su OpenBUGS, un equivalente open source di WinBUGS. Esistono varie versioni di Bugs. --- ## Software https://www.mrc-bsu.cam.ac.uk/software https://drive.google.com/drive/folders/1XzG681AaZpRtPtynkuGUt2BspMTvhYoP **WinBUGS** It is an established and stable, stand-alone version of the software, which will remain available but not further developed. **OpenBUGS** It is an open-source version of the package, on which all future development work will be focused. For reference, see The BUGS project: Evolution, critique and future directions, by David Lunn, David Spiegelhalter, Andrew Thomas and Nicky Best, Statistics in Medicine, 2009. --- **JAGS** JAGS è Just Another Gibbs Sampler. It is a program for analysis of Bayesian hierarchical models using MCMC simulation not wholly unlike BUGS. JAGS was written with three aims in mind: • To have an engine for the BUGS language that runs on Unix • To be extensible, allowing users to write their own functions, distributions and samplers. • To be a plaftorm for experimentation with ideas in Bayesian modelling. To this last end, JAGS is licensed under the GNU General Public License. --- ### Ciò di cui ci dobbiamo solamente preoccupare: il modello 1. L'utente specifica un modello statistico, di (quasi) arbitraria complessità (semplicemente specificando le relazioni tra le variabili coinvolte). 2. Il software include un "**sistema esperto**", che determina uno schema MCMC appropriato (basato in generale sul Gibbs sampling) per analizzare il modello specificato. 3. L'utente controlla quindi l'esecuzione dello schema e può scegliere tra un'ampia gamma di tipi di output. *Tutto ciò di cui ci si deve preoccupare è essenzialmente specificare il modello!* --- class: middle ### Apprendimento bayesiano, inferenza e previsione in due modelli di dati standard: normale e Poisson --- ### Esempi Due esempi BUGS per modelli parametrici uno-dimensionali - Il modello Normale univariato con media `\(\mu\)` incognita e varianza `\(\sigma^2\)` nota. - Il modello Poisson per conteggi di eventi. --- ### Modello Normale univariato Congdon, pag. 17, Example 2.1, Program 2.1 *Systolic Blood Pressure* - Si supponga di estrarre un campione casuale di `\(20\)` letture della pressione sistolica del sangue `\(y_i\)` da una sottopopolazione di uomini adulti, che potrebbe costituire un particolare gruppo diagnostico. - Sappiamo da un'indagine nazionale sulla salute della popolazione che `\(\sigma =13\)`. - Siamo interessati a stimare `\(\mu\)`, la media della pressione sanguigna nel nostro gruppo, oltre a predire il livello della pressione in un tipico nuovo paziente dello stesso gruppo. - Si supponga di specificare un'a priori non informativa per `\(\mu\)`. - Si supponga di sapere che la pressione sanguigna tipica per tutti i maschi adulti sia `\(125\)`, e vogliamo testare se il particolare gruppo diagnostico ha pressione media superiore o inferiore, ossia `\(H_0:\,\mu\geq 125\)` vs. `\(H_1:\,\mu< 125\)`. --- .bold[Modello Normale univariato] - Queste domande potrebbero essere risolte direttamente usando la distribuzione di probabilità Normale. Poichè la verosimiglianza `$$p(y|\theta)\equiv L(\mu;y) \;{\rm è} \; \propto \exp(-1/(2\sigma^2)\,(y-\mu)^2)$$` <br> se l'a priori è della stessa forma, e.g., `$$p(\theta) \;{\rm è} \; \propto exp(-1/(2\tau_0^2)\,(\mu-\mu_0)^2),$$` <br> allora l'a posteriori manterrà anche questa forma. <br> Infatti, `$$p(\theta|y) \;{\rm è} \; \propto \;exp(-1/2(1/\tau_0^2\,(\mu-\mu_0)^2\,+\,1/\sigma^2\,(y-\mu)^2))$$` L'a priori Normale è una famiglia coniugata per il modello Normale - ma una prospettiva campionaria è egualmente percorribile. (fare l'esercizio in BUGS) --- ### Poisson Congdon, pag. 36, Example 2.15, Program 2.15 *Trent Leukaemia Mortality* - Confronto della mortalità tra aree dopo aver standardizzato per l'età (che è un fattore che influenza il rischio di mortalità). - I dati consistono nelle morti (1989) per *myeloid leukaemia* in Derby, che denotiamo `\(y_1\)`, e nel resto della regione Trent dell'Inghilterra, `\(y_2\)`, di cui Derby è una parte. - L'obiettivo è confrontare la mortalità tra Derby, l'area indice, e l'intera regione Trent definita come la popolazione standard. - Siano `\(y_{i,j}\)` le morti osservate, e `\(n_{ij}\)` le popolazioni per i gruppi d'età `\(i\)` nell'area `\(j\)` dove `\(p^*_i=n_{i1}/(n_{i1}+n_{i2})\)` è la quota sul totale. --- #### Confronto della mortalità tra aree: versione 1 del modello Se `\(m_i\)` è il tasso di mortalità per l'età `\(i\)`, allora le morti attese nell'area `\(j\)` `$$E_j=\sum_i\,m_in_{ij}$$` Supponendo che non ci sia interazione tra area e età, le morti osservate nell'area `\(j\)` `$$y_j=\sum_iy_{ij} \sim Pois(E_j\rho_j)$$` dove `\(\rho_j=1\)` indica che i tassi di mortalità dell'area indice e dell'area standard sono uguali. `\(\rho_j\)` è anche noto come SMR dell'area `\(j\)`. -- Nell'esempio, `\(y_1=30\)`, e dopo aver calcolato i tassi di mortalità per età basati sull'intera regione, otteniamo che `\(E_1=22.27\)` --- #### Confronto della mortalità tra aree: versione 2 del modello Se si considerano variabili i tassi di mortalità per età (denotati `\(m_i\)` sopra), allora le morti `\(y_{ij}\)` specifiche per età sono considerate outcomes di poisson `$$y_{ij}\sim\,Pois(\lambda_{ij})$$` con medie `$$\lambda_{ij}=\theta_{ij}n_{ij}$$` dove `\(\theta_{ij}\)` è il sottostante tasso di mortalità dell'area `\(j\)` per età `\(i\)`. Le morti attese per età `\(i\)` in tutta la regione sono quindi `\(\lambda_{i1}+\lambda_{i2}\)` e le morti attese nell'area indice sono `$$E_1=\sum_i\,p_i^{*}(\lambda_{i1}+\lambda_{i2})$$` dove `\(p^*_i\)` è la quota di popolazione di età `\(i\)` nell'area indice sul totale.