High pareto K diagnostics when adding measurement error to a model

Modeling
1

Hi there,

This is an update on a question I posted the other day here.

My data includes the mean and SE of antibody titres measured in several groups of individuals, over time. I don’t have any individual titre measurements. Therefore, I want to fit a model to the mean antibody titre, whilst accounting for the known SE. To do this, I followed the format in the stan guide for a Bayesian Measurement Error Model. The titres (mu_obs) are assumed to follow a model of biphasic decay, and I estimate the true titre values (mu_true).

When I fit the model without the measurement error, everything converges, the pareto k diagnostics are fine etc. When I add in the measurement error, the model still converges but now all the pareto k diagnostics are > 1. This makes sense to me, as I now have more parameters than data, as I estimate a mu_true for each mu_obs, in addition to the parameters estimated as part of the biphasic model.

Is there any way to reduce the number of parameters, whilst still accounting for the uncertainty in the antibody titre measurement?

I have included the model below:

data{
   int N;               // Number of time points data collected 
   int T ;              // Number of time points to solve the biphasic model at
   int B ;              // Number of baseline serostatuses
   int K ;              // Number of serotypes 
   real time[T];        // time 
   real mu_meas[B,K,N]; // measured average titres   
   real tau[B,K,N];     // SE of titres 
}


parameters {
  real<upper = 0> pi_1;     // decay param 1
  real pi_2;                // decay param 2
  real<lower = 0> ts   ;    // time switch decay rate
  real<lower = 0> sigma ;   // sd
  real mu_true[B,K,N] ;     // unknown true titres  
  
}

transformed parameters {
  real neut[B,K,T] ;
  real n_N[B,K,N] ;
  
  // biphasic model 
  for(b in 1:B){
   for(k in 1:K){
    for(t in 1:T)
    neut[b,k,t] = mu_true[b,k,1] * (exp(pi_1 * time[t] + pi_2 * ts) + exp(pi_2 * time[t] + pi_1 * ts)) / (exp(pi_1 * ts) + exp(pi_2 * ts)) ;

   n_N[b,k,1] =    neut[b,k,1] ;
   n_N[b,k,2] =    neut[b,k,6] ;
   n_N[b,k,3] =    neut[b,k,12] ;
   n_N[b,k,4] =    neut[b,k,24] ;
   n_N[b,k,5] =    neut[b,k,36] ;
   n_N[b,k,6] =    neut[b,k,48] ;
   }
    }
}


model {

// likelihood  
  for(b in 1:B)
   for(k in 1:K)
   target += normal_lpdf(mu_meas[b,k,]| mu_true[b,k,], tau[b,k,]);

// priors
  for(b in 1:B)
   for(k in 1:K)
   mu_true[b,k,] ~ normal(n_N[b,k], sigma) ; 
 
  pi_1 ~ normal(0,0.1) ;
  pi_2 ~ normal(0,0.1) ;
  ts ~ normal(1,2) ;
  sigma ~ cauchy(0,5);
}

generated quantities{
 vector[B*K*N] log_lik ;
 
 for(t in 1:N){
   log_lik[t]     = normal_lpdf(mu_meas[1,1,t]| mu_true[1,1,t], tau[1,1,]) ;
   log_lik[t+N]   = normal_lpdf(mu_meas[1,2,t]| mu_true[1,2,t], tau[1,2,]) ;
   log_lik[t+2*N] = normal_lpdf(mu_meas[1,3,t]| mu_true[1,3,t], tau[1,3,]) ;
   log_lik[t+3*N] = normal_lpdf(mu_meas[1,4,t]| mu_true[1,4,t], tau[1,4,]) ;
   log_lik[t+4*N] = normal_lpdf(mu_meas[2,1,t]| mu_true[2,1,t], tau[2,1,]) ;
   log_lik[t+5*N] = normal_lpdf(mu_meas[2,2,t]| mu_true[2,2,t], tau[2,2,]) ;
   log_lik[t+6*N] = normal_lpdf(mu_meas[2,3,t]| mu_true[2,3,t], tau[2,3,]) ;
   log_lik[t+7*N] = normal_lpdf(mu_meas[2,4,t]| mu_true[2,4,t], tau[2,4,]) ;
 }

}

Thank you in advance!

2

You could have also continued in that thread to keep the information in one thread.

The previous model did not include this one term. This makes the model over-parameterised. Each observation has its own parameter, and unless sigma ise very small, then leaving out one observation changes the posterior of the corresponding parameter a lot. See a corresponding example in Roaches cross-validation demo. That demo shows also how it can be possible to marginalize over that one parameter to get importance sampling to work.

3

Hi Aki,

Thank you for your quick response! Sorry, next time I will continue on the same thread.

I am working on updating the model and integrating loo, thanks for sharing the example.