/
Basic_model_fitting_procedure.R
262 lines (214 loc) · 7.92 KB
/
Basic_model_fitting_procedure.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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
#==================================================================
# FULL MODEL
#==================================================================
# Colors
cols <- c("#000000","#999999", "#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00")
# Packages needed
library(smfsb)
# Set working directory
setwd("~/Documents")
# Define functions for fitting ODEs
source(file = "Part3_equations_new_dates_130221.R")
# Read in parameter values and ranges
params <- read.csv("params_enviro_120221.csv",header=F)
# Read in experimental data points to fit to
data <- read.csv("experimental_data_full_pens_WB_180520.csv",header=T)
# Correct days in dataset
data$day[1] <- 0
data$day[2] <- 14
data$day[3] <- 36
data$day[4] <- 84
data$day[5] <- 126
data$day[6] <- 0
data$day[7] <- 14
data$day[8] <- 36
data$day[9] <- 84
data$day[10] <- 126
# read in temp data
temps <- read.csv("Weather_data_cleaned_150921_WB.csv",header=T)
temps <- as.numeric(temps[(1:length(temps[,1])),1])
temps <- (temps-32)*5/9
length(temps)
length(temps[which(temps<(4))])
data_empty <- read.csv("experimental_data_enviro_WB_180520.csv",header=T)
data_empty <- data_empty[c(1:5),]
# Define starting conditions
NE0 <- 10^(data$ML_generic[1])
RE0 <- 10^(data$ML_resistant[1])
SE0 <- NE0-RE0
NC0 <- 10^(data$ML_generic[6])
RC0 <- 10^(data$ML_resistant[6])
SC0 <- NC0-RC0
NE0_empty <- 10^data_empty$ML_generic[1]
RE0_empty <- 10^(data_empty$ML_resistant[1])
SE0_empty <- NE0_empty-RE0_empty
# Define experimental data points
t1_mean <- data$ML_generic[1]
t2_mean <- data$ML_generic[2]
t3_mean <- data$ML_generic[3]
t4_mean <- data$ML_generic[4]
t5_mean <- data$ML_generic[5]
t1_mean_cat <- data$ML_generic[6]
t2_mean_cat <- data$ML_generic[7]
t3_mean_cat <- data$ML_generic[8]
t4_mean_cat <- data$ML_generic[9]
t5_mean_cat <- data$ML_generic[10]
t1_Median_r <- data$ML_resistant[1]
t2_Median_r <- data$ML_resistant[2]
t3_Median_r <- data$ML_resistant[3]
t4_Median_r <- data$ML_resistant[4]
t5_Median_r <- data$ML_resistant[5]
t1_Median_r_cat <- data$ML_resistant[6]
t2_Median_r_cat <- data$ML_resistant[7]
t3_Median_r_cat <- data$ML_resistant[8]
t4_Median_r_cat <- data$ML_resistant[9]
t5_Median_r_cat <- data$ML_resistant[10]
t1_mean_empty <- data_empty$ML_generic[1]
t2_mean_empty <- data_empty$ML_generic[2]
t3_mean_empty <- data_empty$ML_generic[3]
t4_mean_empty <- data_empty$ML_generic[4]
t5_mean_empty <- data_empty$ML_generic[5]
t1_Median_r_empty <- data_empty$ML_resistant[1]
t2_Median_r_empty <- data_empty$ML_resistant[2]
t3_Median_r_empty <- data_empty$ML_resistant[3]
t4_Median_r_empty <- data_empty$ML_resistant[4]
t5_Median_r_empty <- data_empty$ML_resistant[5]
# define time-points
t1 <- data$day[1]
t2 <- data$day[2]
t3 <- data$day[3]
t4 <- data$day[4]
t5 <- data$day[5]
time_max <- (max(data$day))*24 # model runtime in days
param_names <- c("f", "Ke", "mue", "beta_rs","GSE","g","Kc","muc","s","c","Ke_empty","gse_empty","mue_empty","a")
# Define fixed parameters
GSE_min <-params[1,3] # baseline growth rate in environment
GSE_max <- params[2,3]
f_min <- params[3,3] # fitness cost
f_max <- params[4,3]
Ke_min <- params[5,3] # carrying capacity in the environment
Ke_max <- params[6,3]
mue_min <- params[7,3] # death rate in the environment
mue_max <- params[8,3]
beta_rs_min <- params[9,3] # horizantal transmission rate
beta_rs_max <- params[10,3]
g_min <- params[15,3]
g_max <- params[16,3]
Kc_min <- params[17,3]
Kc_max <- params[18,3]
muc_min <- params[19,3]
muc_max <- params[20,3]
s_min <- params[21,3]
s_max <- params[22,3]
c_min <- params[23,3]
c_max <- params[24,3]
I_min <- params[31,3]
I_max <- params[32,3]
n_min <- params[36,3] #remove carrying capacity change over time
n_max <- params[37,3]
rmodel <- function(th){
# parameter values
f <- f_min+th[1]*(f_max-f_min)
Ke <- Ke_min+th[2]*(Ke_max-Ke_min)
mue <- mue_min+th[3]*(mue_max-mue_min)
beta_rs <- beta_rs_min+th[4]*(beta_rs_max-beta_rs_min)
gse <- GSE_min+th[5]*(GSE_max-GSE_min)
g <- g_min+th[6]*(g_max-g_min)
Kc <- Kc_min+th[7]*(Kc_max-Kc_min)
muc <- muc_min+th[8]*(muc_max-muc_min)
s <- s_min+th[9]*(s_max-s_min)
c <- c_min+th[10]*(c_max-c_min)
#I <- I_min+th[11]*(I_max-I_min)
#n <- n_min+th[12]*(n_max-n_min)
#n <- th[12]*0
n<- 0
temperatures1 <- temps[1:(time_max+1)]
#gse_t <- temperatures1*I + gse
gse_t <- temperatures1*0 + gse # remove temperature from fitting procedure
gse_t[which(temperatures1<4)] <- 0
# run model
resist_full.t <- seq(0,time_max,1)
resist_full.init <- c(SE0,RE0,SC0,RC0)
pars <- c(f,beta_rs,g,Kc,muc,s,c,Ke,mue,n,gse_t)
resist_full.sol <<- resist.dyn(resist_full.t, resist_full.init, pars)
# label the output
time <- resist_full.sol[,1]
SE <- (resist_full.sol[,2]) # Take log10 of output
RE <- (resist_full.sol[,3])
NE <<- log10(SE+RE)
RE <<- log10(RE)
SC <- (resist_full.sol[,4]) # Take log10 of output
RC <- (resist_full.sol[,5])
NC <<- log10(SC+RC)
RC <<- log10(RC)
#n <- th[12]*0
s <- 0
c <- 0
Ke_empty <- 1.5+th[11]*(2.5-1.5) # change to having new parameters for empty pen: carrying capacity 1.5 to 2.5, baseline growth rate and mortality rate (same ranges)
gse_empty <- GSE_min+th[12]*(GSE_max-GSE_min)
mue_empty <- mue_min+th[13]*(mue_max-mue_min)
#gse_t <- temperatures1*I + gse_empty
gse_t <- temperatures1*0 + gse_empty
gse_t[which(temperatures1<4)] <- 0
#Ke_empty <- 2
#mue_empty <- mue
# run model
resist_full.t <- seq(0,time_max,1)
resist_full.init <- c(SE0_empty,RE0_empty,0,0)
pars <- c(f,beta_rs,g,Kc,muc,s,c,Ke_empty,mue_empty,n,gse_t)
resist_full_empty.sol <<- resist.dyn(resist_full.t, resist_full.init, pars)
SE_empty <- (resist_full_empty.sol[,2]) # Take log10 of output
RE_empty <- (resist_full_empty.sol[,3])
NE_empty <<- log10(SE_empty+RE_empty)
RE_empty <<- log10(RE_empty)
}
# Difference function - use model function
rdist <- function(th){
rmodel(th)
sum(sqrt((RE[(t2*24+1)]-t2_Median_r)^2)+
sqrt((RE[(t3*24+1)]-t3_Median_r)^2)+
sqrt((RE[(t4*24+1)]-t4_Median_r)^2)+
sqrt((RE[(t5*24+1)]-t5_Median_r)^2)+
sqrt((NE[(t2*24+1)]-t2_mean)^2)+
sqrt((NE[(t3*24+1)]-t3_mean)^2)+
sqrt((NE[(t4*24+1)]-t4_mean)^2)+
sqrt((NE[(t5*24+1)]-t5_mean)^2)+
sqrt((RC[(t2*24+1)]-t2_Median_r_cat)^2)+
sqrt((RC[(t3*24+1)]-t3_Median_r_cat)^2)+
sqrt((RC[(t4*24+1)]-t4_Median_r_cat)^2)+
sqrt((RC[(t5*24+1)]-t5_Median_r_cat)^2)+
sqrt((NC[(t2*24+1)]-t2_mean_cat)^2)+
sqrt((NC[(t3*24+1)]-t3_mean_cat)^2)+
sqrt((NC[(t4*24+1)]-t4_mean_cat)^2)+
sqrt((NC[(t5*24+1)]-t5_mean_cat)^2)+
sqrt((RE_empty[(t2*24+1)]-t2_Median_r_empty)^2)+
sqrt((RE_empty[(t3*24+1)]-t3_Median_r_empty)^2)+
sqrt((RE_empty[(t4*24+1)]-t4_Median_r_empty)^2)+
sqrt((RE_empty[(t5*24+1)]-t5_Median_r_empty)^2)+
sqrt((NE_empty[(t2*24+1)]-t2_mean_empty)^2)+
sqrt((NE_empty[(t3*24+1)]-t3_mean_empty)^2)+
sqrt((NE_empty[(t4*24+1)]-t4_mean_empty)^2)+
sqrt((NE_empty[(t5*24+1)]-t5_mean_empty)^2),na.rm=T)
}
rprior <- function() {
# parameter values
runif(13,0,1)
}
dprior <- function(x, ...) { dunif(x[1], 0,1, ...) +
dunif(x[2], 0,1, ...) + dunif(x[3], 0,1, ...) + dunif(x[4], 0,1, ...)+ dunif(x[5], 0,1, ...)+ dunif(x[6], 0,1, ...)+ dunif(x[7], 0,1, ...)+ dunif(x[8], 0,1, ...)+ dunif(x[9], 0,1, ...)+ dunif(x[10], 0,1, ...)+ dunif(x[11], 0,1, ...)+ dunif(x[12], 0,1, ...)+ dunif(x[13], 0,1, ...)}
rperturb <- function(th){th + rnorm(13, 0, 0.005)}
dperturb <- function(thNew, thOld, ...){sum(dnorm(thNew, thOld, 0.005, ...))}
# Model settings:
total_iterations <- 25000# old dates: 1.5 minutes for 100 iterations; took 49222 seconds (13.6 hours) for 20000; new dates - 30872s for 10000
start <-proc.time()
out = abcSmc(total_iterations, rprior, dprior, rdist, rperturb,
dperturb, verb=T, steps=10, factor=10)
a <- numeric(total_iterations)
for(i in 1:total_iterations){
a[i] <- rdist(out[i,])
}
out <- cbind(out,a)
colnames(out) <- c(param_names)
# save output
#write.csv(out,paste("output_basic.csv"))
#print(proc.time()-start)