Dear all,
I have a question regarding the model fit diagnostics for a NLMM with censoring in the response.
I fitted a NLMM model in stan where there is censoring in the response.
I used an equivalent version of Tobit regression to fit this model.
Now, when I want to assess the model fit by using, for example, the plot between predicted values v.s observed values, I am facing now the difficulty to deal with censored data points.
The plot now I have is something which looks like:
Clearly, I am searching for an appropriate way to be able to assess the fit of censored data points (which based on the plot: the vertical line at around 1.6 since the data are plotted on the log scale, and the response has lower censoring at 5 which means those values less than 5 are not identified to the accurate value of it).
In this link, http://monolix.lixoft.com/data-and-models/censoreddata/, they proposed to sample the censored data points from this distribution:
p(y_{censored} | y_{observed} ,ψ,θ ) where ψ , θ are the estimated population and individual parameters. " This is done by adding a residual error on top of the prediction, using a truncated normal distribution to make sure that the simulated y remains within the censored interval."
Can anyone help to give me a hint how to sample this using Stan, given the fact that I have a NLMM model where the response values are predicted based on the time point they were measured?
Thank you.
The code to fit the model is provided below:
data {
int<lower=1> N; // Total number of data points
real<lower=0> time[N]; // Time points
real<lower=0> y[N]; // Antibody levels on log10 scale
int<lower=1> J; // Number of subjects
int<lower=1, upper=J> subj[N]; // Subject id
int<lower=0, upper=1> cens[N]; // Censoring indication
int<lower=0, upper=1> C2[N]; // Child indication
real ABMum[N]; // AB of pregnant women at birth
}
parameters {
real<lower=0> beta; real<lower=0> gamma; real<lower=4 ,upper=12> h;
real<lower=0> alpha; real<lower=0,upper=100> A0;
real<lower=0> sigma;
vector<lower=0>[J] b_raw; //To apply non-centered parameterization
vector<lower=0>[J] bbeta; vector<lower=0>[J] bgamma; vector<lower=0>[J] balpha;
real<lower=0,upper=100> sigma_A0; real<lower=0,upper=100> sigma_beta;
real<lower=0,upper=100> sigma_gamma; real<lower=0,upper=100> sigma_alpha;
real v1A0 ; real v2A0 ;
real v1beta;
real v1gamma;
real v1alpha;
real v1h;
}
transformed parameters{
real betapart[N]; real gammapart[N]; real alphapart[N]; real<lower=0> A0part[N]; real<lower=4> hpart[N];
real mu[N]; real diff[N];
for (i in 1:N) {
A0part[i] = ( A0 + v1A0*C2[i] + v2A0*ABMum[i] ) * b_raw[subj[i]] ;
betapart[i] = ( beta + v1beta*C2[i] ) * bbeta[subj[i]] ;
gammapart[i] = ( gamma + v1gamma*C2[i] ) * bgamma[subj[i]] ;
alphapart[i] = ( alpha + v1alpha*C2[i] ) * balpha[subj[i]] ;
hpart[i] = h + v1h*C2[i] ;
if ( (time[i] <= 2 ) && (cens[i] ==0) ) { //tvacc[i]
mu[i] = log(A0part[i]) - betapart[i]*time[i] ; }
else if ( time[i] > 2 && time[i] <= hpart[i] && ( cens[i] ==0 ) ) { //tvacc[i]
mu[i] = log(A0part[i]) - betapart[i]*2 + gammapart[i]*(time[i] - 2) ; } //tvacc[i]
else if ( (time[i] > hpart[i] ) && (cens[i] ==0) ) {
mu[i] = log(A0part[i]) - betapart[i]*2 + gammapart[i]*(hpart[i] - 2) - alphapart[i]*(time[i] - hpart[i]); } //vacc[i]
else if ((time[i] <= 2) && (cens[i] ==1)) { //tvacc[i]
mu[i] = log(A0part[i]) - betapart[i]*time[i] ; }
else if ( time[i] > 2 && time[i] <= hpart[i] && (cens[i] ==1) ) { //tvacc[i]
mu[i] = log(A0part[i]) - betapart[i]*2 + gammapart[i]*(time[i] - 2) ; } //tvacc[i]
else if ((time[i] > hpart[i] ) && ( cens[i] == 1 ) ) {
mu[i] = log(A0part[i]) - betapart[i]*2 + gammapart[i]*(hpart[i] - 2 ) - alphapart[i]*(time[i] - hpart[i]);} //2
diff[i] = (y[i] - mu[i])/sigma;
}
}
model {
// Priors
bbeta ~ gamma( sigma_beta , sigma_beta); balpha ~ gamma( sigma_alpha , sigma_alpha ) ;
b_raw ~ gamma( sigma_A0 , sigma_A0 ); bgamma ~ gamma( sigma_gamma , sigma_gamma );
sigma_A0 ~ gamma(0.01, 0.01); sigma_alpha ~ gamma(0.01, 0.01);
sigma_beta ~ gamma(0.01 , 0.01); sigma_gamma ~ gamma(0.01, 0.01);
// Likelihood
for (i in 1:N){
if ( (time[i] <= 2 ) && (cens[i] ==0) ) { //tvacc[i]
target += normal_lpdf(y[i] | mu[i], sigma); }
else if ( time[i] > 2 && time[i] <= hpart[i] && ( cens[i] ==0 ) ) { //tvacc[i]
target += normal_lpdf( y[i] | mu[i], sigma); }
else if ( (time[i] > hpart[i] ) && (cens[i] ==0) ) {
target += normal_lpdf( y[i] | mu[i], sigma); }
else if ((time[i] <= 2) && (cens[i] ==1)) { //tvacc[i]
target += normal_lcdf(y[i] | mu[i], sigma);}
else if ( time[i] > 2 && time[i] <= hpart[i] && (cens[i] ==1) ) { //tvacc[i]
target += normal_lcdf( y[i] | mu[i] , sigma ) ; }
else if ((time[i] > hpart[i] ) && ( cens[i] == 1 ) ) {
target += normal_lcdf(y[i] | mu[i] , sigma ) ; }
}
}