\(Models\ comparison\ trough\ Widely\ applicable\ information\ criterion\ (WAIC),\ using\ the dataset\\ with\ the\ visit\ length\ at\ the\ feeder\ when\ the\ next\ visit\ was\ greater\ than\ or\ equal\ to\ 600sc, from\ 7\ trials,\\ in\ pig\ production\)

\(The\ matrix\ to\ each\ model\ contain:\\ elpd_{waic}:\ expected\ log\ predictive\ pointwise\ density\\p_{waic}:\ effective\ number\ of\ parameters\\ waic:\ information\ ciretion\ converted\ to\ deviance\ scale=\ -2*elpd_{waic}\)

1. Model comparison visit length as the next visit is greater than 600 seconds (6256 records)

1.1. WAIC criterion

print(m1.600swo.waic)
## 
## Computed from 30000 by 6256 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -20763.4  91.5
## p_waic       166.8   5.5
## waic       41526.9 183.0
## Warning: 37 (0.6%) p_waic estimates greater than 0.4. We recommend trying loo
## instead.
print(m2.600swo.waic)
## 
## Computed from 30000 by 6256 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -20764.9  91.5
## p_waic       173.8   5.6
## waic       41529.8 183.0
## Warning: 39 (0.6%) p_waic estimates greater than 0.4. We recommend trying loo
## instead.
print(m3.600swo.waic)
## 
## Computed from 30000 by 6256 log-likelihood matrix
## 
##           Estimate    SE
## elpd_waic -20764.8  91.5
## p_waic       172.8   5.6
## waic       41529.6 183.0
## Warning: 40 (0.6%) p_waic estimates greater than 0.4. We recommend trying loo
## instead.

2. Pareto Smoothed Importance Sampling-Leave one out cross-validation (PSI-LOO)

\(Models\ comparison\ trough\ Leave-one-out,\ cross-validation\ (LOO,\ Vehtari\ and\ Gelman,\ 2015)\ using\ the\ dataset\\ with\ the\ visit\ length\ at\ the\ feeder\ when\ the\ next\ visit\ was\ greater\ than\ or\ equal\ to\ 600sc,\\ from\ 7\ trials,\ in\ pig\ production\)

\(The\ matrix\ to\ each\ model\ contain\ the\ estimates\ and\ their\ standard\ errors\ of:\\ \hat{elpd}_{loo}:\ expected\ log\ predictive\ pointwise\ density\\ \ \hat{p}_{loo}:\ effective\ number\ of\ parameters\\ looic:\ information\ ciretion\ converted\ to\ deviance\ scale=\ -2*\hat{elpd}_{loo}\)

The table show a summary of Pareto \(k\) diagnostics, whis is used to asses the reliability of the estimates. If we have:

\(**\ \mathbf{k\ge 1},\ neither\ variance\ nor\ mean\ the\ raw\ importance\ ratios\ exist.\ we\ are\ not\ able\ to\ compute\ an\ estimate\\ for\ Monte\ Carlo\ standar\ error\ (SE)\ of\ the\ \hat{elpd}_{loo} \\**\ \mathbf{0.5\le\ k\ < 1}\, the\ variance\ the\ raw\ importance\ ratios\ is\ infinite,\ but\ the\ mean\ exist. If\ \mathbf{k}\ is\ between\ \mathbf{0.5}\ and\ \mathbf{0.7}\ then\ observed\\ Monte\ Carlo\ error\ estimates\ with\ PSIS\\ **\ \mathbf{k\ \le 0.5},\ the\ distribution\ of\ raw\ importance\ ratios\ has\ finite\ variance\ and\ Monte\ Carlo\ standar\ error\ is\ small\)

\(Compare\ fitted\ models\ based\ on\ expected\ log\ pointwise\ predictive\ density\ (\hat{elpd}_{loo})\ for\ new\ data,\ the\ matrix\ contain:\\ elpd_{diff}: is\ the\ diference\ in\ elpd\ for\ two\ models,\ if\ more\ than\ two\ models\ are\ compared,\ the\ difference\ is\ computed\ relative\ to\ the\ model\\ with\ highest\ elpd\\ se_{diff}:\ standard\ error\ of\ difference\)

2.1. M1 mixed model

print(m16002.loo,plot_k = T, label_points = TRUE)
## 
## Computed from 30000 by 6256 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -20764.5  91.6
## p_loo       168.0   5.5
## looic     41528.9 183.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

2.2. M2 mixed model

print(m26002.loo,plot_k = T, label_points = TRUE)
## 
## Computed from 30000 by 6256 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -20765.2  91.5
## p_loo       174.5   5.7
## looic     41530.4 183.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

2.3. M3 mixed model

print(m36002.loo,plot_k = T, label_points = TRUE)
## 
## Computed from 30000 by 6256 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo -20765.7  91.5
## p_loo       174.8   5.7
## looic     41531.3 183.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Diagnostics for PSIS

a<-rbind(mcse_loo(m16002.loo),
      mcse_loo(m26002.loo),
      mcse_loo(m36002.loo))
rownames(a)<-c("M1", "M2", "M3")
colnames(a)<-"MCSE"
kable(a,caption = "Monte Carlo Standar Error (MCSE) by each model")%>%kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%column_spec(1, bold = T)
Monte Carlo Standar Error (MCSE) by each model
MCSE
M1 0.0837480
M2 0.0862830
M3 0.0895754 
b<-round(rbind(summary(pareto_k_values(m16002.loo)),
      summary(pareto_k_values(m26002.loo)),
      summary(pareto_k_values(m36002.loo))),4)
rownames(b)<-c("M1", "M2", "M3")
kable(b,caption = "Distribution of Pareto k estimates by each model")%>%
  kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%column_spec(1, bold = T)
Distribution of Pareto k estimates by each model
Min. 1st Qu. Median Mean 3rd Qu. Max.
M1 -0.1598 -0.0154 0.0201 0.0243 0.0587 0.3826
M2 -0.1523 -0.0086 0.0280 0.0314 0.0671 0.3793
M3 -0.1806 -0.0076 0.0275 0.0304 0.0643 0.3870 
a<-round(rbind(summary(psis_n_eff_values(m16002.loo)),
summary(psis_n_eff_values(m26002.loo)),
summary(psis_n_eff_values(m36002.loo))),3)
rownames(a)<-c("M1", "M2", "M3")
kable(a,caption = "Distribution of estimates PSIS effective sample size by each model")%>%
  kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%column_spec(1, bold = T)
Distribution of estimates PSIS effective sample size by each model
Min. 1st Qu. Median Mean 3rd Qu. Max.
M1 4332.437 28720.20 31879.99 30046.14 33098.05 36550.70
M2 4043.915 25462.44 34559.47 31721.25 38466.16 46460.92
M3 4225.012 23058.38 30816.53 29169.09 35374.85 43894.11

3. Model comparison

\(Compare\ fitted\ models\ based\ on\ expected\ log\ pointwise\ predictive\ density\ (\hat{elpd}_{waic},\ \hat{elpd}_{loo}),\\ the\ matrix\ contain:\\ elpd_{diff}: is\ the\ diference\ in\ elpd\ for\ two\ models,\ if\ more\ than\ two\ models\ are\ compared,\ the\ difference\ is\ computed\ relative\ to\ the\ model\\ with\ highest\ elpd\\ se_{diff}:\ standard\ error\ of\ difference\\ se-elpd_{waic},\ se-elpd_{loo} :\ standard \ error \ expected\ log\ predictive\ pointwise\ density\\ se-p_{waic},\ se-p_{loo}:\ standard \ error\ effective\ number\ of\ parameters\)

a<-loo_compare(m16002.loo, m26002.loo, m36002.loo)[,1:8]
#colnames(a)[c(7,8)]<-c("","")
kable(a,caption = "Model comparison LOO-CV subset distant replacement")%>%
  kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%
  column_spec(1, bold = T)%>%
  column_spec(6, background = "#AAA0CB")
Model comparison LOO-CV subset distant replacement
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
model1 0.0000000 0.0000000 -20764.46 91.55110 168.0132 5.528268 41528.91 183.1022
model2 -0.7615038 0.6900596 -20765.22 91.47788 174.4620 5.674473 41530.44 182.9558
model3 -1.2058893 0.6783456 -20765.66 91.51310 174.7821 5.683235 41531.33 183.0262 
a<-loo_compare(m1.600swo.waic,m2.600swo.waic, m3.600swo.waic)[,1:8]
kable(a,caption = "Model comparison WAIC subset distant replacement without obs")%>%
  kable_styling(bootstrap_options = c("striped"),full_width = F, position = "left",font_size = 14)%>%
  column_spec(1, bold = T)%>%
  column_spec(6, background = "#AAA0CB")
Model comparison WAIC subset distant replacement without obs
elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic se_waic
model1 0.000000 0.0000000 -20763.44 91.50668 166.7899 5.459760 41526.88 183.0134
model3 -1.354450 0.5986169 -20764.79 91.48142 172.7695 5.605754 41529.59 182.9628
model2 -1.467862 0.6479409 -20764.91 91.47858 173.8155 5.643163 41529.81 182.9572