R has an inbuilt data set called Orange, which measures the circumference of 5 orange trees (in mm) at various ages (measured in days since 1968). I’m using this data set for this example model.
My model is as follows:
Circumference at time t is c(t) = \dfrac{\beta_1}{1 + \exp(-(t - \beta_2) / \beta_3)}, where all parameters are positive.
\beta_1 is the circumference at maturity (when t is large), \beta_2 is the time point at which half the final circumference is reached (i.e., if t = \beta_2, then c(t) = \dfrac{\beta_1}{2}), and \beta_3 measures the rate of growth of the tree.
In this thread, I managed to develop an unpooled version of this nonlinear regression model:
```{stan output.var="non_pooled"}
data {
int<lower=0> n;
vector<lower=0>[n] t;
int<lower=0> groups;
matrix<lower=0>[7, 5] y;
}
parameters {
matrix<lower=0>[3, groups] beta;
real<lower=0> sigma;
}
transformed parameters {
matrix[7, groups] mu;
for (j in 1:groups) {
for (i in 1:7) {
mu[i, j] = beta[1, j] * inv(1 + exp(-(t[i] - beta[2, j]) / beta[3, j]));
}
}
}
model {
for (j in 1:groups) {
for (i in 1:7) {
y[i, j] ~ normal(mu[i, j], sigma);
}
beta[1, j] ~ normal(170, 5);
beta[2, j] ~ normal(85, 20);
beta[3, j] ~ normal(500, 500);
}
sigma ~ student_t(2, 0, 1);
}
```
```{r}
data.in <- list(t = Orange$age, n = length(Orange$circumference), y = matrix(Orange$circumference, 7, 5), group = as.integer(Orange$Tree), groups <- length(unique(Orange$Tree)))
model.fit2 <- sampling(non_pooled, data=data.in)
```
```{r}
print(model.fit2, digits=4, pars=c("mu", "beta", "sigma"))
```
Explanation
We have groups = 5 groups and n = 35 observations within those groups.
I created the 7 \times 5 matrix y, where the columns denote the groups and the rows denote the observations within those groups. Therefore, we have 7 \times 5 = 35 total observations, as required.
Since this is an unpooled model, we need to generate parameters \beta_1^j, \beta_2^j, \beta_3^j for each group j = 1, \dots, 5. Therefore, beta is a 3 \times 5 matrix.
mu is a 7 \times 5 matrix, where the columns denote the group and the rows denote the observations. This is because each group has 7 observations, and so there are 7 observations within each group. This matrix contains the conditional means at various time points across the groups.
Hierarchical Model
I now want to try and go one step further and build a hierarchical nonlinear regression model that puts hierarchical priors onto the three parameters. So if \beta_i^j , i = 1, 2, 3, denotes the three regression parameters for group/tree j, then we would assume that they come from the following populations:
\beta_1^j \sim \mathcal{N}(\mu_{\beta_1}, \sigma^2_{\beta_1}), \ j = 1, \dots, 5
\beta_2^j \sim \mathcal{N}(\mu_{\beta_2}, \sigma^2_{\beta_2}), \ j = 1, \dots, 5
\beta_3^j \sim \mathcal{N}(\mu_{\beta_3}, \sigma^2_{\beta_3}), \ j = 1, \dots, 5
And with appropriate priors on the hyper-parameters \mu_{\beta_i} and \sigma^2_{\beta_i}, \ i = 1, 2, 3
I’ve tried to build my hierarchical model using (1) Andrew Gelman’s textbook Bayesian Data Analysis, Third Edition and (2) Richard McElreath’s textbook Statistical Rethinking:
(from Bayesian Data Analysis, Third Edition, by Andrew Gelman)
And a couple of other Stan examples that I was able to find online.
Despite this, the atrocious n_eff and Rhat statistics of my model seem to suggest that I’ve done this incorrectly.
My code is as follows:
```{stan output.var="hierarchical"}
data {
int<lower=0> n;
vector<lower=0>[n] t;
int<lower=0> groups;
matrix<lower=0>[7, 5] y;
}
parameters {
matrix<lower=0>[3, groups] beta;
vector<lower=0>[3] mu_beta;
vector<lower=0>[3] sigma_beta;
real<lower=0> sigma;
}
transformed parameters {
matrix[7, groups] mu;
for (j in 1:groups) {
for (i in 1:7) {
mu[i, j] = beta[1, j] * inv(1 + exp(-(t[i] - beta[2, j]) / beta[3, j]));
}
}
}
model {
for (j in 1:groups) {
for (i in 1:7) {
y[i, j] ~ normal(mu[i, j], sigma);
}
beta[1, j] ~ normal(mu_beta[1], sigma_beta[1]);
beta[2, j] ~ normal(mu_beta[2], sigma_beta[2]);
beta[3, j] ~ normal(mu_beta[3], sigma_beta[3]);
}
mu_beta[3] ~ normal(170, 5);
mu_beta[2] ~ normal(85, 20);
mu_beta[1] ~ normal(500, 500);
sigma_beta ~ student_t(2, 0, 1);
sigma ~ student_t(2, 0, 1);
}
```
```{r}
data.in <- list(t = Orange$age, n = length(Orange$circumference), y = matrix(Orange$circumference, 7, 5), group = as.integer(Orange$Tree), groups <- length(unique(Orange$Tree)))
model.fit3 <- sampling(hierarchical, data=data.in)
```
```{r}
print(model.fit3, digits=4, pars=c("mu", "beta", "sigma", "mu_beta"))
```
As can be seen, I attempted to adapt my unpooled model to the hierarchical model.
I would greatly appreciate it if people could please take the time to show the corrected code, with perhaps, if they’d be so kind, some accompanying explanation of what the errors were, so that I may have a clear understanding of how to do this myself in the future.
