Positive_ordered constraint in brms

Interfaces brms
1

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:

y_i=dmax-2( dmax-dtheta)/(exp(s_0 (age_i-theta) )+exp(s_1 (age_i-theta) ) )+e_i

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.

1 Like
2

I think you are stretching brms to its limits - it was never primarily intended for these hugely non-linear models, so you are at the edge of what the abstraction brms uses let’s you do easily… I think that you are most likely better served by writing this model in Stan.

With that said, I think you might be able to hack around a bit. I see 3 ways:

You can do what Stan does for positive_ordered under the hood - have s1 with 0 lower bound, define s2_raw with lower bound of 0 and then use (s1 + s2_raw) whenever you would use s2.

Alternatively, it might be possible to pass lb = "s1" in the prior for s2… But not completely sure brms will eat this up.

Finally, the most hacky way would be by using stanvars to define a new parameter with arbitrary Stan code and then somehow convince brms to use… Maybe by defining a custom likelihood that takes the postive_ordered parameter… But we are getting into really hacky territory here.

Best of luck with your model!

3

Thanks for the reply. The following solution worked “…what Stan does for positive_ordered under the hood - have s1 with 0 lower bound, define s2_raw with lower bound of 0 and then use (s1 + s2_raw) whenever you would use s2.” In the example shown below, parameters s1 and s2 are labeled as s0 and s1.

The only thing is that it requires editing of brms generated stan code to avoid duplication error for the b_s1 parameter. Is there a way to comment out a portion of the stan code or to redefine a parameter within the brms?

library(tidyverse)
library(rstan)
library(brms)
library(stringi)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE, javascript = FALSE) 

# 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)


# 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)
)

# Stan will set uniform priors on the s1_raw parameter. 
# Alternatively, prior on s1_raw can be set as described here https://rdrr.io/cran/brms/man/stanvar.html

stanvars = stanvar(scode = "vector<lower=0>[1] s1_raw;", block = "parameters")+
           stanvar(scode = "vector<lower=0>[1] b_s1=b_s0+s1_raw;", block = "tparameter")

# Data 
sdata <- make_standata(bf, data=df, prior = bp, stanvars= stanvars)

# Code
scode <- make_stancode(bf, data=df, prior = bp , stanvars= stanvars)

# We now need to edit the scode otherwise brms will throw the following error 
# "Duplicate declaration of variable, name=b_s1;...."


# Edit scode to comment out the originally defined parameter b_se i.e., vector<lower=0>[K_s1] b_s1;
# library(stringi)
edited_scode <- stri_replace_all_fixed(scode, 
                "vector<lower=0>[K_s1] b_s1", 
                "// vector<lower=0>[K_s1] b_s1;")

# write edited_scode as .stan file 
writeLines(edited_scode, "edited_scode.stan")

# if want to check the edited_scode, use cat(edited_scode).
# Alternatively, write it in the directory and then view and then delete it 
# sink("edited_scode.stan", type="output")
# writeLines(edited_scode)
# sink()
# file.show("edited_scode" , delete.file=TRUE)


# Re-compile the model using edited_scode via rstan
rstan_fit_edited_scode <- stan(file = "edited_scode.stan",data = sdata, chains = 0)

# Create an empty brms object -> Set empty = TRUE
brms_fit_edited_scode <- brm(bf, prior = bp, data=df, stanvars= stanvars, 
                         seed = 123, chains = 2, iter = 2000, cores = 2, 
                         empty = TRUE)

# Now replace scode with edited_scode in brms 
brms_fit_edited_scode$fit = rstan_fit_edited_scode

# Finally, fit the brms model with edited_scode -> Set recompile = FALSE
brms_fit_edited_scode = update(brms_fit_edited_scode, stanvars= stanvars,
                        seed = 123, chains = 2, iter = 2000, cores = 2,
                        recompile = FALSE)


summary(brms_fit_edited_scode)

Thanks

1 Like
4

I don’t think that’s possible, but note that my first proposal was actually simpler (sorry, maybe the additional options were just confusing). Just to make things clear - we start with a simpler model:

bf <- bf(height ~ exp(s0) + exp(s1) + s1,
         s0     ~ 1,
         s1     ~ 1,
         nl = TRUE)

To force s0 and s1 to be positive ordered (with s0 smaller) you could IMHO (not tested) replace all instances of s1 in the formula with s0 + s1_raw and add lower bound of zero for both s0 and s1_raw:

bf <- bf(height ~ exp(s0) + exp(s0 + s1_raw) + (s0 + s1_raw),
         s0     ~ 1,
         s1_raw     ~ 1,
         nl = TRUE)

bp = c(
  prior(cauchy(0, 1),           nlpar = "s0",     lb=0),
  prior(cauchy(0, 1),           nlpar = "s1_raw",     lb=0),
)

Does that make sense?

5

Worked perfectly! Many thanks.