# Multivariate nonlinear mixed effect model using brms

**URL:** <https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621>\
**Category:** brms\
**Created:** [March 10, 2020, 12:35am UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621 "2020-03-10T00:35:44Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![oluwaseun](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/oluwaseun/32/7056_2.png) [@oluwaseun](https://discourse.mc-stan.org/u/oluwaseun)\
**Post date:** [March 10, 2020, 12:35am UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/1 "2020-03-10T00:35:44Z")

</div>

Hi all,

I have this nonlinear model:  
Yijk = gij(tik; Xi;Wij) + eijk;  
where gij(t; x;w ; n) = (Pj \* exp(Wij)) /(1+ exp(-(log(tik)-Uij(Xi))/(Gj))),

Uij(Xi) = MUj +B’j\*Xi.  
B’j = b1j,b2j,b3j

Wi~MVN(0,Sig) and Sig is J\*J dimensional variance-covariance matrix for the random effects. J= 3, j= 1 to J. i.e. j = DAMAG125, DAMAG151 and DAMAG292.  
eijk ~ N(0, (sigma^2)j).

My aim is to get Pj (P1,P2,P3), Sig, MUj(MU1,MU2,MU3), Gj (G1,G2,G3), (b1j,b2j,b3j), (sigma^2j).

Here is my code:

f1 = DAMAGE\_Y ~ (P_exp(w1+w2+w3)) / (1 + exp(-(log(t)-(MU+(b1_x1+b2_x2+b3_x3)))/G))  
prior\_11 = c(set\_prior(“normal(0, 1)”,nlpar = “b1”, coef = “DAMAG125”),  
set\_prior(“normal(0, 1)”,nlpar = “b1”, coef = “DAMAG151”),  
set\_prior(“normal(0, 1)”,nlpar = “b1”, coef = “DAMAG292”),  
set\_prior(“normal(0, 1)”, nlpar = “b2”, coef = “DAMAG125”),  
set\_prior(“normal(0, 1)”, nlpar = “b2”, coef = “DAMAG151”),  
set\_prior(“normal(0, 1)”, nlpar = “b2”, coef = “DAMAG292”),  
set\_prior(“normal(0, 1)”, nlpar = “b3”,coef = “DAMAG125”),  
set\_prior(“normal(0, 1)”, nlpar = “b3”,coef = “DAMAG151”),  
set\_prior(“normal(0, 1)”, nlpar = “b3”,coef = “DAMAG292”),  
set\_prior(“normal(4, 3)”, nlpar = “P”, coef = “DAMAG125”),  
set\_prior(“normal(4, 3)”, nlpar = “P”, coef = “DAMAG151”),  
set\_prior(“normal(4, 3)”, nlpar = “P”, coef = “DAMAG292”),  
set\_prior(“normal(0, 1)”, nlpar = “MU”, coef = “DAMAG292”),  
set\_prior(“normal(0, 1)”, nlpar = “MU”, coef = “DAMAG125”),  
set\_prior(“normal(0, 1)”, nlpar = “MU”, coef = “DAMAG151”),  
set\_prior(“normal(-2, 2)” , class = “sd”, nlpar = “W1”),  
set\_prior(“normal(-2, 2)” , class = “sd”, nlpar = “W2”),  
set\_prior(“normal(-2, 2)” , class = “sd”, nlpar = “W3”),  
set\_prior(“normal(0, 1)”, nlpar = “G”,coef = “DAMAG125”),  
set\_prior(“normal(0, 1)”, nlpar = “G”,coef = “DAMAG151”),  
set\_prior(“normal(0, 1)”, nlpar = “G”,coef = “DAMAG292”))

form = bf(f1,nl = TRUE)+list(phi1~0+DAMAG,phi3~0+DAMAG,MUJ~0+DAMAG, w1~1+(1|DAMAG),w2~1+(1|DAMAG),w3~1+(1|DAMAG), b1~0+DAMAG,b2~0+DAMAG,b3~0+DAMAG)  
n1\_b \<- brm(form,data = rdat,prior = prior\_11,chains=1,iter=1000)

Here is my result:

```
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS

```

sd(w1\_Intercept) 0.86 0.71 0.10 2.68 1.24 8 29  
sd(w2\_Intercept) 0.62 0.48 0.05 1.69 1.06 20 60  
sd(w3\_Intercept) 0.62 0.49 0.11 1.79 1.15 17 27

Population-Level Effects:  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS  
phi1\_DAMAG125 -1.85 1.17 -4.43 -0.21 1.09 15  
phi1\_DAMAG151 -1.07 0.74 -2.85 -0.09 1.14 8  
phi1\_DAMAG292 -1.37 0.98 -3.79 -0.21 1.13 7  
phi3\_DAMAG125 0.66 0.04 0.58 0.71 1.17 5  
phi3\_DAMAG151 -1.83 0.13 -2.07 -1.59 1.32 3  
phi3\_DAMAG292 0.55 0.04 0.46 0.63 1.16 9  
MUJ\_DAMAG125 4.50 0.09 4.33 4.68 1.06 11  
MUJ\_DAMAG151 13.29 0.33 12.75 13.86 1.53 2  
MUJ\_DAMAG292 4.15 0.09 3.97 4.32 1.00 13  
w1\_Intercept -14.19 13.41 -33.85 11.11 1.01 8  
w2\_Intercept 36.57 13.36 12.59 59.02 1.14 8  
w3\_Intercept -22.86 7.69 -32.59 -7.11 1.07 9  
b1\_DAMAG125 -0.49 0.03 -0.54 -0.45 1.51 2  
b1\_DAMAG151 4.96 0.18 4.56 5.25 1.11 8  
b1\_DAMAG292 -0.47 0.03 -0.52 -0.41 1.37 3  
b2\_DAMAG125 -0.97 0.11 -1.14 -0.76 1.05 15  
b2\_DAMAG151 -1.60 0.86 -3.01 0.59 1.24 3  
b2\_DAMAG292 -1.32 0.16 -1.63 -0.99 1.01 20  
b3\_DAMAG125 0.32 0.02 0.29 0.36 1.00 16  
b3\_DAMAG151 0.19 0.28 -0.34 0.70 1.05 14  
b3\_DAMAG292 0.33 0.03 0.27 0.38 1.32 3

Family Specific Parameters:  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS Tail\_ESS  
sigma 0.13 0.00 0.12 0.13 1.06 21 90

How can I split the sigma? Instead of having sigma = 0.13,  
I want to have sigma for DAMAG125, DAMAGE151 and DAMAGE252.

Can anyone help me out?.

---

<div class="post-metadata">

**Author:** ![Guido\_Biele](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/guido_biele/32/3414_2.png) [@Guido\_Biele](https://discourse.mc-stan.org/u/Guido_Biele)\
**Post date:** [March 10, 2020, 7:50am UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/2 "2020-03-10T07:50:58Z")

</div>

If I understand correctly, you want the residual variance to vary by DAMAG125, DAMAGE151 and DAMAGE252.

See the section **More Complex Multivariate Models** [here](https://cran.r-project.org/web/packages/brms/vignettes/brms_multivariate.html) for how to estimate varying residual variance.

Also not that you can use latex here (put the equation between $$), which makes equation more accessible. (`$Y_{ijk} = g_{ij}(t_{ik}, X_{i},W_{ij}) + e_{ijk}$` makes Y\_{ijk} = g\_{ij}(t\_{ik}, X\_{i},W\_{ij}) + e\_{ijk})

---

<div class="post-metadata">

**Author:** ![oluwaseun](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/oluwaseun/32/7056_2.png) [@oluwaseun](https://discourse.mc-stan.org/u/oluwaseun)\
**Post date:** [March 10, 2020, 12:43pm UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/3 "2020-03-10T12:43:26Z")

</div>

Thank you so much for your reply. I checked More Complex Multivariate Models you suggested. They have two response variables (tarsus and back) but in my own case, I have only one response variable and also a random effect (Damage number). The random effect has 3 levels (125, 151, 292). I want to get the residual variance for each of the random effects.  
How can I do that please?

---

<div class="post-metadata">

**Author:** ![Guido\_Biele](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/guido_biele/32/3414_2.png) [@Guido\_Biele](https://discourse.mc-stan.org/u/Guido_Biele)\
**Post date:** [March 10, 2020, 2:31pm UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/4 "2020-03-10T14:31:45Z")

</div>

I think you can still use the same formula approach. Did you.try that?

---

<div class="post-metadata">

**Author:** ![oluwaseun](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/oluwaseun/32/7056_2.png) [@oluwaseun](https://discourse.mc-stan.org/u/oluwaseun)\
**Post date:** [March 12, 2020, 11:49am UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/5 "2020-03-12T11:49:26Z")

</div>

How will I do that with one response variable? Please guide me.

---

<div class="post-metadata">

**Author:** ![Guido\_Biele](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/guido_biele/32/3414_2.png) [@Guido\_Biele](https://discourse.mc-stan.org/u/Guido_Biele)\
**Post date:** [March 13, 2020, 1:33pm UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/6 "2020-03-13T13:33:38Z")

</div>

The first example in the brms vignette [Estimating Distributional Models with brms](https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html) is close to what (I think) you want to do.

You have to specify two brms formulas, one for you outcome and one for the residual variance, for example:

```
outcome_model = bf(y ~ x)
resvar_model = bf(sigma ~ group)
fit = brms(outcome_model + resvar_model, ...)

```

Here, the residual variance would vary for individuals with different values in the variable `group`.

`resvar_model = bf(sigma ~ group)` is specific to a gaussian model, and would need to change for a different `family`. The models for outcome and residual terms can be more complex or hierarchical.

---

<div class="post-metadata">

**Author:** ![oluwaseun](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/oluwaseun/32/7056_2.png) [@oluwaseun](https://discourse.mc-stan.org/u/oluwaseun)\
**Post date:** [April 8, 2020, 1:55pm UTC](https://discourse.mc-stan.org/t/multivariate-nonlinear-mixed-effect-model-using-brms/13621/7 "2020-04-08T13:55:06Z")

</div>

Thank you. I am still not getting the result.