Having trouble with convergence of shared parameters in my Joint model

[Olive Buhule]

I am running simulations to test my joint model with the codes below. All parameters are converging to the true values except for the shared parameters (gama and r). I have tried different priors but it's not helping.  Any help will be highly appreciated.

data{
int<lower=1> Nobs; // number of observations
real time[Nobs];
real y1[Nobs];// Continuous longitudinal outcome
real K;
int<lower=1> N; // number of subjects
int<lower=1,upper=N> id[Nobs]; // subject ID
int<lower=0,upper=1>  y2[N];// Binary outcome
vector[4] zero; // a vector of zeros
real survt[N]; // survival time

}

parameters{
real delta; // fixeffs
vector[3] beta;// fixeffs
real alpha;// fixeffs
vector[3] gama;// shared parameter for the survival outcome
vector[3] r;// shared parameter for the binary outcome
vector[4] b[N]; // subj raneffs
real pa; // weibull para
real<lower=0> sigma_e; // residual variance
vector<lower=0>[4] sigma_b; // subj variance
corr_matrix[4] Omega_b; // corr matrix for
// subj raneffs

}

transformed parameters{
cov_matrix[4] Sigma_b; // varcov matrix for subj ranefs

for(l in 1:4){
for(c in 1:4){
Sigma_b[l,c] <- sigma_b[l] * sigma_b[c] * Omega_b[l,c];

}
}
}

model{
real mu[Nobs];
real A[Nobs];
real mu3[N];

for(i in 1:Nobs){
A[i]<- time[i]-(delta + b[id[i],4]);

mu[i] <- K+ if_else(A[i]>=0,(beta[1]*(time[i]-(delta + b[id[i],4])) + beta[2]*((time[i]-(delta + b[id[i],4]))*(time[i]-(delta + b[id[i],4]))) + b[id[i],1]*(time[i]-(delta + b[id[i],4]))+ b[id[i],2]*((time[i]-(delta + b[id[i],4]))*(time[i]-(delta + b[id[i],4])))),
                      (beta[3]*(time[i]-(delta + b[id[i],4])) + b[id[i],3]*(time[i]-(delta + b[id[i],4]))));

 }
y1 ~ normal(mu,sigma_e);

for(j in 1:N){
mu3[j] <- exp(-(gama[1]*b[j,1]+ gama[2]*b[j,2]+gama[3]*b[j,3])/pa);
}

increment_log_prob(weibull_log(survt,pa,mu3)); // Jacobian adjustment

for(j in 1:N){
y2 ~ bernoulli_logit(alpha + r[1]*b[j,1]+ r[2]*b[j,2]+r[3]*b[j,3]);

}

b ~ multi_normal(zero,Sigma_b);
 

pa ~ gamma(1, 1.0E-3);
sigma_b ~ gamma(1.5,1.0E-4);
sigma_e ~ gamma(1.5,1.0E-4);
Omega_b ~ lkj_corr(2.0);
delta ~ normal(0, 100);
beta ~ normal(0, 100);

}

Thanks in advance.

Olive

[Bob Carpenter]

I can't quite decompose that big fat expression.

If parameters don't converge, it's usually because
of identifiability issues. You want to look at what they
are correlated with.

Also, you probably pretty much never want to use if_else()
as it evaluates both expressions (with all the derivative
calculation overhead) but only returns one.

So rather than

x <- a + if_else(c, d, e);

you want to write

if (c)
x <- a + d;
else
x <- a + e;

Also, those long-tailed gamma priors can be problematic both
statistically and computationally.

I would strongly urge you to use the correlation matrix Cholesky
factor parameterization --- it's much much more efficient. There's
a description in the manual chapter on regression.

And there are a few loops that could be vectorized by transposing
your data representations.

None of the efficiency concerns will help with convergence, though.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Olive Buhule]

Thank you, Bob! I think the problem is with  the random effects  not being estimated  correctly, because it is their coefficients that are not converging. I am fitting a shared random effects model.   

Thanks, again!

Olive