Dynamically created model formula

Interfaces brms
1

Short summary of the problem
trying to model through a for loop and change model formula dynamically, stuck at the model formula throwing error, wonder if you have input on how to create the formula allowing for a change dynamically

code_to_run_your_model(if_applicable)
allnames<-c(“BRADY”,“ADVANCEDAGE”,“SEPSIS”,“COPD”,“CKD”,“LIVER”)

library(gtools)
library(brms)
permutated<-data.frame(permutations(length(allnames), length(allnames), v=allnames, set=F, repeats.allowed=T))
permutated<-permutated[!duplicated(permutated[c(1:length(allnames))]),]
permutated$fineq<-apply(permutated,1,function(x) gsub("._", " “, paste0(”+",tolower(sort(unique(x))),“int”,"",sort(unique(x)),collapse = “”)))
permutated$fingroups<-apply(permutated[,1:length(allnames)],1,function(x) gsub("._", " “, paste0(”,",tolower(sort(unique(x))),“int”,"~1+ (1|ID/BOV)",collapse = “”)))
permutated$finpri<-apply(permutated[,1:length(allnames)],1,function(x) gsub("._", " “, paste0(”,",“prior(normal(0.001, 100), nlpar = '”,tolower(sort(unique(x))),“int’,lb=0.0001)”,collapse = “”)))
permutated<-permutated[!duplicated(permutated$fineq),]
permutatedfin<-permutated[,(length(allnames)+1):(length(allnames)+3),drop=F]

i<-2
mod_formula<-bf(print(as.formula(paste0(“QT ~ (((RR^alpha)QT0+genSEX)+Acos(3.1416/12(CYRCADIAN-pha)))”,permutatedfin[i,][,1]))),
for (j in 1:length(strsplit(paste0(",alpha ~ 1+ (1|ID/BOV),gen ~ 1,QT0 ~ 1+ (1|ID/BOV),A ~ 1+ (1|ID/BOV),pha~1+ (1|ID/BOV)",permutatedfin[i,][,2]),",")[1][[1]][-1])){
print(as.formula(strsplit(paste0(",alpha ~ 1+ (1|ID/BOV),gen ~ 1,QT0 ~ 1+ (1|ID/BOV),A ~ 1+ (1|ID/BOV),pha~1+ (1|ID/BOV)",permutatedfin[i,][,2]),",")[1][[1]][-1][j]))
},nl = TRUE )

If possible, add also code to simulate data or attach a (subset of) the dataset you work with.

Please also provide the following information in addition to your question:

  • Operating System: Windows
  • brms Version: 2.16.1

Don’t forget to add relevant tags to your topic (top right of this form) especially for application area. Delete this text before posting your question :-) Thx!

Looking forward to your topic!

Thanks,

Andras

2

just in case the code does not run when pasted here it is again:

allnames<-c(“BRADY”,“ADVANCEDAGE”,“SEPSIS”,“COPD”,“CKD”,“LIVER”)

library(gtools)
permutated<-data.frame(permutations(length(allnames), length(allnames), v=allnames, set=F, repeats.allowed=T))
permutated<-permutated[!duplicated(permutated[c(1:length(allnames))]),]
permutated$fineq<-apply(permutated,1,function(x) gsub(“._“, " “, paste0(”+”,tolower(sort(unique(x))),“int”,"”,sort(unique(x)),collapse = “”)))
permutated$fingroups<-apply(permutated[,1:length(allnames)],1,function(x) gsub(“._“, " “, paste0(”,”,tolower(sort(unique(x))),“int”,“~1+ (1|ID/BOV)”,collapse = “”)))
permutated$finpri<-apply(permutated[,1:length(allnames)],1,function(x) gsub("._”, " “, paste0(”,",“prior(normal(0.001, 100), nlpar = '”,tolower(sort(unique(x))),“int’,lb=0.0001)”,collapse = “”)))
permutated<-permutated[!duplicated(permutated$fineq),]
permutatedfin<-permutated[,(length(allnames)+1):(length(allnames)+3),drop=F]

i<-2
library(brms)
mod_formula<-bf(print(as.formula(paste0(“QT ~ (((RR^alpha)QT0+genSEX)+Acos(3.1416/12(CYRCADIAN-pha)))”,permutatedfin[i,][,1]))),
for (j in 1:length(strsplit(paste0(“,alpha ~ 1+ (1|ID/BOV),gen ~ 1,QT0 ~ 1+ (1|ID/BOV),A ~ 1+ (1|ID/BOV),pha~1+ (1|ID/BOV)”,permutatedfin[i,][,2]),“,”)[1][[1]][-1])){
print(as.formula(strsplit(paste0(“,alpha ~ 1+ (1|ID/BOV),gen ~ 1,QT0 ~ 1+ (1|ID/BOV),A ~ 1+ (1|ID/BOV),pha~1+ (1|ID/BOV)”,permutatedfin[i,][,2]),“,”)[1][[1]][-1][j]))
},nl = TRUE )

3

Did you find a solution to this? I’m wondering also about how to create dynamic formulas.

4

I’m not sure I understand the specifics of what you’re trying to do. Can you provide a simpler example that captures the essence of your question and maybe provide a bit more explanation of what you expect the output to look like? For now, here’s an example of building several formulas and priors programmatically:

library(tidyverse)
library(brms)

allnames = names(mtcars)[1:5]

# Get all formulas with mpg as the outcome and three predictor variables
formulas = apply(combn(allnames[-1], 3), 2, function(x) as.formula(paste(allnames[1], " ~", paste(x, collapse=" + "))))

# Set priors separately for each predictor
sd = c(1, 0.5, 0.8, 0.2) %>% set_names(allnames[-1])
priors = map(formulas, ~{
  vars = all.vars(.x)[-1]
  map_df(vars, ~prior_string(prior=paste0("normal(0,", sd[match(.x, names(sd))], ")"), 
                             coef=names(sd)[match(.x, names(sd))]))
  })

formulas  
#> [[1]]
#> mpg ~ cyl + disp + hp
#> <environment: 0x7fa14536fba0>
#> 
#> [[2]]
#> mpg ~ cyl + disp + drat
#> <environment: 0x7fa145376e48>
#> 
#> [[3]]
#> mpg ~ cyl + hp + drat
#> <environment: 0x7fa16390d048>
#> 
#> [[4]]
#> mpg ~ disp + hp + drat
#> <environment: 0x7fa163910ef8>

priors
#> [[1]]
#>          prior class coef group resp dpar nlpar   lb   ub source
#>    normal(0,1)     b  cyl                       <NA> <NA>   user
#>  normal(0,0.5)     b disp                       <NA> <NA>   user
#>  normal(0,0.8)     b   hp                       <NA> <NA>   user
#> 
#> [[2]]
#>          prior class coef group resp dpar nlpar   lb   ub source
#>    normal(0,1)     b  cyl                       <NA> <NA>   user
#>  normal(0,0.5)     b disp                       <NA> <NA>   user
#>  normal(0,0.2)     b drat                       <NA> <NA>   user
#> 
#> [[3]]
#>          prior class coef group resp dpar nlpar   lb   ub source
#>    normal(0,1)     b  cyl                       <NA> <NA>   user
#>  normal(0,0.8)     b   hp                       <NA> <NA>   user
#>  normal(0,0.2)     b drat                       <NA> <NA>   user
#> 
#> [[4]]
#>          prior class coef group resp dpar nlpar   lb   ub source
#>  normal(0,0.5)     b disp                       <NA> <NA>   user
#>  normal(0,0.8)     b   hp                       <NA> <NA>   user
#>  normal(0,0.2)     b drat                       <NA> <NA>   user

Then you could fit a model for each formula:

models = map(seq_along(formulas), 
             ~brm(prior=priors[[.x]], formula=formulas[[.x]], data=mtcars, 
                  family=gaussian(), backend="cmdstanr"))
5

Thanks Joel, your solution will get me the results I was looking for… Andras