Growth Curve Model Compilation Error

Modeling
1

I have data the height of a human (y_i) at time points (t_i), i = 1, \dots, n, with t_1 < t_2 < \dots < t_{n - 1} < t_n. The model is

Y_i | \mu_i, \sigma^2 \sim N(\mu_i, \sigma^2), i = 1, \dots, n

\mu_i = h_1 - \dfrac{2(h_1 - h_{\theta})}{\exp(s_0(t_i - \theta)) + \exp(s_1(t_i - \theta))}

The model has the following 6 parameters:

  1. The height h_1 is the adult height of the individual (i.e. when t becomes large).
  2. The time point \theta, and the (expected) height h_\theta at that time point. Note, we would expect h_\theta \le h_1 or, in other words, 0 \le \theta < t_n.
  3. Two rate parameters s_0 and s_1, which should be positive. Note that these are actually not distinguishable, i.e. (s_0, s_1) = (0.4, 0.8) and (s_0, s_1) = (0.8, 0.4) will yield the same fits (if all other parameters are set to the same value). To make these rate parameters identifiable, it would be natural to impose constraints of the form 0 < s_0 < s_1.
  4. The error variability \sigma^2.

My Stan code is as follows:

 data {
  int<lower=1> n;
  ordered[n] t; // age
  ordered[n] y; // height of human
 }
 
 parameters {
  positive_ordered[2] h;
  real<lower=0, upper=t[n-1]>theta;
  positive_ordered[2] s; // we will put flat priors on these rate parameters
  real<lower=0>sigma;
 }
 
 model {
 
  theta ~ uniform(0, 1);
  h[1] ~ uniform(0, y[n]);
  h[2] ~ normal(180, 20);
  sigma ~ cauchy(0, 5);
  
  y ~ normal(h[2] - (2*(h[2] - h[1]) / (exp(s[1]*(t - theta)) + exp(s[2]*(t - theta)))), sigma);
 }

I am getting the following error when I try to compile this model:

20%20am
Screen Shot 2018-05-21 at 11.21.20 am.png1258×1046 79.7 KB

I’ve reviewed the code, but I’m unable to find where the error is. I would appreciate it if people could please take the time to help me with this error.

2

It tells you the error is on line 21 and it tells you the problem is that there is no implementation of a real number divided by a vector. Do

h[2] - (2*(h[2] - h[1]) * inv(exp(s[1]*(t - theta)) + exp(s[2]*(t - theta)))
1 Like
3

Thank you for the assistance. This seems to have worked, but I’m now getting another error when sampling:

```{r}
child = read.csv("child.csv")
t = child$Age
y = child$Height
n = length(t)
data.in <- list(t = t, y = y, n = n)
model.fit99 <- sampling(mod94, data=data.in)
```

SAMPLING FOR MODEL ‘stan-15d9195d6b44’ NOW (CHAIN 1).
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can’t start sampling from this initial value.
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can’t start sampling from this initial value.
…
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can’t start sampling from this initial value.
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can’t start sampling from this initial value.
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
error occurred during calling the sampler; sampling not done

Any idea what I’ve done wrong here?

4

Try calling it with better initial values by specifying init_r = 0.5 or so.

5

Thanks for the advice.

I tried init_r = 0.5, which failed, and I tried with init_r = 3, which worked, but the n_eff values are very low:

child = read.csv("child.csv")
t = child$Age
y = child$Height
n = length(t)
data.in <- list(t = t, y = y, n = n)
model.fit99 <- sampling(mod94, data=data.in, init_r = 3)

There were 2620 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
Runtime warnings and convergence problems the pairs() plot to diagnose sampling problems

print(model.fit99, pars=c("theta", "h[1]", "h[2]", "sigma", "s[1]", "s[2]"), probs=c(0.1,0.5,0.9), digits=5)

17%20pm
17%20pm1130×564 73.3 KB

I then took the advice of the warning message and set adapt_delta = 0.99:

child = read.csv("child.csv")
t = child$Age
y = child$Height
n = length(t)
data.in <- list(t = t, y = y, n = n)
stan(control = list(adapt_delta = 0.99))
model.fit99 <- sampling(mod94, data=data.in, init_r = 3)

print(model.fit99, pars=c("theta", "h[1]", "h[2]", "sigma", "s[1]", "s[2]"), probs=c(0.1,0.5,0.9), digits=5)

Which still produces very low n_eff values:

52%20pm
52%20pm1126×578 73.9 KB

I have attached my data:

Child.csv (301 Bytes)

6

Your model has characteristics that make it hard for NUTS to sample efficiently. It is a bit of a surprise it works at all. I would do

parameters {
  real<lower=0,upper=y[n]> h_theta;
  real<lower=h_theta> h_1;
  real<lower=0, upper=t[n-1]>theta;
  positive_ordered[2] s; // it s a bad idea to put flat priors on these rate parameters
  real<lower=0>sigma;
}
model {
  target += normal_lpdf(h1 | 180, 20);
  target += cauchy_lpdf(sigma | 0, 5);
  target += normal_lpdf(y | h[2] - (2*(h[2] - h[1]) * inv(exp(s[1]*(t - theta)) + exp(s[2]*(t - theta))), sigma)
}
1 Like
7

Thank you for the assistance.

I used your feedback and changed a couple of extra things, and everything seems to be fine now.

Thanks again!