How to handle high Pareto k when observations are meaningful?

Modeling
, ,
1

Hello,

I am attempting to fit discrete time event models to some frog tadpole metamorphosis data.

I have been able to successfully fit the model, but LOO for model comparisons demonstrate high pareto K values for roughly 0.3% of observations (10 out of n = 2879).

I have tried moment-matching, using non-random intercepts, longer chain iterations and predictor transformations but that only really removes on k > 0.7 observation if at all. I know the standard is to remove outliers, but the outliers correspond to metamorphosis observations (specifically high abundance of competitors), which I think is what is driving the event to occur so removing the observations seems ill-advised.

I know there are some alternatives, such as Mixture Importance sampling or using exact leave-one-out, a la reloo in brms, but since I am computing directly through rstan I am lost on how to achieve this.

I know that with high k values, that the output of resulting LOO comparisons are unreliable, but I don’t have the know-how to determine if the 10 observations are “that bad” and can be “ignored” or implement mixis or reloo methods through rstan.

I have attached my model and resulting outputs. Any help is appreciated!

Stan Model Code

data {
  int<lower=1> N;                   // Number of data rows
  int<lower=1> T;                   // Number of time steps
  int<lower=1> K;                   // Number of predictors
  int<lower = 1> Ta;                // Number of tanks
  int<lower = 1, upper = Ta> Tank[N];// Tank ID
    int<lower = 1> En;                // Number of enclosures
  int<lower = 1, upper = En> Enclosure[N];// Enclosure ID
  int<lower=1> Time[N];             // Time interval (e.g., 1, 2, 3,...)
  int<lower=0> M[N];                // Metamorphosis (0 or 1)
  matrix[N,K] X;                    // Predictors
}

parameters {
  vector[T] alpha;                 // Time-specific intercepts
  vector[K] beta;                  // Coefficients for predictors
  vector[Ta] gamma;                // random effects for tank
  vector[En] tau;                // random effects for enclosure
  real <lower=0> sigma_gamma;      // standard deviation of tank random effects
  real <lower=0> sigma_alpha;      // standard deviation of temporal random effects
  real <lower=0> sigma_tau;      // standard deviation of enclosure random effects
}

transformed parameters {
  vector[N] h;                     // Hazard probabilities
  for (i in 1:N) {
    int t = Time[i];
    int ta = Tank[i];
    int en = Enclosure[i];
    h[i] = 1 - exp(-exp(alpha[t] + gamma[ta] + tau[en] + X[i]* beta));
  }
}


model {
  // Priors (modify these as appropriate)
  beta ~ normal(0, 1);
  
  // Random effects
  gamma ~ normal(0, sigma_gamma);
  sigma_gamma ~ exponential(1);
  
  alpha ~ normal(0, sigma_alpha);
  sigma_alpha ~ exponential(1);
  
  tau ~ normal(0, sigma_tau);
  sigma_alpha ~ exponential(1);
  // Likelihood
  M ~ bernoulli(h);
}

generated quantities {
  vector[N] log_lik;
  int M_rep[N];  // simulated metamorphosis outcomes
  for (i in 1:N) {
    log_lik[i] = bernoulli_lpmf(M[i] | h[i]);
    M_rep[i] = bernoulli_rng(h[i]);  // simulate new dataset from the model
  }
}

An example of LOO output is (which is identical to moment matched LOO):

Computed from 8000 by 2879 log-likelihood matrix.

         Estimate   SE
elpd_loo   -222.6 21.7
p_loo        43.4  6.5
looic       445.2 43.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     2869  99.7%   365     
   (0.7, 1]   (bad)        10   0.3%   <NA>    
   (1, Inf)   (very bad)    0   0.0%   <NA>    
See help('pareto-k-diagnostic') for details.```

with the LOO compare table being:

   elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo  looic se_looic
g1      0.00    0.00  -222.60       21.72 43.35     6.50 445.20    43.44
g2     -2.33    1.87  -224.93       22.04 45.20     6.84 449.86    44.07
2

Generically, what reloo does is simply refit the model J times (once for each problematic \hat{k} value) and obtain the exact held out log_lik draws. So you’ll probably need to modify your Stan model to take in training (i.e., all the data points except the held out one) data and test (the held out point) data, then compute the log_lik draws in generated quantities just for this one point. Then you can loop over the 10 problematic observations and store them as a S \times J matrix (or array).

You should be able to substitute these log_lik columns for the PSIS
draws within the regular loo functionality but I am less familiar with the internal machinery and can’t comment exactly on these part.

3

High Pareto-k observations are not necessarily outliers, see the discussion in CV-FAQ. They can be due to model being flexible or due to random chance with finite sample size. As you had already tried running more iterations, it is likely that it’s not just due to random chance. I have been recently been experimenting more with mixis and based on what you tell, it is likely that it could work well in your case. loo package has Mixture importance sampling vignette which is using rstan.

4

Thank you for the direction! I followed the vignette and made the following code:

data {
  int<lower=1> N;                   // Number of data rows
  int<lower=1> T;                   // Number of time steps
  int<lower=1> K;                   // Number of predictors (includes intercept)
  int<lower = 1> Ta;                // Number of tanks
  int<lower = 1, upper = Ta> Tank[N];// Tank ID
    int<lower = 1> En;                // Number of enclosures
  int<lower = 1, upper = En> Enclosure[N];// Enclosure ID
  int<lower=1> Time[N];             // Time interval (e.g., 1, 2, 3,...)
  int<lower=0> M[N];                // Metamorphosis (0 or 1)
  matrix[N,K] X;                    // Predictors
  int <lower=0> mixis;
}

parameters {
  vector[T] alpha;                 // Time-specific intercepts
  vector[K] beta;                  // Coefficients for predictors
  vector[Ta] gamma;                // random effects for tank
  vector[En] tau;                // random effects for enclosure
  real <lower=0> sigma_gamma;      // standard deviation of tank random effects
  real <lower=0> sigma_alpha;      // standard deviation of temporal random effects
  real <lower=0> sigma_tau;      // standard deviation of enclosure random effects
}

transformed parameters {
  vector[N] h;                     // Hazard probabilities
  for (i in 1:N) {
    int t = Time[i];
    int ta = Tank[i];
    int en = Enclosure[i];
    h[i] = 1 - exp(-exp(alpha[t] + gamma[ta] + tau[en] + X[i]* beta));
  }
}


model {
  // Priors (modify these as appropriate)
  beta ~ normal(0, 1);
  
  // Random effects
  gamma ~ normal(0, sigma_gamma);
  sigma_gamma ~ exponential(1);
  
  alpha ~ normal(0, sigma_alpha);
  sigma_alpha ~ exponential(1);
  
  tau ~ normal(0, sigma_tau);
  sigma_tau ~ exponential(1);
  // Likelihood
  M ~ bernoulli(h);
  
  vector[N] log_lik;
  for(n in 1:N){
    log_lik[n] = bernoulli_lpmf(M[n] | h[n]);
  }
  target += sum(log_lik);
  if (mixis){
    target += log_sum_exp(-log_lik);
  }
}

generated quantities {
  vector[N] log_lik;
  int M_rep[N];  // simulated metamorphosis outcomes
  for (i in 1:N) {
    log_lik[i] = bernoulli_lpmf(M[i] | h[i]);
    M_rep[i] = bernoulli_rng(h[i]);  // simulate new dataset from the model
  }
}


llm<- extract(testfit, pars = "log_lik")$log_lik
library(matrixStats)
lcm<- rowLogSumExps(-llm)
lw<- -llm-lcm
elpdmixix<- logSumExp(-lcm) - rowLogSumExps(t(lw))   
print(paste("ELPD (MixIS)=",round(sum(elpdmixix),2)))

ELPD (MixIS)= -190.58

If I am understanding this correctly, I set all mixis values in the data equal to one? Then from there use the matrix stats functions to generate an ELPD? I am interpreting and thus implementing this corrrectly?

Thank you for your help thus far!

5

You should drop M ~ bernoulli(h); as otherwise you have the likelihood contribution twice (also when you write M ~ bernoulli(h); it would be better to call that term “data model” as you are saying that M has Bernoulli distribution)

As you had only 0.7 < k < 1, the PSIS estimate and mixis estimate should be quite close to each other

6

Thank you for such an expedient reply! Per your suggestions the model block would simply look like this?

model {
  // Priors (modify these as appropriate)
  beta ~ normal(0, 1);
  
  // Random effects
  gamma ~ normal(0, sigma_gamma);
  sigma_gamma ~ exponential(1);
  
  alpha ~ normal(0, sigma_alpha);
  sigma_alpha ~ exponential(1);
  
  tau ~ normal(0, sigma_tau);
  sigma_tau ~ exponential(1);
  // Likelihood
  
  vector[N] log_lik;
  for(n in 1:N){
    log_lik[n] = bernoulli_lpmf(M[n] | h[n]);
  }
  target += sum(log_lik);
  if (mixis){
    target += log_sum_exp(-log_lik);
  }

For clarification, I am correct in saying that the mixis variable is just a dummy variable of 1 as it is passed to the stan data?

7

Yes (or you can create s separate Stan file which doesn’t use that variable and that if statement)

8

Fantastic! Thank you so much for your help!

9

Please post here when you have some results or paper or something you can share. Always interesting to see what people do