Non-linear hierarchical fit: First order compartment model in cross-over design

[Dieter Menne]

We have recorded hormone concentration time series from 17 subjects, each under 4 treatments in a cross-over design. See HFit.pdf for the raw data and a first order compartment model fit similar to SSfol in R-package nlme. Treatment em (1..4) for each of the coefficients is the parameter of interest:

y[n] <- D_1 * Ke_1 * Ka_1 * (exp(-Ke_1*minut) - exp(-Ka_1*minut)) / (Ka_1 - Ke_1);

A fit similar to that in 17.3 of ARM (treatment and airport flight simulator) gives reasonable results

lKe_1_emn <- lKe_1_em[emn]  + lKe_1_subj[subjn]; # and similar for the other parameters

but this does not reflect the nested structure of the design; I could not find a nested formulation that converges comparably.

Any ideas are welcome. The example files attached are self-contained, all data are in .Rdata

[Bob Carpenter]

I think this may have received no attention because everyone's
too lazy to unzip an archive. I'm attaching the model below a repeat
of the original question.

I think the main problem may be in these lines from the model block:

D_1_subj ~ normal(0, sigmaD_1_subj);
D_1_em ~ normal(100,100);

emn <- em[n];
subjn <- subj[n];

D_1_emn <- D_1_em[emn] + D_1_subj[subjn]; # D_1_subj not required

D_1_s[n] ~ normal(D_1_emn, sigmaD_1_s);

This is going to cause problems because you can add to D_1_em and
subtract from D_1_subj and get the same answer. So you have to make
sure your priors are strong enough to control inference. Making models
like these hierarchical can often help (as in the IRT case). That 100
you have for the std dev in D_1_em may be problematic depending on the scale
of the problem. Making the prior for D_1_em hierarchical for both mean and
sd may help with fitting.

- Bob

data {
int<lower=0> N; # number of data points
int<lower=0> n_em; # number of treatments, fixed to 4
int<lower=0> n_series; # number of time series
int<lower=0> n_subj; # number of subjects
vector[N] conc;
vector[N] minute;
int subj[n_series];
int em[n_series]; # Mapping of series to treatments
int series[N] ; # time series
}

parameters {

# Per series
vector<lower=0, upper=0.01>[n_series] Ke_1_s;
vector<lower=0, upper=0.08>[n_series] Ka_1_s;
vector<lower=0, upper=1000>[n_series] D_1_s;

real<lower = 0> sigmalKe_1_s;
real<lower = 0> sigmalKa_1_s;
real<lower = 0> sigmaD_1_s;

vector<lower = -3, upper = 3>[n_subj] lKe_1_subj;
vector<lower = -3, upper = 3>[n_subj] lKa_1_subj;
vector<lower = -20, upper = 20>[n_subj] D_1_subj; # not required

vector<lower = -7.5>[n_em] lKe_1_em;
vector<lower = -6.5>[n_em] lKa_1_em;
vector<lower = 0>[n_em] D_1_em;

real<lower=0, upper = 2> sigmalKe_1_subj;
real<lower=0, upper = 2> sigmalKa_1_subj;
real<lower=0, upper = 20> sigmaD_1_subj;

real<lower = 0> sigma;
}

transformed parameters {
vector[N] y;
real<lower=0> constr;

for (n in 1:N) {
real minut;
real Ke_1;
real Ka_1;
real D_1;
int is;

minut <- minute[n];
is <- series[n];

Ke_1 <- Ke_1_s[is];
Ka_1 <- Ka_1_s[is];
D_1 <- D_1_s[is];
# occasionally, a chain gets stuck without this
constr <- Ka_1 - Ke_1; # Constraint to avoid flips

y[n] <- D_1 * Ke_1 * Ka_1 * (exp(-Ke_1*minut) - exp(-Ka_1*minut)) / (Ka_1 - Ke_1);
}
}

model {

lKe_1_em ~ normal(-6,2);
lKa_1_em ~ normal(-5,2);
D_1_em ~ normal(100,100);

# sigmalKe_1_s ~ cauchy(0,.5);
# sigmalKa_1_s ~ cauchy(0,.5);
# sigmaD_1_s ~ cauchy(100,100);

# sigmalKe_1_subj ~ cauchy(0,.5); # Use uniform, faster
# sigmalKa_1_subj ~ cauchy(0,.5);
# sigmaD_1_subj ~ cauchy(0,20);

lKe_1_subj ~ normal(0, sigmalKe_1_subj) ; # 0.8
lKa_1_subj ~ normal(0, sigmalKa_1_subj) ; # 0.5
D_1_subj ~ normal(0, sigmaD_1_subj) ; # 5

for (n in 1:n_series){
real lKe_1_emn;
real lKa_1_emn;
real D_1_emn;
int subjn;
int emn;

emn <- em[n];
subjn <- subj[n];

lKe_1_emn <- lKe_1_em[emn] + lKe_1_subj[subjn];

lKa_1_emn <- lKa_1_em[emn] + lKa_1_subj[subjn];
D_1_emn <- D_1_em[emn] + D_1_subj[subjn]; # D_1_subj not required
# D_1_emn <- D_1_em[emn];

Ke_1_s[n] ~ lognormal(lKe_1_emn, sigmalKe_1_s);
Ka_1_s[n] ~ lognormal(lKa_1_emn, sigmalKa_1_s);
D_1_s[n] ~ normal(D_1_emn, sigmaD_1_s);
}
conc ~ normal(y, sigma);
}

[Dieter Menne]

Thanks, Bob,

and apologies for not replying earlier. I was dragged into another project and escaped only yesterday.

Re Zip: I thought it was more polite when everything is zipped in one file, but it looks like that was the wrong idea.

Dieter

[Bob Carpenter]

Sorry --- I didn't mean to give you the impression that zip's
the wrong idea for packaging (gzipped tarballs are also OK for the
devs, but please please, no other oddball formats, especially
Windows ones because none of the devs use Windows).

The problem is more that it seems to entail a lot of work and we
don't have a dedicated staff. So like on most mailing lists like
this, difficult (in terms of time required) questions languish and
the easy ones get answered right away.

If you have further questions, please let us know.

- 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.
> For more options, visit https://groups.google.com/d/optout.