Hi,

I am implementing the Preece-Baines growth model using brms package. The stan implementation of this model has been discussed earlier (Growth Curve Model Compilation Error). The model is defined as:

Where y_i is height at age_i and dmax represents the maximum size (i.e., height). The dtheta is size during the peak growth period (theta) and s_0 and s_1 denote pre-peak and peak growth rates. The e_i are residuals.

The key features of this model are: 1) all five growth parameters (dmax, dtheta, s_0 , s_1 and theta) are positive with dmax > dtheat, and 2) the peak growth rate is higher than the pre-peak growth rate i.e., 0 <s_0<s_1

Therefore, it is natural to put these constraints while fitting the model. For the first condition, it is straight forward to put lower bound as <lower=0>. However, I could not figure out how to use positive_ordered in brms which is otherwise easy to implement in rstan.

Below I present a simple example of fitting the above model using rstan (satisfying both conditions 1 and 2) and brms (satisfying condition 1).

```
library(tidyverse)
library(rstan)
library(brms)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
# data
df <- read.table(text = "
age boy01 boy02 boy03 boy04 boy05 boy06 boy07 boy08 boy09 boy10 boy11 boy12 boy13 boy14 boy15 boy16 boy17 boy18
1 1 81.3 76.2 76.8 74.1 74.2 76.8 72.4 73.8 75.4 78.8 76.9 81.6 78 76.4 76.4 76.2 75 79.7
2 1.25 84.2 80.4 79.8 78.4 76.3 79.1 76 78.7 81 83.3 79.9 83.7 81.8 79.4 81.2 79.2 78.4 81.3
3 1.5 86.4 83.2 82.6 82.6 78.3 81.1 79.4 83 84.9 87 84.1 86.3 85 83.4 86 82.3 82 83.3
4 1.75 88.9 85.4 84.7 85.4 80.3 84.4 82 85.8 87.9 89.6 88.5 88.8 86.4 87.6 89.2 85.4 84 86.5
5 2 91.4 87.6 86.7 88.1 82.2 87.4 84.2 88.4 90 91.4 90.6 92.2 87.1 91.4 92.2 88.4 85.9 88.9
6 3 101. 97 94.2 98.6 89.4 94 93.2 97.3 97.3 100. 96.6 99.3 96.2 101. 101. 101 95.6 99.4
7 4 110. 105. 100. 104. 96.9 102. 102. 107. 103. 111 105. 106. 104 106. 110. 107. 102. 104.
8 5 116. 112. 107. 111 104. 109. 109 113. 108. 118. 112 113. 111 113. 117. 115. 109. 112.
9 6 122. 119. 112. 116. 111. 116. 117. 119. 114. 126. 119. 120. 117. 120. 122. 121. 118. 119
10 7 130 125 119. 123. 116. 122. 123. 126. 120. 131. 125. 127. 124. 129. 130. 128 125. 128", header = T, row.names = 1);
df <- df %>%
gather(boy, height, -age)
# head(df)
# rstan
rstan_code <- "
data {
int N;
vector[N] y;
vector[N] age;
}
parameters {
real<lower=0> dmax;
real<lower=0> dtheta;
real<lower=0> theta;
positive_ordered[2] s;
real<lower=0> sigma;
}
transformed parameters {
real<lower=0> s0;
real<lower=0> s1;
s0 = s[1];
s1 = s[2];
}
model {
vector[N] mu;
for (i in 1:N)
mu[i] = dmax - 2 * (dmax - dtheta) / (exp(s0 * (age[i] - theta)) + (exp(s1 * (age[i] - theta))));
dmax ~ normal(max(y), 5);
dtheta ~ normal(max(y)*.95, 5);
theta ~ normal(7, 2);
s ~ cauchy(0, 1);
y ~ normal(mu, sigma);
}
"
# Fit model
rstan_fit <- stan(
model_code = rstan_code,
data = list(
N = nrow(df),
y = df$height,
age = df$age),
seed = 123, chains = 2, iter = 2000, cores = 2)
summary(rstan_fit, pars = c("dmax", "dtheta", "theta", "s0", "s1", "sigma"))$summary
# brms
bf <- bf(height ~ dmax - 2 * (dmax - dtheta) / (exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))),
dmax ~ 1,
dtheta ~ 1,
s0 ~ 1,
s1 ~ 1,
theta ~ 1,
nl = TRUE)
bp = c(
prior(normal(max(Y), 5), nlpar = "dmax", lb=0),
prior(normal(max(Y)*0.95, 5), nlpar = "dtheta", lb=0),
prior(cauchy(0, 1), nlpar = "s0", lb=0),
prior(cauchy(0, 1), nlpar = "s1", lb=0),
prior(normal(7, 2), nlpar = "theta", lb=0)
)
scode <- make_stancode(bf, data=df, prior = bp)
brms_fit <- brm(bf, prior = bp, data=df,
seed = 123, chains = 2, iter = 2000, cores = 2)
Warning: The largest R-hat is 1.83, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess
summary(brms_fit)
```

The first model (rstan) did not produce any warning.

Clearly it seems that putting positive_ordered constraints helps in identifying the model.