Nonlinear model with grouped data and a single control, model specification problem in BRMS

Modeling
, ,
13

Yes deleting out the control group works fine and there are vignettes out there on constructing these very simple data structures and running a simple single nonlinear relationship in BRMS.

No, all the data set of a DOE should go into the regression analysis, but it does have to be coded correctly and have the correct attribute columns in the dataset to specify the grouping. Besides, why shouldn’t it work?

Here is an interesting discussion on a very similar problem:

Declaration Nonlinear mixed effects model in BRMS

Funnily enough there is a recommendation that a fixed effect can just get inserted into the nonlinear formula. Works as you showed in your equations and my initial attempt. Quite valid if this is just a constant getting added dependent on some grouping factor….and yes this does work for this synthetic dataset, but I can definitely see the advantages of adding this in its own formula (like what if the rate also is affected, dumping this into one huge formula would make it unreadable). I just cant seem to get BRMS to accept the coding. I think its something to do with the specification of the priors.

In the above reference mmihaylova added a fixed effect “Cat” as a test out…but there must have been an addition to the priors section for “Cat” that is not mentioned.

Here is the brm coding so far:

I’ve simplified the formula (same synthetic data though):

# taking linear terms out of nonlinear formula

# nonlinear, continuous scale for response
# Initial kicks
init_func <- function(chain_id=1) {
  list ( b_Asym  =  array(14) ,
         b_A  =  array(6) ,
         b_rate  =  array(0.2),
         #b_Group = array(0.5),
         b_Control = array(2)
         )
}

# setting up chains
init_list <- list(
  init_func(chain_id = 1),
  init_func(chain_id = 2),
  init_func(chain_id = 3),
  init_func(chain_id = 4)
)

# NUTS configuration
control <- list(
  adapt_engaged = TRUE,
  adapt_delta = 0.9, #increased from default of 0.8
  stepsize = 0.08, # 0.05 default
  max_treedepth = 15 # was 15
)

prior1 <- prior(normal(14, 1), nlpar = "Asym") + 
  prior(normal(6, 1), nlpar = "A")+         
  prior(normal(0.2, 0.1), nlpar = "rate")+
  prior(normal(0.4,0.03 ), class="b", group="Control", nlpar="Asym")
  

# prior(normal(0.5, 1), nlpar = "Group"  )+

fit1 <- brm(bf(Response ~ Asym - ((A)*exp(-rate * Treatment)), 
               Asym ~ 1 + Control, # + Group,
               A ~ 1,
               rate ~ 1,
               nl = TRUE),
            data = SynthData, 
            prior = prior1,
            #family = cumulative(link="logit", threshold="flexible"),
            iter = 4000 , warmup = 2000, thin = 2,
            chains = 4, cores=6,
            init = init_list, 
            control = control,
            #file = "./model/fit14"
)

summary(fit1)

# fit process visualization
ggs(fit1, burnin = TRUE) %>%
  filter(Parameter == "b_KG5_Intercept") %>% 
  mutate(chain = factor(Chain),
         intercept = value) %>% 
  
  ggplot(aes(x = Iteration, y = intercept, color = chain)) +
  geom_line()+
  ggtitle("default settings")+
  theme_bw()
```

I’ve commented out “Group” out of the model and priors & initial values. Its only a small effect in the synthetic data. I’m making some guesses here. The grouping has to be declared in the prior, so that presumably “group=”Control”, and I am assuming this has to be tied to a parameter in the nonlinear formula declaration and this is the nlpar=”Asym”.

image
image675×108 11 KB

There is definitely something wrong in the prior:

image

This ref too:

Specifying a nonlinear mixed effects model in BRMS

looking at adding random effects to a nonlinear model, no fixed effect though.