Modélisation avec \(\texttt{Stan}\)

C’est le printemps!

Et l’occasion, comme la délicate fleur de cerisier, de s’ouvrir au monde et de faire la connaissance de \(\texttt{Stan}\). Nous allons brièvement voir comment coder des GLM avec \(\texttt{Stan}\), et progressivement complexifier notre modèle pour entrevoir l’ensemble des possibilités offertes avec ce language. On va commencer par une analyse classique, l’ANOVA, mais l’envisager comme un modèle hiérarchique sous le paradigme bayésien.

On va donc s’intéresser, puisque c’est de saison, à la fleuraison des cerisiers… sakura en japonais.

Description du jeu de données

### clean-up
rm(list = ls())

lapply(c("tidyverse", "readxl", "ggthemes", "sp", "sf", "rstan", "loo"), library, character.only=TRUE)
## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0       v purrr   0.3.1  
## v tibble  2.0.1       v dplyr   0.8.0.1
## v tidyr   0.8.3       v stringr 1.4.0  
## v readr   1.3.1       v forcats 0.4.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For improved execution time, we recommend calling
## Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
## although this causes Stan to throw an error on a few processors.
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
## This is loo version 2.0.0.
## **NOTE: As of version 2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. Visit mc-stan.org/loo/news for details on other changes.
## 
## Attaching package: 'loo'
## The following object is masked from 'package:rstan':
## 
##     loo
## [[1]]
##  [1] "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
##  [6] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"    
## [11] "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [16] "base"     
## 
## [[2]]
##  [1] "readxl"    "forcats"   "stringr"   "dplyr"     "purrr"    
##  [6] "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse"
## [11] "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [16] "methods"   "base"     
## 
## [[3]]
##  [1] "ggthemes"  "readxl"    "forcats"   "stringr"   "dplyr"    
##  [6] "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [11] "tidyverse" "stats"     "graphics"  "grDevices" "utils"    
## [16] "datasets"  "methods"   "base"     
## 
## [[4]]
##  [1] "sp"        "ggthemes"  "readxl"    "forcats"   "stringr"  
##  [6] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"   
## [11] "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices"
## [16] "utils"     "datasets"  "methods"   "base"     
## 
## [[5]]
##  [1] "sf"        "sp"        "ggthemes"  "readxl"    "forcats"  
##  [6] "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"    
## [11] "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics" 
## [16] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[6]]
##  [1] "rstan"       "StanHeaders" "sf"          "sp"          "ggthemes"   
##  [6] "readxl"      "forcats"     "stringr"     "dplyr"       "purrr"      
## [11] "readr"       "tidyr"       "tibble"      "ggplot2"     "tidyverse"  
## [16] "stats"       "graphics"    "grDevices"   "utils"       "datasets"   
## [21] "methods"     "base"       
## 
## [[7]]
##  [1] "loo"         "rstan"       "StanHeaders" "sf"          "sp"         
##  [6] "ggthemes"    "readxl"      "forcats"     "stringr"     "dplyr"      
## [11] "purrr"       "readr"       "tidyr"       "tibble"      "ggplot2"    
## [16] "tidyverse"   "stats"       "graphics"    "grDevices"   "utils"      
## [21] "datasets"    "methods"     "base"
theme_set(theme_bw(base_size = 20))


### set paths
setwd(WorkDir <- getwd())
DataDir <- paste(WorkDir, "data", sep = "/")
OutDir <- paste(WorkDir, "output", sep = "/")

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

datafile <- "Sweet_cherry_phenology_data_1978-2015.xlsx"

sakura <- read_excel(path = paste(DataDir, datafile, sep = "/"), sheet = 1)

### remove missing values
sakura <- sakura %>% 
  filter(!is.na(Full_Flowering),
         !is.na(Plantation)
         )

Nous travaillons avec un jeu de données , collectées entre 1978 et 2016. (et publiées dans l’article de Bénédicte Wendel de l’INRA Bordeaux et al., paru dans Scientific Data, 2016). Ces données correspondent aux dates de pleine floraison de \(9691\) individus/clones, chacun étant un cerisier Prunus avium en Europe. Pour chaque individu, il est possible de calculer son âge:

sakura <- sakura %>% 
  mutate(age = Year - Plantation,
         age = ifelse(age > 14, 14, age)
         )

## quick summary by age
sakura %>% 
  group_by(age) %>% 
  summarize(effectif = n(),
            flowering = round(mean(Full_Flowering, na.rm = TRUE), 1)
            )
## # A tibble: 14 x 3
##      age effectif flowering
##    <dbl>    <int>     <dbl>
##  1     1       44      89.2
##  2     2     1041      94  
##  3     3     1484      93.6
##  4     4     1455      94.8
##  5     5     1466      93.7
##  6     6     1298      93.7
##  7     7     1021      91.9
##  8     8      842      92.3
##  9     9      434      94.1
## 10    10      263      92.5
## 11    11      116      93.7
## 12    12       89      96.8
## 13    13       82      96  
## 14    14       56      97.1

On dispose également du site où se trouve cet individu:

## quick summary by Site (12 sites)
sakura %>% 
  group_by(Site) %>% 
  summarize(effectif = n(),
            flowering = round(mean(Full_Flowering, na.rm = TRUE), 1)
            )
## # A tibble: 12 x 3
##    Site              effectif flowering
##    <chr>                <int>     <dbl>
##  1 Balandran             2714      91.5
##  2 Bozas                  100     103. 
##  3 Carpentras             552      92.3
##  4 Etoile                1686      96.2
##  5 Montauban               79      88.9
##  6 Obernai                110     100. 
##  7 St Epain               221      97.5
##  8 St Gilles              401      95.4
##  9 St Laurent d'Agny      236      99.4
## 10 Torreilles             128      92  
## 11 Toulenne              3454      93  
## 12 Vienna                  10      99.8

ou de l’année d’observation :

## quick summary by decade (4)
sakura %>% 
  mutate(period = as.integer(cut(Year, breaks = seq(1980, 2020, 10)))) %>% 
  group_by(period) %>% 
  summarize(effectif = n(),
            flowering = round(mean(Full_Flowering, na.rm = TRUE), 1)
            ) 
## # A tibble: 4 x 3
##   period effectif flowering
##    <int>    <int>     <dbl>
## 1      1      769      93.7
## 2      2     3476      92.5
## 3      3     4328      94.3
## 4      4     1118      94.2

Un modèle d’ANOVA

La période ne semble pas avoir un effet marqué et on l’ignora dans la suite de l’analyse. Nous allons encore une fois regarder les dates de floraisons de ces cerisier (\(i \in [1:9691]\)), en tenant compte notamment de l’âge (\(\texttt{j}\)) et de la localisation des individus (\(\texttt{k}\)):

\[y_{ijk}|\mu_{ijk}, \sigma \sim \mathcal{N}(\mu_{ijk}, \sigma)\]

\[\mu_{ijk} = \beta_0 + \alpha_j + \delta_k\]

Les paramètres \(\alpha_j\) quantifient l’effet de l’âge des arbres sur la date de floraison observée, alors que les paramètres \(\delta_k\) quantifient l’effet de la localisation. Ici, deux options sont possibles pour le modélisateur:

  • On peut choisir de modéliser chaque paramètre \(\alpha_j\) (ou \(\delta_k\)) de manière indépendante, et donc de leur leur assigner à chacun un prior qui leur est propre (fixed effect), ou
  • on peut faire l’hypothèse moins forte que les \(\alpha_j\) (ou \(\delta_k\)) ne sont pas indépendants (au sens statistique du terme) mais qu’ils sont échangeables au sens de deFinetti, c’est-à-dire que notre jugement probabiliste (leur loi conjointe a priori) ne dépend pas de l’ordre où ils sont présentés: c’est dire qu’ils se ressemblent mais qu’il n’est pas possible de les classer a priori (ils sont invariants par permutation des labels \(j\) pour les \(\alpha_j\), ou \(k\) pour les \(\delta_k\)).

Si l’on munit les \(\alpha_j\) (ou les \(\delta_k\)) d’une distribution commune, cette hypothèse d’échangeabilité est réalisée.

\[\forall j, \alpha_j|\sigma_{\mathrm{age}} \sim \mathcal{N}(0, \sigma_{\mathrm{age}})\]

\[\forall k, \delta_k|\sigma_{\mathrm{site}} \sim \mathcal{N}(0, \sigma_{\mathrm{site}})\]

L’échangeabilité traduit une corrélation due à une ressemblance entre les arbres d’un même age et/ou d’un même site. Cette construction hiérarchique permet à l’inférence d’un des paramètres \(\alpha_{j_0}\) de s’appuyer –en plus des données concernant directement le groupe \(j_0\)– également et indirectement sur l’estimation des autres \(\alpha_{j}\) (l’expression anglaise consacrée est to borrow strength from neighbours)

Il ne reste plus qu’à spécifier des priors sur les paramètres \(\sigma\), \(\sigma_{\mathrm{age}}\) et \(\sigma_{\mathrm{site}}\). Ici, nous allons choisir une paramétrisation pour partitioner la variance totale \(\sigma_{\mathrm{tot}}^2\) :

\[\sigma_{\mathrm{tot}}^2 = \sigma^2 + \sigma_{\mathrm{age}}^2 + \sigma_{\mathrm{site}}^2\]

Soit

  • \(\sigma = \sigma_{\mathrm{tot}} \times \sqrt(\pi_1)\),
  • \(\sigma_{\mathrm{age}} = \sigma_{\mathrm{tot}} \times \sqrt(\pi_2)\), et
  • \(\sigma_{\mathrm{site}} = \sigma_{\mathrm{tot}} \times \sqrt(\pi_3)\)

avec \(\sum_{l=1}^{3} \pi_l = 1\) et \(\sigma_{\mathrm{tot}} \sim \mathcal{S}^+(0, 1, 3)\).

Ce modèle est une ANOVA à deux facteurs (\(\mathrm{age}\) et \(\mathrm{site}\)).

Ecritude du modèle d’ANOVA sous \(\texttt{Stan}\)

On va écrire ce modèle en \(\texttt{Stan}\), un language de programmation qui permet d’ajuster des modèles à des données en utilisant par défaut un algorithme Monte Carlo Hamiltonien. La première étape pour utiliser \(\texttt{Stan}\) est d’écrire soi même le modèle que l’on souhaite ajuster. Le code est très structuré et comporte plusieurs blocs :

functions {
    
}

data {
    
}

transformed data {
    
}

parameters {
    
}

transformed parameters {
    
}

model {
    
}

generated quantities {
    
}

On va ici écrire notre modèle pour estimer la date moyenne de floraison. Il faut commencer par définir les données que l’on va utiliser.

data {
  int<lower = 1> n_obs;
  int<lower = 1> n_age;
  int<lower = 1> n_site;
  vector[n_obs] FLOWERING;
  int<lower = 1, upper = n_age> AGE[n_obs];
  int<lower = 1, upper = n_site> SITE[n_obs];
  real prior_location_beta0; // hyper-parameter
  real<lower = 0.0> prior_scale_beta0; // hyper-parameter
}

Puis, il faut définir les paramètres (ou inconnues) du modèle.

parameters {
  simplex[3] pi;
  real beta0;
  real<lower = 0.0> sigma_tot;
  vector[n_age + 1] alpha;
  vector[n_site + 1] delta;
}

Il est possible de définir autant de quantités dérivées de ces paramètres (et des données) dans le bloc

transformed parameters {
  real R_squared;
  real sigma;
  real sigma_age;
  real sigma_site;
  vector[n_obs] mu; // linear predictor
  sigma = sqrt(pi[1]) * sigma_tot;
  sigma_age = sqrt(pi[2]) * sigma_tot;
  sigma_site = sqrt(pi[3]) * sigma_tot;
  mu = rep_vector(beta0, n_obs) + alpha[AGE] + delta[SITE]; // vectorized form
  R_squared = 1 - variance(FLOWERING - mu) / variance(FLOWERING);
}

Enfin, il faut spécifier les priors et la vraisemblance des données. Dans le bloc \(\texttt{model}\), il est attendu d’utiliser le symbole \(\sim\). C’est le seul bloc où on pourra d’ailleurs utiliser ce symbole signifiant tirage aléatoire.

model {
  beta0 ~ normal(prior_location_beta0, prior_scale_beta0);
  sigma_tot ~ student_t(3, 0.0, 1.0);
  alpha ~ normal(0.0, sigma_age);
  delta ~ normal(0.0, sigma_site);
  FLOWERING ~ normal(mu, sigma);
}

Enfin, il est possible de travailler avec la distribution a posteriori directement avec \(\texttt{Stan}\) pour faire par exemple de la prédiction.

generated quantities {
  vector[n_obs] log_lik; // for WAIC/loo computations
  real flower_new; // prediction
  vector[n_obs] y_rep; // posterior predictive check
  real earliest;
  real latest;
  flower_new = normal_rng(beta0 + alpha[n_age + 1] + delta[n_site + 1], sigma);
  for(i in 1:n_obs) {
   log_lik[i] = normal_lpdf(FLOWERING[i]| mu[i], sigma); // log probability density function
   y_rep[i] = normal_rng(mu[i], sigma); // prediction from posterior
  }
  earliest = min(y_rep); // ppc
  latest = max(y_rep); // ppc
}

Une fois ce modèle écrit (en \(\texttt{C++}\)), \(\texttt{Stan}\) va le compiler en un exécutable \(\texttt{.exe}\) que l’on va pouvoir appeler depuis notre environnement \(\texttt{R}\).

stan_anova <- '
 data {
  int<lower = 1> n_obs;
  int<lower = 1> n_age;
  int<lower = 1> n_site;
  vector[n_obs] FLOWERING;
  int<lower = 1, upper = n_age> AGE[n_obs];
  int<lower = 1, upper = n_site> SITE[n_obs];
  real prior_location_beta0;
  real<lower = 0.0> prior_scale_beta0;
 }
 
 parameters {
  simplex[3] pi;
  real beta0;
  real<lower = 0.0> sigma_tot;
  vector[n_age + 1] alpha;
  vector[n_site + 1] delta;
 }

 transformed parameters {
  real R_squared;
  real sigma;
  real sigma_age;
  real sigma_site;
  vector[n_obs] mu;
  sigma = sqrt(pi[1]) * sigma_tot;
  sigma_age = sqrt(pi[2]) * sigma_tot;
  sigma_site = sqrt(pi[3]) * sigma_tot;
  mu = rep_vector(beta0, n_obs) + alpha[AGE] + delta[SITE];
  //for (i in 1:n_obs) {
  // mu[i] = beta0 + alpha[AGE[i]] + delta[SITE[i]];
  //}
  R_squared = 1 - variance(FLOWERING - mu) / variance(FLOWERING);
 }

 model {
  beta0 ~ normal(prior_location_beta0, prior_scale_beta0);
  sigma_tot ~ student_t(3, 0.0, 1.0);
  alpha ~ normal(0.0, sigma_age);
  delta ~ normal(0.0, sigma_site);
  FLOWERING ~ normal(mu, sigma);
 }

 generated quantities {
  vector[n_obs] log_lik;
  vector[n_obs] y_rep; // posterior predictive check
  real flower_new;
  real earliest;
  real latest;
  flower_new = normal_rng(beta0 + alpha[n_age + 1] + delta[n_site + 1], sigma);
  for(i in 1:n_obs) {
   log_lik[i] = normal_lpdf(FLOWERING[i]| mu[i], sigma);
   y_rep[i] = normal_rng(mu[i], sigma);
  }
  earliest = min(y_rep);
  latest = max(y_rep);
 }
'

## Compilation
model_anova <- stan_model(model_code = stan_anova,
                          model_name = "Ma premiere ANOVA avec Stan"
                          )

La commande \(\texttt{rstan::stan_model}\) permet de compiler un modèle. Ce temps de compilation peut-être long et de peu d’intérêt ici : une fois le modèle compilé, on peut l’appeler grâe à la commande \(\texttt{rstan::sampling}\) sans avoir besoin de recompiler. La commande \(\texttt{rstan::stan}\) encapsule les commandes \(\texttt{rstan::stan_model}\) et \(\texttt{rstan::sampling}\) en une seule, mais l’inconvénient est que l’on va recompiler à chaque appel de cette fonction.

Ajustons le modèle aux données

n_chains <- 4
n_iter <- 1000
n_warm <- 500
n_thin <- 1

fit_1 <- sampling(object = model_anova,
                  data = list(n_obs = nrow(sakura),
                              n_age = length(unique(sakura$age)),
                              n_site = length(unique(sakura$Site)),
                              FLOWERING = sakura$Full_Flowering,
                              AGE = sakura$age,
                              SITE = as.numeric(factor(sakura$Site, levels = unique(sakura$Site))),
                              prior_location_beta0 = 90,
                              prior_scale_beta0 = 10
                              ), 
                  pars = c("beta0", "alpha", "delta", "sigma", "sigma_age", "sigma_site", "sigma_tot", "pi", "flower_new", "earliest", "latest", "R_squared", "log_lik"),
                  chains = n_chains, 
                  iter = n_iter, 
                  warmup = n_warm, 
                  thin = n_thin
                  )

get_elapsed_time(fit_1)
##          warmup  sample
## chain:1 129.618  94.972
## chain:2 153.414  80.935
## chain:3 170.641  89.257
## chain:4 127.589 114.477

La première chose à faire est de vérifier la convergence des paramètres.

stan_rhat(fit_1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Inférence bayésienne : examen de la loi conjointe a posteriori

On peut également regarder rapidement les estimations :

print(fit_1, digits = 3, pars = c("beta0", "alpha", "delta", "sigma", "sigma_age", "sigma_site", "sigma_tot", "pi", "flower_new", "earliest", "latest", "R_squared"))
## Inference for Stan model: Ma premiere ANOVA avec Stan.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##               mean se_mean    sd    2.5%     25%     50%     75%   97.5%
## beta0       95.558   0.075 1.481  92.551  94.569  95.495  96.604  98.378
## alpha[1]    -1.232   0.030 0.941  -3.061  -1.860  -1.202  -0.583   0.587
## alpha[2]     1.045   0.025 0.523   0.024   0.703   1.031   1.375   2.140
## alpha[3]     0.459   0.024 0.497  -0.474   0.132   0.449   0.800   1.442
## alpha[4]     1.609   0.025 0.501   0.654   1.289   1.596   1.931   2.589
## alpha[5]     0.498   0.025 0.503  -0.458   0.173   0.492   0.821   1.521
## alpha[6]     0.520   0.025 0.508  -0.435   0.170   0.503   0.871   1.534
## alpha[7]    -1.303   0.025 0.515  -2.273  -1.647  -1.301  -0.963  -0.301
## alpha[8]    -1.162   0.025 0.532  -2.156  -1.517  -1.171  -0.813  -0.108
## alpha[9]     0.020   0.024 0.581  -1.155  -0.366   0.033   0.397   1.192
## alpha[10]   -2.295   0.025 0.635  -3.616  -2.707  -2.293  -1.857  -1.111
## alpha[11]   -1.359   0.025 0.753  -2.818  -1.873  -1.361  -0.854   0.167
## alpha[12]    1.748   0.029 0.859   0.139   1.163   1.691   2.316   3.460
## alpha[13]    1.012   0.030 0.834  -0.657   0.468   0.995   1.557   2.729
## alpha[14]    0.765   0.027 0.922  -1.041   0.143   0.794   1.369   2.596
## alpha[15]    0.117   0.043 1.644  -3.034  -0.921   0.128   1.149   3.362
## delta[1]    -6.182   0.072 1.579  -9.218  -7.244  -6.162  -5.113  -3.131
## delta[2]    -4.377   0.071 1.404  -7.032  -5.331  -4.342  -3.451  -1.498
## delta[3]    -2.916   0.071 1.400  -5.520  -3.868  -2.885  -1.970  -0.071
## delta[4]     1.961   0.072 1.486  -0.815   0.931   1.963   2.953   4.885
## delta[5]    -3.155   0.071 1.422  -5.899  -4.122  -3.142  -2.219  -0.266
## delta[6]     0.204   0.072 1.411  -2.444  -0.743   0.231   1.145   2.946
## delta[7]     4.216   0.075 1.465   1.440   3.210   4.237   5.163   7.161
## delta[8]     7.288   0.073 1.556   4.397   6.199   7.288   8.339  10.511
## delta[9]    -0.174   0.070 1.430  -2.882  -1.143  -0.167   0.789   2.696
## delta[10]    4.648   0.074 1.529   1.845   3.591   4.655   5.634   7.733
## delta[11]    2.530   0.076 2.334  -1.800   0.940   2.490   4.028   7.284
## delta[12]   -3.124   0.070 1.527  -6.085  -4.142  -3.074  -2.105  -0.076
## delta[13]    0.122   0.121 4.296  -8.300  -2.675   0.199   2.776   8.839
## sigma        7.693   0.001 0.055   7.588   7.656   7.693   7.729   7.800
## sigma_age    1.570   0.014 0.419   0.964   1.276   1.498   1.787   2.625
## sigma_site   4.445   0.038 0.957   2.949   3.763   4.310   4.965   6.817
## sigma_tot    9.068   0.021 0.510   8.373   8.724   8.959   9.308  10.446
## pi[1]        0.726   0.003 0.075   0.544   0.685   0.737   0.779   0.842
## pi[2]        0.032   0.001 0.018   0.011   0.020   0.028   0.039   0.082
## pi[3]        0.242   0.003 0.076   0.124   0.187   0.231   0.286   0.424
## flower_new  95.753   0.228 9.068  78.243  89.610  95.726 102.097 112.598
## earliest    62.859   0.057 2.420  56.969  61.549  63.139  64.626  66.566
## latest     125.526   0.059 2.609 121.371 123.711 125.162 126.985 131.679
## R_squared    0.099   0.000 0.001   0.098   0.099   0.099   0.100   0.100
##            n_eff  Rhat
## beta0        392 1.008
## alpha[1]    1010 1.001
## alpha[2]     455 1.001
## alpha[3]     440 1.000
## alpha[4]     410 1.001
## alpha[5]     414 1.002
## alpha[6]     410 1.002
## alpha[7]     430 1.001
## alpha[8]     447 1.002
## alpha[9]     587 1.000
## alpha[10]    645 1.000
## alpha[11]    933 1.000
## alpha[12]    885 0.999
## alpha[13]    762 0.999
## alpha[14]   1139 1.001
## alpha[15]   1451 1.002
## delta[1]     480 1.009
## delta[2]     392 1.010
## delta[3]     390 1.010
## delta[4]     423 1.008
## delta[5]     407 1.010
## delta[6]     384 1.010
## delta[7]     381 1.010
## delta[8]     449 1.008
## delta[9]     415 1.009
## delta[10]    424 1.007
## delta[11]    939 1.002
## delta[12]    473 1.008
## delta[13]   1266 1.001
## sigma       2162 0.999
## sigma_age    836 1.003
## sigma_site   622 1.002
## sigma_tot    596 1.001
## pi[1]        637 1.001
## pi[2]        794 1.005
## pi[3]        626 1.002
## flower_new  1581 1.001
## earliest    1789 0.999
## latest      1987 1.000
## R_squared    856 1.001
## 
## Samples were drawn using NUTS(diag_e) at Mon Apr 15 15:03:15 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Avec le modèle qu’il est écrit, on a directement accès aux fractions respectives de la variance totale que l’on peut attribuer à l’âge, la zone ou la part résiduelle via les estimation des paramètres \(\hat{\pi}\).

apply(rstan::extract(fit_1, "pi")$pi, 2, mean)
## [1] 0.72588664 0.03200054 0.24211282
apply(rstan::extract(fit_1, "pi")$pi, 2, sd)
## [1] 0.07456004 0.01787135 0.07645566

Cela peut se représenter aussi graphiquement :

lower <- function(x, alpha = 0.8) { coda::HPDinterval(coda::as.mcmc(x), prob = alpha)[1] }
upper <- function(x, alpha = 0.8) { coda::HPDinterval(coda::as.mcmc(x), prob = alpha)[2] }
get_summary <- function(x, alpha = 0.8) { c(mean(x), sd(x), coda::HPDinterval(coda::as.mcmc(x), prob = alpha)) }

summary_anova <- as.data.frame(
  do.call('rbind', lapply(c("sigma", "sigma_age", "sigma_site"),
                          function(param) {
                            get_summary(as.numeric(rstan::extract(fit_1, param)[[1]]))
                            }
                          )
          )
)

names(summary_anova) <- c("mean", "se", "lower", "upper")
summary_anova$component <- c("residual", "age", "site")

summary_anova %>% 
  mutate(component = factor(component, levels = c("residual", "age", "site")[order(mean)])) %>% 
  ggplot(aes(x = component, y = mean)) +
  geom_linerange(aes(x = component, ymin = lower, ymax = upper)) +
  geom_point(size = 2) +
  ylab("Estimate") + xlab("Source of variation") +
  coord_flip() +
  theme_bw()

Pour tester formellement les différentes composantes de variation, il faudrait ici ajuster d’autres modèles sans ces composantes, puis comparer ces modèles entre eux via le \(\texttt{WAIC}\).

Pour le moment, on peut néanmoins conclure avec confiance qu’il y a une part de la variation observée qui est liée au site (on pourrait éventuellement considérer un effet spatial plus explicite). On peut regarder les estimations \(\hat{\delta}_k\).

freq <- sakura %>% 
  group_by(Site) %>% 
  summarize(effectif = n(),
            flowering = mean(Full_Flowering, na.rm = TRUE)
            ) 

moyenne = with(freq, sum(effectif * flowering) / sum(effectif))
freq[13, 1] = "Average"
freq[13, 2] = mean(freq$effectif)
freq[13, 3] = moyenne

post_site <- as.data.frame(t(apply(matrix(rep(rstan::extract(fit_1, "beta0")$beta0, each = length(unique(sakura$Site)) + 1), ncol = length(unique(sakura$Site)) + 1, byrow = TRUE) + rstan::extract(fit_1, "delta")$delta, 2, get_summary)))
names(post_site) <- c("mean", "se", "lower", "upper")
post_site$where <- c(unique(sakura$Site), "Average")
post_site <- cbind(post_site,
                   freq[match(post_site$where, freq$Site), c('flowering', 'effectif')]
                   )

post_site %>% 
  mutate(where = factor(where, levels = c(unique(sakura$Site), "Average")[order(mean)])) %>% 
  ggplot(aes(x = where, y = mean)) +
  geom_linerange(aes(x = where, ymin = lower, ymax = upper)) +
  geom_point(size = 2) +
  geom_point(aes(x = where, y = flowering, size = effectif), color = 'red', alpha = 0.3) +
  scale_y_continuous(name = "Estimate (days)", breaks = 95 -10:10) + 
  xlab("Site") +
  coord_flip() +
  theme_bw()

Accéder directement à la vraisemblance avec \(\texttt{Stan}\)

Mais nous allons maintenant nous découvrir un avantage bien particulier à \(\texttt{Stan}\) : la fonction \(\texttt{target +=}\). Cette fonction permet en effet d’accèder directement à l’incrémentation de la logvraisemblance, et donc de pouvoir écrire en théorie n’importe quel modèle dont on connaît la vraisemblance. En d’autres termes, à partir du moment où l’on sait écrire la vraisemblance du modèle qui nous intéresse, on sait l’ajuster avec \(\texttt{Stan}\)!

Supposons qu’il existe un potentiel génétique pour une floraison précoce chez les cérisiers. Cela se voit uniquement au moment de la floraison, et sur aucun autre aspect phénotypique. L’idée ici est de chercher à identifier les individus avec ce potentiel.

Nous allons supposer qu’il existe un mélange de deux sortes d’individus, mais nous ne savons pas vraiment qui est qui. Cela signifie qu’il existe dans notre jeu de données une proportion \(\mathrm{proba}\) de cerisiers à floraison tardive, et son complément \(1 - \mathrm{proba}\)) à floraison précoce. Donc un modèle plus juste pour les dates de floraison \(\mathbb{y}\) serait :

\[y_{ijk}| \beta_1, \sigma, \alpha_j, \delta_k \sim \mathcal{N}(\beta_1 + \alpha_j +\delta_k, \sigma)\] si l’individu \(i\) est précoce (avec une probabilité \(1 - \mathrm{proba}\)); et

\[y_{ijk}| \beta_2, \sigma, \alpha_j, \delta_k \sim \mathcal{N}(\beta_2 + \alpha_j +\delta_k, \sigma)\] si l’individu \(i\) est tardif (avec une probabilité \(\mathrm{proba}\)).

On peut donc voir le problème comme celui de données manquantes : on a perdu la variable \(z_{ijk}\) (i.e. un label latent) qui indique si un cerisier \(i\) d’âge \(j\) et situé à \(k\) est précoce (\(z_{ijk} = 0\)) ou tardif (\(z_{ijk} = 1\)). Mais l’algorithme de Stan ne sait pas travailler avec des variables discrètes rangées dans le bloc des inconnues (parameters) car pour résoudre les équations hamiltoniennes qui lui permettent de visiter l’espace des inconnues, il doit opérer un calcul de dérivées. Il faudra donc marginaliser l’effet des \(z_{ijk} = 1\) lors du calcul de la vraisemblance.

On peut écrire le modèle suivant :

\(z_{ijk} \sim \mathcal{B}(\mathrm{proba})\)

\(y_{ijk}|z_{ijk}, \alpha_j, \delta_k, \beta_1, \beta_2, \sigma \sim \mathcal{N}(\beta_{z_i} + \alpha_j + \delta_k, \sigma)\)

Cela est un modèle hiérarchique : il y a deux niveaux dans le modèle : celui des observations \(y_{ijk}\) et celui de la variable latente \(z_{ijk}\). Cela peut tout à fait s’ajuste avec \(\texttt{Stan}\). Sous ce modèle à mélange, la vraisemblance de la donnée \(y_{ijk}\) est (avec quelques abus de notation):

\[\mathcal{L}(y_{ijk}) = (1 - \mathrm{proba}) \times \mathcal{N}(\beta_1 + \alpha_j + \delta_k, \sigma) + \mathrm{proba} \times \mathcal{N}(\beta_2 + \alpha_j + \delta_k, \sigma)\]

La log-vraisemblance est donc (toujours avec des abus de notations):

\[\ell(y_{ijk}) = \log((1 - \mathrm{proba}) \times \mathcal{N}(\beta_1 + \alpha_j + \delta_k, \sigma) + \mathrm{proba} \times \mathcal{N}(\beta_2 + \alpha_j + \delta_k, \sigma))\]

Et elle s’écrit de cette manière là avec \(\texttt{Stan}\) :

  for(i in 1:n_obs) {
   target += log_sum_exp(log(proba) + normal_lpdf(FLOWERING[i] | mu[1, i], sigma),
   log1m(proba) + normal_lpdf(FLOWERING[i] | mu[2, i], sigma));
  }

On peut aussi utiliser le bloc \(\texttt{functions}\) pour définir de nouvelles distributions, ou fonctions etc..

functions {
  real TwoGaussianMixture_lpdf(real y, real prob, vector location, real scale) {
        real log_pdf[2];
    log_pdf[1] = log1m(prob) + normal_lpdf(y| location[1], scale);
    log_pdf[2] = log(prob) + normal_lpdf(y| location[2], scale);
    return log_sum_exp(log_pdf);
    }
}

D’où le modèle :

stan_mix1 <- '
 functions {
  real TwoGaussianMixture_lpdf(real y, real prob, vector location, real scale) {
        real log_pdf[2];
    log_pdf[1] = log1m(prob) + normal_lpdf(y| location[1], scale);
    log_pdf[2] = log(prob) + normal_lpdf(y| location[2], scale);
    return log_sum_exp(log_pdf);
  }
  real TwoGaussianMixture_rng(real prob, vector location, real scale) {
   int z;
   z = bernoulli_rng(prob);
   return  z ? normal_rng(location[2], scale) : normal_rng(location[1], scale); 
  }
 }

 data {
  int<lower = 1> n_obs;
  int<lower = 1> n_age;
  int<lower = 1> n_site;
  vector[n_obs] FLOWERING;
  int<lower = 1, upper = n_age> AGE[n_obs];
  int<lower = 1, upper = n_site> SITE[n_obs];
  real prior_location_beta0;
  real<lower = 0.0> prior_scale_beta0;
  real prior_location_diff;
  real<lower = 0.0> prior_scale_diff;
 }
 
 parameters {
  real<lower = 0.0, upper = 1.0> proba;
  simplex[3] pi;
  real beta0;
  real<lower = 0> difference;
  real<lower = 0.0> sigma_tot;
  vector[n_age + 1] alpha;
  vector[n_site + 1] delta;
 }

 transformed parameters {
  real R_squared;
  real sigma;
  real sigma_age;
  real sigma_site;
  vector[n_obs] mu[2];
  vector[2] beta;
  beta[1] = prior_location_beta0 + beta0 * prior_scale_beta0;
  beta[2] = beta[1] + difference;
  sigma = sqrt(pi[1]) * sigma_tot;
  sigma_age = sqrt(pi[2]) * sigma_tot;
  sigma_site = sqrt(pi[3]) * sigma_tot;
  //for (i in 1:n_obs) {
  // mu[1, i] = beta[1] + alpha[AGE[i]] + delta[ZONE[i]];
  // mu[2, i] = beta[2] + alpha[AGE[i]] + delta[ZONE[i]];
  //}
  mu[1] = rep_vector(beta[1], n_obs) + alpha[AGE] + delta[SITE];
  mu[2] = rep_vector(beta[2], n_obs) + alpha[AGE] + delta[SITE];
  R_squared =  1 - variance((1 - proba) * mu[1] + proba * mu[2] - FLOWERING) / variance(FLOWERING);
 }

 model {
  sigma_tot ~ student_t(3, 0.0, 1.0);
  beta0 ~ normal(0.0, 1.0);
  difference ~ normal(prior_location_diff, prior_scale_diff);
  alpha ~ normal(0.0, sigma_age);
  delta ~ normal(0.0, sigma_site);
  for(i in 1:n_obs) {
   //target += log_sum_exp(log1m(proba) + normal_lpdf(FLOWERING[i] | mu[1, i], sigma),
   //log(proba) + normal_lpdf(FLOWERING[i] | mu[2, i], sigma));
   //target += TwoGaussianMixture_lpdf(FLOWERING[i]| proba, to_vector(mu[1:2, i]), sigma);
   FLOWERING[i] ~ TwoGaussianMixture(proba, to_vector(mu[1:2, i]), sigma);
  }
 }

 generated quantities {
  vector[n_obs] log_lik;
  vector[n_obs] y_rep;
  int<lower = 0, upper = 1> type;
  real flower_new;
  real earliest;
  real latest;
  type = bernoulli_rng(proba);
  flower_new = normal_rng(beta[1 + type] + alpha[n_age + 1] + delta[n_site + 1], sigma);
  for(i in 1:n_obs) {
   log_lik[i] = log_sum_exp(log1m(proba) + normal_lpdf(FLOWERING[i] | mu[1, i], sigma),
   log(proba) + normal_lpdf(FLOWERING[i] | mu[2, i], sigma));
   y_rep[i] = TwoGaussianMixture_rng(proba, to_vector(mu[1:2, i]), sigma);
  }
  earliest = min(y_rep);
  latest = max(y_rep);
 }
'

## Compilation
mix1 <- stan_model(model_code = stan_mix1,
                   model_name = "Ma premiere ANOVA a melange avec Stan"
                   )

Ajustons le modèle aux données

fit_2 <- sampling(object = mix1,
                  data = list(n_obs = nrow(sakura),
                              n_age = length(unique(sakura$age)),
                              n_site = length(unique(sakura$Site)),
                              FLOWERING = sakura$Full_Flowering,
                              AGE = sakura$age,
                              SITE = as.numeric(factor(sakura$Site, levels = unique(sakura$Site))),
                              prior_location_beta0 = 90,
                              prior_scale_beta0 = 10,
                              prior_location_diff = 7,
                              prior_scale_diff = 3
                              ), 
                  pars = c("proba", "beta", "alpha", "delta", "sigma", "sigma_age", "sigma_site", "sigma_tot", "pi", "difference", "flower_new", "earliest", "latest", "R_squared", "log_lik"),
                  chains = n_chains, 
                  iter = n_iter, 
                  warmup = n_warm, 
                  thin = n_thin
                  )

get_elapsed_time(fit_2)
##          warmup  sample
## chain:1 1549.68 637.458
## chain:2 1597.79 606.092
## chain:3 1398.40 575.447
## chain:4 1500.39 896.093

Comme toujours, on commence par vérifier la convergence des paramètres.

stan_rhat(fit_2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Tous les paramètres semblent avoir convergé!

print(fit_2, digits = 3, pars = c("proba", "beta", "alpha", "delta", "sigma", "sigma_age", "sigma_site", "sigma_tot", "pi", "difference", "flower_new", "earliest", "latest", "R_squared"))
## Inference for Stan model: Ma premiere ANOVA a melange avec Stan.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##               mean se_mean    sd    2.5%     25%     50%     75%   97.5%
## proba        0.606   0.000 0.020   0.568   0.593   0.606   0.620   0.644
## beta[1]     89.384   0.069 1.398  86.573  88.502  89.384  90.270  92.118
## beta[2]     99.447   0.071 1.395  96.716  98.534  99.443 100.322 102.248
## alpha[1]    -1.602   0.029 1.034  -3.773  -2.252  -1.554  -0.907   0.360
## alpha[2]     1.043   0.027 0.582  -0.070   0.647   1.032   1.423   2.276
## alpha[3]     0.527   0.027 0.564  -0.623   0.153   0.519   0.886   1.653
## alpha[4]     1.618   0.027 0.565   0.516   1.239   1.600   1.994   2.766
## alpha[5]     0.571   0.028 0.567  -0.520   0.209   0.589   0.935   1.694
## alpha[6]     0.423   0.028 0.572  -0.686   0.058   0.427   0.799   1.605
## alpha[7]    -1.441   0.027 0.571  -2.594  -1.831  -1.439  -1.076  -0.310
## alpha[8]    -1.341   0.028 0.583  -2.525  -1.719  -1.339  -0.965  -0.208
## alpha[9]     0.027   0.027 0.613  -1.186  -0.365   0.028   0.432   1.232
## alpha[10]   -2.735   0.028 0.690  -4.212  -3.177  -2.725  -2.270  -1.408
## alpha[11]   -1.631   0.028 0.825  -3.353  -2.149  -1.613  -1.086  -0.028
## alpha[12]    2.328   0.029 0.856   0.762   1.738   2.305   2.894   4.029
## alpha[13]    1.201   0.028 0.911  -0.433   0.547   1.198   1.792   3.057
## alpha[14]    1.017   0.031 0.978  -0.882   0.376   0.989   1.660   2.964
## alpha[15]    0.058   0.044 1.760  -3.534  -1.067   0.102   1.143   3.572
## delta[1]    -6.521   0.066 1.501  -9.547  -7.461  -6.471  -5.523  -3.635
## delta[2]    -4.214   0.067 1.278  -6.853  -4.983  -4.240  -3.439  -1.594
## delta[3]    -2.835   0.067 1.278  -5.503  -3.607  -2.843  -2.051  -0.203
## delta[4]     2.082   0.067 1.340  -0.681   1.230   2.022   2.922   4.875
## delta[5]    -3.179   0.067 1.306  -5.749  -3.966  -3.208  -2.353  -0.535
## delta[6]     0.394   0.067 1.280  -2.246  -0.373   0.367   1.178   3.042
## delta[7]     4.074   0.068 1.341   1.423   3.221   4.055   4.931   6.698
## delta[8]     6.889   0.069 1.450   3.903   5.941   6.922   7.811   9.814
## delta[9]    -0.297   0.066 1.323  -2.994  -1.088  -0.315   0.517   2.556
## delta[10]    5.016   0.068 1.432   2.208   4.081   4.994   5.943   7.924
## delta[11]    1.868   0.068 2.320  -2.500   0.268   1.837   3.377   6.458
## delta[12]   -2.856   0.070 1.420  -5.714  -3.784  -2.857  -1.918  -0.057
## delta[13]   -0.033   0.098 4.334  -8.904  -2.712  -0.021   2.862   8.438
## sigma        5.919   0.002 0.097   5.731   5.854   5.917   5.982   6.113
## sigma_age    1.775   0.013 0.444   1.129   1.453   1.692   2.032   2.812
## sigma_site   4.226   0.035 0.893   2.878   3.598   4.077   4.704   6.399
## sigma_tot    7.533   0.022 0.559   6.753   7.147   7.429   7.806   8.953
## pi[1]        0.627   0.003 0.083   0.440   0.575   0.637   0.686   0.763
## pi[2]        0.059   0.001 0.029   0.023   0.038   0.051   0.072   0.128
## pi[3]        0.315   0.003 0.086   0.176   0.254   0.302   0.367   0.515
## difference  10.063   0.006 0.247   9.557   9.899  10.071  10.236  10.510
## flower_new  95.413   0.221 9.271  77.224  89.225  95.456 101.818 113.336
## earliest    64.725   0.044 2.005  60.002  63.653  64.964  66.117  67.886
## latest     122.567   0.051 2.207 119.121 121.007 122.298 123.805 127.743
## R_squared    0.099   0.000 0.001   0.097   0.098   0.099   0.099   0.100
##            n_eff  Rhat
## proba       1779 1.002
## beta[1]      408 1.007
## beta[2]      390 1.008
## alpha[1]    1260 1.001
## alpha[2]     458 1.005
## alpha[3]     427 1.004
## alpha[4]     444 1.005
## alpha[5]     423 1.003
## alpha[6]     429 1.006
## alpha[7]     440 1.004
## alpha[8]     439 1.004
## alpha[9]     518 1.006
## alpha[10]    622 1.005
## alpha[11]    846 1.001
## alpha[12]    871 1.001
## alpha[13]   1068 1.003
## alpha[14]    991 1.001
## alpha[15]   1624 1.001
## delta[1]     512 1.005
## delta[2]     366 1.008
## delta[3]     361 1.008
## delta[4]     397 1.008
## delta[5]     376 1.009
## delta[6]     361 1.009
## delta[7]     387 1.009
## delta[8]     447 1.005
## delta[9]     406 1.008
## delta[10]    439 1.006
## delta[11]   1156 1.001
## delta[12]    408 1.008
## delta[13]   1959 1.001
## sigma       2166 0.999
## sigma_age   1198 1.000
## sigma_site   653 1.003
## sigma_tot    638 1.003
## pi[1]        700 1.002
## pi[2]       1285 1.000
## pi[3]        693 1.003
## difference  1660 1.000
## flower_new  1765 1.001
## earliest    2112 0.999
## latest      1884 1.000
## R_squared    975 1.002
## 
## Samples were drawn using NUTS(diag_e) at Mon Apr 15 15:44:43 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Est-ce que ce modèle est meilleur que le précédent?

loo::compare(loo::waic(rstan::extract(fit_1, "log_lik")$log_lik),
             loo::waic(rstan::extract(fit_2, "log_lik")$log_lik)
             )
## elpd_diff        se 
##      46.3       9.7
  • Que remarquez-vous sur l’estimation de \(\hat{\sigma}\)?
  • Que remarquez-vous sur l’estimation de \(\hat{\sigma}_{\mathrm{site}}\)?
summary_anova <- as.data.frame(
  do.call('rbind', lapply(c("sigma", "sigma_age", "sigma_site"),
                          function(param) {
                            get_summary(as.numeric(rstan::extract(fit_2, param)[[1]]))
                            }
                          )
          )
)

names(summary_anova) <- c("mean", "se", "lower", "upper")
summary_anova$component <- c("residual", "age", "site")

summary_anova %>% 
  mutate(component = factor(component, levels = c("residual", "age", "site")[order(mean)])) %>% 
  ggplot(aes(x = component, y = mean)) +
  geom_linerange(aes(x = component, ymin = lower, ymax = upper)) +
  geom_point(size = 2) +
  ylab("Estimate") + xlab("Source of variation") +
  coord_flip() +
  theme_bw()

Inclure des covariables

Il est également possible d’inclure des prédicteurs comme le cultivar (\(190\) différents). Cela revient à faire une régression logistique sur la variable latente \(z_{ijk}\) :

\[z_{ijk} \sim \mathcal{B}(\mathrm{proba}_i)\]

\[\mathrm{logit}(\mathrm{proba}_i) = \eta_0 + \eta_{\mathrm{Cultivar}_i}\]

Cela revient à ajuster simulatément deux GLMs.

stan_mix2 <- '
 functions {
  real TwoGaussianMixture_lpdf(real y, real prob, vector location, real scale) {
        real log_pdf[2];
    log_pdf[1] = log1m(prob) + normal_lpdf(y| location[1], scale);
    log_pdf[2] = log(prob) + normal_lpdf(y| location[2], scale);
    return log_sum_exp(log_pdf);
  }
  real TwoGaussianMixture_rng(real prob, vector location, real scale) {
   int z;
   z = bernoulli_rng(prob);
   return  z ? normal_rng(location[2], scale) : normal_rng(location[1], scale); 
  }
 }

 data {
  int<lower = 1> n_obs;
  int<lower = 1> n_age;
  int<lower = 1> n_site;
  int<lower = 1> n_cultivar;
  vector[n_obs] FLOWERING;
  int<lower = 1, upper = n_age> AGE[n_obs];
  int<lower = 1, upper = n_site> SITE[n_obs];
  int<lower = 1, upper = n_cultivar> CULTIVAR[n_obs];
  real prior_location_beta0;
  real<lower = 0.0> prior_scale_beta0;
  real prior_location_diff;
  real<lower = 0.0> prior_scale_diff;
  real prior_location_eta0;
  real<lower = 0.0> prior_scale_eta0;
 }
 
 parameters {
  simplex[3] pi;
  real beta0;
  real<lower = 0> difference;
  real<lower = 0.0> sigma_tot;
  real<lower = 0.0> sigma_cultivar;
  vector[n_age + 1] alpha;
  vector[n_site + 1] delta;
  real eta0;
  vector[n_cultivar] eta;
 }

 transformed parameters {
  real R_squared;
  real sigma;
  real sigma_age;
  real sigma_site;
  vector[n_obs] mu[2];
  vector[n_obs] proba;
  vector[2] beta;
  beta[1] = prior_location_beta0 + beta0 * prior_scale_beta0;
  beta[2] = beta[1] + difference;
  sigma = sqrt(pi[1]) * sigma_tot;
  sigma_age = sqrt(pi[2]) * sigma_tot;
  sigma_site = sqrt(pi[3]) * sigma_tot;
  mu[1] = beta[1] + alpha[AGE] + delta[SITE];
  mu[2] = beta[2] + alpha[AGE] + delta[SITE];
  proba = inv_logit(rep_vector(eta0, n_obs) + eta[CULTIVAR]);
  R_squared =  1 - variance((1 - proba) .* mu[1] + proba .* mu[2] - FLOWERING) / variance(FLOWERING);
 }

 model {
  sigma_tot ~ student_t(3, 0.0, 1.0);
  sigma_cultivar ~ student_t(3, 0.0, 1.0);
  beta0 ~ normal(0.0, 1.0);
  difference ~ normal(prior_location_diff, prior_scale_diff);
  alpha ~ normal(0.0, sigma_age);
  delta ~ normal(0.0, sigma_site);
  eta0 ~ normal(prior_location_eta0, prior_scale_eta0);
  eta ~ normal(0.0, sigma_cultivar);
  for(i in 1:n_obs) {
   target += log_sum_exp(log1m(proba[i]) + normal_lpdf(FLOWERING[i] | mu[1, i], sigma),
   log(proba[i]) + normal_lpdf(FLOWERING[i] | mu[2, i], sigma));
  }
 }

 generated quantities {
  vector[n_obs] log_lik;
  vector[n_obs] y_rep;
  int<lower = 0, upper = 1> type;
  real flower_new;
  real new_cultivar;
  real earliest;
  real latest;
  new_cultivar = normal_rng(0.0, 1.0);
  type = bernoulli_logit_rng(eta0 + sigma_cultivar * new_cultivar);
  flower_new = normal_rng(beta[1 + type] + alpha[n_age + 1] + delta[n_site + 1], sigma);
  for(i in 1:n_obs) {
   log_lik[i] = log_sum_exp(log1m(proba[i]) + normal_lpdf(FLOWERING[i] | mu[1, i], sigma),
   log(proba[i]) + normal_lpdf(FLOWERING[i] | mu[2, i], sigma));
   y_rep[i] = TwoGaussianMixture_rng(proba[i], to_vector(mu[1:2, i]), sigma);
  }
  earliest = min(y_rep);
  latest = max(y_rep);
 }
'

## Compilation
mix2 <- stan_model(model_code = stan_mix2,
                   model_name = "Ma deuxieme ANOVA a melange avec Stan"
                   )

Ajustons ce nouveau modèle aux données :

fit_3 <- sampling(object = mix2,
                  data = list(n_obs = nrow(sakura),
                              n_age = length(unique(sakura$age)),
                              n_site = length(unique(sakura$Site)),
                              n_cultivar = length(unique(sakura$Cultivar)),
                              FLOWERING = sakura$Full_Flowering,
                              AGE = sakura$age,
                              SITE = as.numeric(factor(sakura$Site, levels = unique(sakura$Site))),
                              CULTIVAR = as.numeric(factor(sakura$Cultivar, levels = unique(sakura$Cultivar))),
                              prior_location_beta0 = 90,
                              prior_scale_beta0 = 10,
                              prior_location_diff = 7,
                              prior_scale_diff = 3,
                              prior_location_eta0 = 0.0,
                              prior_scale_eta0 = 1.5
                              ), 
                  pars = c("eta0", "eta", "sigma_cultivar", "beta", "alpha", "delta", "sigma", "sigma_age", "sigma_site", "sigma_tot", "pi", "difference", "flower_new", "earliest", "latest", "R_squared", "log_lik"),
                  chains = n_chains, 
                  iter = n_iter, 
                  warmup = n_warm, 
                  thin = n_thin
                  )

get_elapsed_time(fit_3)
##          warmup  sample
## chain:1 1394.99 558.599
## chain:2 1316.44 399.479
## chain:3 1319.65 428.503
## chain:4 1424.87 441.154

Vérifions la convergence des paramètres …

stan_rhat(fit_3)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(fit_3, diguts = 3, pars = c("eta0", "sigma_cultivar", "beta", "alpha", "delta", "sigma", "sigma_age", "sigma_site", "sigma_tot", "pi", "difference", "flower_new", "earliest", "latest", "R_squared"))
## Inference for Stan model: Ma deuxieme ANOVA a melange avec Stan.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##                  mean se_mean   sd   2.5%    25%    50%    75%  97.5%
## eta0            -0.60    0.02 0.28  -1.18  -0.78  -0.59  -0.41  -0.07
## sigma_cultivar   2.99    0.03 0.37   2.36   2.74   2.96   3.21   3.83
## beta[1]         91.93    0.11 1.36  89.21  91.02  91.90  92.81  94.59
## beta[2]        100.77    0.11 1.36  98.16  99.85 100.78 101.66 103.44
## alpha[1]        -0.78    0.03 0.90  -2.66  -1.33  -0.77  -0.20   0.89
## alpha[2]         1.22    0.04 0.51   0.25   0.88   1.22   1.54   2.23
## alpha[3]         0.55    0.04 0.49  -0.34   0.23   0.53   0.86   1.55
## alpha[4]         1.53    0.03 0.49   0.60   1.22   1.51   1.85   2.54
## alpha[5]         0.43    0.03 0.49  -0.47   0.11   0.41   0.76   1.43
## alpha[6]         0.54    0.03 0.50  -0.38   0.20   0.54   0.87   1.56
## alpha[7]        -1.28    0.03 0.50  -2.21  -1.62  -1.29  -0.96  -0.24
## alpha[8]        -1.13    0.03 0.51  -2.09  -1.48  -1.15  -0.79  -0.11
## alpha[9]         0.10    0.03 0.56  -0.94  -0.28   0.09   0.46   1.25
## alpha[10]       -2.53    0.03 0.61  -3.74  -2.93  -2.54  -2.12  -1.32
## alpha[11]       -0.99    0.03 0.74  -2.47  -1.46  -0.99  -0.50   0.47
## alpha[12]        1.77    0.03 0.80   0.28   1.22   1.76   2.29   3.42
## alpha[13]        0.76    0.03 0.77  -0.65   0.25   0.74   1.26   2.35
## alpha[14]        0.57    0.03 0.84  -1.00   0.01   0.56   1.13   2.32
## alpha[15]        0.03    0.03 1.58  -2.97  -0.96   0.00   1.05   3.18
## delta[1]        -5.91    0.10 1.46  -8.91  -6.87  -5.90  -4.93  -3.17
## delta[2]        -4.29    0.10 1.29  -6.96  -5.09  -4.25  -3.44  -1.78
## delta[3]        -2.85    0.10 1.29  -5.58  -3.68  -2.83  -1.98  -0.37
## delta[4]         2.13    0.10 1.35  -0.69   1.30   2.16   3.05   4.70
## delta[5]        -2.56    0.10 1.32  -5.34  -3.39  -2.53  -1.69  -0.01
## delta[6]         0.04    0.10 1.29  -2.64  -0.76   0.06   0.89   2.54
## delta[7]         4.39    0.10 1.37   1.69   3.48   4.45   5.27   7.10
## delta[8]         6.80    0.10 1.44   3.95   5.86   6.87   7.75   9.45
## delta[9]        -0.41    0.10 1.33  -3.25  -1.26  -0.37   0.49   2.17
## delta[10]        5.08    0.10 1.41   2.17   4.17   5.10   6.01   7.86
## delta[11]        1.71    0.09 2.47  -2.93  -0.01   1.64   3.37   6.62
## delta[12]       -2.72    0.10 1.39  -5.56  -3.64  -2.65  -1.76  -0.05
## delta[13]        0.00    0.09 4.37  -8.38  -2.77  -0.02   2.64   8.72
## sigma            6.36    0.00 0.08   6.18   6.30   6.36   6.42   6.52
## sigma_age        1.51    0.01 0.38   0.94   1.23   1.45   1.72   2.36
## sigma_site       4.26    0.05 0.98   2.77   3.59   4.11   4.71   6.44
## sigma_tot        7.85    0.03 0.59   7.07   7.46   7.74   8.10   9.21
## pi[1]            0.67    0.00 0.08   0.48   0.62   0.67   0.73   0.80
## pi[2]            0.04    0.00 0.02   0.01   0.03   0.03   0.05   0.09
## pi[3]            0.29    0.00 0.09   0.15   0.23   0.28   0.34   0.49
## difference       8.84    0.01 0.23   8.37   8.68   8.85   9.00   9.28
## flower_new      95.93    0.20 8.99  79.03  89.47  95.91 102.38 113.46
## earliest        65.26    0.05 2.03  60.53  64.02  65.49  66.69  68.57
## latest         124.50    0.05 2.32 120.80 122.86 124.16 125.84 129.91
## R_squared        0.27    0.00 0.00   0.27   0.27   0.27   0.27   0.28
##                n_eff Rhat
## eta0             285 1.01
## sigma_cultivar   198 1.01
## beta[1]          165 1.03
## beta[2]          164 1.03
## alpha[1]        1020 1.00
## alpha[2]         207 1.03
## alpha[3]         195 1.03
## alpha[4]         212 1.03
## alpha[5]         211 1.03
## alpha[6]         211 1.03
## alpha[7]         229 1.03
## alpha[8]         284 1.02
## alpha[9]         325 1.02
## alpha[10]        312 1.02
## alpha[11]        744 1.01
## alpha[12]        664 1.01
## alpha[13]        700 1.01
## alpha[14]        892 1.01
## alpha[15]       2307 1.00
## delta[1]         222 1.01
## delta[2]         173 1.01
## delta[3]         171 1.02
## delta[4]         188 1.01
## delta[5]         182 1.01
## delta[6]         173 1.01
## delta[7]         190 1.01
## delta[8]         210 1.01
## delta[9]         182 1.01
## delta[10]        214 1.01
## delta[11]        783 1.00
## delta[12]        205 1.01
## delta[13]       2331 1.00
## sigma            333 1.01
## sigma_age        802 1.00
## sigma_site       402 1.01
## sigma_tot        377 1.01
## pi[1]            445 1.01
## pi[2]            690 1.00
## pi[3]            427 1.01
## difference       355 1.00
## flower_new      2069 1.00
## earliest        2017 1.00
## latest          1792 1.00
## R_squared        536 1.00
## 
## Samples were drawn using NUTS(diag_e) at Mon Apr 15 16:18:59 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

… avant de regarder si ce modèle est meilleur

loo::compare(loo::waic(rstan::extract(fit_2, "log_lik")$log_lik),
             loo::waic(rstan::extract(fit_3, "log_lik")$log_lik)
             )
## elpd_diff        se 
##     884.6      37.7

Maintenant, tracer la relation entre cultivar et floraison précoce.

proba_late <- plogis(matrix(rep(extract(fit_3, "eta0")$eta0, each = length(unique(sakura$Cultivar))), byrow = FALSE, ncol = length(unique(sakura$Cultivar))) + extract(fit_3, "eta")$eta)

data.frame(id = unique(sakura$Cultivar),
                       proba = apply(proba_late, 2, mean),
                       lower = apply(proba_late, 2, lower),
                       upper = apply(proba_late, 2, upper)
                       ) %>% 
  mutate(id = factor(id, levels = id[order(proba)])) %>% 
  ggplot(aes(x = id, y = proba)) +
  geom_linerange(aes(x = id, ymin = lower, ymax = upper)) +
  geom_point() +
  xlab("cultivar") + ylab("Pr(Late Flowering)") +
  coord_flip() +
  theme(axis.text.y = element_text(size = 6))