Brms vs. nlme

Interfaces brms
1

I have a 3 non linear models below 1 fit using nlme, 1 fit using brms that is not multi-level and one that is fit using brms that is multilevel. The goodness of fit, S statistic, is 180, 231 and 245 respectively so the nlme model has the best fit. You can see this if you look at the plot.

I think that is because the prior on the “a” coefficient is prior(normal(0,2),nlpar=“a”) where the variance is fixed. I’d like to add a hyperparameter: prior(normal(0,sigma_group),nlpar=“a”) so that each group could have different variance for the a coefficient. How can that be done in brms? In Stan it is typically real<lower=0> sigma_group; declared as a parameter and then the prior is declared normal(0,sigma_group);

Also any other suggestions to get the brms models to get near the nlme fit?

library(nlme)
library(brms)
library(ggplot2)
set.seed(1)
x = 1:20
y = x^2+rnorm(20,0,50)
x2 = 1:20
y2 = x2^2+rnorm(20,0,80)
x3 = 1:20
y3 = x3^2+rnorm(20,10,200)
x4 = 1:20
y4 = x4^2+rnorm(20,20,300)
dat1 = data.frame(x =c(x,x2,x3,x4), y=c(y,y2,y3,y4), group = c(rep("A",20),rep("B",20),rep("C",20),rep("D",20) )  )
f1<-formula( y ~ a*x^b | group)
#### using nlme
m <- nlsList(f1,data=dat1,start=list(a=1,b=1)) 
m
summary(m)
dat1$fitted_nlme = predict(m)
#### using brms multilevel model
fit1 =brm(   bf( y~ a*x^b,   
                 nonlinear =  list(
                   a ~  1 + (1|group),
                   b ~  1 + (1|group)),
                 nl= TRUE),
             data = dat1,
             prior= c(
               prior(normal(0,2),nlpar="a"),
               prior(normal(0,2),nlpar="b")
             ) ,control = list(adapt_delta = .99) )
summary(fit1)
dat1 = cbind(dat1, fitted(fit1)[,-2])

fixed.effects(fit1)[1,1]
random.effects(fit1)$group[,,1]

fixed.effects(fit1)[2,1]
random.effects(fit1)$group[,,2]

a_pop=    fixed.effects(fit1)[1,1]
b_pop=    fixed.effects(fit1)[2,1]

#calculate goodnes of fit to compare to nlme model
for(i in 1:nrow(dat1)){
  #i = 1
  if(dat1$group[i]=="A"){
    a_group = random.effects(fit1)$group[,,1][1]
    b_group = random.effects(fit1)$group[,,2][1]
  }else if(dat1$group[i]=="B"){
    a_group = random.effects(fit1)$group[,,1][2]
    b_group = random.effects(fit1)$group[,,2][2]
  }else if(dat1$group[i]=="C"){
    a_group = random.effects(fit1)$group[,,1][3]
    b_group = random.effects(fit1)$group[,,2][3]
  }else{
    a_group = random.effects(fit1)$group[,,1][4]
    b_group = random.effects(fit1)$group[,,2][4]
  }
  dat1$fitted_brms2[i] = (a_pop+a_group)*dat1$x[i]^(b_pop+b_group) 
}

#### using brms not multilevel
fit2  =brm(   bf( y~ a*x^b,   
             nonlinear =  list(
               a ~  0 + group,
               b ~  0 + group), 
             nl= TRUE),
             data = dat1,
             prior= c(
               prior(normal(0,2),nlpar="a"),
               prior(normal(0,2),nlpar="b")
             ) ,control = list(adapt_delta = .99) 
             )
summary(fit2)
coef =  posterior_samples(fit2, pars = c("a_groupA","b_groupA",
                                 "a_groupB","b_groupB",
                                 "a_groupC","b_groupC",
                                 "a_groupD","b_groupD"
                                 ))  
mean_brms_coef = rbind(apply(coef,2,mean) )
#calculate goodnes of fit to compare to nlme model
for(i in 1:nrow(dat1)){
  #i = 1
  if(dat1$group[i]=="A"){
    index=1
  }else if(dat1$group[i]=="B"){
    index = 2
  }else if(dat1$group[i]=="C"){
    index=3
    }else(index =4)
  dat1$fitted_brms[i] = mean_brms_coef[index]*dat1$x[i]^mean_brms_coef[index+1]
}

## compare models
sqrt( sum((dat1$y -dat1$fitted_nlme   )^2)/72)
sqrt( sum((dat1$y -dat1$fitted_brms  )^2)/72)
sqrt( sum((dat1$y -dat1$fitted_brms2   )^2)/72)


ggplot(dat1, aes(x = x, y = y, color = group))+geom_point()+
  geom_line(data = dat1, aes(x = x, y = fitted_nlme, color ="nlme"),inherit.aes=FALSE  )+
  geom_line(data = dat1, aes(x = x, y = Estimate, color ="brms_multilevel_fitted"),inherit.aes=FALSE  )+
  geom_line(data = dat1, aes(x = x, y = fitted_brms, color ="brms_not_multilevel"),inherit.aes=FALSE  )+
  facet_wrap(~group)
1 Like
2

Try out

bf(
  y ~ a*x^b,
  a ~ 0 + group,
  b ~ 0 + group
)

This doesn’t use mutlilevel structure, which does not make much sense if you only have two groups.

3

Thank you for the response. In my real data I have more than 10 groups and I’d like to use a multilevel structure. If I did that would there be a population level a and population level b coefficients-similar to a linear multilevel model in STAN? Can you advise how to make it a multilevel mode and do have any suggestion on improving the fit?

Note: I updated the problem to include more groups (4). I also added a comparison of fit to the nlme model. The nlme model has a much lower residual standard error.

4

To fit a multilevel model, go for

bf(
  y ~ a*x^b,
  a ~ 1 + (1|group),
  b ~ 1 + (1|group)
)

After model fitting, you can extract the (posterior distributions of) the coefficients for each group via the coef() method,

1 Like
5

Please see edit which compares 3 models and the nlme is winning.

6

I wouldn’t call it “winning” since you are just evaluating in-sample fit.

Also, the computation of the fitted values for the brms model is not quite correct.
Summary statistics should only be computed after all other transformations have been done,
not the other way round. To extract fitted values for brms model, please use the fitted() method.

dat1$fitted_brms_new <- fitted(fit1)[, 1]
dat1$fitted_brms_new2 <- fitted(fit2)[, 1]
sqrt( sum((dat1$y -dat1$fitted_brms_new   )^2)/72)
sqrt( sum((dat1$y -dat1$fitted_brms_new2   )^2)/72)

looks much closer to what we find with nlme.

7

Thank you for the correction and the brms package. Is there a way to add hyperparameters for the variance of the “a” parameter similar to declaring real<lower=0> sigma_group and then having a prior for a be prior(normal(0,sigma_group),nlpar=“a”)?

8

See ?set_prior

9

Hi @paul.buerkner
I don’t have access to R so I am primarily using Pystan. Is Brms available from python? I am really interested in trying it out. Thanks.

10

Nope, you have to use R.

1 Like
11

Hi @Amit_Kushwaha and @paul.buerkner

You can use the library bambi in python. It was developed, partially inspired by brms and lme4, to provide a formulaic interface to build Bayesian mixed effects models.

You can access both pystan and pymc3 as the back-end engines. I find it incredibly helpful to have a counterpart to brms in the python ecosystem. This way I test my models in Rand Python. Here is a link

bambi

12

Paul:

But couldn’t a Python user make use of brms in an indirect way?

  1. Set up R on your machine.
  2. Write a model in brms, which then produces a Stan program.
  3. Save the Stan program as a file and then run it using PyStan or CmdStanPy.

This looks clunky but it’s not so bad, as you’re only using R as a device to create your Stan program. You can do all your data processing in Python. And, once you have your Stan program and you’re running it from Python, you never have to look back: just alter the Stan program as needed.

13

Or use the brms wrapper


As a python user R syntax doesn’t always make sense, so this will enable user to skip writing any R.

2 Likes
14

Andrew: This is certainly an option if you only use brms to generate Stan code. However, a lot of brms’ merits are its flexible post-processing which you loose as soon as you port only the Stan code over to Python.