I have a longitudinal dataset for plant growth, where each row represents an individual plant's measurements over multiple seasons. The variables in the dataset include
Here's an example of the dataset:
library(lme4)
df <- data.frame(
id = c("a", "a", "a", "a", "a", "a", "b", "b", "b", "b", "b", "c", "c", "c", "c", "c", "d", "d", "d", "d", "d", "d", "e", "e", "e", "e", "e", "e"),
season = c("s1", "s1", "s2", "s3", "s3", "s4", "s1", "s1", "s2", "s3", "s3", "s1", "s1", "s2", "s3", "s3", "s1", "s1", "s2", "s3", "s3", "s4", "s1", "s1", "s2", "s3", "s3", "s4"),
td = c(59, 95, 210, 134, 62, 169, 60, 96, 208, 134, 62, 60, 94, 210, 133, 63, 59, 99, 206, 134, 62, 169, 59, 122, 182, 134, 62, 171),
Lt = c(40, 44, 72, 84, 87, 101, 40, 44, 66, 75, 77, 43, 47, 67, 79, 80, 38, 42, 64, 75, 76, 80, 44, 47, 69, 79, 80, 99),
Lt_1 = c(35.5, 40, 44.5, 72, 84.2, 87, 35, 40.1, 43.8, 65.9, 75.4, 40.8, 43, 47, 66.8, 78.7, 34, 37.5, 41.9, 63.6, 75.2, 76.5, 39.2, 43.5, 46.9, 68.7, 79.4, 80.5)
)
head(df)
> head(df)
id season td Lt Lt_1
1 a s1 59 40.0 35.5
2 a s1 95 44.5 40.0
3 a s2 210 72.0 44.5
4 a s3 134 84.2 72.0
5 a s3 62 87.0 84.2
6 a s4 169 101.4 87.0
Now, I am attempting to fit a model to estimate the plant growth using the previous measurement (Lt_1) and the duration of time (td) as predictors. The model formula is as follows:
formula = Lt ~ 150*((1 + ((150/Lt_1)^(1/p)-1)*exp(-k*td/365))^(-p))
Question #1: how to specify the parameters to be estimated in the lmer function?
I want to include a random effect in the parameter "k" for "id" in the model. In the nlme function from the nlme package, I can explicitly specify the parameters to be estimated like this.
library(nlme)
model <- nlme(formula,
fixed = k + p ~ 1,
random = k ~ 1,
groups = ~ id,
data = df,
start = c(k = 0.2,
p = 1))
However, I have no idea how to achieve the same with the lme4::lmer function. A lot of examples with lmer use a generic function to estimate the parameter like [this][1]. Can anybody help me how to do this using a specific function?
Question #2: how to estimate "k" for each season?
I want to estimate the parameter "k" separately for each season in the dataset. I assume that each season may have different "k" values, such as "k1" for "s1" and "k2" for "s2". How can I specify this in the code using the lmer function? I tried something like this but did not work even with nlme function.
formula2 = Lobs ~ 150*((1 + ((150/Lt_1)^(1/p)-1)*exp(-k*season*td/365))^(-p))
model <- nlme(formula2,
fixed = k + p ~ 1,
random = k ~ 1,
groups = ~ id,
data = df,
start = c(k = 0.2,
p = 1))
Thank you for your advice.
lmer simply won't work for what you want to do, and nlmer is unfortunately rather immature — if you can make it work, nlme is probably still the best choice. I think you want something like this, but the data set you've given is a bit too small to figure out if the problems I'm having are due to the setup or to the limitations of the data.
It is almost always worthwhile to set up a function that can provide the gradient as well as the objective function:
dfun <- deriv(formula[-2], ## RHS of formula
namevec = c("p", "k"),
function.arg = c("k", "p", "Lt_1", "td"))
## test
dfun(p = 1, k = 0.2, Lt_1 = 35, Lt = 40)
## [1] 35.74052
## attr(,"gradient")
## p k
## [1,] 0.4133755 3.7294
Something like this (I prefer to specify grouping variables on the right side of a bar (|) rather than via the groups argument):
library(nlme)
model <- nlme(Lt ~ dfun(k, p, Lt_1, td),
fixed = list(p ~ 1, k ~ season),
random = k ~ 1| id,
data = df,
start = list(fixed = c(1, 0.2, 0, 0, 0))
)
This gets to
Error in nlme.formula(Lt ~ dfun(k, p, Lt_1, td), fixed = list(p ~ 1, k ~ : step halving factor reduced below minimum in PNLS step
but it's hard to know in general whether this means there's something wrong with the model setup/starting values/etc. or whether the example data set is too small to work — could be either.