# Fitting a multilevel GPD to negated data

I want to fit a multi-level generalised pareto distribution (GPD) to negated data. I have used stan code from here. While running the model, some chains are not initializing at all but other few are running fine and converging properly. I get the following error for those chains which are not initializing

More specifically, in my study, I have 9 subgroups, for which TG depends upon SG. As single threshold cannot classify extreme in all groups, I am using 9 different thresholds obtained from threshold stability plots. As I am interested in minimum (Prob (TG<=0)), I have used negated value of TG to fit a GPD. Therefore, I am trying to compute the probability (neg_TG>=0).
I’m wondering whether it’s possible to fit a multilevel GPD in STAN to negated data with varying threshold. If yes, why all the chains are not initializing? Also, is it possible to introduce the threshold as a parameter in parameter block along with scale and shape parameter instead of as data? Data has a ragged structure. I am attaching the stan code, r script and data.
data_SG.csv (14.4 KB)
multilevel_model.R (1.9 KB)
stan_code.stan (1.6 KB)

``````functions {
real gpareto_lpdf(vector y, real ymin, real xi, real sigma, real s) {
// generalised Pareto log pdf
real N = s;
real inv_xi = inv(xi);
if (xi<0 && max(y-ymin)/sigma > -inv_xi)
reject("xi<0 and max(y-ymin)/sigma > -1/xi; found xi, sigma =", xi, sigma);
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(xi) > 1e-15)
return -(1+inv_xi)*sum(log1p((y-ymin) * (xi/sigma))) -N*log(sigma);
else
return -sum(y-ymin)/sigma -N*log(sigma); // limit xi->0
}

real gpareto_rng(real ymin, real xi, real sigma) {
// generalised Pareto rng
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(xi) > 1e-15)
return ymin + (uniform_rng(0,1)^-xi -1) * sigma / xi;
else
return ymin - sigma*log(uniform_rng(0,1)); // limit xi->0
}
}

data {
// number of groups
int<lower=1> K;
//threshold
vector<lower=-10>[K] ymin;
// number of observations
int<lower=1> N;
//observations
vector [N] y;
// group sizes
int<lower=0> s[K];
// ymax
vector<lower=ymin>[K] ymax;
}

parameters {
//sigma mu
real<lower=0>sigma_mu;
//sigma sd
real<lower=0>sigma_sigma;
//y-specific
vector<lower=0>[K] sigma;
vector<lower=0>[K] xi_offset;
}

transformed parameters{
vector[K] xi = -sigma ./ (ymax-ymin) + xi_offset;
}

model {
int pos;
pos = 1;

//priors
sigma_mu ~ normal(0,3);
sigma_sigma ~ normal(0,3);
xi ~ uniform(-1,1);
sigma ~ normal(sigma_mu, sigma_sigma);

for(i in 1:K){
segment(y, pos, s[i]) ~ gpareto(ymin[i], xi[i], sigma[i], s[i]);
pos = pos + s[i];
}
}

****``````

You should be getting more debug information out of this, like the line of the failure. I’m not sure why you’re not seeing that.

You need to initialize chains at a value of the parameters where the density is positive (i.e., log density is finite). Your model’s too complicated for me to understand easily, but you’re presumably rejecting because of illegal values.

One thing you can do if you can’t get more debug info out of the command (I don’t know how to do that), then you can write print() statements into your code where you can do `print("place marker", target())` and see when `target()` becomes -infinity to try to diagnose where things are going wrong.

1 Like

Sir, thank you very much for taking the time to read my question.

I am a beginner in stan and do not have much coding knowledge. In my data set, 3 subgroups have fewer observations (less than 10). After removing those, the problem got fixed, and all the chains are running and mixing well. Is that might have been the issue?

Further, I would like your comment on my second question “Is it possible to introduce the threshold as a parameter in the parameter block along with scale and shape parameters instead of as data?