I am trying to build a somewhat complex model using nls that has a lot of precedent in my field, but I repeatedly get the error " Missing value or an infinity produced when evaluating the model". From previous posts I know this is likely driven by bad starting values, but I am having trouble generating higher quality starting values. My model takes the general structure:
efflux = a * e(b * temperature) + ln(c + d*moisture) .
I was able to generate strong starting values for parameters a and b, because all of the equation left of the plus sign can be separated into its own simple equation, the log can be taken of both sides, and then the linear model of these log transformed data can easily have their starting parameters extracted from them. However, I have trouble extracting parameters for c and d, as that same procedure doesn't work nearly as well for the full equation. Here is my code, attempting to haphazardly guess at some c and d starting parameters:
library(dplyr)
model_efflux <- d_resp_data_2020 %>%
group_by(stand, plot) %>%
do(model = nls(avefflux ~ a * exp(b * avtemp_data_2020) + log(c + d*(moisture_data_2020)), start = list(a = 0.8, b = 0.1, c = -1, d = 1), data = .)) %>% ungroup()
And here is some data to work with:
,timestamp,timestampfull,date,stand,plot,position,depth,moisture_data_2020,avefflux,avtemp_data_2020
1,9:30,2020-11-16 09:30:00,11/16/2020,ac1911,2,9,d,11.8,0.691,5.85
2,10:00,2020-11-16 10:00:00,11/16/2020,ac1911,2,15,d,14.5,0.455,6.11
3,10:00,2020-11-16 10:00:00,11/16/2020,ac1911,1,3,d,16.7,0.345,5.34
4,10:30,2020-11-16 10:30:00,11/16/2020,ac1911,1,9,d,14.5,1.11,6.44
5,10:30,2020-11-16 10:30:00,11/16/2020,ac1911,1,15,d,9.2,0.861,5.72
6,11:30,2020-11-16 11:30:00,11/16/2020,ac1911,1,3,d,11.4,0.236,5.61
7,12:00,2020-11-16 12:00:00,11/16/2020,ac1911,1,9,d,15.4,1.28,5.53
8,12:00,2020-11-16 12:00:00,11/16/2020,ac1911,1,15,d,12.9,0.651,5.88
9,12:00,2020-11-16 12:00:00,11/16/2020,ac1911,2,3,d,16.1,0.781,5.78
10,12:30,2020-11-16 12:30:00,11/16/2020,ac1911,2,9,d,12.4,0.839,5.71
11,13:00,2020-11-16 13:00:00,11/16/2020,ac1911,2,15,d,15.9,1.295,6.02
12,14:00,2020-11-16 14:00:00,11/16/2020,ac1911,7,3,d,10.1,0.831,5.46
13,14:30,2020-11-16 14:30:00,11/16/2020,ac1911,7,9,d,11.4,0.626,5.75
14,14:30,2020-11-16 14:30:00,11/16/2020,ac1911,7,15,d,14.2,0.686,5.22
15,15:30,2020-11-16 15:30:00,11/16/2020,ac1911,8,3,d,16.7,0.611,5.77
16,15:30,2020-11-16 15:30:00,11/16/2020,ac1911,8,9,d,12.1,0.954,5.11
17,16:00,2020-11-16 16:00:00,11/16/2020,ac1911,8,15,d,7.6,0.709,5.91
18,16:30,2020-11-16 16:30:00,11/16/2020,ac1911,1,3,d,15.9,0.614,6.06
19,16:30,2020-11-16 16:30:00,11/16/2020,ac1911,1,9,d,13.6,0.971,5.96
20,16:30,2020-11-16 16:30:00,11/16/2020,ac1911,1,15,d,15.2,1.235,5.84
21,11:30,2020-11-17 11:30:00,11/17/2020,ac1911,104,3,d,11.8,0.544,4.7
22,11:30,2020-11-17 11:30:00,11/17/2020,ac1911,104,9,d,11.4,0.228,3.68
23,12:00,2020-11-17 12:00:00,11/17/2020,ac1911,104,15,d,10.5,0.657,4.38
24,12:30,2020-11-17 12:30:00,11/17/2020,ac1911,103,3,d,8.8,0.581,3.99
25,12:30,2020-11-17 12:30:00,11/17/2020,ac1911,103,9,d,11.4,0.518,3.99
26,13:00,2020-11-17 13:00:00,11/17/2020,ac1911,103,15,d,10.9,0.418,3.53
27,15:30,2020-11-17 15:30:00,11/17/2020,ac1911,2,3,d,17.2,0.709,3.3
28,15:30,2020-11-17 15:30:00,11/17/2020,ac1911,2,9,d,10.7,0.568,3.47
29,16:00,2020-11-17 16:00:00,11/17/2020,ac1911,2,15,d,9.5,0.949,3.77
30,16:00,2020-11-17 16:00:00,11/17/2020,ac1911,3,3,d,11.9,0.76,4.92
31,16:00,2020-11-17 16:00:00,11/17/2020,ac1911,3,9,d,12.9,0.794,5.17
32,16:30,2020-11-17 16:30:00,11/17/2020,ac1911,3,15,d,8.9,0.602,4.3
33,16:30,2020-11-17 16:30:00,11/17/2020,ac1911,1,3,d,9.6,0.762,3.51
34,16:30,2020-11-17 16:30:00,11/17/2020,ac1911,1,9,d,11.3,0.435,3.16
35,17:00,2020-11-17 17:00:00,11/17/2020,ac1911,1,15,d,11.2,0.631,3.79
36,16:00,2020-11-18 16:00:00,11/18/2020,ac1911,e,3,d,11.2,1.33,4.45
37,16:00,2020-11-18 16:00:00,11/18/2020,ac1911,e,9,d,14.1,2.135,5.21
38,16:00,2020-11-18 16:00:00,11/18/2020,ac1911,e,15,d,10.9,0.949,4.48
39,16:30,2020-11-18 16:30:00,11/18/2020,ac1911,x,3,d,16.9,2.565,5.08
40,16:30,2020-11-18 16:30:00,11/18/2020,ac1911,x,9,d,11.6,0.891,4.6
41,17:00,2020-11-18 17:00:00,11/18/2020,ac1911,x,15,d,10.9,1.145,4.71
42,17:00,2020-11-18 17:00:00,11/18/2020,ac1911,f,3,d,14.4,0.778,5.02
43,17:00,2020-11-18 17:00:00,11/18/2020,ac1911,f,9,d,10.1,1.165,5.87
44,17:30,2020-11-18 17:30:00,11/18/2020,ac1911,f,15,d,10.7,0.726,4.52
45,9:00,2020-11-19 09:00:00,11/19/2020,ac1911,3,3,d,10.6,0.292,2.53
46,9:30,2020-11-19 09:30:00,11/19/2020,ac1911,3,9,d,9.3,0.613,3.46
47,9:30,2020-11-19 09:30:00,11/19/2020,ac1911,3,15,d,9.9,0.438,2.43
48,9:30,2020-11-19 09:30:00,11/19/2020,ac1911,2,3,d,12.4,0.502,3.4
49,10:00,2020-11-19 10:00:00,11/19/2020,ac1911,2,9,d,10.6,0.33,3.08
50,10:00,2020-11-19 10:00:00,11/19/2020,ac1911,2,15,d,9.9,0.326,2.89
51,10:00,2020-11-19 10:00:00,11/19/2020,ac1911,1,3,d,8.8,0.645,3.78
52,10:30,2020-11-19 10:30:00,11/19/2020,ac1911,1,15,d,9,0.609,3.5
53,11:00,2020-11-19 11:00:00,11/19/2020,ac1911,2,3,d,8.8,0.74,4.14
54,11:00,2020-11-19 11:00:00,11/19/2020,ac1911,2,9,d,10.3,0.782,4.58
55,11:00,2020-11-19 11:00:00,11/19/2020,ac1911,2,15,d,10.5,0.735,4.49
56,11:30,2020-11-19 11:30:00,11/19/2020,ac1911,1,3,d,11.9,0.869,4.6
57,11:30,2020-11-19 11:30:00,11/19/2020,ac1911,1,15,d,10.4,1.024,5.09
58,12:00,2020-11-19 12:00:00,11/19/2020,ac1911,3,3,d,11.2,1.205,5.13
59,12:00,2020-11-19 12:00:00,11/19/2020,ac1911,3,9,d,12.6,1.048,4.7
60,12:30,2020-11-19 12:30:00,11/19/2020,ac1911,3,15,d,16.2,0.663,5.21
61,12:30,2020-11-19 12:30:00,11/19/2020,ac1911,1,3,d,14.8,0.928,5.2
62,13:00,2020-11-19 13:00:00,11/19/2020,ac1911,1,9,d,11.9,0.923,5.69
63,13:00,2020-11-19 13:00:00,11/19/2020,bf52,1,15,d,12.9,1.62,5.52
64,13:00,2020-11-19 13:00:00,11/19/2020,bf52,2,3,d,10.8,0.582,5.61
65,13:00,2020-11-19 13:00:00,11/19/2020,bf52,2,9,d,7.8,1.255,5.63
66,13:30,2020-11-19 13:30:00,11/19/2020,bf52,2,15,d,8.3,1.11,5.76
67,13:30,2020-11-19 13:30:00,11/19/2020,bf52,3,3,d,9.7,0.821,5.99
68,13:30,2020-11-19 13:30:00,11/19/2020,bf52,3,9,d,8.8,1.35,5.68
69,15:00,2020-11-19 15:00:00,11/19/2020,bf52,3,15,d,10.3,1.215,5.71
70,15:00,2020-11-19 15:00:00,11/19/2020,bf52,1,3,d,11.8,0.703,6.21
71,15:00,2020-11-19 15:00:00,11/19/2020,bf52,1,9,d,8.6,0.821,6.68
72,15:30,2020-11-19 15:30:00,11/19/2020,bf52,1,15,d,7.3,0.629,6.64
73,15:30,2020-11-19 15:30:00,11/19/2020,bf52,2,3,d,10.5,0.45,5.12
74,15:30,2020-11-19 15:30:00,11/19/2020,bf52,2,9,d,11.2,0.788,5.34
75,16:00,2020-11-19 16:00:00,11/19/2020,bf52,2,15,d,12.9,1.135,5.6
76,16:00,2020-11-19 16:00:00,11/19/2020,bf52,107,3,d,7.7,0.884,6.69
77,16:00,2020-11-19 16:00:00,11/19/2020,bf52,107,9,d,10.5,1.095,6.1
78,16:30,2020-11-19 16:30:00,11/19/2020,bf52,107,15,d,10.5,1.195,6.16
79,10:00,2020-11-20 10:00:00,11/20/2020,bf52,1,3,d,12.9,1.59,6.36
80,10:00,2020-11-20 10:00:00,11/20/2020,bf52,1,9,d,13.2,1.73,6.61
81,10:00,2020-11-20 10:00:00,11/20/2020,bf52,1,15,d,11,2.37,6.46
82,10:30,2020-11-20 10:30:00,11/20/2020,bf52,2,3,d,9.8,1.51,6.76
83,10:30,2020-11-20 10:30:00,11/20/2020,bf52,2,9,d,NA,NA,NA
84,11:00,2020-11-20 11:00:00,11/20/2020,bf52,2,15,d,14.4,0.778,6.53
85,11:30,2020-11-20 11:30:00,11/20/2020,bf52,3,3,d,9.8,0.987,6.02
86,11:30,2020-11-20 11:30:00,11/20/2020,bf52,3,9,d,10.1,1.775,6.54
87,11:30,2020-11-20 11:30:00,11/20/2020,bf52,3,15,d,NA,NA,NA
88,12:30,2020-11-20 12:30:00,11/20/2020,bf52,a,3,d,12.3,0.802,7.45
89,12:30,2020-11-20 12:30:00,11/20/2020,bf52,a,9,d,12.4,1.27,7.36
90,12:30,2020-11-20 12:30:00,11/20/2020,bf52,a,15,d,7.4,1.24,7.46
91,13:00,2020-11-20 13:00:00,11/20/2020,bf52,c,3,d,13.4,1.25,7.75
92,13:00,2020-11-20 13:00:00,11/20/2020,bf52,c,9,d,7.9,1.375,7.86
93,13:00,2020-11-20 13:00:00,11/20/2020,bf52,c,15,d,10.1,0.976,7.3
94,13:30,2020-11-20 13:30:00,11/20/2020,bf52,b,3,d,10.2,1.255,7.28
95,13:30,2020-11-20 13:30:00,11/20/2020,bf52,b,9,d,10.7,1.4,7.26
96,13:30,2020-11-20 13:30:00,11/20/2020,bf52,b,15,d,8.4,1.55,7.5
97,10:10,2020-07-28 10:00:00,20200728,bf52,1,3,d,6.3,3.335,21.33
98,10:24,2020-07-28 10:30:00,20200728,bf52,1,9,d,5.8,6.325,20.82
99,10:31,2020-07-28 10:30:00,20200728,bf52,1,15,d,4.9,3.445,21.23
100,10:42,2020-07-28 10:30:00,20200728,bf52,2,3,d,9.3,2.12,22.55
101,10:52,2020-07-28 11:00:00,20200728,bf52,2,9,d,6.4,4.155,22.52
102,10:59,2020-07-28 11:00:00,20200728,bf52,2,15,d,11.2,6.135,22.08
103,11:16,2020-07-28 11:30:00,20200728,bf52,1,3,d,10.3,4.965,21.89
104,11:57,2020-07-28 12:00:00,20200728,bf52,1,15,d,7.7,4.52,20.99
Can anyone advise me on why this error might be coming up, and a better way to maybe get some starting values?
Using nlsLM in minpack.lm we get the following. Note that while this function is more likely to provide an answer than nls it is often at the expense of premature convergence.
Note that we negated b in the formula.
Below, we first fit the combined model and then use the result from that as starting values to fit the individual groups.
We added an n column to show the number of rows in each group.
library(dplyr)
library(minpack.lm)
fo <- avefflux ~ a * exp(-b * moisture_data_2020) + log(c + d*avtemp_data_2020)
st <- list(a = 0.8, b = 0.1, c = 0.1, d = 1)
fm <- nlsLM(fo, d_resp_data_2020, start = st)
model_efflux <- d_resp_data_2020 %>%
group_by(stand, plot) %>%
mutate(n = n()) %>%
group_by(n, .add = TRUE) %>%
do(model = nlsLM(fo, data = ., start = coef(fm))) %>%
ungroup %>%
rowwise %>%
mutate(a = coef(model)[1], b = coef(model)[2], c = coef(model)[3],
d = coef(model)[4]) %>%
ungroup
giving the following. It seems it used the overall fit in many cases and in other cases the fit seems highly suspect.
> model_efflux
# A tibble: 17 x 8
stand plot n model a b c d
<chr> <chr> <int> <list> <dbl> <dbl> <dbl> <dbl>
1 ac1911 1 18 <nls> 311. 1.02 0.977 0.241
2 ac1911 103 3 <nls> 30.2 0.563 -0.930 0.677
3 ac1911 104 3 <nls> 30.2 0.563 -0.930 0.677
4 ac1911 2 14 <nls> 2642. 1.09 0.838 0.251
5 ac1911 3 9 <nls> -38121775. 2.19 0.570 0.371
6 ac1911 7 3 <nls> 30.2 0.563 -0.930 0.677
7 ac1911 8 3 <nls> 30.2 0.563 -0.930 0.677
8 ac1911 e 3 <nls> 30.2 0.563 -0.930 0.677
9 ac1911 f 3 <nls> 30.2 0.563 -0.930 0.677
10 ac1911 x 3 <nls> 30.2 0.563 -0.930 0.677
11 bf52 1 12 <nls> -239464524. 4.02 -21.6 4.29
12 bf52 107 3 <nls> 30.2 0.563 -0.930 0.677
13 bf52 2 12 <nls> 24349. 1.73 -11.3 2.48
14 bf52 3 6 <nls> 196010890. 2.28 -9.55 2.12
15 bf52 a 3 <nls> 30.2 0.563 -0.930 0.677
16 bf52 b 3 <nls> 30.2 0.563 -0.930 0.677
17 bf52 c 3 <nls> 30.2 0.563 -0.930 0.677
The overall fit is:
> fm
Nonlinear regression model
model: avefflux ~ a * exp(-b * moisture_data_2020) + log(c + d * avtemp_data_2020)
data: d_resp_data_2020
a b c d
30.1573 0.5629 -0.9301 0.6775
residual sum-of-squares: 45.99
Number of iterations to convergence: 15
Achieved convergence tolerance: 1.49e-08
We redo this using nls. Because the a values above varied wildly we fixed the a value at the overall fit's value for a and then fit the others separately. As expected nls fails for the small groups where it would have otherwise given an answer that is likely invalid.
fo <- avefflux ~ a * exp(-b * moisture_data_2020) + log(c + d*avtemp_data_2020)
st <- list(a = 0.8, b = 0.1, c = 0.1, d = 1)
fm <- nls(fo, d_resp_data_2020, start = st)
a <- coef(fm)[1]
model_efflux <- d_resp_data_2020 %>%
group_by(stand, plot) %>%
mutate(n = n()) %>%
group_by(n, .add = TRUE) %>%
do(model = try(nls(fo, data = ., start = coef(fm)[-1]))) %>%
ungroup %>%
rowwise %>%
mutate(a = if (inherits(model, "nls")) a else NA,
b = if (inherits(model, "nls")) coef(model)[["b"]] else NA,
c = if (inherits(model, "nls")) coef(model)[["c"]] else NA,
d = if (inherits(model, "nls")) coef(model)[["d"]] else NA) %>%
ungroup
> model_efflux
# A tibble: 17 x 8
stand plot n model a b c d
<chr> <chr> <int> <list> <dbl> <dbl> <dbl> <dbl>
1 ac1911 1 18 <nls> 30.2 0.774 0.984 0.239
2 ac1911 103 3 <try-errr [1]> NA NA NA NA
3 ac1911 104 3 <try-errr [1]> NA NA NA NA
4 ac1911 2 14 <nls> 30.2 0.617 0.822 0.252
5 ac1911 3 9 <try-errr [1]> NA NA NA NA
6 ac1911 7 3 <try-errr [1]> NA NA NA NA
7 ac1911 8 3 <try-errr [1]> NA NA NA NA
8 ac1911 e 3 <try-errr [1]> NA NA NA NA
9 ac1911 f 3 <try-errr [1]> NA NA NA NA
10 ac1911 x 3 <try-errr [1]> NA NA NA NA
11 bf52 1 12 <try-errr [1]> NA NA NA NA
12 bf52 107 3 <try-errr [1]> NA NA NA NA
13 bf52 2 12 <nls> 30.2 0.768 -11.6 2.54
14 bf52 3 6 <nls> 30.2 0.524 -6.22 1.52
15 bf52 a 3 <try-errr [1]> NA NA NA NA
16 bf52 b 3 <try-errr [1]> NA NA NA NA
17 bf52 c 3 <try-errr [1]> NA NA NA NA
This is the same as the last section except we fixed both a and b at the value given by the overall fit. This time all groups provide a result.
fo <- avefflux ~ a * exp(-b * moisture_data_2020) + log(c + d*avtemp_data_2020)
st <- list(a = 0.8, b = 0.1, c = 0.1, d = 1)
fm <- nls(fo, d_resp_data_2020, start = st)
a <- coef(fm)[1]
b <- coef(fm)[2]
model_efflux <- d_resp_data_2020 %>%
group_by(stand, plot) %>%
mutate(n = n()) %>%
group_by(n, .add = TRUE) %>%
do(model = try(nls(fo, data = ., start = coef(fm)[3:4]))) %>%
ungroup %>%
rowwise %>%
mutate(a = if (inherits(model, "nls")) a else NA,
b = if (inherits(model, "nls")) b else NA,
c = if (inherits(model, "nls")) coef(model)[["c"]] else NA,
d = if (inherits(model, "nls")) coef(model)[["d"]] else NA) %>%
ungroup
> model_efflux
# A tibble: 17 x 8
stand plot n model a b c d
<chr> <chr> <int> <list> <dbl> <dbl> <dbl> <dbl>
1 ac1911 1 18 <nls> 30.2 0.563 0.620 0.291
2 ac1911 103 3 <nls> 30.2 0.563 0.681 0.210
3 ac1911 104 3 <nls> 30.2 0.563 -0.803 0.551
4 ac1911 2 14 <nls> 30.2 0.563 0.727 0.262
5 ac1911 3 9 <nls> 30.2 0.563 0.330 0.382
6 ac1911 7 3 <nls> 30.2 0.563 4.10 -0.396
7 ac1911 8 3 <nls> 30.2 0.563 9.87 -1.43
8 ac1911 e 3 <nls> 30.2 0.563 -28.0 6.95
9 ac1911 f 3 <nls> 30.2 0.563 -1.33 0.711
10 ac1911 x 3 <nls> 30.2 0.563 -80.8 18.0
11 bf52 1 12 <nls> 30.2 0.563 -8.21 1.90
12 bf52 107 3 <nls> 30.2 0.563 16.6 -2.24
13 bf52 2 12 <nls> 30.2 0.563 -9.17 2.04
14 bf52 3 6 <nls> 30.2 0.563 -6.24 1.55
15 bf52 a 3 <nls> 30.2 0.563 102. -13.4
16 bf52 b 3 <nls> 30.2 0.563 0.208 0.450
17 bf52 c 3 <nls> 30.2 0.563 -6.31 1.20