# User defined function to produce samples in mixed effect models and produce posterior predicted probability plot

Hello,
I am using user defined function for mixed model. I would like to generate posterior predicted probability, but it seems like my function is not working. Instead of writing the function in the function block I used increment log probability. I also wrote log probability on the function block but couldn’t be able to generate the samples in generated blocks. …rng function did not work. Anyone please helps me to generate the samples by using my function. The code I have created is as follows. My code is working perfectly other than generating the samples. I also attached my data so that it will be easier for your to run the code. Please help me to produce posterior predicted plot of response varaible y and y_rep.

``````data {
int<lower=0> n; // number of observations
int<lower=0> M; // number of subjects
real<lower=0,upper=1> y[n]; // response variable
int K; // number of covariate
matrix [n,K] X; // covariate
int id[n]; // id variable
vector [K] b0;
matrix[K,K] B0;
}

parameters {
vector [K] beta;
real <lower=0> sigma2;
vector [M] bi_r;
}
transformed parameters {
real <lower=0> sigma;     //sigma in original bugs model
real <lower=0> mu[n];
vector [M] bi;
vector[n] eta;
sigma=sqrt(sigma2);
eta=X*beta;
bi=sigma*bi_r;
for(i in 1:n){
mu[i] = inv_logit(eta[i] + bi[id[i]]);
}
}

model {
for (j in 1:M){
bi[j] ~ normal(0, sigma);
bi_r[j] ~ normal(0, 1);
}
beta~multi_normal(b0, B0);
sigma2 ~ uniform(0, 10);
for (i in 1:n){
increment_log_prob(2*log1m(mu[i])
- log(mu[i])
- 3*log1m(y[i])
- (y[i]*(1 - mu[i]))/(mu[i]*(1 - y[i])));
}
}
generated quantities{
vector[n] y_pred;
real dev;
vector[n] log_lik;
dev = 0;
for (i in 1:n) {
log_lik[i] = (2*log1m(mu[i])
- log(mu[i])
- 3*log1m(y[i])
- ((y[i]*(1 - mu[i]))/(mu[i]*(1 - y[i]))));
y_pred[i] = mu[i];
}
dev = dev+(-2)*sum(log_lik);
}
``````

mydata.csv (9.7 KB)

Nirajan Bam

I reformatted your code - code blocks need to be enclosed in triple back-quotes.

It looks like you’re trying to translate some BUGS model directly into Stan. I suggest you start by reading the appendix to the User’s Guide on this - 36 Transitioning from BUGS | Stan User’s Guide. Also, see chapters on mixture models - 5 Finite Mixtures | Stan User’s Guide - and section on vectorization - 25.8 Vectorization | Stan User’s Guide

Hello Mitzi,
Thank you for your response. I even don’t know Bugs. I did not translate anything from Bugs. My code is working fine from my end. I am able to recover my parameters. Of course, there are some cleaning issues and vectoring the code. My only question is to generate response varaible y by using the log likelihood function:
(2log1m(mu[i])
- log(mu[i])
- 3
log1m(y[i])
- (y[i](1 - mu[i]))/(mu[i](1 - y[i])))
In my model mu is only on parameter. The model is not a mixture model.
You said you reformatted the code. Would you please send me your reformatted codes. Thanking you,

Nirajan Bam

I edited your post and added the backticks so that the code displays as code.
I also ran your code through an editor to get proper indenting, eliminated use of now-deprecated assignment operator `->` - changed to `=`, and added whitespace for readability. again, I would refer you to the Stan User’s Guide, which has a section on coding standards.

there’s a comment in your code `sigma in original bugs model`, which is why I thought that you’re translating a BUGS model. perhaps you meant something different.

I realize that you are very busy and you just want to make your model work. however, I see no obvious fix to the code you’ve shared - the way forward it to start with a very simple model - perhaps not a mixture model, just a simple regression model. or maybe consider an off-the-shelf solution using RStanArm or BRMS. then add the posterior predictive check - 27 Posterior Predictive Sampling | Stan User’s Guide - if you can get one RNG function in the generated quantities block working, you’ll be in better shape to deal with more complicated likelihoods.

I was going to try running the code, but the CSV file doesn’t seem to include all of the variables defined in the code. Am I missing something?

Hello Khol,
You need do create the data in the list form like this:
“X<-model.matrix(~x1i+x2i+x3i+x4i, data=mydata)
K= ncol(X)
mydata<-list(n=length(mydata\$y),K= ncol(X), M=max(mydata\$id), y=mydata\$y, X=X, id=mydata\$id,b0=rep(0,K),B0=diag(100,K))”
Thanks,

Thank you Mitzi !!

Do you get the solution to your problem? I am also having a similar problem with the new probability distribution model, whose quantile function is not in closed form so could not generate rng. Do you have any ideas? Please share them with me.

Hello Buddy,
Glad to hear from you. I started simulation without PPC. My problem has not been solved yet. Surely, I will share with you if solved. STAN manual is not much clear about PRNG functions. Thanks!!!