I'd like to build a nonlinear mixed effect model that describes the relationship between two variables, "x" and "y", which vary randomly by a third variable "r" using an exponential rise to a maximum as described by the equation: y = theta(1-exp(-beta*x)).
I've been able to create the nonlinear model for x and y using nls(), but I have not been successful in incorporating a random effect into nlme().
When I build the model using nlme() I end up with an error message: "Error in eval(predvars, data, env) : object 'theta' not found". This error is unexpected to me since the nls() model ran without issue using the same dataframe.
To build the dataset:
x = c(33,35,16,8,31,31,31,23,7,7,7,7,11,11,3,3,6,6,32,32,1,17,17,17,25,40,40,6,6,29,29,13,23,23,44,44,43,43,13,4,6,15,17,22,28,8,11,22,32,6,12,20,27,15,29,29,29,29,29,12,12,16,16,12,12,2,49,49,14,14,14,37,2.87,4.86,7.90,11.95,16.90,16.90,18.90,18.89,22.00,24.08,27.14,30.25,31.22,32.26,7,14,19,31,36,7,14,19,31,36,7,16,16,16,16,16,16,32,32,32,32,32,32,11,11,11,13,13,13,13,13,13,13,13,13,13,13,13,9,9)
y = c(39.61,32.66,27.06,21.74,22.18,38.19,35.02,23.13,9.70,14.20,13.40,15.30,18.80,19.00,3.80,4.90,15.00,14.20,24.90,16.56,1.76,29.29,28.49,18.64,27.10,9.47,14.14,10.27,8.44,26.15,25.43,22.00,19.00,13.00,73.19,67.76,32.34,36.86,8.00,1.57,8.33,16.20,14.69,18.95,20.52,4.92,8.28,15.27,18.37,6.60,10.98,12.56,19.04,5.49,21.00,12.90,17.30,11.40,12.20,15.63,15.22,33.80,17.78,19.33,3.86,8.57,30.40,13.39,11.93,4.55,6.18,12.70,2.71,7.23,5.61,22.74,15.71,16.95,18.31,20.78,17.64,20.00,19.52,24.86,30.06,24.92,4.17,11.02,10.08,14.94,25.98,0.00,3.67,3.67,6.69,11.90,5.06,13.21,10.33,0.00,0.00,6.47,8.38,28.57,25.26,28.67,27.92,33.69,29.61,6.11,7.13,6.93,4.81,15.34,4.90,14.94,8.88,10.24,8.80,10.46,10.48,9.19,9.67,9.40,24.98,50.79)
r = c("A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","C","C","D","E","E","E","F","G","G","H","H","H","H","I","I","I","J","J","J","J","K","L","L","L","L","L","L","L","L","L","L","L","L","L","L","M","N","N","N","N","N","O","P","P","P","P","P","Q","R","R","S","S","S","T","U","U","U","U","U","U","U","U","U","U","U","U","U","U","V","V","V","V","V","V","V","V","V","V","W","X","X","X","X","Y","Y","Z","Z","Z","Z","Z","Z","AA","AA","AA","AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","AC","AC")
df = data.frame(x,y,r)
To build the nonlinear model without "r" as a random effect.
nls_test = nls(y~theta*(1-exp(-beta*x)),
data = df,
start = list(beta = 0.2, theta = 38),
trace = TRUE)
In my model, the only fixed effect is x and the only random effect is r. I've tried building an nlme() model that reflects this, based on the nlme package documentation (https://cran.r-project.org/web/packages/nlme/nlme.pdf),more specifically these lines of code found on page 186 of the documentation linked above.
The nlme() object that I've tried to create with my data is as follows:
nlme_test = nlme(y ~ theta*(1-exp(-beta*x)),
fixed = x~1,
random = r~1,
data = df,
start = c(theta = 38,
beta = 0.2))
And results in the following error.
Error in eval(predvars, data, env) : object 'theta' not found
From what I gather, this is related to 'theta' not being included in the dataframe ("df") used to build the nlme object, but it is unclear to me why this occurs as most examples that I have found for this error are related to the use of the predict() function and missing column or disagreement between column names.
Also, since the nls() model (nls_Test) worked fine using the same start = c(theta = 38, beta = 0.2) and without a 'theta' or 'beta' data column in df, I'm a bit confused as to why I'm receiving this error about column name error.
Does anyone have suggestions or references to help me incorporate the random effect into my nlme model?
Thanks!
Expanding on my (now deleted, because incomplete) comment, I assume this is what you want to do. Please confirm carefully by reading the help page about nlme (i.e. ?nlme::nlme).
nlme_test <- nlme(y ~ theta*(1-exp(-beta*x)),
fixed = theta + beta ~ 1,
random = theta + beta ~ 1,
groups = ~ r,
data = df,
start = c(theta = 38,
beta = 0.2))
The fixed and random arguments should not name the variables in your model formula but the regression parameters. This way the function knows which parts of the model are variables (to be found in data) and which parts are parameters. Also, you missed the groups argument in order to specify how the data is clustered.
Output:
summary(nlme_test)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: y ~ theta * (1 - exp(-beta * x))
## Data: df
## AIC BIC logLik
## 887.6224 904.6401 -437.8112
##
## Random effects:
## Formula: list(theta ~ 1, beta ~ 1)
## Level: r
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## theta 1.145839e+01 theta
## beta 1.061366e-05 0.01
## Residual 6.215030e+00
##
## Fixed effects: theta + beta ~ 1
## Value Std.Error DF t-value p-value
## theta 21.532188 2.8853414 96 7.462614 0e+00
## beta 0.104404 0.0251567 96 4.150144 1e-04
## Correlation:
## theta
## beta -0.548
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.89510795 -0.51882772 -0.09466037 0.34471808 3.66855121
##
## Number of Observations: 126
## Number of Groups: 29