# Checking my setup on multivariate nonlinear model using brms

Hi all,

I have this multivariate nonlinear model that I want to set up on brms. Can you please check my setup on brms and correct me? I think I am not getting the setup.

Here is my data structure:

Let n = 48 test units, i = 1,…,n, J = 3 wavelenghts(DAMAGE_NUM1250, DAMAGE_NUM1510, DAMAGE_NUM2925), j = 1,…,J,
p=3 covariates(Uv, Rh,Temp), x_i is covariate for units i, x_i = (x_i1,…x_ip), k = 1,…, m_i, where m_i is the total number of measured time points for units i. We also have time t_ik (at unit i and k).

My model is this:

y_ijk = g_ij(t_ik, x_i, w_ij) + e_ijk

e_ijk ~ N(0, r_ j)

w_i ~ MVN(0,R), where R is J by J (3 by 3) dimensional variance-covariance matrix for the random effects on J.

w_i = (w_i1,....w_ij)',  w_i is a vector of the random effecs associated with unit i for J degradation.

g_ij(t_ik, x_i, w_ij)  = (k_ j*exp⁡(w_ij))/(1+exp⁡(-(log⁡(t_ik - u_ij (**x**_i))/γ_ j )

u_ij (**x**_i) = u_ j + β'_ j * **x**_i


This is my setup on brms:

f1 = y ~ (k*exp(w)) / (1 + exp(-(log(t)-(u +(β1*x1+β2*x2+β3*x3)))/ γ))

prior_11 = c(set_prior("normal(0, 1)",nlpar = "β1", coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 1)",nlpar = "β1", coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 1)",nlpar = "β1", coef = "DAMAGE_NUM2925"),
set_prior("normal(0, 1)", nlpar = "β2", coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 1)", nlpar = "β2", coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 1)", nlpar = "β2", coef = "DAMAGE_NUM2925"),
set_prior("normal(0, 1)", nlpar = "β3",coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 1)", nlpar = "β3",coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 1)", nlpar = "β3",coef = "DAMAGE_NUM2925"),
set_prior("normal(0, 1)", nlpar = "k", coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 1)", nlpar = "k", coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 1)", nlpar = "k", coef = "DAMAGE_NUM2925"),
set_prior("normal(0, 1)", nlpar = "u", coef = "DAMAGE_NUM2925"),
set_prior("normal(0, 1)", nlpar = "u", coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 1)", nlpar = "u", coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 1)", nlpar = "γ",coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 1)", nlpar = "γ",coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 1)", nlpar = "γ",coef = "DAMAGE_NUM2925"),
set_prior("normal(0, 0.002)" , class = "sd",group="SPEC_NUM", nlpar = "w",coef = "DAMAGE_NUM1250"),
set_prior("normal(0, 0.001)" , class = "sd",group="SPEC_NUM", nlpar = "w",coef = "DAMAGE_NUM1510"),
set_prior("normal(0, 0.006)" , class = "sd",group="SPEC_NUM", nlpar = "w",coef = "DAMAGE_NUM2925"),
set_prior("lkj(1)", class ="cor"))

form = bf(f1,nl = TRUE)+list(k~0+DAMAGE_NUM, γ~0+DAMAGE_NUM, u~0+DAMAGE_NUM, w~0+DAMAGE_NUM|SPEC_NUM, β1~0+DAMAGE_NUM,β2~0+DAMAGE_NUM,β3~0+DAMAGE_NUM,sigma~0+DAMAGE_NUM)
n2_b <- brm(form,data = rdat,prior = prior_11,chains=1,iter=5000, warmup=2000,  control = list(adapt_delta = 0.999999999))


I used SPEC_NUM (48 levels) to serve as my “n”.

Note : k must be negative, u_j and γ_j must be positive.

Thank you all.

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

• Operating System:
• brms Version:

I made some edits so your code is easier to read. It’s best if you can use the  to enclose chunks of code.

Thank you.

Can anyone help please?

The question might be a little too specific for a confident answer of mine but your model looks ok to me. You could simplify it by using a linear formula for **x**_i` as explaind in the brms doc and brms_nonlinear vignette instead of manually writing out \beta_1*x_i … manually.

Okay. I will try that.

Thank you.