paul-buerkner/brms

Error in fcor() for nonlinear model

Closed this issue · 1 comments

Dear brms experts,

I get the following error "Error in get(covars[i], data) : object 'R' not found" when fitting a nonlinear model with known residual covariance matrix. Below, I have created a small simulated example that recreates the error. Note, that brms() runs smoothly for a linear model, using the sampling error covariance matrix and fits the nonlinear model using the sampling standard error.

I used brms_2.21.1 under R version 4.3.1.

Any ideas?

Many thanks in advance,
Anders Strathe


Simulated example
library(brms)

Simulated aggregated data
sim = data.frame(y=c(-2.56,-4.573,-7.622,-9.656,-10.792,-11.705,-12.312,-12.645,-12.977,-12.77),
se=c(0.201,0.26,0.388,0.48,0.61,0.715,0.853,1.033,1.155,1.218),
time=c(2,4,8,12,16,20,28,36,44,52) )
str(sim)
sampling error covariance matrix
R = matrix(c(0.04,0.042,0.055,0.062,0.065,0.075,0.061,0.071,0.083,0.085,0.042,0.067,0.087,0.095,0.11,0.126,0.126,0.138,0.152,0.156,0.055,0.087,0.151,0.172,0.202,0.236,0.249,0.272,0.286,0.283,0.062,0.095,0.172,0.231,0.275,0.314,0.334,0.372,0.394,0.394,0.065,0.11,0.202,0.275,0.373,0.42,0.463,0.532,0.569,0.565,0.075,0.126,0.236,0.314,0.42,0.511,0.576,0.664,0.713,0.705,0.061,0.126,0.249,0.334,0.463,0.576,0.727,0.842,0.897,0.904,0.071,0.138,0.272,0.372,0.532,0.664,0.842,1.067,1.165,1.194,0.083,0.152,0.286,0.394,0.569,0.713,0.897,1.165,1.335,1.382,0.085,0.156,0.283,0.394,0.565,0.705,0.904,1.194,1.382,1.483), 10, 10)

Fit linear discrete time-course model
fit <- brm(y ~ 0 + factor(time) + fcor(R),
data = sim,
data2 = list(R = R),
prior = prior(constant(1), class="sigma"),
control = list(adapt_delta = 0.95, max_treedepth = 20),
iter = 2000, chains = 2, cores = 2)

fit

Fit nonlinear time-course model, using sampling standard error
nonlin <- brm(bf(y | se(se) ~ emaxT*(1 - exp( -exp(kT)*time )),
emaxT ~ 1, kT ~ 1,
nl=TRUE),
data = sim,
prior = prior(normal(-10, 100), nlpar = "emaxT") +
prior(normal(-2.5,100), nlpar = "kT"),
control = list(adapt_delta = 0.95, max_treedepth = 20),
iter = 2000, chains = 2, cores = 2)
nonlin

Fit nonlinear time-course model, including sampling covariance matrix
nonlin2 <- brm(bf(y ~ emaxT*(1 - exp( -exp(kT)*time )) + fcor(R),
emaxT ~ 1, kT ~ 1,
nl=TRUE),
data = sim,
data2 = list(R = R),
prior = prior(normal(-10, 100), nlpar = "emaxT") +
prior(normal(-2.5,100), nlpar = "kT") +
prior(constant(1), class="sigma"),
control = list(adapt_delta = 0.95, max_treedepth = 20),
iter = 2000, chains = 2, cores = 2)
nonlin2

Oh, I do see the issue. You cannot add fcor() (or any other autocor term) to the non-linear formula. in this case, use the autocor argument