-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfit.response.surface.data.R
176 lines (158 loc) · 4 KB
/
fit.response.surface.data.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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
# libraries required to run the code
library(bbmle)
library(dotwhisker)
library(dplyr)
library(odeintr)
# to use/run the 2sp model
source('lib/model.2sp.R')
# create the simulated dataset
# see file below for response surface "exerimental design"
source('lib/generate.response.surface.data.R')
# load the functions used for model fitting
source('lib/model.fitting.functions.R')
#################################
# example model fitting approach
#################################
# first estimate the best-fit model parameters for i only data
# starting guesses for i only parameters
params <- list(
log_r_i = log(10),
log_phi_i = log(max(subset(simulated.data, focal=="i")$fecundity))
)
# use mle2 to infer unknown values for i only data
# NOTE: unidentifiable parameters are held at fixed values
model.fit.i <- bbmle::mle2(
minuslogl = fecundity.model.NLL,
start = unlist(params),
fixed = list(
gamma_i = 0,
gamma_j = 0,
mu_i = 0,
mu_j = 0,
nu_i = 0,
nu_j = 0,
beta_i = 0.001,
beta_j = 0.001,
K_i = 1,
K_j = 1,
log_r_j = 0,
log_phi_j = 0,
alpha_ij = 0,
alpha_ji = 0
),
data=subset(simulated.data, plants.j == 0),
vecpar = TRUE,
control = list(
maxit=100000
),
skip.hessian=FALSE
)
# second estimate the best-fit model parameters for j only data
# starting guesses for j only parameters
params$log_r_j <- log(10)
params$log_phi_j = log(max(subset(simulated.data, focal=="j")$fecundity))
# use mle2 to infer unknown values for j only data
# NOTE: unidentifiable parameters are held at fixed values
model.fit.j <- bbmle::mle2(
minuslogl = fecundity.model.NLL,
start = unlist(params),
fixed = list(
gamma_i = 0,
gamma_j = 0,
mu_i = 0,
mu_j = 0,
nu_i = 0,
nu_j = 0,
beta_i = 0.001,
beta_j = 0.001,
K_i = 1,
K_j = 1,
log_r_i = 0,
log_phi_i = 0,
alpha_ij = 0,
alpha_ji = 0
),
data=subset(simulated.data, plants.i == 0),
vecpar = TRUE,
control = list(
maxit=100000
),
skip.hessian=FALSE
)
# lastly we will use the estimated parameters above as starting values when fitting ALL data
params$log_r_i <- as.numeric(coef(model.fit.i)["log_r_i"])
params$log_r_j <- as.numeric(coef(model.fit.j)["log_r_j"])
params$log_phi_i <- as.numeric(coef(model.fit.i)["log_phi_i"])
params$log_phi_j <- as.numeric(coef(model.fit.j)["log_phi_j"])
# assume no heterospecific interactions as starting guess
# (admittedly, this is likely a poor guess)
params$alpha_ij <- 0
params$alpha_ji <- 0
# use mle2 to estimate ALL best-fit model parameters and associated SEs
# NOTE: unidentifiable parameters are held at fixed values
model.fit <- bbmle::mle2(
minuslogl = fecundity.model.NLL,
start = unlist(params),
fixed = list(
gamma_i = 0,
gamma_j = 0,
mu_i = 0,
mu_j = 0,
nu_i = 0,
nu_j = 0,
K_i = 1,
K_j = 1,
beta_i = 0.001,
beta_j = 0.001
),
data=simulated.data,
vecpar = TRUE,
control = list(
maxit=100000
),
skip.hessian=FALSE
)
# print out the coefficient summary table
print(summary(model.fit))
# evaluate a model with the parameters used for the simulations
known.params <- list()
known.params$log_r_i <- log(simulated.params$r[1])
known.params$log_r_j <- log(simulated.params$r[2])
known.params$log_phi_i <- log(simulated.params$phi[1])
known.params$log_phi_j <- log(simulated.params$phi[2])
known.params$alpha_ij <- simulated.params$alpha_ij
known.params$alpha_ji <- simulated.params$alpha_ji
# use mle2 to calculate loglikelihood, but mostly to generate a coefficient table
model.known <- bbmle::mle2(
minuslogl = fecundity.model.NLL,
start = unlist(known.params),
fixed = list(
gamma_i = 0,
gamma_j = 0,
mu_i = 0,
mu_j = 0,
nu_i = 0,
nu_j = 0,
K_i = 1,
K_j = 1,
beta_i = 0.001,
beta_j = 0.001
),
data=simulated.data,
vecpar = TRUE,
eval.only=TRUE
)
# plot known parameters versus inferred parameters
dwplot(list(model.known, model.fit)) %>%
relabel_predictors(
c(
log_r_i = "ln(r_i)",
log_r_j = "ln(r_j)",
log_phi_i = "ln(phi_i)",
log_phi_j = "ln(phi_j)"
)
) +
scale_colour_discrete(
name = "Parameters",
labels = c("Inferred", "Known")
)