-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathtest.R
159 lines (135 loc) · 7.93 KB
/
test.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
# load in packages
library(purrr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(nls.multstart)
library(broom)
library(rTPC)
library(MuMIn)
# write function to label ggplot2 panels
label_facets_num <- function(string){
len <- length(string)
string = paste('(', 1:len, ') ', string, sep = '')
return(string)
}
# write function to convert label text size to points
pts <- function(x){
as.numeric(grid::convertX(grid::unit(x, 'points'), 'mm'))
}
# load in data
data("Chlorella_TRC")
# change rate to be non-log transformed
d <- mutate(Chlorella_TRC, rate = exp(ln.rate))
chlorella_tpc <- select(d, curve_id, growth_temp = growth.temp, process, flux, temp, rate)
usethis::use_data(chlorella_tpc, rTPC)
# filter for just the first curve
d_1 <- filter(d, curve_id == 1)
# plot
ggplot(d_1, aes(temp, rate)) +
geom_point() +
theme_bw()
get_start_vals(d_1$temp, d_1$rate, model_name = 'sharpeschoolhigh_1981')
get_start_vals(d_1$temp, d_1$rate, model_name = 'sharpeschoolfull_1981')
get_start_vals(d_1$temp, d_1$rate, model_name = 'sharpeschoollow_1981')
get_start_vals(d_1$temp, d_1$rate, model_name = 'briere2_1999')
get_start_vals(d_1$temp, d_1$rate, model_name = 'thomas_2012')
get_start_vals(d_1$temp, d_1$rate, model_name = 'lactin2_1995')
get_start_vals(d_1$temp, d_1$rate, model_name = 'quadratic_2008')
get_start_vals(d_1$temp, d_1$rate, model_name = 'ratkowsky_1983')
get_start_vals(d_1$temp, d_1$rate, model_name = 'gaussian_1987')
get_start_vals(d_1$temp, d_1$rate, model_name = 'rezende_2019')
get_start_vals(d_1$K, d_1$rate, model_name = 'delong_2017')
get_lower_lims(d_1$temp, d_1$rate, model_name = 'rezende_2019')
get_upper_lims(d_1$temp, d_1$rate, model_name = 'rezende_2019')
get_lower_lims(d_1$temp, d_1$rate, model_name = 'sharpeschoolhigh_1981')
get_upper_lims(d_1$temp, d_1$rate, model_name = 'sharpeschoollow_1981')
mod <- nls.multstart::nls_multstart(rate ~ hinshelwood_1947(temp = temp, a, e, b, eh),
data = d_1,
iter = 500,
start_lower = get_start_vals(d_1$temp, d_1$rate, model_name = 'hinshelwood_1947') * 0.8,
start_upper = get_start_vals(d_1$temp, d_1$rate, model_name = 'hinshelwood_1947') * 1.1,
supp_errors = 'Y')
mod <- nls.multstart::nls_multstart(rate ~ delong_2017(temp_K = K, c, eb, ef, tm, ehc),
data = d_1,
iter = 500,
start_lower = get_start_vals(d_1$K,d_1$rate, model_name = 'delong_2017') -1,
start_upper = get_start_vals(d_1$K,d_1$rate, model_name = 'delong_2017') +1,
supp_errors = 'Y')
mod <- nls.multstart::nls_multstart(rate ~ thomas_2017(temp = temp, a, b, c, d, e),
data = d_1,
iter = 500,
start_lower = get_start_vals(d_1$temp,d_1$rate, model_name = 'thomas_2017') -1,
start_upper = get_start_vals(d_1$temp,d_1$rate, model_name = 'thomas_2017') +1,
supp_errors = 'Y')
y = thomas_2017(d_1$temp, a = 1.174, b = 0.64, c = 1.119, d = 0.267, e = 0.103)
plot(y ~ d_1$temp)
preds <- augment(mod)
ggplot(preds) +
geom_point(aes(temp, rate)) +
geom_line(aes(temp, .fitted))
number_of_models <- 5
d_models <- group_by(d_1, curve_id, growth.temp, process, flux) %>%
nest() %>%
mutate(., sharpeschoolhigh = map(data, ~nls_multstart(rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 20),
data = .x,
iter = 500,
start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') - 10,
start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') + 10,
supp_errors = 'Y')),
sharpeschoolfull = map(data, ~nls_multstart(rate ~ sharpeschoolfull_1981(temp = temp, r_tref, e, el, tl, eh, th, tref = 20),
data = .x,
iter = 500,
start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolfull_1981') - 10,
start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolfull_1981') + 10,
supp_errors = 'Y',
lower = c(r_tref = 0, e = 0, el = 0, tl = 0, eh = 0, th = 0),
upper = c(r_tref = 10, e = 10, el = 10, tl = 293.15, eh = 10, th = 400))),
briere2 = map(data, ~nls_multstart(rate ~ briere2_1999(temp = temp, tmin, tmax, a, b),
data = .x,
iter = 500,
start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'briere2_1999') - 1,
start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'briere2_1999') + 2,
supp_errors = 'Y',
lower = c(tmin = -10, tmax = 20, a = -10, b = -10),
upper = c(tmin = 20, tmax = 80, a = 10, b = 10))),
thomas = map(data, ~nls_multstart(rate ~ thomas_2012(temp = temp, a, b, c, topt),
data = .x,
iter = 500,
start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'thomas_2012') - 1,
start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'thomas_2012') + 1,
supp_errors = 'Y',
lower = c(a= 0, b = -10, c = 0, topt = 0))),
rezende = map(data, ~nls_multstart(rate ~ rezende_2019(temp = temp, q10, a, b, c),
data = .x,
iter = 500,
start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'rezende_2019') -1,
start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'rezende_2019') +1,
upper = get_upper_lims(.x$temp, .x$rate, model_name = 'rezende_2019'),
lower = get_lower_lims(.x$temp, .x$rate, model_name = 'rezende_2019'),
supp_errors = 'Y')))
d_stack <- gather(d_models, 'model', 'output', 6:ncol(d_models))
newdata <- tibble(temp = seq(min(d_1$temp), max(d_1$temp), length.out = 100),
K = seq(min(d_1$K), max(d_1$K), length.out = 100))
d_preds <- d_stack %>%
mutate(., pred = map(output, augment, newdata = newdata)) %>%
unnest(., pred)
extra_params <- d_stack %>%
mutate(., est = map(output, est_params)) %>%
select(., -c(data, output)) %>%
unnest(., est) %>%
mutate(., rmax = round(rmax, 2))
ggplot(filter(d_preds))
# plot
ggplot(d_preds, aes(temp, rate)) +
geom_point(aes(temp, rate), d_1) +
geom_text(aes(-Inf, Inf, label = paste('Topt =', topt, 'ºC\n', 'rmax = ', rmax)), hjust = -0.1, vjust = 1.5, extra_params, size = pts(10)) +
geom_line(aes(temp, .fitted, col = model)) +
facet_wrap(~model, labeller = labeller(model = label_facets_num)) +
theme_bw(base_size = 16) +
theme(legend.position = 'none',
strip.text = element_text(hjust = 0),
strip.background = element_blank()) +
xlab('Temperature (ºC)') +
ylab('rate') +
geom_hline(aes(yintercept = 0), linetype = 2)