I’m working with a set of nested models for largeish amount of data (~5000 observations), and when estimating PSIS-LOO-CV with the loo package, I typically get a small number of observations flagged as bad and sometimes an even smaller number as as very bad. How worried should I be about this? Looking closer at the flagged observations, it’s clear that they are observations with an unexpectedly high count value, even when the model uses a negative binomial to account for over-dispersion/variability. The main purpose of the analysis is prediction (not hypothesis testing), so I don’t want to remove the flagged observations and pretend they don’t exist; there is no reason to suspect the data is wrong. Also, the number of flagged observations differs among models, and for poor models (mainly those that assumes a Poisson distribution for the observation model). There are multiple data sets I’m running the models for, but here is the diagnostics from the best and worst fit models (based on loo) for one of them:
Best:
Estimate SE
elpd_loo -9758.8 122.2
p_loo 173.5 9.5
looic 19517.7 244.4
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 4646 99.0% 426
(0.5, 0.7] (ok) 30 0.6% 108
(0.7, 1] (bad) 14 0.3% 25
(1, Inf) (very bad) 1 0.0% 10
See help(‘pareto-k-diagnostic’) for details.
Worst (has a lot more problematic observations):
Estimate SE
elpd_loo -20913.2 789.2
p_loo 2054.9 154.7
looic 41826.4 1578.4
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 4424 94.3% 177
(0.5, 0.7] (ok) 123 2.6% 64
(0.7, 1] (bad) 68 1.4% 14
(1, Inf) (very bad) 76 1.6% 1
The Stan code:
data {
int<lower=1> NrOfLans;
int<lower=1> NrOfKretsar;
int<lower=1> Reports;
int<lower=0> K[Reports];
real<lower=0> A[Reports];
real logrelA[Reports];
int<lower=1> KretsNrSeqID[Reports];
int<lower=1> LansSeqID[NrOfLans];
int<lower=1> LanSeqID[Reports];
int<lower=1> KretsarSeqID[NrOfKretsar];
int<lower=1> KretsListLanSeqID[NrOfKretsar];
int<lower=0, upper=1> includealpha;//Boolean declaration
int<lower=0, upper=1> includeKeffectonalpha;//
int<lower=0, upper=1> includeLeffectonalpha;//
int<lower=0, upper=1> includephi;//
int<lower=0, upper=1> includeKeffectonphi;//
int<lower=0, upper=1> includeLeffectonphi;//
int<lower=0, upper=1> includeLdifferencevaraibilityonm;
}
parameters {
real logLm[NrOfLans];
real logKm[NrOfKretsar];
real logLalpha[includealpha && includeLeffectonalpha ? NrOfLans : 0];
real logLphi[includephi && includeLeffectonphi ? NrOfLans : 0];
real logSm;
real logSalpha[includealpha ? 1 : 0];
real logSphi[includephi ? 1 : 0];
real<lower=0> taoLm[includeLdifferencevaraibilityonm ? NrOfLans : 0];
real<lower=0> Tm;
real<lower=0> Talpha[includealpha && includeLeffectonalpha ? 1 : 0];
real<lower=0> Tphi[includephi && includeLeffectonphi ? 1 : 0];
real<lower=0> sigmataoLm[includeLdifferencevaraibilityonm ? 1 :0];
real<lower=0> taoSm;
real logKalpha[includealpha && includeKeffectonalpha ? NrOfKretsar : 0];// size of array is zero if includeKeffectonalpha is false (0)
real logKphi[includephi && includeKeffectonphi ? NrOfKretsar : 0];
real<lower=0> taoLalpha[includealpha && includeKeffectonalpha ? NrOfLans : 0];
real<lower=0> taoLphi[includephi && includeKeffectonphi ? NrOfLans : 0];
real<lower=0> sigmataoLalpha[includealpha && includeKeffectonalpha ? 1 : 0];
real<lower=0> sigmataoLphi[includephi && includeKeffectonphi ? 1 : 0];
real<lower=0> taoSalpha[includealpha && includeKeffectonalpha ? 1 : 0];
real<lower=0> taoSphi[includephi && includeKeffectonphi ? 1 : 0];
}
transformed parameters {
real<lower=0> alpha[includealpha ? NrOfKretsar : 0];
real phi[NrOfKretsar];
real<lower=0> m[NrOfKretsar];
real<lower=0> tm[NrOfLans];
real<lower=0> talpha[includealpha && includeKeffectonalpha ? NrOfLans : 0];
real<lower=0> tphi[includephi && includeKeffectonphi ? NrOfLans : 0];
for (ik in 1:NrOfKretsar){
m[ik] = exp( logSm + logLm[KretsListLanSeqID[ik]] + logKm[ik] );
if (includealpha ){
if (includeKeffectonalpha) {
alpha[ik] = exp(logSalpha[1] + logLalpha[KretsListLanSeqID[ik]] + logKalpha[ik] );
} else {
if (includeLeffectonalpha ) {
alpha[ik] = exp(logSalpha[1] + logLalpha[KretsListLanSeqID[ik]]);
} else {
alpha[ik] = exp(logSalpha[1]);
}
}
}
if (includephi){
if (includeKeffectonphi) {
phi[ik] = logSphi[1] + logLphi[KretsListLanSeqID[ik]] + logKphi[ik];
} else {
if (includeLeffectonphi){
phi[ik] = logSphi[1] + logLphi[KretsListLanSeqID[ik]];
} else {
phi[ik] = logSphi[1];
}
}
} else {
phi[ik]=0;
}
}
for (il in 1:NrOfLans){
if (includeLdifferencevaraibilityonm){
tm[il] = taoSm * taoLm[il]; //tm is the precission of between circuit variability of logK for each county
}else{
tm[il] = taoSm;
}
if (includealpha && includeKeffectonalpha){
talpha[il] = taoSalpha[1] * taoLalpha[il]; //
}
if (includephi && includeKeffectonphi){
tphi[il] = taoSphi[1] * taoLphi[il]; //
}
}
}
model{
for (r in 1:Reports){
// Number of killed animals in report r is negative binomial distributed,
// but defined from the mean (m[KretsNrSeqID[r]] * A[r]) and shape (c[Lan[r]]) of the gamma distribution of the equivalent gamma-poisson mixture
if (includealpha ){
K[r] ~ neg_binomial_2(m[KretsNrSeqID[r]] * A[r] * exp(logrelA[r] * phi[KretsNrSeqID[r]]) ,alpha[KretsNrSeqID[r]]);
} else {
K[r] ~ poisson(m[KretsNrSeqID[r]] * A[r] * exp(logrelA[r] * phi[KretsNrSeqID[r]]));
}
}
for(ik in 1:NrOfKretsar){
logKm[ik] ~ normal( 0 , tm[KretsListLanSeqID[ik]]^-.5 ); // logKm[Kretsar[ik]] is the circuit effect on logm
if (includealpha && includeKeffectonalpha){
logKalpha[ik] ~ normal( 0 , talpha[KretsListLanSeqID[ik]]^-.5 ); // logKalpha[Kretsar[ik]] is the circuit effect on logalpha
}
if (includephi && includeKeffectonphi){
logKphi[ik] ~ normal( 0 , tphi[KretsListLanSeqID[ik]]^-.5 ); // logKphi[Kretsar[ik]] is the circuit effect on logphi
}
}
for (il in 1:NrOfLans){
logLm[il] ~ normal( 0 , Tm^-.5 );// Tm is the std of between county variability for logm
if (includealpha && includeLeffectonalpha){
logLalpha[il] ~ normal( 0 , Talpha[1]^-.5 );// Talpha is the std of between county variability for logalpha
}
if (includephi && includeLeffectonphi){
logLphi[il] ~ normal( 0 , Tphi[1]^-.5 );// Tphi is the std of between county variability for logphi
}
if(includeLdifferencevaraibilityonm){
taoLm[il] ~ lognormal( 0 , sigmataoLm[1]^-.5 ); // taoLm[il] is the std of variability of logm between circuits in county il
}
if (includealpha && includeKeffectonalpha){
taoLalpha[il] ~ lognormal( 0 , sigmataoLalpha[1]^-.5 ); //taoLalpha [il] is the std of variability of logalpha between circuits in county il and is log-normally distributed with std sigmataoLalpha
}
if (includephi && includeKeffectonphi){
taoLphi[il] ~ lognormal( 0 , sigmataoLphi[1]^-.5 ); //taoLphi [il] is the precision of variability of logphi between circuits in county il and is log-normally distributed with std sigmataoLphi
}
}
logSm ~ normal(0,100 );
if (includealpha){
logSalpha[1] ~ normal(0,2.3496 );//based on 95% coef var between 0.1 and 10
if (includeLeffectonalpha) {
Talpha[1] ~ gamma(1.9178,0.2171); // based on 95% between 1 and 25
}
if (includeKeffectonalpha){
sigmataoLalpha[1] ~ lognormal(0,1 );
taoSalpha[1] ~ lognormal(0,2 );
}
}
if (includephi){
logSphi[1] ~ normal(0,0.2984503);// based on 95% between -0.585 and 0.585, which corresponds to maximum 300% increas in hunting rate with twice the area
if (includeLeffectonphi){
Tphi[1] ~ gamma(1.9178,0.2171); // based on 95% between 1 and 25
}
if (includeKeffectonphi){
sigmataoLphi[1] ~ lognormal(0,1 );
taoSphi[1] ~ lognormal(0,2 );
}
}
Tm ~ gamma(1,.1 );
if(includeLdifferencevaraibilityonm){
sigmataoLm[1] ~ gamma(1,.1 );
}
taoSm ~ gamma(1,.1 );
}
generated quantities {
vector[Reports] log_lik;
for (r in 1:Reports){
if (includealpha){
log_lik[r] = neg_binomial_2_lpmf(K[r] | m[KretsNrSeqID[r]] * A[r] * exp(logrelA[r] * phi[KretsNrSeqID[r]]), alpha[KretsNrSeqID[r]]);
} else {
log_lik[r] = poisson_lpmf(K[r] | m[KretsNrSeqID[r]] * A[r] * exp(logrelA[r] * phi[KretsNrSeqID[r]]));
}
}
}