Looic SE yielding to NA value

Modeling
1

When running the following Stan and R codes, I am getting an estimate for loo but the estimate of its SE yields an NA value. However, model diagnoses seem to be ok and fit the observed data. How could I interpret this?

R code:
datalist2 ← list(N = nobs, nsubjects = indivs,
subj_ids = new_ids,
y = antibodies, ts = ts, t0 = 0.0, ysex = sex, yAV = studygroup, yage = age)
rAB0 = rep(0, length(unique(ids)))
ruAB = rep(0, length(unique(ids)))
r = cbind(rAB0, ruAB)
str(r)
beta3 = c(-0.1,0.2,0.0001)
init2 = function(){
list(logpABS = 0, logpABL = 0, loguAB = log(0.05), loguSASC = log(0.01),
logAB0 = 7, logsigma = 0, logsigmaAB0 = -3, logsigma_uAB = -2,
rho_rhobit = 0.5, r = r, beta = beta3)
}

Model6cov ← stan_model(“C:/Users/IGarciaFogeda/Documents/Stan/Rstan/HIERARCHICAL/Longitudinal/longitudinal_model6.stan”)
fit_model6cov ← sampling(Model6cov,
data = datalist2, init = init2,
iter = 1000,
chains = 2, warmup = 500,
seed = 0, control = list (stepsize = 0.1))

functions {
  real[] model6(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[3];
    dydt[1] = theta[1]*y[2] + theta[2]*y[3] - theta[3]*y[1];
    dydt[2] = -theta[4]*y[2];
    dydt[3] = 0;
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  real ts[N];
  real t0;
  real ysex[N];
  real yAV[N];
  real yage[N];
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logpABS;
  real logpABL;
  real loguAB;
  real loguSASC;
  
  real logAB0;
  real logsigma;
  real logsigmaAB0;
  real logsigma_uAB;
  real rho_rhobit;
  
  vector[2] r[nsubjects];
  real beta[3];
}
model {
  real new_mu[N,3];
  real mu[1,3];
  real eval_time[1];
  real theta[4];
  vector[nsubjects] rAB0;
  vector[nsubjects] ruAB;
  
  real inits[3]; 
  
  real pABS = exp(logpABS);
  real pABL = exp(logpABL);
  real uAB = exp(loguAB);
  real uSASC = exp(loguSASC);
  
  real AB0 = exp(logAB0);

  real sigma = exp(logsigma);
  real sigmaAB0 = exp(logsigmaAB0);
  real sigma_uAB = exp(logsigma_uAB);
  
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  real rho = (exp(rho_rhobit) - 1)/(exp(rho_rhobit) + 1);
  
  theta[1] = pABS;
  theta[2] = pABL;
  theta[4] = uSASC;
  
  //prior distributions
  logpABS ~ normal(0, 10);
  logpABL ~ normal(0, 10);
  loguSASC ~ normal(0, 10);
  loguAB ~ normal(0, 10);
  
  logAB0 ~ normal(7, 2);

  logsigmaAB0 ~ normal(0, 1);
  logsigma_uAB ~ normal(0, 1);
  rho_rhobit ~ uniform(-10,10);
  
  mean_vec[1] = -pow(sigmaAB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_uAB, 2.0)/2;
  
  Sigma[1,1] = pow(sigmaAB0, 2.0);
  Sigma[2,2] = pow(sigma_uAB, 2.0);
  Sigma[1,2] = rho*sigmaAB0*sigma_uAB;
  Sigma[2,1] = rho*sigmaAB0*sigma_uAB;
  
  for (subj_id in 1:nsubjects) {
        r[subj_id] ~ multi_normal(mean_vec, Sigma);
        rAB0[subj_id] = r[subj_id,1];
        ruAB[subj_id] = r[subj_id,2];
  }
  
  for (obs_id in 1:N){
    inits[1] = AB0*exp(rAB0[subj_ids[obs_id]]);
    inits[2] = inits[1]/pABS;
    inits[3] = inits[1]/pABL;
    
    theta[3] = uAB*exp(ruAB[subj_ids[obs_id]]  + beta[1]*ysex[obs_id]+ beta[2]*yAV[obs_id] + beta[3]*yage[obs_id]); 
     
    eval_time[1] = ts[obs_id];
    if (eval_time[1] == 0){
      new_mu[obs_id,1] = log(inits[1]) - pow(sigma, 2.0)/2;
      new_mu[obs_id,2] = inits[2];
      new_mu[obs_id,3] = inits[3];}
    else{
      mu = integrate_ode_rk45(model6, inits, t0, eval_time, theta, x_r, x_i);
      new_mu[obs_id,1] = log(mu[1,1]) - pow(sigma, 2.0)/2;
      new_mu[obs_id,2] = mu[1,2];
      new_mu[obs_id,3] = mu[1,3];
    }  
    //likelihood function for AB levels
    y[obs_id] ~ lognormal(new_mu[obs_id,1], sigma); 
  }
}
generated quantities {
  real z_pred[N];
  real sigma2 = exp(logsigma);
  for (t in 1:N){
   z_pred[t] = lognormal_rng(log(y[t] + 1e-5), sigma2); 
  }
}

paretodiagnosis
paretodiagnosis909×483 19.2 KB

2

How are you computing loo? Your Stan code doesn’t include computation for log_lik which is the usual way. The figure showing the loo output is missing some lines from top, but as there is only one k-value, it indicates you are providing for loo function only on observation, and standard error cannot be computed with one observation producing NA.

3

Hi,

Thanks for the message. I have gotten the loo in the following way through calling ‘lp__’, which I thought it would yield the log_lik?

log_lik_1 ← extract_log_lik(fit_model6cov, parameter_name = “lp__”, merge_chains = TRUE)
loo_1 ← loo(log_lik_1)
print(loo_1)
waic(log_lik_1)

4

This includes the log prior density and all the pointwise log likelihood components summed together. You need to compute pointwise log likelihood terms separately in the generated quantities. See Writing Stan programs for use with the loo package • loo

1 Like
5

Thanks!