# Cargar los paquetes
library(bayesrules)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rstan)
## Loading required package: StanHeaders
##
## rstan version 2.32.6 (Stan version 2.32.2)
##
## 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 within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
##
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
##
## Attaching package: 'rstan'
##
## The following object is masked from 'package:tidyr':
##
## extract
library(rstanarm)
## Loading required package: Rcpp
## This is rstanarm version 2.32.1
## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
## - For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores())
##
## Attaching package: 'rstanarm'
##
## The following object is masked from 'package:rstan':
##
## loo
library(bayesplot)
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
library(tidybayes)
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(broom.mixed)
library(ggplot2)
library(GGally)
## Warning: package 'GGally' was built under R version 4.4.2
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(corrplot)
## corrplot 0.95 loaded
library(patchwork)
#Cargar los datos
data(weather_perth)
#nombre de las variables
weather_perth %>%
names()
## [1] "mintemp" "maxtemp" "rainfall" "windgustdir"
## [5] "windgustspeed" "winddir9am" "winddir3pm" "windspeed9am"
## [9] "windspeed3pm" "humidity9am" "humidity3pm" "pressure9am"
## [13] "pressure3pm" "temp9am" "temp3pm" "raintoday"
## [17] "risk_mm" "raintomorrow" "year" "month"
## [21] "day_of_year"
#Consultamos descripción de la base de datos
?weather_perth
## starting httpd help server ... done
#Estructura interna del objeto "weather_perth"
str(weather_perth)
## tibble [1,000 × 21] (S3: tbl_df/tbl/data.frame)
## $ mintemp : num [1:1000] 19.6 18.4 11.7 13 17.8 9.7 17.9 9.7 11.2 12.9 ...
## $ maxtemp : num [1:1000] 32.3 34.7 19.8 24.9 28.5 21.6 31 27.4 17.9 19.6 ...
## $ rainfall : num [1:1000] 0 0 0.6 0 0 2.2 0 0 9 13.6 ...
## $ windgustdir : Ord.factor w/ 16 levels "N"<"NNE"<"NE"<..: 5 4 14 11 10 11 10 11 14 13 ...
## $ windgustspeed: num [1:1000] 31 37 48 44 39 24 30 33 52 31 ...
## $ winddir9am : Ord.factor w/ 16 levels "N"<"NNE"<"NE"<..: 5 5 14 9 9 3 5 4 13 2 ...
## $ winddir3pm : Ord.factor w/ 16 levels "N"<"NNE"<"NE"<..: 11 7 14 10 11 13 12 11 13 13 ...
## $ windspeed9am : num [1:1000] 7 17 15 9 15 7 7 9 19 9 ...
## $ windspeed3pm : num [1:1000] 15 15 15 24 20 11 13 17 20 15 ...
## $ humidity9am : int [1:1000] 55 43 62 53 65 84 48 51 47 90 ...
## $ humidity3pm : int [1:1000] 44 16 85 46 51 59 43 34 47 57 ...
## $ pressure9am : num [1:1000] 1010 1011 1014 1019 1013 ...
## $ pressure3pm : num [1:1000] 1007 1006 1012 1016 1011 ...
## $ temp9am : num [1:1000] 25.3 24.6 17.4 19.7 24.2 15.4 25.3 21 16.7 15.1 ...
## $ temp3pm : num [1:1000] 30 33.7 16.5 22.8 26.8 20.3 28.1 25.8 17.5 18.1 ...
## $ raintoday : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 2 1 1 2 2 ...
## $ risk_mm : num [1:1000] 0 0 3.4 0 0 0 0 0 5 4.8 ...
## $ raintomorrow : Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 2 2 ...
## $ year : num [1:1000] 2014 2011 2013 2014 2011 ...
## $ month : num [1:1000] 2 1 9 12 12 9 12 10 10 8 ...
## $ day_of_year : num [1:1000] 45 11 261 347 357 254 364 293 304 231 ...
## - attr(*, "na.action")= 'omit' Named int [1:211] 1 3 9 19 37 38 51 55 65 87 ...
## ..- attr(*, "names")= chr [1:211] "1" "3" "9" "19" ...
#Observamos las primeras 3 filas de los datos
head(weather_perth,3)
## # A tibble: 3 × 21
## mintemp maxtemp rainfall windgustdir windgustspeed winddir9am winddir3pm
## <dbl> <dbl> <dbl> <ord> <dbl> <ord> <ord>
## 1 19.6 32.3 0 E 31 E SW
## 2 18.4 34.7 0 ENE 37 E SE
## 3 11.7 19.8 0.6 WNW 48 WNW WNW
## # ℹ 14 more variables: windspeed9am <dbl>, windspeed3pm <dbl>,
## # humidity9am <int>, humidity3pm <int>, pressure9am <dbl>, pressure3pm <dbl>,
## # temp9am <dbl>, temp3pm <dbl>, raintoday <fct>, risk_mm <dbl>,
## # raintomorrow <fct>, year <dbl>, month <dbl>, day_of_year <dbl>
#Summary de los datos
summary(weather_perth)
## mintemp maxtemp rainfall windgustdir windgustspeed
## Min. : 0.3 Min. :13.90 Min. : 0.000 SW :241 Min. :13.00
## 1st Qu.: 9.7 1st Qu.:20.70 1st Qu.: 0.000 SSW :150 1st Qu.:30.00
## Median :13.7 Median :24.60 Median : 0.000 NE : 96 Median :35.00
## Mean :13.3 Mean :25.58 Mean : 1.495 WSW : 88 Mean :35.07
## 3rd Qu.:17.0 3rd Qu.:30.12 3rd Qu.: 0.200 ENE : 50 3rd Qu.:41.00
## Max. :29.7 Max. :42.70 Max. :51.200 ESE : 44 Max. :76.00
## (Other):331
## winddir9am winddir3pm windspeed9am windspeed3pm humidity9am
## E :151 SW :242 Min. : 2.00 Min. : 2.00 Min. :13.00
## NE :121 WSW :135 1st Qu.: 7.00 1st Qu.:11.00 1st Qu.:49.00
## ENE :112 SSW :113 Median :11.00 Median :15.00 Median :59.00
## SSE : 86 W : 99 Mean :11.22 Mean :14.59 Mean :60.87
## ESE : 79 E : 51 3rd Qu.:15.00 3rd Qu.:19.00 3rd Qu.:73.00
## NNE : 73 SE : 49 Max. :28.00 Max. :28.00 Max. :99.00
## (Other):378 (Other):311
## humidity3pm pressure9am pressure3pm temp9am temp3pm
## Min. : 8.0 Min. :1001 Min. :1000 Min. : 7.30 Min. :11.2
## 1st Qu.:35.0 1st Qu.:1013 1st Qu.:1011 1st Qu.:15.20 1st Qu.:19.5
## Median :46.0 Median :1017 Median :1015 Median :19.10 Median :23.1
## Mean :45.8 Mean :1018 Mean :1015 Mean :19.09 Mean :24.0
## 3rd Qu.:56.0 3rd Qu.:1022 3rd Qu.:1019 3rd Qu.:22.80 3rd Qu.:28.2
## Max. :99.0 Max. :1038 Max. :1034 Max. :39.00 Max. :41.8
##
## raintoday risk_mm raintomorrow year month
## No :816 Min. : 0.000 No :814 Min. :2008 Min. : 1.000
## Yes:184 1st Qu.: 0.000 Yes:186 1st Qu.:2011 1st Qu.: 3.000
## Median : 0.000 Median :2013 Median : 7.000
## Mean : 1.634 Mean :2014 Mean : 6.601
## 3rd Qu.: 0.200 3rd Qu.:2017 3rd Qu.:10.000
## Max. :47.200 Max. :2020 Max. :12.000
##
## day_of_year
## Min. : 1.0
## 1st Qu.: 86.0
## Median :194.5
## Mean :185.6
## 3rd Qu.:283.0
## Max. :365.0
##
#seleccionamos solo variables numéricas para ver correlaciones
weather_perth_varnum <- select(weather_perth, -windgustdir,-winddir9am,
-winddir3pm,-raintoday,-raintomorrow)
#visualización de correlaciones
M <- cor(weather_perth_varnum )
corrplot(M, method = "number")
#Matriz de correlaciones y visualizaciones conjuntamente con variable categórica raintoday
ggpairs(weather_perth, columns = c("day_of_year","rainfall" ,"windgustspeed",
"windspeed9am", "humidity9am","pressure9am", "temp9am","risk_mm",
"temp3pm"),
title = "Relación entre variables numéricas y conjuntamente con la variable categórica raintoday (si ha llovido o no)",
upper = list(continuous = wrap("cor",
size = 3)),
lower = list(continuous = wrap("smooth",
alpha = 0.3,
size = 0.1)),
mapping = aes(color = raintoday))
#Observamos boxplots de temp3pm condicionado a las variables categóricas
#raintomorrow y winddir9am respecto a raintoday
ggplot(weather_perth, aes(y = temp3pm, x = raintomorrow, fill = raintoday )) +
geom_boxplot() +
labs(title = "Boxplots de temp3pm por categorías de raintoday y raintomorrow")
ggplot(weather_perth, aes(y = temp3pm, x = winddir9am, fill = raintoday )) +
geom_boxplot() +
labs(title = "Boxplots de temp3pm por categorías de raintoday y winddir9am")
#Nuestra variable independiente será temp3pm. Propondremos modelos para explicarla.
#Luego simulamos los modelos mediante MCMC
##Establecemos los parámetros del modelo normal para β0c y visualizamos
#media y desviación estándar para los datos de la variable temp3pm
mean(weather_perth$temp3pm)
## [1] 24.0009
sd(weather_perth$temp3pm)
## [1] 5.829782
#visualización de densidad de temp3pm y modelo normal propuesto para β0c
curve(dnorm(x, 24, 5.8),5,45, # density plot
xlab = "temp3pm",
ylab = "Density",
main = "Densidad de temp3pm vs Normal utilizada para β0c",
lwd = 4, # thickness of line
col = "orange",
ylim = c(0,0.07))
# Densidad streams
lines(density(weather_perth$temp3pm), # density plot
lwd = 4,
col = "black")
# Legend
legend(x ="topright", # Ubicación de la leyenda
c("Densidad de temp3pm","Normal utilizada para β0c"),
col = c("black","orange"),
lwd = c(2,2,2,2),
bty = "n")
############# MODELO 1 PROPUESTO ######################
##Observamos relación entre temp9am y temp3pm e incidencia de
#las variables categóricas raintoday y raintomorrow
#ggplot para temp9am y temp3pm (Regresión lineal simple)
#y = temp3pm, x = temp9am
ggplot(weather_perth, aes(x = temp9am, y = temp3pm)) +
geom_point(size = 1.5 ,color = "black") +
geom_smooth(formula = 'y ~ x',method = "lm", se = TRUE) +
labs(title = "temp9am vs temp3pm")+
#y = temp3pm, x = temp9am, color = raintoday (2 rectas de regresión lineal simple)
ggplot(weather_perth, aes(y = temp3pm, x = temp9am, color = raintoday)) +
geom_point(size = 1.5) +
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE) +
labs(title = "temp9am vs temp3pm por categorías de raintoday")
#Densidades de "temp3pm" según variables categórica raintoday
ggplot(weather_perth, aes(x = temp3pm, fill = raintoday)) +
geom_density(alpha = 0.5) +
labs(title = "Densidades de temp3pm por categorías de raintoday")
#Modelo 1: Elegimos quedarnos con temp9am y raintoday para explicar temp3pm
#Simulación de modelo 1 a priori
weather_model_1_prior <- stan_glm(
temp3pm ~ temp9am + raintoday,
data = weather_WU, family = gaussian,
prior_intercept = normal(24, 5.8),
prior = normal(0, 2.5, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4, iter = 5000*2, seed = 131184,
prior_PD = TRUE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001722 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 17.22 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.137 seconds (Warm-up)
## Chain 1: 0.144 seconds (Sampling)
## Chain 1: 0.281 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.14 seconds (Warm-up)
## Chain 2: 0.138 seconds (Sampling)
## Chain 2: 0.278 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.6e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.167 seconds (Warm-up)
## Chain 3: 0.173 seconds (Sampling)
## Chain 3: 0.34 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.141 seconds (Warm-up)
## Chain 4: 0.208 seconds (Sampling)
## Chain 4: 0.349 seconds (Total)
## Chain 4:
prior_summary(weather_model_1_prior)
## Priors for model 'weather_model_1_prior'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 24, scale = 5.8)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0], scale = [2.5,2.5])
## Adjusted prior:
## ~ normal(location = [0,0], scale = [ 3.09,57.28])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.13)
## ------
## See help('prior_summary.stanreg') for more details
set.seed(131184)
weather_perth %>%
add_predicted_draws(weather_model_1_prior, ndraws = 100) %>%
ggplot(aes(x = .prediction, group = .draw)) +
geom_density() +
xlab("temp3pm") +
lims(x = c(-250, 250)) +
labs(title = "M1: Simulación a priori 100 densidades para temp3pm") +
weather_perth %>%
add_epred_draws(weather_model_1_prior, ndraws = 100) %>%
ggplot(aes(x = temp9am, y = temp3pm, color = raintoday)) +
geom_line(aes(y = .epred, group = paste(raintoday, .draw))) +
labs(title = "M1: Simulación a priori 100 rectas para temp9am vs temp3pm según raintoday")
#Simulación de modelo 1 a posteriori
weather_model_1_posterior <- update(weather_model_1_prior, prior_PD = FALSE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000604 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.04 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.191 seconds (Warm-up)
## Chain 1: 0.261 seconds (Sampling)
## Chain 1: 0.452 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.2e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.191 seconds (Warm-up)
## Chain 2: 0.25 seconds (Sampling)
## Chain 2: 0.441 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.184 seconds (Warm-up)
## Chain 3: 0.251 seconds (Sampling)
## Chain 3: 0.435 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.185 seconds (Warm-up)
## Chain 4: 0.253 seconds (Sampling)
## Chain 4: 0.438 seconds (Total)
## Chain 4:
#3 primeras filas del dataframe de valores simulados
weather_model_1_posterior_df <- as.data.frame(weather_model_1_posterior)
head(weather_model_1_posterior_df , 3)
## (Intercept) temp9am raintodayYes sigma
## 1 6.333133 0.9506423 -2.404892 3.853199
## 2 3.264444 1.0805059 -2.478352 4.550963
## 3 4.204468 1.0390663 -2.212055 4.195689
#Diagnósticos de la simulación MCMC del modelo 1 a posteriori
# Effective sample size ratio and Rhat
neff_ratio(weather_model_1_posterior)
## (Intercept) temp9am raintodayYes sigma
## 1.22960 1.24525 1.20235 1.18640
rhat(weather_model_1_posterior)
## (Intercept) temp9am raintodayYes sigma
## 0.9999531 0.9999279 1.0002026 0.9999720
#Autocorrelation
mcmc_acf(weather_model_1_posterior)
# Trace plots of parallel chains
mcmc_trace(weather_model_1_posterior, size = 0.1)
# Density plots of parallel chains
mcmc_dens_overlay(weather_model_1_posterior)
#Porcentaje de valores de β1>0
weather_model_1_posterior_df %>%
mutate(exceeds_0 = temp9am > 0) %>%
tabyl(exceeds_0)
## exceeds_0 n percent
## TRUE 20000 1
#simulación de densidades y rectas de regresión lineal (modelo 1 a posteriori)
set.seed(131184)
weather_perth %>%
add_predicted_draws(weather_model_1_posterior, ndraws = 100) %>%
ggplot(aes(x = .prediction, group = .draw)) +
geom_density() +
xlab("temp3pm") +
lims(x = c(-10, 60)) +
labs(title = "M1: Simulación a posteriori 100 densidades para temp3pm") +
weather_perth %>%
add_epred_draws(weather_model_1_posterior, ndraws = 100) %>%
ggplot(aes(x = temp9am, y = temp3pm, color = raintoday)) +
geom_line(aes(y = .epred, group = paste(raintoday, .draw))) +
labs(title = "M1: Simulación a posteriori 100 rectas para temp9am vs temp3pm según raintoday")
#Summary del modelo 1 a posteriori
tidy(weather_model_1_posterior, effects = c("fixed", "aux"))
## # A tibble: 5 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 (Intercept) 4.62 0.996
## 2 temp9am 1.03 0.0476
## 3 raintodayYes -2.01 0.891
## 4 sigma 4.07 0.210
## 5 mean_PPD 24.5 0.411
#Intervalos de credibilidad a posteriori del 80% para el modelo 1
posterior_interval(weather_model_1_posterior, prob = 0.80,
pars = c("(Intercept)","temp9am", "raintodayYes","sigma"))
## 10% 90%
## (Intercept) 3.3460457 5.8800773
## temp9am 0.9689882 1.0908342
## raintodayYes -3.1368544 -0.8502729
## sigma 3.8176291 4.3511361
#Simulación para modelo 1 de conjunto de predicciones para 15 grados en temp9am
set.seed(131184)
temp3pm_prediction <- posterior_predict(
weather_model_1_posterior,
newdata = data.frame(temp9am = c(15, 15),
raintoday = c("No", "Yes")))
#Plot de modelo 1 de predicción a posteriori
mcmc_areas(temp3pm_prediction) +
ggplot2::scale_y_discrete(labels = c("Rain Today: No", "Rain Today: Yes")) +
xlab("temp3pm")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
############# MODELO 2 PROPUESTO ######################
#ggplot para humidity9am y temp3pm (Regresión lineal simple)
#y = temp3pm, x = humidity9am
ggplot(weather_perth, aes(x = humidity9am, y = temp3pm)) +
geom_point(size = 1.5 ,color = "black") +
geom_smooth(formula = 'y ~ x',method = "lm", se = TRUE) +
labs(title = "humidity9am vs temp3pm") +
#y = temp3pm, x = humidity9am, color = raintoday (2 rectas de regresión lineal simple)
ggplot(weather_perth, aes(y = temp3pm, x = humidity9am, color = raintoday)) +
geom_point(size = 1.5) +
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE) +
labs(title = "humidity9am vs temp3pm por categorías de raintoday")
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE)
## geom_smooth: na.rm = FALSE, orientation = NA, se = FALSE
## stat_smooth: na.rm = FALSE, orientation = NA, se = FALSE, method = lm, formula = y ~ x
## position_identity
#Scatterplot y = temp3pm, x = temp9am con variación de color según humidity9am
ggplot(weather_perth, aes(y = temp3pm, x = temp9am, color = humidity9am)) +
geom_point(size = 4) +
labs(title = "temp9am vs temp3pm según el nivel de humidity9am") +
#ggplot y = temp3pm, x = temp9am con variable humidity9am discretizada en 2 niveles ("low","high")
ggplot(weather_perth,
aes(y = temp3pm, x = temp9am,
color = cut(humidity9am, 2, labels = c("low", "high")))) +
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE, lwd = 1.5) +
labs(color = "humidity level 9 am") +
lims(y = c(10, 45)) +
labs(title = "temp9am vs temp3pm según 2 niveles para humidity9am")
#Modelo 2: Elegimos quedarnos con temp9am, raintoday, humidity9am
#y la interacción raintoday:humidity9am para explicar temp3pm
#Simulación de modelo 2 a priori
weather_model_2_prior <- stan_glm(
temp3pm ~ temp9am + raintoday + humidity9am + raintoday:humidity9am,
data = weather_WU, family = gaussian,
prior_intercept = normal(24, 5.8),
prior = normal(0, 2.5, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4, iter = 5000*2, seed = 131184,
prior_PD = TRUE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.194 seconds (Warm-up)
## Chain 1: 0.258 seconds (Sampling)
## Chain 1: 0.452 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.361 seconds (Warm-up)
## Chain 2: 0.354 seconds (Sampling)
## Chain 2: 0.715 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.242 seconds (Warm-up)
## Chain 3: 0.35 seconds (Sampling)
## Chain 3: 0.592 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.223 seconds (Warm-up)
## Chain 4: 0.184 seconds (Sampling)
## Chain 4: 0.407 seconds (Total)
## Chain 4:
prior_summary(weather_model_2_prior)
## Priors for model 'weather_model_2_prior'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 24, scale = 5.8)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
## Adjusted prior:
## ~ normal(location = [0,0,0,...], scale = [ 3.09,57.28, 0.82,...])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.13)
## ------
## See help('prior_summary.stanreg') for more details
#Simulación de modelo 2 a posteriori
weather_model_2_posterior <- update(weather_model_2_prior, prior_PD = FALSE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.404 seconds (Warm-up)
## Chain 1: 1.367 seconds (Sampling)
## Chain 1: 2.771 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.7e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.241 seconds (Warm-up)
## Chain 2: 1.419 seconds (Sampling)
## Chain 2: 2.66 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.198 seconds (Warm-up)
## Chain 3: 1.37 seconds (Sampling)
## Chain 3: 2.568 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.174 seconds (Warm-up)
## Chain 4: 1.852 seconds (Sampling)
## Chain 4: 3.026 seconds (Total)
## Chain 4:
#Diagnósticos de la simulación MCMC del modelo 2 a posteriori
# Effective sample size ratio and Rhat
neff_ratio(weather_model_2_posterior)
## (Intercept) temp9am raintodayYes
## 0.55760 0.66330 0.43425
## humidity9am raintodayYes:humidity9am sigma
## 0.59920 0.43110 0.78045
rhat(weather_model_2_posterior)
## (Intercept) temp9am raintodayYes
## 1.0001829 1.0001164 1.0002012
## humidity9am raintodayYes:humidity9am sigma
## 1.0000427 1.0002395 0.9998676
#Autocorrelation
mcmc_acf(weather_model_2_posterior)
# Trace plots of parallel chains
mcmc_trace(weather_model_2_posterior, size = 0.1)
# Density plots of parallel chains
mcmc_dens_overlay(weather_model_2_posterior)
#3 primeras filas del dataframe de valores simulados
weather_model_2_posterior_df <- as.data.frame(weather_model_2_posterior)
head(weather_model_2_posterior_df , 3)
## (Intercept) temp9am raintodayYes humidity9am raintodayYes:humidity9am
## 1 18.05686 0.7052910 -5.421261 -0.1439764 0.1038593
## 2 17.78329 0.7276913 -7.185396 -0.1448841 0.1288900
## 3 12.28486 0.8869308 -17.501799 -0.1035953 0.2457589
## sigma
## 1 3.564038
## 2 3.550448
## 3 3.555492
#Summary del modelo 2 a posteriori
tidy(weather_model_2_posterior, effects = c("fixed", "aux"))
## # A tibble: 7 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 (Intercept) 13.8 1.47
## 2 temp9am 0.850 0.0479
## 3 raintodayYes -9.94 4.85
## 4 humidity9am -0.113 0.0142
## 5 raintodayYes:humidity9am 0.138 0.0617
## 6 sigma 3.54 0.181
## 7 mean_PPD 24.5 0.357
#Intervalos de credibilidad a posteriori del 80% para el modelo 2
posterior_interval(weather_model_2_posterior, prob = 0.80,
pars = c("(Intercept)","temp9am", "raintodayYes","humidity9am","raintodayYes:humidity9am","sigma"))
## 10% 90%
## (Intercept) 11.97515100 15.72236655
## temp9am 0.78795287 0.91033545
## raintodayYes -16.07013318 -3.83367370
## humidity9am -0.13066425 -0.09458304
## raintodayYes:humidity9am 0.06029086 0.21623869
## sigma 3.32319987 3.78649264
#Simulación,a priori y a posteriori 100 densidades para temp3pm con el modelo 2
set.seed(131184)
weather_perth %>%
add_predicted_draws(weather_model_2_prior, ndraws = 100) %>%
ggplot(aes(x = .prediction, group = .draw)) +
geom_density() +
xlab("temp3pm") +
lims(x = c(-250, 250)) +
labs(title = "M2: Simulación a priori 100 densidades para temp3pm")+
weather_perth %>%
add_predicted_draws(weather_model_2_posterior, ndraws = 100) %>%
ggplot(aes(x = .prediction, group = .draw)) +
geom_density() +
xlab("temp3pm") +
lims(x = c(-10, 60)) +
labs(title = "M2: Simulación a posteriori 100 densidades para temp3pm")
## Warning: Removed 60 rows containing non-finite outside the scale range
## (`stat_density()`).
############# MODELO 3 PROPUESTO ######################
#ggplot para humidity9am y temp3pm (Regresión lineal simple)
#y = temp3pm, x = pressure9am
ggplot(weather_perth, aes(x = pressure9am, y = temp3pm)) +
geom_point(size = 1.5 ,color = "black") +
geom_smooth(formula = 'y ~ x',method = "lm", se = TRUE) +
labs(title = "pressure9am vs temp3pm") +
#y = temp3pm, x = humidity9am, color = raintoday (2 rectas de regresión lineal simple)
ggplot(weather_perth, aes(y = temp3pm, x = pressure9am, color = raintoday)) +
geom_point(size = 1.5) +
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE) +
labs(title = "pressure9am vs temp3pm por categorías de raintoday")+
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE)
#ggplot para day_of_yeary temp3pm (Regresión lineal simple)
#y = temp3pm, x = day_of_year
ggplot(weather_perth, aes(x = day_of_year, y = temp3pm)) +
geom_point(size = 1.5 ,color = "black") +
geom_smooth(formula = 'y ~ poly(x, 3)',method = "lm", se = TRUE) +
labs(title = "day_of_year vs temp3pm") +
#y = temp3pm, x = day_of_year, color = raintoday (2 rectas de regresión lineal simple)
ggplot(weather_perth, aes(y = temp3pm, x = day_of_year, color = raintoday)) +
geom_point(size = 1.5) +
geom_smooth(formula = 'y ~ poly(x, 3)',method = "lm", se = FALSE) +
labs(title = "day_of_year vs temp3pm por categorías de raintoday")
#ggplot para pressure9am y humidity9am (Regresión lineal simple)
#y = humidity9am, x = pressure9am
ggplot(weather_perth, aes(x = pressure9am, y = humidity9am)) +
geom_point(size = 1.5 ,color = "black") +
geom_smooth(formula = 'y ~ x',method = "lm", se = TRUE) +
labs(title = "pressure9am vs humidity9am") +
#y = temp3pm, x = humidity9am, color = raintoday (2 rectas de regresión lineal simple)
ggplot(weather_perth, aes(y = humidity9am, x = pressure9am, color = raintoday)) +
geom_point(size = 1.5) +
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE) +
labs(title = "pressure9am vs humidity9am por categorías de raintoday")+
geom_smooth(formula = 'y ~ x',method = "lm", se = FALSE)
#Simulación de modelo 3 a priori
weather_model_3_prior <- stan_glm(
temp3pm ~ temp9am + raintoday + humidity9am + raintoday:humidity9am
+pressure9am+poly(day_of_year, 3),
data = weather_WU, family = gaussian,
prior_intercept = normal(24, 5.8),
prior = normal(0, 2.5, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4, iter = 5000*2, seed = 131184,
prior_PD = TRUE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.231 seconds (Warm-up)
## Chain 1: 0.339 seconds (Sampling)
## Chain 1: 0.57 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.229 seconds (Warm-up)
## Chain 2: 0.214 seconds (Sampling)
## Chain 2: 0.443 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.229 seconds (Warm-up)
## Chain 3: 0.225 seconds (Sampling)
## Chain 3: 0.454 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.22 seconds (Warm-up)
## Chain 4: 0.199 seconds (Sampling)
## Chain 4: 0.419 seconds (Total)
## Chain 4:
prior_summary(weather_model_3_prior)
## Priors for model 'weather_model_3_prior'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 24, scale = 5.8)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
## Adjusted prior:
## ~ normal(location = [0,0,0,...], scale = [ 3.09,57.28, 0.82,...])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.13)
## ------
## See help('prior_summary.stanreg') for more details
#Simulación de modelo 3 a posteriori
weather_model_3_posterior <- update(weather_model_3_prior, prior_PD = FALSE)
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000131 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.379 seconds (Warm-up)
## Chain 1: 1.407 seconds (Sampling)
## Chain 1: 2.786 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.371 seconds (Warm-up)
## Chain 2: 1.658 seconds (Sampling)
## Chain 2: 3.029 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.355 seconds (Warm-up)
## Chain 3: 1.78 seconds (Sampling)
## Chain 3: 3.135 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.254 seconds (Warm-up)
## Chain 4: 1.153 seconds (Sampling)
## Chain 4: 2.407 seconds (Total)
## Chain 4:
#3 primeras filas del dataframe de valores simulados
weather_model_3_posterior_df <- as.data.frame(weather_model_3_posterior)
head(weather_model_3_posterior_df , 3)
## (Intercept) temp9am raintodayYes humidity9am pressure9am
## 1 63.667415 0.6489775 -14.43381 -0.1154947 -0.04444561
## 2 72.142314 0.7542307 -12.36045 -0.1231294 -0.05513229
## 3 -6.901602 0.7821147 -17.78366 -0.1376189 0.02283258
## poly(day_of_year, 3)1 poly(day_of_year, 3)2 poly(day_of_year, 3)3
## 1 -1.0807479 11.505310 3.898003
## 2 -0.4938564 5.344119 1.528438
## 3 -11.7968398 3.204313 3.844550
## raintodayYes:humidity9am sigma
## 1 0.1744064 3.428140
## 2 0.1858683 3.641258
## 3 0.2383241 3.557035
#Diagnósticos de la simulación MCMC del modelo 3 a posteriori
# Effective sample size ratio and Rhat
neff_ratio(weather_model_3_posterior)
## (Intercept) temp9am raintodayYes
## 0.95280 0.54760 0.56475
## humidity9am pressure9am poly(day_of_year, 3)1
## 0.69340 0.94945 0.91835
## poly(day_of_year, 3)2 poly(day_of_year, 3)3 raintodayYes:humidity9am
## 0.59875 0.99430 0.55870
## sigma
## 0.94110
rhat(weather_model_3_posterior)
## (Intercept) temp9am raintodayYes
## 0.9999159 1.0004724 1.0002745
## humidity9am pressure9am poly(day_of_year, 3)1
## 1.0001886 0.9999050 1.0001660
## poly(day_of_year, 3)2 poly(day_of_year, 3)3 raintodayYes:humidity9am
## 1.0002714 1.0000876 1.0002039
## sigma
## 0.9999059
#Autocorrelation
mcmc_acf(weather_model_3_posterior)
# Trace plots of parallel chains
mcmc_trace(weather_model_3_posterior, size = 0.1)
# Density plots of parallel chains
mcmc_dens_overlay(weather_model_3_posterior)
#Summary del modelo 3 a posteriori
tidy(weather_model_3_posterior, effects = c("fixed", "aux"))
## # A tibble: 11 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 (Intercept) 8.96 49.2
## 2 temp9am 0.840 0.0794
## 3 raintodayYes -10.4 4.79
## 4 humidity9am -0.116 0.0162
## 5 pressure9am 0.00512 0.0483
## 6 poly(day_of_year, 3)1 -2.98 3.80
## 7 poly(day_of_year, 3)2 -0.161 5.83
## 8 poly(day_of_year, 3)3 3.85 3.63
## 9 raintodayYes:humidity9am 0.144 0.0613
## 10 sigma 3.57 0.185
## 11 mean_PPD 24.5 0.359
#Simulación,a priori y a posteriori 100 densidades para temp3pm con el modelo 2
set.seed(131184)
weather_perth %>%
add_predicted_draws(weather_model_3_prior, ndraws = 100) %>%
ggplot(aes(x = .prediction, group = .draw)) +
geom_density() +
xlab("temp3pm") +
lims(x = c(-250, 250)) +
labs(title = "M3: Simulación a priori 100 densidades para temp3pm")+
weather_perth %>%
add_predicted_draws(weather_model_2_posterior, ndraws = 100) %>%
ggplot(aes(x = .prediction, group = .draw)) +
geom_density() +
xlab("temp3pm") +
lims(x = c(-10, 60)) +
labs(title = "M3: Simulación a posteriori 100 densidades para temp3pm")
## Warning: Removed 56 rows containing non-finite outside the scale range
## (`stat_density()`).
#3 primeras filas del dataframe de valores simulados
weather_model_3_posterior_df <- as.data.frame(weather_model_3_posterior)
head(weather_model_3_posterior_df , 3)
## (Intercept) temp9am raintodayYes humidity9am pressure9am
## 1 63.667415 0.6489775 -14.43381 -0.1154947 -0.04444561
## 2 72.142314 0.7542307 -12.36045 -0.1231294 -0.05513229
## 3 -6.901602 0.7821147 -17.78366 -0.1376189 0.02283258
## poly(day_of_year, 3)1 poly(day_of_year, 3)2 poly(day_of_year, 3)3
## 1 -1.0807479 11.505310 3.898003
## 2 -0.4938564 5.344119 1.528438
## 3 -11.7968398 3.204313 3.844550
## raintodayYes:humidity9am sigma
## 1 0.1744064 3.428140
## 2 0.1858683 3.641258
## 3 0.2383241 3.557035
###### Evaluación y comparación de los 3 modelos ############
#comparamos 50 simulaciones del modelo 2 a posteriori con la densidad de los datos de temp3pm
set.seed(131184)
pp_check(weather_model_1_posterior, nreps = 50) +
xlab("temp3pm") +
labs(title = "Ajuste del Modelo 1")+
pp_check(weather_model_2_posterior, nreps = 50) +
xlab("temp3pm") +
labs(title = "Ajuste del Modelo 2")+
pp_check(weather_model_3_posterior, nreps = 50) +
xlab("temp3pm") +
labs(title = "Ajuste del Modelo 3")
#Evaluación de la precisión predictiva con cross-validation
set.seed(131184)
prediction_summary_cv(model = weather_model_1_posterior, data = weather_perth, k = 10)
## $folds
## fold mae mae_scaled within_50 within_95
## 1 1 1.709597 0.5939473 0.54 0.96
## 2 2 2.101302 0.7468403 0.46 0.94
## 3 3 1.868849 0.6560954 0.53 0.95
## 4 4 2.029431 0.7094259 0.49 0.99
## 5 5 1.861783 0.6490289 0.52 0.98
## 6 6 1.940690 0.6780777 0.49 0.95
## 7 7 1.801927 0.6297113 0.52 0.96
## 8 8 2.077115 0.7321259 0.47 0.93
## 9 9 1.897787 0.6643686 0.51 0.94
## 10 10 1.910047 0.6678968 0.50 0.94
##
## $cv
## mae mae_scaled within_50 within_95
## 1 1.919853 0.6727518 0.503 0.954
prediction_summary_cv(model = weather_model_2_posterior, data = weather_perth, k = 10)
## $folds
## fold mae mae_scaled within_50 within_95
## 1 1 1.808025 0.6504433 0.53 0.97
## 2 2 1.949178 0.7116478 0.46 0.95
## 3 3 1.725510 0.6262860 0.51 0.95
## 4 4 1.887759 0.6898302 0.50 0.94
## 5 5 1.870561 0.6811117 0.50 0.94
## 6 6 1.812362 0.6515825 0.52 0.98
## 7 7 1.810221 0.6589992 0.51 0.93
## 8 8 1.838694 0.6632064 0.50 0.96
## 9 9 1.771596 0.6421914 0.51 0.95
## 10 10 1.786782 0.6401028 0.53 0.97
##
## $cv
## mae mae_scaled within_50 within_95
## 1 1.826069 0.6615401 0.507 0.954
prediction_summary_cv(model = weather_model_3_posterior, data = weather_perth, k = 10)
## $folds
## fold mae mae_scaled within_50 within_95
## 1 1 1.720374 0.6950441 0.46 0.95
## 2 2 1.819989 0.7268127 0.45 0.97
## 3 3 1.760549 0.7017057 0.49 0.96
## 4 4 1.800401 0.7140948 0.47 0.96
## 5 5 1.775009 0.7006110 0.48 0.98
## 6 6 1.570487 0.6222231 0.53 0.96
## 7 7 2.045346 0.8247590 0.40 0.93
## 8 8 1.467977 0.5965239 0.55 0.93
## 9 9 1.873665 0.7497035 0.44 0.97
## 10 10 1.633280 0.6589296 0.51 0.95
##
## $cv
## mae mae_scaled within_50 within_95
## 1 1.746708 0.6990407 0.478 0.956
#Evaluación de la precisión predictiva con expected log-predictive densities (ELPD).
set.seed(131184)
loo_1 <- loo(weather_model_1_posterior)
loo_2 <- loo(weather_model_2_posterior)
loo_3 <- loo(weather_model_3_posterior)
loo_compare(loo_1, loo_2, loo_3)
## elpd_diff se_diff
## weather_model_2_posterior 0.0 0.0
## weather_model_3_posterior -3.2 1.3
## weather_model_1_posterior -26.0 8.1