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!