Modèle Hiérarchique avec Stan
Matthieu Authier & Eric Parent
19 avril 2019
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.1On 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.8ou 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.2Un 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.477La 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.24211282apply(rstan::extract(fit_1, "pi")$pi, 2, sd)## [1] 0.07456004 0.01787135 0.07645566Cela 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.154Vé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))