Checking my setup on multivariate nonlinear model using brms

Interfaces brms
1

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:
2

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

3

Thank you.

Can anyone help please?

4

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.

5

Okay. I will try that.

Thank you.