This book is in Open Review. I want your feedback to make the book better for you and other readers. To add your annotation, select some text and then click the on the pop-up menu. To see the annotations of others, click the in the upper right hand corner of the page

Cours 27 Introduction à Stan

library(readr)
library(tidyverse)
library(rstan)
rstan_options(auto_write = TRUE) #option pour ne pas recompiler à chaque fois !!! gain de temps
options(mc.cores = parallel::detectCores()) #option pour ajouter des coeurs au calcul
library(bayesplot) #visualiser la chaine de Markov
library(shinystan) #for interactive stan output visualization
library(rstanarm) #for Bayesian automatic regression modelling using stan
library(brms) #Bayesian generalized multivariate non-linear multilevel models using stan

27.1 Pourquoi les stats bayesiennes ?

  • on peut exprimer nos croyances/expertises sur les paramètres (prior)
  • prend en compte l’incertitude proprement
  • permet de prendre tout niveau de complexité de modele

Theoreme de Bayes * Prior : proba des paramètres apriori * Vraissemblance : proba des données sachant les paramètres * Posterior : distribution des paramètres sachant les données

Stan = un language : - entre crochet - un point-virgule à la fin de chaque ligne - // pour commencer un commentaire - <lower=0> : pour borner la variable

27.2 Stan program :

Ôuvrir un fichier stan. aucun # n’est toléré !! 3 blocks de commande :

data block on déclare : - la taille du jeu de donnée - les differentes variables

parameters block On déclare le nom des paramètres et leurs bornes si on le souhaite

model block - on déclare le prior (loi que suit le paramètre) (permet d’augmenter la vitesse d’analyse) - sinon prior non-informatif - écriture du modèle (vraissemenblance)

Les priors doivent être définis sous sens biologique/écologique

GLOPNET <- read_csv("GLOPNET.csv", skip=10)
LES <- GLOPNET %>% 
  filter(BIOME=="TROP_RF", GF == "T") %>% 
  select(Dataset, Species, "log LL", "log LMA") %>% 
  na.omit() #le bayesien ne supporte pas les NA
LES %>% 
  ggplot(aes(`log LMA`, `log LL`, col=Dataset)) +
  geom_point() +
  xlab("Logarithm of Leaf Mass per Area (LMA)") +
  ylab("Logarythm of Leaf Lifespan (LL)")

Modèle proposé :

_log LL ~ N(alpha + beta *log LMA, sigma^2)_ -> régression linéaire à mettre dans un fichier stan

data <- list(
  N = dim(LES) [1], #les noms des var doivent etre les memes que dans le fichier stan
  logLMA = LES$"log LMA",
  logLL = LES$"log LL"
  )

fit1 <- stan("stan.stan", data = data) #LET'S GOOOO! Données injectées, compilation lancée !
  • Nombre de chaine par défaut mais on peut choisir, 4 c’est le minimum pour interpréter les graphes.
  • thin = période d’amaigrissement = pas de l’itération
  • warmup = periode de chauffe : petite balade aléatoire
  • lp_ = log de la vraisemblance
  • n_eff = nbr d’itérations effectives : qui sont relevantes
  • Rhat doit être égal à 1 (ou 1.1): ça veut dire que l’estimation des paramètres a réussi à converger

Ici les chaînes n’ont pas réusi à converger, et ce sont perdues (tres peu de chaînes effectives).

on peut rajouter des lignes ggplot pour modifier les plots

mcmc_trace(as.array(fit1), #as.array = comme vecteurs
           facet_args=list(labeller=label_parsed), #pour mettre en lettre greques
           np = nuts_params(fit1)) # np pour afficher la divergeance

MCMC : diagnostic des chaines PPC : comparaison de modèles

Il faut une exploration indépendante des paramètres pour trouver le point où ça cohabite = maximum de vraissemblance (voir dessin carnet)

#les paramètres doivent être indépendants et former une patate à leur valeur correspondant aux données.
mcmc_pairs(as.array(fit1))

On voit que alpha et lp_ sont liés par une relation. On voit sur les graphes des chaines, que alpha veut aller dans le négatif mais ne peut pas car défini sur R+ (loi gamma)

Logique biologique : alpha doit être égal à 0 quand LMA est égal à 0 puisque s’il n’y a pas de masse il n’y pas de feuille dont pas de duree de vie !!! Donc on vire alpha du modèle.

Nouveau modele : _log LL ~ N( beta *log LMA, sigma^2)_ A mettre dans un fichier stan On peut faire une copie de l’ancien fichier : file > coche le stan file > more > copy

fit2 <- stan("LLLMA.stan", data = data)
mcmc_trace(as.array(fit2), #as.array : comme vecteurs
           facet_args=list(labeller=label_parsed), np = nuts_params(fit2))#pour mettre en lettre greques
mcmc_pairs(as.array(fit2))

Ca fait des belles patates donc c’est bon !!!! Et de beaux histos !

On veut mtn connaitre les distributions des paramètres à posteriori: on veux des paramètres : - différents de 0 (pour une relation entre var expli et var réponse) - et un sigma petit (pour un bon fit)

mcmc_areas(as.array(fit2), prob=0.95, pars = c("beta", "sigma")) # pars pour n'afficher que les paramètres que je veux, pas la vraissemblance

Posteriors sont normaux et significatifs.

mcmc_intervals(as.array(fit2), prob=0.95, pars = c("beta", "sigma")) #vu du dessus
launch_shinystan(fit2)

Conclusion : - se renseigner sur les formes de lois, de modèle - la définition des lois - centrer-réduire les variables pour faire converger plus vite (meme echelle) - borner les paramètres si possible