Multivariate nonlinear mixed effect model using brms

Interfaces brms
1

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 ~ (Pexp(w1+w2+w3)) / (1 + exp(-(log(t)-(MU+(b1x1+b2x2+b3x3)))/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?.

2

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

See the section More Complex Multivariate Models here 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})

3

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?

4

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

5

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

6

The first example in the brms vignette Estimating Distributional Models with brms 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.

1 Like
7

Thank you. I am still not getting the result.