/
Validation_model_markdown_13June2019.Rmd
196 lines (158 loc) · 6.59 KB
/
Validation_model_markdown_13June2019.Rmd
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
---
title: "Validating against ODE"
---
## Validation of Gillespie algorithm and disease transmission code in C
1. Comparing results from pure Stochastic Gillespie Algorithm (pure SSA) i.e. all disease status transition occurs probalistically with ignoring the history of their status (e.g. how long being already being infectious doesn ot affect the timing of an individual to become immune) to ODE
2. And comparing results from modified Gillespie Algorithm (modified SSA) i.e. disease status transition occurs considering their disease status history
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Read simulated output from C
```{r}
#=========================================================================================
# make sure the path is correct: you need to put simulated output from modified_SSA.c in "E:/tem/analysis" and outout from SSA.c in "E:/tem/analysis/SSA", or alternatively anywhere that suits you. If latter, make sure you change
# location variable to a proper name
# Read data for modified SSA
location = "E:/tem/analysis"
setwd(location)
lf_r<-list.files(path=location,full.names=F,pattern="csv")
data_list<-lapply(lf_r,function(i){
read.table(i,header = FALSE, sep = ",")
})
for(i in seq(data_list))
{
num<-gsub("\\D","",substring(lf_r[i],21,22))
assign(paste0("df",num),data_list[[i]])
}
# Read SSA data
location = "E:/tem/analysis/SSA"
setwd(location)
lf_ssa<-list.files(path=location,full.names=F,pattern="csv")
data_ssa_list<-lapply(lf_ssa,function(i){
read.table(i,header = FALSE, sep = ",")
})
for(i in seq(data_ssa_list))
{
num<-gsub("\\D","",substring(lf_ssa[i],24,25))
assign(paste0("df_ssa",num),data_ssa_list[[i]])
}
```
## Preprare Plot for 100 simulated results from C code both pure SSA and modified SSA
```{r}
#=========================================================================================
library(ggplot2)
p_s <-p_e<-p_i<-p_r<- ggplot()
for(i in seq(data_list))
{
tem<-get(paste0("df",i-1))
tem <- data.frame(t(tem))
tem <- tem[tem$X1>0 & !is.na(tem$X1),]
colnames(tem)[1] <- "TIME"
colnames(tem)[2] <- "S"
colnames(tem)[3] <- "E"
colnames(tem)[4] <- "I"
colnames(tem)[5] <- "R"
tem <- tem[tem$TIME<=100,]
loop_input <- "geom_line(data=tem,aes(x = TIME, y = S),color='lightblue')"
p_s <- p_s + eval(parse(text=loop_input))
loop_input <- "geom_line(data=tem,aes(x = TIME, y = E),color='lightblue')"
p_e <- p_e + eval(parse(text=loop_input))
loop_input <- "geom_line(data=tem,aes(x = TIME, y = I),color='lightblue')"
p_i <- p_i + eval(parse(text=loop_input))
loop_input <- "geom_line(data=tem,aes(x = TIME, y = R),color='lightblue')"
p_r <- p_r + eval(parse(text=loop_input))
}
#Do same for SSA
p_s_ssa <-p_e_ssa <-p_i_ssa <-p_r_ssa <- ggplot()
for(i in seq(data_ssa_list))
{
tem<-get(paste0("df_ssa",i-1))
tem <- data.frame(t(tem))
tem <- tem[tem$X1>0 & !is.na(tem$X1),]
colnames(tem)[1] <- "TIME"
colnames(tem)[2] <- "S"
colnames(tem)[3] <- "E"
colnames(tem)[4] <- "I"
colnames(tem)[5] <- "R"
tem <- tem[tem$TIME<=100,]
loop_input <- "geom_line(data=tem,aes(x = TIME, y = S),color='lightblue')"
p_s_ssa <- p_s_ssa + eval(parse(text=loop_input))
loop_input <- "geom_line(data=tem,aes(x = TIME, y = E),color='lightblue')"
p_e_ssa <- p_e_ssa + eval(parse(text=loop_input))
loop_input <- "geom_line(data=tem,aes(x = TIME, y = I),color='lightblue')"
p_i_ssa <- p_i_ssa + eval(parse(text=loop_input))
loop_input <- "geom_line(data=tem,aes(x = TIME, y = R),color='lightblue')"
p_r_ssa <- p_r_ssa + eval(parse(text=loop_input))
}
```
## ODE for a within-farm transmission
```{r}
#===========================================================================
library(deSolve)
#setup disease dynamic function
SEIR.dyn <- function(t,var,par){
S <- var[1]
E <- var[2]
I <- var[3]
R <- var[4]
N <- S + E + I + R
beta <- par[1] #S->E
gamma <- par[2] # E->I
delta <- par[3] #I->R
# Derivatives
dS <- -beta*S*I/N
dE <- beta*S*I/N - gamma*E
dI <- gamma*E - delta*I
dR <- delta*I
# Return the derivatives
list(c(dS,dE,dI,dR))
}
# Set up parameter
beta <- 0.5
gamma <- 1/25 #expected value for E->I in IBM
delta <- 1/50 #expected value for I->R in IBM
SEIR.par <- c(beta,gamma,delta)
SEIR.init <- c(499,0,1,0) #one infectious
SEIR.t <- seq(0,100,by=0.1)
# Run ODE
SEIR.sol <- lsoda(SEIR.init,SEIR.t, SEIR.dyn, SEIR.par)
#Rename results
SEIR.sol <-data.frame(SEIR.sol)
colnames(SEIR.sol)[1] <- "TIME"
colnames(SEIR.sol)[2] <-"S"
colnames(SEIR.sol)[3] <-"E"
colnames(SEIR.sol)[4] <-"I"
colnames(SEIR.sol)[5] <-"R"
```
## Plot and compare results from C (pure SSA and modified SSA) and R (ODE)
```{r}
#========================================================================
## Pure SSA
# Susceptible: blue(stochastic output from C) and black(ODE from R)
p_s_ssa_combined <- p_s_ssa + geom_line(data=SEIR.sol,aes(x=TIME,y=S),color="black",size=1)
# Exposed: blue(stochastic output from C) and black(ODE from R)
p_e_ssa_combined <- p_e_ssa + geom_line(data=SEIR.sol,aes(x=TIME,y=E),color="black",size=1)
# Infectious: blue(stochastic output from C) and black(ODE from R)
p_i_ssa_combined <- p_i_ssa + geom_line(data=SEIR.sol,aes(x=TIME,y=I),color="black",size=1)
# Immune: blue(stochastic output from C) and black(ODE from R)
p_r_ssa_combined <- p_r_ssa + geom_line(data=SEIR.sol,aes(x=TIME,y=R),color="black",size=1)
p_s_ssa_combined
p_e_ssa_combined
p_i_ssa_combined
p_r_ssa_combined
## Stochastic Gillespie results are surrounding the result from R ODE system
## Modified SSA
# Susceptible: blue(stochastic output from C) and black(ODE from R)
p_s_combined <- p_s + geom_line(data=SEIR.sol,aes(x=TIME,y=S),color="black",size=1)
# Exposed: blue(stochastic output from C) and black(ODE from R)
p_e_combined <- p_e + geom_line(data=SEIR.sol,aes(x=TIME,y=E),color="black",size=1)
# Infectious: blue(stochastic output from C) and black(ODE from R)
p_i_combined <- p_i + geom_line(data=SEIR.sol,aes(x=TIME,y=I),color="black",size=1)
# Immune: blue(stochastic output from C) and black(ODE from R)
p_r_combined <- p_r + geom_line(data=SEIR.sol,aes(x=TIME,y=R),color="black",size=1)
p_s_combined
p_e_combined
p_i_combined
p_r_combined
## Modified Stochastic Gillespie results are similar to R ODE system but not identical, especially for R group. This is because modified algorithm takes into account the duration which an individual spent in each disease status, which is more realistic than assuming that disease status change happens anytime (i.e. the probability of an individual moving to the next disease status is the same for an individual who spent 100 days in the status and an individual who spent only one day in the status)
```