User defined non linear model in brms

Interfaces brms
1

Dear All,
I am pretty new with brms and I was wondering if it is possible to specify a user defined function for parameters estimate of a non linear model.

In my case the dependent variables are a summation of other variables. A reproducible example could be something below:

### create data
set.seed(123)
class <-c(rep(3,10),rep(5,10),rep(6,10),rep(2,10),rep(8,10),rep(1,10),rep(7,10))
x=runif(length(class), 1,10)

### create the true response
allclasses=(unique(class))
dependent=rep(NA,length(allclasses))
for(k in 1:length(allclasses)){
  deptemp=x[which(class==allclasses[k])]
  sumvalt=sum(deptemp,na.rm=TRUE)
  dependent[k]=sumvalt
}
response=2+5*dependent

### define the model
model<-function(vect1,class,param){
  
  allclasses=(unique(class))
  dependent=rep(NA,length(allclasses))
  for(k in 1:length(allclasses)){
    deptemp=vect1[which(class==allclasses[k])]
    sumvalt=sum(deptemp,na.rm=TRUE)
    dependent[k]=sumvalt
  }
  output=param[1]+param[2]*dependent
  return(output)
}
2

Yes, this is possible via the stanvars argument of brm.

3

Dear Paul thanks a lot.
As far I understood stanvar seems to modify the prior limits. In my case the response variable has different length of the predictor and I need to apply a function to the predictors (specifically a sum by classes) to connect them to the response variable.

Any hints?

Thanks again and best regards
giuseppe

4

You may also use stanvar to pass self-defined Stan functions via stanvar(..., block = "functions"). So if you can define your function in Stan language, you can pass it via stanvar and then use it in the predictor term.

5

Hi Paul,

looking at some previous post and trying with my previous code in this post I was doing something like this below. I have problem to define the data in brm because the response has different length of the predictors.

### create data
set.seed(123)
class <-c(rep(3,10),rep(5,10),rep(6,10),rep(2,10),rep(8,10),rep(1,10),rep(7,10))
x=runif(length(class), 1,10)


### create the true response
allclasses=sort(unique(class))
dependent=rep(NA,length(allclasses))
for(k in 1:length(allclasses)){
  deptemp=x[which(class==allclasses[k])]
  sumvalt=sum(deptemp,na.rm=TRUE)
  dependent[k]=sumvalt
}
response=2+5*dependent



# Prior distribution
prior <- 
  prior(normal(0.5, 0.01), class = "b", nlpar = "theta1", lb = 0) +
  prior(normal(0.5, 0.1), class = "b",  nlpar = "theta2", lb = 0) 
 



# Custom code for the model
nlfun <- stanvar(
  scode = " vector nlfun(vector allValues, vector allClasses,vector uniqueClasses,theta1,theta2) {
  int numClass = rows(uniqueClass);
  
  vector[numClass] output;
  for (k in 1:numClass) {
    for(c in 1:length(allClasses)){
      if(allClasses[c]==uniqueClasses[k]){
        output[k] += theta1+ theta2* allValues[c];
      }
    }
  }
  return mu;
}",
    block = "functions"
)

# BRMS formula for the model
formula <- 
  bf(response ~ nlfun(x, classes, sort(unique(classes)), theta1, theta2),
     theta1 + theta1 ~ 1,
     nl = TRUE,
     loop = FALSE)


# Running brms on the model
mod <- 
    brm(formula = formula,
        prior   = prior,
        data=
        warmup  = 500,
        iter    = 1500,
        thin    = 1)
```
6

Hmm in that case (response and predictor having different length), I would think it’s actually much easier to just use Stan directly instead of brms.