exponential decay - from NLS to nlme

Question

I am somewhat familiar with lme4, but not with the nlme package, and I struggling to fit the mixed-effects that I desire.

I ran the following model with NLS. (Comments about the starting values are welcome too)

m3 <- nls(s  ~ f0 + f1 * 2^(-MINUTES / K),
          data = dat ,
          start = list(f0 = 117, f1 = 16, K = 2.5))
summary(m3)

Q1 - What's the code to run the previous model with nlme, but with a random intercept for "ID"?

I tried

m1 <- nlme(s  ~ f0 + f1 * 2^(-MINUTES / K),
    data = dat ,
    fixed = f0 + f1 + MINUTES + K ~ 1,
    random = ID ~ 1 , 
    start = list(f0 = 117, f1 = 16, K = 2.5))

But gave error

Error in getGroups.data.frame(dataMix, eval(parse(text = paste("~1", deparse(groups[[2]]), : invalid formula for groups

Q2 - Same as before, but also with MINUTES*group interaction?

dataframe:

  dat <-  structure(list(ID = c(1, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 
10, 10, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 
13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 15, 15, 
15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 
17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 2, 2, 2, 2, 2, 2, 2, 
3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 40, 40, 40, 40, 40, 
40, 40, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 
43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 45, 45, 
45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 46, 46, 49, 49, 49, 49, 
49, 49, 49, 5, 5, 5, 5, 5, 5, 5, 52, 52, 52, 52, 52, 52, 52, 
53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 54, 55, 55, 55, 
55, 55, 55, 55, 56, 56, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 
57, 57, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 
8, 8, 8, 9, 9, 9, 9, 9, 9, 9), MINUTES = c(0, 5, 10, 15, 20, 
25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 
5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 
20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 
0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 
15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 
30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 
10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 
25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 
5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 
20, 25, 30, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 
0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 
15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 
30, 0, 5, 10, 15, 20, 25, 30, 0, 5, 10, 15, 20, 25, 30), group = c("fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff", "fff", "fff", "fff", "fff", "fff", 
"fff", "fff"), s = c(148, 151, 142, 147, 141, 142, 145, 
131.5, 124, 114, 117, 119, 127, 119, 155.5, 110, 112, 111, 110, 
108, 110, 149.5, 128, 117, 123, 125, 124, 130, 131, 117, 119, 
125, 129, 125, 125, 147.5, 132, 117, 114, 116, 117, 115, 133, 
125, 122, 125, 127, 124, 128, 106, 120, 115, 116, 119, 113, 115, 
121, 97, 114, 99, 97, 102, 102, 133.5, 105, 104, 107, 103, 104, 
108, 131, 118, 122, 121, 127, 126, 119, 133, 129, 130, 130, 129, 
124, 126, 125, 117, 115, 116, 114, 122, 119, 137.5, 116, 111, 
117, 111, 117, 108, 124, 124, 125, 123, 119, 121, 112, 136.5, 
123, 126, 123, 119, 127, 112, 120.5, 104, 100, 107, 110, 106, 
104, 116.5, 117, 116, 117, 120, 113, 121, 149, 134, 117, 115, 
123, 117, 110, 146, 129, 131, 126, 127, 137, 128, 129.5, 105, 
119, 125, 107, 120, 120, 113, 118, 114, 109, 109, 117, 112, 120, 
110, 116, 107, 112, 111, 113, 108, 103, 105, 106, 104, 111, 123, 
84, 84, 87, 88, 81, 84, 116, 97, 107, 99, 96, 94, 104, 133, 127, 
118, 127, 123, 119, 121, 140.5, 116, 118, 116, 121, 121, 116, 
143, 138, 137, 130, 129, 119, 133, 120, 108, 98, 103, 105, 100, 
97, 130, 107, 113, 114, 111, 106, 108, 156, 137, 138, 134, 142, 
135, 145)), row.names = c(8L, 9L, 10L, 11L, 12L, 13L, 14L, 22L, 
23L, 24L, 25L, 26L, 27L, 28L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 
43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 
56L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 78L, 79L, 80L, 81L, 82L, 
83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 99L, 100L, 101L, 
102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 112L, 
120L, 121L, 122L, 123L, 124L, 125L, 126L, 134L, 135L, 136L, 137L, 
138L, 139L, 140L, 148L, 149L, 150L, 151L, 152L, 153L, 154L, 155L, 
156L, 157L, 158L, 159L, 160L, 161L, 162L, 163L, 164L, 165L, 166L, 
167L, 168L, 169L, 170L, 171L, 172L, 173L, 174L, 175L, 176L, 177L, 
178L, 179L, 180L, 181L, 182L, 183L, 184L, 185L, 186L, 187L, 188L, 
189L, 190L, 191L, 192L, 193L, 194L, 195L, 196L, 197L, 198L, 199L, 
200L, 201L, 202L, 203L, 204L, 205L, 206L, 207L, 208L, 209L, 210L, 
218L, 219L, 220L, 221L, 222L, 223L, 224L, 225L, 226L, 227L, 228L, 
229L, 230L, 231L, 232L, 233L, 234L, 235L, 236L, 237L, 238L, 239L, 
240L, 241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 249L, 250L, 
251L, 252L, 253L, 254L, 255L, 256L, 257L, 258L, 259L, 260L, 261L, 
262L, 263L, 264L, 265L, 273L, 274L, 275L, 276L, 277L, 278L, 279L, 
287L, 288L, 289L, 290L, 291L, 292L, 293L, 294L, 295L, 296L, 297L, 
298L, 299L, 300L, 308L, 309L, 310L, 311L, 312L, 313L, 314L), class = "data.frame")

score 1 · Answer 1 · answered Oct 25 '22 at 19:01

I believe these would be:

library(nlme)
m1 <- nlme(s  ~ f0 + f1 * 2^(-MINUTES / K),
    data = dat ,
    fixed = f0 + f1 + K ~ 1,
    random = f0 ~ 1 | ID,
    start = c(f0 = 117, f1 = 16, K = 2.5))


m2 <- nlme(s  ~ f0 + f1 * 2^(-MINUTES / K),
    data = dat ,
    fixed = f0 + f1 + K ~ 1,
    random = f0 + f1 ~ 1 | ID,
    start = c(f0 = 117, f1 = 16, K = 2.5))

By "MINUTES*group interaction" I am assuming you only care about the magnitude (f1) and not the scale (K) of the MINUTES effect to vary by group: in principle random = f0 + f1 + K ~ 1 | ID is possible, but I ran into a "maximum number of iterations reached without convergence" error. You can extend that (e.g. control = nlmeControl(maxIter = 1000)), but I ran into more problems down the line ("step halving factor reduced below minimum"). You're likely to run into more and more computational issues as you make the model that complicated, unless you have very clean data ...

Thanks a lot once again, Ben. It works! It was surprising to not see "MINUTES" listed in "fixed =....". If I may, do you have any thoughts about weather this model or splines with two linear equations or splines with 1 linear + 1 quadratic would be a better fit? Thank you. — jtjtjtjt, Oct 26 '22 at 08:27
`MINUTES` is a covariate, not a parameter. `fixed` is about the variation of *parameters* across groups/with respect to covariates (i.e., linear submodels for the parameters). I don't have strong feelings about the model - seems context-dependent. But this seems reasonable. — Ben Bolker, Oct 26 '22 at 12:48

exponential decay - from NLS to nlme

1 Answers1