Growth models with many groups

Hello,

I am interested in fitting growth models to some data with a variety of treatments/groups in the data. Right now my project has 15 groups (genotypes) but I am hoping to get some advice on general scalability. In the example below the smallData.fit model will run and give good results but the completeData.fit model will fail to initialize.

My initial values are confined to be positive by rgamma since I am using lognormal priors, I have also tried inits=0 and setting inits as a list of lists containing values for a, b, and c for each chain (idea taken from here), but the model has consistently failed to initialize when I have more than 2 groups.

Is there a set of guidelines on how many groups a brms model should be able to handle? Has anyone had a similar situation and found a solution besides running separate models on subsets of the data (that option does work but it seems like a poor workaround)?

library(brms)
library(tidyverse)

## Setup simulated data
growthSim <- function(x,a,b,c){ 
  a_r <- a+rnorm(1,mean = 0,sd=10)
  b_r <- b+rnorm(1,mean=0,sd=2)
  c_r <- c+rnorm(1,mean=0,sd=.035)
  return(a_r*exp(-b_r*exp(-c_r*x)))}
x <- 1:25
smallDf <- rbind(do.call(rbind,lapply(letters[1:2], function(L) do.call(rbind, lapply(1:20,function(i) data.frame("sample"=paste0("sample_",i),"treatment"=L,"time"=x,"y"=growthSim(x,200-(5*which(letters==L)),13-(.2*which(letters==L)),.2+(0.001*which(letters==L))),stringsAsFactors = F))))))

completeDf <- rbind(do.call(rbind,lapply(letters[1:15], function(L) do.call(rbind, lapply(1:20,function(i) data.frame("sample"=paste0("sample_",i),"treatment"=L,"time"=x,"y"=growthSim(x,200-(5*which(letters==L)),13-(.2*which(letters==L)),.2+(0.001*which(letters==L))),stringsAsFactors = F))))))

completeDf%>%
  ggplot(aes(time,y,group=interaction(treatment,sample)))+
  geom_line(aes(color=treatment), show.legend=F)

## Define Priors

prior1 <- prior(lognormal(log(130), .25),nlpar = "a") +
  prior(lognormal(log(12), .25), nlpar = "b") + 
  prior(lognormal(log(1.2), .25), nlpar = "c") + 
  prior(student_t(3,0,5), dpar="sigma") +
  prior(gamma(2,0.1), class="nu")

## Run models

smallData.fit <- brm(bf(y ~ a*exp(-b*exp(-c*time)), 
               sigma~s(time,by=treatment), 
               a + b + c ~ 0+treatment,
               autocor = ~arma(~time|sample:treatment,1,1),nl = TRUE),
            family = student, prior = prior1, data = smallDf, iter = 2000, 
            cores = 4, chains = 4, backend = "cmdstanr", 
            control = list(adapt_delta = 0.999,max_treedepth = 20),
            inits = function(){list(b_a=rgamma(2,1),b_b=rgamma(2,1),b_c=rgamma(2,1))})

completeData.fit <- brm(bf(y ~ a*exp(-b*exp(-c*time)), 
               sigma~s(time,by=treatment), 
               a + b + c ~ 0+treatment,
               autocor = ~arma(~time|sample:treatment,1,1),nl = TRUE),
            family = student, prior = prior1, data = completeDf, iter = 2000, 
            cores = 4, chains = 4, backend = "cmdstanr", 
            control = list(adapt_delta = 0.999,max_treedepth = 20),
            inits = function(){list(b_a=rgamma(2,1),b_b=rgamma(2,1),b_c=rgamma(2,1))})

Error messages from trying to make completeData.fit for all chains are:

Chain 4 Rejecting initial value:
Chain 4   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4   Stan can't start sampling from this initial value.
Chain 4 Initialization between (-2, 2) failed after 100 attempts.
Chain 4  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
Warning: Chain 4 finished unexpectedly!

Also this is my first post so if I should provide more information or I have tagged this topic incorrectly please let me know! Thank you!

So does your model run without the grouping/treatments? And what does the simple model’s diagnostics look like (are there errors or warnings?).

The model runs with only one genotype and I have used a lot of very similar models comparing two treatment groups and those work well. The amount of time the models take is pretty variable even with really similar/identical data.

Recently I started seeing these warnings when the Stan code compiles:

Warning in '/var/folders/t6/mhdppcf93w5_21y8nv0s7m_40000gn/T/RtmpODvZUa/model-d7f16a179603.stan', line 43, column 2: Declaration of arrays by placing brackets after a variable name is deprecated and will be removed in Stan 2.32.0. Instead use the array keyword before the type. This can be changed automatically using the auto-format flag to stanc

And this warning has been a constant when using lognormal priors.

Warning message: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound. If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately. Warning occurred for prior b_a ~ lognormal(log(130), 0.25) b_b ~ lognormal(log(12), 0.25) b_c ~ lognormal(log(1.2), 0.25)

Rhats are all 1 to 1.01 with the version of the model using smaller data that will initialize.

I went ahead and ran that example model again just now with 2 chains/cores (computer is doing other things at the moment), I’ll attach the full summary of that model at the end of this comment.

The completeData.fit model always fails to initialize. The only thing I have thought of is that there is some prior on an intercept term not shown by get_prior() and that part of the model is using default inits on that variable which then causes the failure. But I would have thought inits=0 would solve the problem if that were the case.

> summary(smallData.fit)
 Family: student 
  Links: mu = identity; sigma = log; nu = identity 
Formula: y ~ a * exp(-b * exp(-c * time)) 
         autocor ~ tructure(list(), class = "formula", .Environment = <environment>)
         sigma ~ s(time, by = treatment)
         a ~ 0 + treatment
         b ~ 0 + treatment
         c ~ 0 + treatment
   Data: smallDf (Number of observations: 1000) 
  Draws: 2 chains, each with iter = 1000; warmup = 0; thin = 1;
         total post-warmup draws = 2000

Smooth Terms: 
                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sds(sigma_stimetreatmenta_1)     6.81      2.39     3.83    12.64 1.00      617
sds(sigma_stimetreatmentb_1)     5.30      1.84     2.90    10.10 1.00      630
                             Tail_ESS
sds(sigma_stimetreatmenta_1)     1061
sds(sigma_stimetreatmentb_1)     1116

Population-Level Effects: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma_Intercept              1.91      0.04     1.83     1.99 1.00     1799
a_treatmenta               185.62      3.44   178.97   192.58 1.00     1809
a_treatmentb               185.35      2.77   180.04   190.95 1.00     1725
b_treatmenta                13.20      0.29    12.66    13.82 1.00     1524
b_treatmentb                11.83      0.34    11.24    12.54 1.00     1481
c_treatmenta                 0.20      0.00     0.20     0.21 1.00     1462
c_treatmentb                 0.20      0.00     0.20     0.21 1.00     1329
sigma_stime:treatmenta_1    31.97      4.74    23.00    41.27 1.00     1431
sigma_stime:treatmentb_1    24.07      4.23    16.03    32.59 1.00     1065
                         Tail_ESS
sigma_Intercept              1527
a_treatmenta                 1306
a_treatmentb                 1495
b_treatmenta                 1225
b_treatmentb                 1406
c_treatmenta                 1182
c_treatmentb                 1357
sigma_stime:treatmenta_1     1547
sigma_stime:treatmentb_1     1133

Family Specific Parameters: 
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu     4.80      0.76     3.58     6.53 1.00     1655     1420

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
1 Like

In a new version of brms I now have a new warning message which makes me think this could be related to the splines initialization.

Missing init values for the following parameters:
 - chain 1: Intercept_sigma, bs_sigma, zs_sigma_1_1, sds_sigma_1_1, zs_sigma_2_1, sds_sigma_2_1, zs_sigma_3_1, sds_sigma_3_1, nu
 - chain 2: Intercept_sigma, bs_sigma, zs_sigma_1_1, sds_sigma_1_1, zs_sigma_2_1, sds_sigma_2_1, zs_sigma_3_1, sds_sigma_3_1, nu
 - chain 3: Intercept_sigma, bs_sigma, zs_sigma_1_1, sds_sigma_1_1, zs_sigma_2_1, sds_sigma_2_1, zs_sigma_3_1, sds_sigma_3_1, nu
 - chain 4: Intercept_sigma, bs_sigma, zs_sigma_1_1, sds_sigma_1_1, zs_sigma_2_1, sds_sigma_2_1, zs_sigma_3_1, sds_sigma_3_1, nu

So I tried setting the inits for those variables explicitly with the list below but that did not change anything.

init = function(){list(b_a=rgamma(2,1),b_b=rgamma(2,1),b_c=rgamma(2,1), 
                                   Intercept_sigma=runif(2,1,3), bs_sigma=runif(2,1,3),
                                   zs_sigma_1_1=runif(2,1,3),sds_sigma_1_1=runif(2,1,3),
                                   zs_sigma_2_1=runif(2,1,3), sds_sigma_2_1=runif(2,1,3), nu=runif(2,1,3))}

I also stumbled into a comment on a post about change-point/piecewise models suggesting recoding categorical variables as dummy variables (although that was in response to a different error than I have had), I gave that a shot with a subset of the data and did not see any improvements.

dummyCompleteDf<-completeDf%>%
  mutate(sample2=sample, treatment2=treatment)%>%
  pivot_wider(names_from = treatment2, values_from=sample2, names_prefix = "trt_")%>%
  mutate(across(.cols = contains("trt_"), .fns = ~ifelse(is.na(.x),0,1)))
dummyData.fit <- brm(bf(y ~ a*exp(-b*exp(-c*time)), 
                              sigma~s(time,by=treatment), 
                              a + b + c ~ 0+trt_a+trt_b+trt_c,
                              autocor = ~arma(~time|sample:(treatment),1,1),
                           nl = TRUE),
                           family = student, prior = prior1, data = dummyCompleteDf%>%filter(treatment %in% c("a", "b", "c")), iter = 1000,  # I started by trying with just a few groups 
                           cores = 4, chains = 4, backend = "cmdstanr", 
                           control = list(adapt_delta = 0.999,max_treedepth = 20),
                           init = 0) # tried with variety of inits

Ended up solving this while working with initialization on an entirely different model. Incorrect length of initial values was leaving the remaining levels of treatment initialized randomly instead of recycling values. Solution was to make enough initial values with rgamma(15,1) for each nlpar. This particular model is very slow with 15 groups, but it does initialize and start to run. It is more than 15 times slower (took ~10 hours) than fitting a model to a single condition (6-10 minutes per) but it does work.

length(unique(completeDf$treatment))

completeData.fit <- brm(bf(y ~ a*exp(-b*exp(-c*time)), 
                           sigma~s(time,by=treatment), 
                           a + b + c ~ 0+treatment,
                           autocor = ~arma(~time|sample:treatment,1,1),nl = TRUE),
                        family = student, prior = prior1, data = completeDf, iter = 2000, 
                        cores = 4, chains = 4, backend = "cmdstanr", 
                        control = list(adapt_delta = 0.999,max_treedepth = 20),
                        inits = function(){list(b_a=rgamma(15,1),b_b=rgamma(15,1),b_c=rgamma(15,1))})

completeData.fit
Family: student 
  Links: mu = identity; sigma = log; nu = identity 
Formula: y ~ a * exp(-b * exp(-c * time)) 
         autocor ~ tructure(list(), class = "formula", .Environment = <environment>)
         sigma ~ s(time, by = treatment)
         a ~ 0 + treatment
         b ~ 0 + treatment
         c ~ 0 + treatment
   Data: completeDf (Number of observations: 7500) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Smooth Terms: 
                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sds(sigma_stimetreatmenta_1)     5.34      1.98     2.83    10.33 1.00     1828
sds(sigma_stimetreatmentb_1)     5.95      2.09     3.24    11.12 1.00     1696
sds(sigma_stimetreatmentc_1)     5.54      1.90     3.00    10.27 1.00     1956
sds(sigma_stimetreatmentd_1)     6.95      2.43     3.80    12.89 1.00     1546
sds(sigma_stimetreatmente_1)     6.21      2.22     3.35    11.78 1.00     1713
sds(sigma_stimetreatmentf_1)     5.23      1.90     2.76    10.09 1.00     2042
sds(sigma_stimetreatmentg_1)     5.68      1.99     3.06    10.79 1.00     1941
sds(sigma_stimetreatmenth_1)     5.71      1.98     3.12    10.55 1.00     1952
sds(sigma_stimetreatmenti_1)     5.34      2.08     2.81    10.55 1.00     1764
sds(sigma_stimetreatmentj_1)     5.31      2.00     2.78    10.28 1.00     1715
sds(sigma_stimetreatmentk_1)     4.19      1.56     2.20     8.13 1.00     2235
sds(sigma_stimetreatmentl_1)     4.86      1.93     2.49     9.34 1.00     1711
sds(sigma_stimetreatmentm_1)     5.05      1.79     2.74     9.33 1.00     1717
sds(sigma_stimetreatmentn_1)     4.82      1.87     2.48     9.43 1.00     1903
sds(sigma_stimetreatmento_1)     4.84      1.76     2.55     9.32 1.00     2226
                             Tail_ESS
sds(sigma_stimetreatmenta_1)     2142
sds(sigma_stimetreatmentb_1)     2391
sds(sigma_stimetreatmentc_1)     2737
sds(sigma_stimetreatmentd_1)     2390
sds(sigma_stimetreatmente_1)     2165
sds(sigma_stimetreatmentf_1)     2749
sds(sigma_stimetreatmentg_1)     2625
sds(sigma_stimetreatmenth_1)     2550
sds(sigma_stimetreatmenti_1)     2365
sds(sigma_stimetreatmentj_1)     2273
sds(sigma_stimetreatmentk_1)     2501
sds(sigma_stimetreatmentl_1)     2065
sds(sigma_stimetreatmentm_1)     2565
sds(sigma_stimetreatmentn_1)     2375
sds(sigma_stimetreatmento_1)     2638

Population-Level Effects: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma_Intercept              1.80      0.02     1.76     1.83 1.00     4957
a_treatmenta               190.53      3.41   184.11   197.33 1.00     4993
a_treatmentb               190.41      1.83   186.92   194.13 1.00     5864
a_treatmentc               178.83      2.25   174.53   183.30 1.00     5327
a_treatmentd               177.29      2.73   171.97   182.68 1.00     5276
a_treatmente               173.74      2.25   169.30   178.15 1.00     5924
a_treatmentf               167.96      2.14   163.86   172.30 1.00     5554
a_treatmentg               162.33      2.00   158.47   166.20 1.00     5333
a_treatmenth               155.09      2.18   150.74   159.46 1.00     4986
a_treatmenti               152.76      1.89   149.13   156.46 1.00     5208
a_treatmentj               147.19      1.85   143.65   151.03 1.00     5312
a_treatmentk               142.77      2.06   138.87   147.01 1.00     5769
a_treatmentl               138.30      2.22   133.77   142.61 1.00     5581
a_treatmentm               134.28      1.34   131.71   136.94 1.00     5673
a_treatmentn               130.13      1.57   127.01   133.16 1.00     4558
a_treatmento               121.21      1.63   118.11   124.45 1.00     6151
b_treatmenta                12.67      0.33    12.10    13.36 1.00     5014
b_treatmentb                11.90      0.26    11.43    12.44 1.00     5280
b_treatmentc                12.09      0.26    11.61    12.63 1.00     5465
b_treatmentd                12.74      0.25    12.28    13.27 1.00     5456
b_treatmente                12.07      0.33    11.48    12.77 1.00     5702
b_treatmentf                10.88      0.35    10.26    11.62 1.00     5328
b_treatmentg                11.74      0.34    11.12    12.48 1.00     5088
b_treatmenth                12.10      0.30    11.53    12.73 1.00     4557
b_treatmenti                11.24      0.29    10.70    11.85 1.00     4516
b_treatmentj                11.74      0.33    11.12    12.44 1.00     4727
b_treatmentk                10.60      0.29    10.07    11.22 1.00     5148
b_treatmentl                11.01      0.37    10.31    11.78 1.00     4404
b_treatmentm                 9.70      0.20     9.34    10.12 1.00     4855
b_treatmentn                 9.83      0.30     9.29    10.47 1.00     4689
b_treatmento                 9.26      0.27     8.77     9.87 1.00     4546
c_treatmenta                 0.20      0.00     0.19     0.21 1.00     4483
c_treatmentb                 0.22      0.00     0.21     0.23 1.00     5208
c_treatmentc                 0.22      0.00     0.21     0.22 1.00     5014
c_treatmentd                 0.21      0.00     0.20     0.21 1.00     5037
c_treatmente                 0.22      0.00     0.21     0.23 1.00     5280
c_treatmentf                 0.21      0.00     0.20     0.22 1.00     5530
c_treatmentg                 0.21      0.00     0.21     0.22 1.00     4831
c_treatmenth                 0.23      0.00     0.22     0.24 1.00     4411
c_treatmenti                 0.22      0.00     0.21     0.23 1.00     4508
c_treatmentj                 0.22      0.00     0.21     0.23 1.00     4305
c_treatmentk                 0.22      0.00     0.21     0.23 1.00     5305
c_treatmentl                 0.21      0.01     0.20     0.23 1.00     4514
c_treatmentm                 0.24      0.00     0.23     0.25 1.00     4512
c_treatmentn                 0.22      0.01     0.21     0.23 1.00     4308
c_treatmento                 0.21      0.01     0.20     0.22 1.00     5121
sigma_stime:treatmenta_1    24.37      4.54    15.57    33.61 1.00     4368
sigma_stime:treatmentb_1    26.53      4.58    17.56    35.79 1.00     4852
sigma_stime:treatmentc_1    25.15      4.13    17.32    33.39 1.00     4972
sigma_stime:treatmentd_1    31.90      4.64    23.16    41.34 1.00     4692
sigma_stime:treatmente_1    27.64      4.48    19.02    36.54 1.00     4259
sigma_stime:treatmentf_1    22.81      4.19    14.92    31.29 1.00     4331
sigma_stime:treatmentg_1    24.59      4.73    15.62    33.72 1.00     4791
sigma_stime:treatmenth_1    26.33      4.26    18.03    34.77 1.00     4871
sigma_stime:treatmenti_1    22.29      4.68    13.18    31.72 1.00     4229
sigma_stime:treatmentj_1    21.80      4.76    12.76    31.50 1.00     4676
sigma_stime:treatmentk_1    17.58      4.02     9.92    25.63 1.00     4442
sigma_stime:treatmentl_1    20.94      4.32    12.94    29.71 1.00     4771
sigma_stime:treatmentm_1    20.95      4.28    12.66    29.63 1.00     4672
sigma_stime:treatmentn_1    19.85      4.12    11.93    28.13 1.00     4420
sigma_stime:treatmento_1    20.28      4.17    12.14    28.51 1.00     4180
                         Tail_ESS
sigma_Intercept              3334
a_treatmenta                 3108
a_treatmentb                 3050
a_treatmentc                 2935
a_treatmentd                 3244
a_treatmente                 2994
a_treatmentf                 3248
a_treatmentg                 3107
a_treatmenth                 3030
a_treatmenti                 3101
a_treatmentj                 3375
a_treatmentk                 3320
a_treatmentl                 3471
a_treatmentm                 2489
a_treatmentn                 3519
a_treatmento                 3356
b_treatmenta                 3227
b_treatmentb                 3141
b_treatmentc                 3204
b_treatmentd                 3095
b_treatmente                 3350
b_treatmentf                 3249
b_treatmentg                 2828
b_treatmenth                 3119
b_treatmenti                 3097
b_treatmentj                 3267
b_treatmentk                 3028
b_treatmentl                 2889
b_treatmentm                 2964
b_treatmentn                 3244
b_treatmento                 2776
c_treatmenta                 2522
c_treatmentb                 2744
c_treatmentc                 2957
c_treatmentd                 3299
c_treatmente                 2835
c_treatmentf                 3075
c_treatmentg                 2929
c_treatmenth                 2984
c_treatmenti                 2475
c_treatmentj                 3025
c_treatmentk                 2814
c_treatmentl                 2936
c_treatmentm                 2713
c_treatmentn                 2800
c_treatmento                 3013
sigma_stime:treatmenta_1     2728
sigma_stime:treatmentb_1     3058
sigma_stime:treatmentc_1     3407
sigma_stime:treatmentd_1     3219
sigma_stime:treatmente_1     3389
sigma_stime:treatmentf_1     2931
sigma_stime:treatmentg_1     3111
sigma_stime:treatmenth_1     2605
sigma_stime:treatmenti_1     3036
sigma_stime:treatmentj_1     3203
sigma_stime:treatmentk_1     3102
sigma_stime:treatmentl_1     3280
sigma_stime:treatmentm_1     2958
sigma_stime:treatmentn_1     2903
sigma_stime:treatmento_1     2921

Family Specific Parameters: 
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu     3.43      0.16     3.13     3.76 1.00     5077     3273

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).