Error in Stan code when modelling a weighted logistic regression model

I am reading this paper, which have discussed about how to use survey weights in Bayesian modeling using pseudo‐maximum likelihood method. I just ran the following stan model they have provided in the paper.


real wt_bin_lpmf(int[] y, vector mu, vector weights, int n){
real check_term;
check_term = 0.0;
for( i in 1:n )
{
check_term = check_term +
weights[i] * bernoulli_logit_lpmf(y[i] | mu[i]);
}
return check_term;
}}
data{
int<lower=1> n; // number of observations
int<lower=1> k; // number of linear predictors
int<lower=0, upper = 1> y[n]; // Response variable
vector<lower=0>[n] weights; // observation-indexed weights
matrix[n, k] X; // coefficient matrix
}
parameters{
vector[k] beta; /* regression coefficients from linear predictor */
}
transformed parameters{
vector[n] mu;
mu = X * beta; /* linear predictor */
} /* end transformed parameters block */
model{
/*improper prior on beta in (-inf,inf)*/
/* directly update the log-probability for sampling */
target += wt_bin_lpmf(y | mu, weights, n);
}

In this code, they have created a R function to calculate the pseudo‐maximum likelihood. When I am trying to run this stan file, I am getting this error.


PARSER FAILED TO PARSE INPUT COMPLETELY
STOPPED AT LINE 1:

It has stopped from the first line. What may be the reason for this?

Thank you

You’re missing the functions{ statement above the wt_bin_lpmf declaration, it should look like this:

functions{
real wt_bin_lpmf(int[] y, vector mu, vector weights, int n){
real check_term;
check_term = 0.0;
for( i in 1:n )
{
check_term = check_term +
weights[i] * bernoulli_logit_lpmf(y[i] | mu[i]);
}
return check_term;
}}
1 Like

@andrjohns
I appreciate your reply. I have a related question about this post. The model based on above stan code is very slow and had many divergent transitions and higher R-hat values. Do you have any suggestions to improve it ? I used normal priors for beta.

Thank you

Nothing that will make a large difference, the main thing to check is that your priors and model weights are reasonable.

Also, it might be helpful to know that the model above can be specified in brms very simply:

brm(y | weights(w) ~ x, family=bernoulli("logit"), data=data)

Where w is the vector of weights.

1 Like