|
| 1 | +--- |
| 2 | +interact_link: notebooks/lloydjansen2004/r_odin.ipynb |
| 3 | +title: 'R using odin' |
| 4 | +permalink: 'chapters/lloydjansen2004/r_odin' |
| 5 | +previouschapter: |
| 6 | + url: chapters/lloydjansen2004/intro |
| 7 | + title: 'Deterministic SEIR' |
| 8 | +nextchapter: |
| 9 | + url: chapters/ross2010/intro |
| 10 | + title: 'SIRS dynamics in a large population of households' |
| 11 | +redirect_from: |
| 12 | + - 'chapters/lloydjansen2004/r-odin' |
| 13 | +--- |
| 14 | + |
| 15 | +## Metapopulation SEIR model with cross-coupling and/or migration |
| 16 | + |
| 17 | +*Author*: Constanze Ciavarella @ConniCia |
| 18 | + |
| 19 | +*Date*: 2018-10-02 |
| 20 | + |
| 21 | +## Requirements |
| 22 | +To use odin, we need to install the package along with its dependencies. This is done using the following commands: |
| 23 | + |
| 24 | + |
| 25 | +{:.input_area} |
| 26 | +```R |
| 27 | +if (!require("drat")) install.packages("drat") |
| 28 | +drat:::add("mrc-ide") |
| 29 | +install.packages("dde") |
| 30 | +install.packages("odin") |
| 31 | +``` |
| 32 | + |
| 33 | +We also define the following helper functions: |
| 34 | + |
| 35 | + |
| 36 | +{:.input_area} |
| 37 | +```R |
| 38 | +# beta matrix: random initialisation ----------------------------------------------------- |
| 39 | +beta.mat <- function (nr_patches) { |
| 40 | + ## beta: positive, non-symmetric matrix of effective contact rates |
| 41 | + ## values are higher on the diagonal (transmission is higher *within*-patch), |
| 42 | + ## values decrease gradually further away from the diagonal (*between*-patch transmission) |
| 43 | + beta <- matrix(0, nrow=nr_patches, ncol=nr_patches) |
| 44 | + for (i in 0:(nr_patches-1)) { |
| 45 | + ## superdiagonal: scale vector s.t. numbers decrease further away from diagonal |
| 46 | + rand.vec <- rnorm(nr_patches-i, 10 * (nr_patches-i)/nr_patches, 2) |
| 47 | + beta[row(beta) == col(beta) - i] <- rand.vec |
| 48 | + ## subdiagonal |
| 49 | + rand.vec <- rnorm(nr_patches-i, 10 * (nr_patches-i)/nr_patches, 2) |
| 50 | + beta[row(beta) == col(beta) - i] <- rand.vec |
| 51 | + } |
| 52 | + return(beta) |
| 53 | +} |
| 54 | + |
| 55 | +# mobility matrix: random initialisation ----------------------------------------------------- |
| 56 | +mob.mat <- function (nr_patches) { |
| 57 | + ## mobility matrix (origin-destination matrix of proportion of population that travels) #TODO trip counts or relative? |
| 58 | + C <- matrix(0, nrow=nr_patches, ncol=nr_patches) |
| 59 | + if (nr_patches > 1) { |
| 60 | + diag(C) <- -sample(1:25, nr_patches, replace=TRUE) |
| 61 | + for (i in 1:nr_patches) { |
| 62 | + C[i, which(C[i, ] == 0)] <- rmultinom(n = 1, size = -C[i,i], prob = rep(1, nr_patches-1)) |
| 63 | + } |
| 64 | + C <- C/100 |
| 65 | + } |
| 66 | + return(C) |
| 67 | +} |
| 68 | + |
| 69 | +# plotting --------------------------------------------------------------------# |
| 70 | +plot.pretty <- function(out, nr_patches, what) { |
| 71 | + |
| 72 | + # plot total densities by disease status --------------------------------------# |
| 73 | + if (what == "total") { |
| 74 | + ## compute total densities by disease status |
| 75 | + S_tot <- rowSums(out$S) / nr_patches |
| 76 | + E_tot <- rowSums(out$E) / nr_patches |
| 77 | + I_tot <- rowSums(out$I) / nr_patches |
| 78 | + R_tot <- rowSums(out$R) / nr_patches |
| 79 | + |
| 80 | + ## plot total densities by disease status |
| 81 | + par(mfrow=c(1, 1), las=1, omi=c(1,0,0,0), xpd=NA) |
| 82 | + plot( t, S_tot, col="green", |
| 83 | + type="l", xlab="Days", ylab="Densities", |
| 84 | + main="Total densities by disease status") |
| 85 | + lines(t, E_tot, col="orange") |
| 86 | + lines(t, I_tot, col="red") |
| 87 | + lines(t, R_tot, col="blue") |
| 88 | + legend(-8.5, -0.3, title="Disease statuses", horiz=TRUE, |
| 89 | + legend=c("Susceptible", "Exposed", "Infectious", "Recovered"), |
| 90 | + col=c("green", "orange", "red", "blue"), lty=1) |
| 91 | + } |
| 92 | + |
| 93 | + # plot densities of some patches by disease status ----------------------------# |
| 94 | + if (what == "panels") { |
| 95 | + ## define the plot panels |
| 96 | + if (nr_patches >= 6) { |
| 97 | + mfrow=c(2, 3) |
| 98 | + panels <- as.integer( seq(1, nr_patches, length.out=6)) |
| 99 | + } else if (nr_patches >= 4) { |
| 100 | + mfrow=c(2, 2) |
| 101 | + panels <- as.integer( seq(1, nr_patches, length.out=4)) |
| 102 | + } else { |
| 103 | + mfrow=c(1, nr_patches) |
| 104 | + panels <- 1:nr_patches |
| 105 | + } |
| 106 | + par(mfrow=mfrow, las=1, omi=c(1,0,0.3,0), xpd=NA) |
| 107 | + |
| 108 | + ## plot the disease statuses of some patches |
| 109 | + ymax <- max(out$N[, panels]) |
| 110 | + for (i in panels) { |
| 111 | + plot (t, out$S[, i], col="green", ylim=c(0, ymax), |
| 112 | + type="l", xlab="Days", ylab="Densities", |
| 113 | + main=paste("Patch ", i)) |
| 114 | + lines(t, out$E[, i], col="orange") |
| 115 | + lines(t, out$I[, i], col="red") |
| 116 | + lines(t, out$R[, i], col="blue") |
| 117 | + } |
| 118 | + title("Densities by disease status", outer=TRUE) |
| 119 | + legend(-130, -1.2, title="Disease statuses", horiz=TRUE, |
| 120 | + legend=c("Susceptible", "Exposed", "Infectious", "Recovered"), |
| 121 | + col=c("green", "orange", "red", "blue"), lty=1) |
| 122 | + } |
| 123 | + |
| 124 | +} |
| 125 | +``` |
| 126 | + |
| 127 | +The model is specified in odin as: |
| 128 | + |
| 129 | + |
| 130 | +{:.input_area} |
| 131 | +```R |
| 132 | +SEIR_cont <- odin::odin({ |
| 133 | + nr_patches <- user() |
| 134 | + n <- nr_patches |
| 135 | + |
| 136 | + ## Params |
| 137 | + lambda_prod[ , ] <- beta[i, j] * I[j] |
| 138 | + lambda[] <- sum(lambda_prod[i, ]) # rowSums |
| 139 | + |
| 140 | + mob_prod[ , ] <- S[i] * C[i, j] |
| 141 | + mob_S[] <- sum(mob_prod[, i]) # colSums |
| 142 | + mob_prod[ , ] <- E[i] * C[i, j] |
| 143 | + mob_E[] <- sum(mob_prod[, i]) |
| 144 | + mob_prod[ , ] <- I[i] * C[i, j] |
| 145 | + mob_I[] <- sum(mob_prod[, i]) |
| 146 | + mob_prod[ , ] <- R[i] * C[i, j] |
| 147 | + mob_R[] <- sum(mob_prod[, i]) |
| 148 | + |
| 149 | + N[] <- S[i] + E[i] + I[i] + R[i] |
| 150 | + output(N[]) <- TRUE |
| 151 | + |
| 152 | + ## Derivatives |
| 153 | + deriv(S[]) <- mu - mu*S[i] - S[i]*lambda[i] + M[1] * mob_S[i] |
| 154 | + deriv(E[]) <- S[i]*lambda[i] - (mu + sigma) * E[i] + M[2] * mob_E[i] |
| 155 | + deriv(I[]) <- sigma*E[i] - (mu + gamma)*I[i] + M[3] * mob_I[i] |
| 156 | + deriv(R[]) <- gamma*I[i] - mu*R[i] + M[4] * mob_R[i] |
| 157 | + |
| 158 | + ## Initial conditions |
| 159 | + initial(S[]) <- 1.0 - 1E-6 |
| 160 | + initial(E[]) <- 0.0 |
| 161 | + initial(I[]) <- 1E-6 |
| 162 | + initial(R[]) <- 0.0 |
| 163 | + |
| 164 | + ## parameters |
| 165 | + beta[,] <- user() # effective contact rate |
| 166 | + sigma <- 1/3 # rate of breakdown to active disease |
| 167 | + gamma <- 1/3 # rate of recovery from active disease |
| 168 | + mu <- 1/10 # background mortality |
| 169 | + C[,] <- user() # origin-destination matrix of proportion of population that travels |
| 170 | + M[] <- user() # relative migration propensity by disease status |
| 171 | + |
| 172 | + ## dimensions |
| 173 | + dim(beta) <- c(n, n) |
| 174 | + dim(C) <- c(n, n) |
| 175 | + dim(M) <- 4 |
| 176 | + dim(lambda_prod) <- c(n, n) |
| 177 | + dim(lambda) <- n |
| 178 | + dim(mob_prod) <- c(n, n) |
| 179 | + dim(mob_S) <- n |
| 180 | + dim(mob_E) <- n |
| 181 | + dim(mob_I) <- n |
| 182 | + dim(mob_R) <- n |
| 183 | + dim(S) <- n |
| 184 | + dim(E) <- n |
| 185 | + dim(I) <- n |
| 186 | + dim(R) <- n |
| 187 | + dim(N) <- n |
| 188 | + |
| 189 | +}) |
| 190 | +``` |
| 191 | + |
| 192 | +{:.output_stream} |
| 193 | +``` |
| 194 | +gcc -I/usr/local/lib/R/include -DNDEBUG -I/usr/local/include -fpic -g -O2 -c odin.c -o odin.o |
| 195 | +g++ -shared -L/usr/local/lib/R/lib -L/usr/local/lib -o odin_b0e443ee.so odin.o -L/usr/local/lib/R/lib -lR |
| 196 | +
|
| 197 | +``` |
| 198 | + |
| 199 | +## Running the model |
| 200 | + |
| 201 | +We first define the user-specified parameters. The default values provided below all satisfy model requirements, but can be changed if needed. |
| 202 | + |
| 203 | + |
| 204 | +{:.input_area} |
| 205 | +```R |
| 206 | +# set seed |
| 207 | +set.seed(1) |
| 208 | + |
| 209 | +# SEIR free parameters --------------------------------------------------------# |
| 210 | +## total number of patches in the model |
| 211 | +nr_patches = 20 |
| 212 | +## relative migration propensity by disease status (S, E, I, R) |
| 213 | +M <- c(1, 0.5, 1, 1) |
| 214 | +## matrix of effective contact rates |
| 215 | +beta <- beta.mat(nr_patches) |
| 216 | +## mobility matrix |
| 217 | +C <- mob.mat(nr_patches) |
| 218 | +``` |
| 219 | + |
| 220 | +Now we run the model. Note that each patch is initialised with the same population size. We also run a test to make sure that the total population size has remained constant throughout the simulation. |
| 221 | + |
| 222 | + |
| 223 | +{:.input_area} |
| 224 | +```R |
| 225 | +# run SEIR model --------------------------------------------------------------# |
| 226 | +mod <- SEIR_cont(nr_patches=nr_patches, beta=beta, C=C, M=M) |
| 227 | +t <- seq(0, 50, length.out=50000) |
| 228 | +out <- mod$run(t) |
| 229 | +out <- mod$transform_variables(out) |
| 230 | +# error check -----------------------------------------------------------------# |
| 231 | +if ( ! all( abs(rowSums(out$N) - nr_patches) < 1E-10 ) ) |
| 232 | + warning("Something went wrong, density is increasing/decreasing!\n") |
| 233 | +``` |
| 234 | + |
| 235 | +Now we can plot the disease dynamics in the total population. |
| 236 | + |
| 237 | + |
| 238 | +{:.input_area} |
| 239 | +```R |
| 240 | +# plot total densities by disease status --------------------------------------# |
| 241 | +plot.pretty(out, nr_patches, "total") |
| 242 | +``` |
| 243 | + |
| 244 | + |
| 245 | + |
| 246 | + |
| 247 | + |
| 248 | +We can also plot the disease dynamics of some of the patches. Dynmics between patches may vary due to cross-coupling and migration rates. |
| 249 | + |
| 250 | + |
| 251 | +{:.input_area} |
| 252 | +```R |
| 253 | +# plot densities of some patches by disease status ----------------------------# |
| 254 | +plot.pretty(out, nr_patches, "panels") |
| 255 | +``` |
| 256 | + |
| 257 | + |
| 258 | + |
| 259 | + |
0 commit comments