Declaration of Nonlinear Mixed-Effects Model in brms

Interfaces brms
1
  • Operating System: Windows 10
  • brms Version: 2.8

I am fitting a nonlinear mixed-effects model to longitudinal data.
Here are some useful details about the dataset:
Data has been obtained from N patients before and after an operation.
For each patient, M measurements of the dependent variable have been obtained both before and after the operation.
I am modelling my data using a slightly modified Richard’s curve, the formula given as:

LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta))

where LDL is my dependent variable, Time is an independent variable, and alpha, gamma, delta - parameters to be estimated.
What I would like to do is include a variable Category - with “Before” and “After” values, which would denote whether a measurement from a patient has been taken before or after the operation. I think Category should be a fixed effect, as I believe that there is some consistency in the measurements that have been taken before/after the operation. I would like to declare Category as be a fixed effect, however, I am not quite sure how to declare it in my model. Am I supposed to add it somewhere in the formula given above, or how exactly would it work?
I did go through some examples of declaration of mixed effects models, but nowhere did I find something which answers my question - the examples were either for linear models or all of the fixed effects were included in the formula for the model, which is not my case.
I provide some of my code below:

fit_model ← function(Data){

fit<- brm(
bf(LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta)),
alpha ~ 1 + (1|Sample) ,
gamma ~ 1 + (1|Sample),
delta ~ 1,
nl = TRUE),
prior = c(
prior(normal(1500, 30), class=“b”,nlpar = “alpha”),
prior(normal(2.5,0.07), class=“b”, nlpar = “gamma”),
prior(normal(0.2,0.05), class=“b”, lb=0, ub=0.5, nlpar = “delta”),
prior(student_t(3, 0, 10), class=“sigma”),
prior(student_t(3, 0, 10), class=“sd”, nlpar=“alpha”),
prior(student_t(3, 0, 10), class=“sd”, nlpar=“gamma”),
prior(normal(0.4,0.03 ), class=“sd”, group=“Sample”, coef=“Intercept”, nlpar=“gamma”),
prior(normal(300,10), class=“sd”, group=“Sample”, coef=“Intercept”, nlpar=“alpha”)),
data = Data, family = gaussian(),

I would really appreciate it if someone could give me some insight.

2

I’m not an expert on modeling or brms but it sort of seems like it’s up to you where in the model Category belongs. If you thought it had a simple additive effect on top of the non-linear model you could write

LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta)) + Category

(not sure if brms would make you convert Category to a binary 0/1 indicator variable or not).

Having its effect be additive is probably not a bad place to start your investigation if you have no other conception of where in the data generating process it belongs.

1 Like
3

In additon to what Sean said, you may include them in the formulas for the non-linear parameters if you like, for instance,

gamma ~ Category + (1 | Sample)
1 Like
4

Thank you both for your replies!

I tried both recommendations.

First, I tried only adding Category to the formula:

LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta)) + Category

But the output did not in any way reflect that:

Group-Level Effects:
~Sample (Number of levels: 56)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(alpha_Intercept) 347.44 11.15 325.73 369.06 5991 1.00
sd(gamma_Intercept) 0.38 0.03 0.32 0.43 5800 1.00

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
alpha_Intercept 1692.08 30.76 1632.99 1752.57 8715 1.00
gamma_Intercept 2.53 0.05 2.44 2.62 4954 1.00
delta_Intercept 0.22 0.02 0.19 0.26 12934 1.00

Then I tried the other option:

bf(LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta)),
alpha ~ 1 + Cat + (1|Sample) ,
gamma ~ 1 + Cat + (1|Sample),
delta ~ 1,
nl = TRUE),

This time it got reflected in the output:

Group-Level Effects:
~Sample (Number of levels: 56)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(alpha_Intercept) 337.55 9.18 319.82 355.63 16107 1.00
sd(gamma_Intercept) 0.61 0.02 0.57 0.65 6288 1.00

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
alpha_Intercept 1476.75 28.12 1421.97 1532.18 14621 1.00
alpha_Cat 1444.90 30.06 1386.09 1503.85 26930 1.00
gamma_Intercept 2.12 0.06 2.00 2.23 2019 1.00
gamma_Cat 1.97 0.07 1.84 2.10 4719 1.00
delta_Intercept 0.24 0.00 0.23 0.25 18006 1.00

Now, considering that my Category values are 0/1, does this output mean that the values of alpha in Category 0 are, on average, greater than the alpha-values in Category 1, with 1444.90?

5

Yes, this seems correct to me.

6

Hi,
I have a follow up question which I can’t seem to be able to answer myself:

Since my categorical variable has only two values - B and F,
I decided to try out and use it once as a fixed effect:

Model 1:

bf(LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta)),
alpha ~ 1 + Cat + (1|Sample) ,
gamma ~ 1 + Cat + (1|Sample),
delta ~ 1,
nl = TRUE),

, and once as a random effect:

fit<- brm(
bf(LDL ~ 15 + (alpha-15) / (1 + exp ((gamma-Time) / delta)),
alpha ~ 1 + (1|Cat) + (1|Sample) ,
gamma ~ 1 + (1|Cat) + (1|Sample),
delta ~ 1,
nl = TRUE),

Now, the values I get for the estimates for alpha and gamma for each individual sample (summing up the fixed effects for alpha and the random effects for alpha, and the same for gamma), I seem to be getting quite different results, throughout the samples.

For example, if I calculate the estimates for alpha, gamma, and delta, generated by the second model, and I plug these into the formula:

LDL ~ 15 + (alpha - 15) / (1 + exp ((gamma - Time) / delta))

I get curves which are much closer to the original data points, which is not the case for Model1, where the curves for the “Before” samples, are way off (just estimated by looking at the plot, and not in any more formal way).

So I was curious why this is happening?
Am I expected to always be getting estimates for my parameters which are quite similar for both models, are these models supposed to, effectively, be the same ?

Or is it possible for these to be behaving differently?

Now, it is worth mentioning that Model1 has 81 data points w
ith pareto-k > 0.7,and Model2 - 100, so I assume this could be the reason.

But, in general, if both models were okay (that is, no influential points), are these models supposed to behave in the exact same way - intuitively, the answer seems to be yes, but is it really?

7

For the “fixed effects” formulation, brms will use dummy coding by default that is one level is the reference category (and hence the intercept represents this category) and the other implies a contrast between the two categories. This works the same way as in simple linear regression. In the “random effects” formulation, you will predict the mean of the two categories and then model the effects of the two categories as deviations from the mean. Since you only have two groups and no really informative priors, I would recommend the “fixed effects” version in which you can of course change the contrast from dummy coding to something else (see ?contr.treatment for documation).