Data
I have a longitudinal dataset for plant growth that includes repeated measurements in different seasons. It is a large dataset with over 100 individuals,so I can't put all the data here but the data look like this
> df
head(dat)
# A tibble: 6 x 5
id season Lt_1 td Lobs
<fct> <fct> <dbl> <dbl> <dbl>
1 2020-15 s1 35 60 40.2
2 2020-15 s1 40.2 94 41.7
3 2020-20 s1 37.3 60 41.7
4 2020-21 s1 35.5 59 39.4
5 2020-24 s1 36.6 60 41.8
6 2020-24 s1 41.8 80 45.1
Model to estimate parameters
I am trying to estimate the parameters for a growth equation using this dataset. The growth equation is as follows:
Lobs = 150*((1 + ((150/Lt_1)^(1/p)-1)exp(-k*td/365))^(-p))
Lt represents the plant height at a measurement, Lt_1 represents the plant height at the previous measurement, and td represents the duration of time between measurements. The parameters I want to estimate are k and p. However, I believe that the parameter k may vary depending on the season, so I want to have specific k values for each season. Additionally, I would like to include individual variation as a random effect in the model.
Issue
Now, I am trying to estimate the parameters using the nlme function from the nlme package. However, I am getting the error below.
library(nlme)
# Model formula
formula = Lobs ~ 150*((1 + ((150/Lt_1)^(1/p)-1)*exp(-k*td/365))^(-p))
# specify starting values for the parameters
p = 1.2
k = c(0.5, 1.09, 0.75, 1.09)
# fitting
model <- nlme(formula,
fixed = c(p ~ 1, k ~ 1 + season),
random = k ~ 1 | id,
data = df_r,
start = list(fixed = c(p, k)),
control=list(maxIter=1e6, msMaxIter = 1e6, msVerbose = TRUE)
)
0: 1150.3155: 2.85771
1: 1150.3141: 2.86659
2: 1150.3141: 2.86654
0: 1103.1191: 4.57637
1: 1103.1180: 4.65814
2: 1103.1172: 4.62838
3: 1103.1172: 4.62391
4: 1103.1172: 4.62430
5: 1103.1172: 4.62430
Error in nlme.formula(formula, fixed = c(p ~ 1, k ~ 1 + season), random = k ~ :
step halving factor reduced below minimum in PNLS step
What I tried
I tried with different starting values and the number of iterations, but none of them worked out.
Only when I remove the season variable from the fixed effect, it succeeded to fit.
p = 1.2
k = 0.2
model <- nlme(formula,
fixed = c(p ~ 1, k ~ 1),
random = k ~ 1 | id,
data = dat,
start = list(fixed = c(p, k)),
control=list(maxIter=1e6, msMaxIter = 1e6, msVerbose = TRUE)
)
> summary(model)
Nonlinear mixed-effects model fit by maximum likelihood
Model: formula
Data: dat
AIC BIC logLik
1631.918 1646.773 -811.9589
Random effects:
Formula: k ~ 1 | id
k Residual
StdDev: 0.02607074 3.500357
Fixed effects: c(p ~ 1, k ~ 1)
Value Std.Error DF t-value p-value
p 4.154633 3.263363 181 1.273114 0.2046
k 0.666089 0.055864 181 11.923354 0.0000
Correlation:
p
k -0.976
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.6905297 -0.8627678 -0.2492133 0.4294071 2.4403191
Number of Observations: 303
Number of Groups: 121
Does anyone have an idea how to solve this problem?
