210 likes | 351 Views
Spatial Modeling Using R: A case study with three landscape types and a disturbance agent. NEMO 2010 Joseph Pekol Dr. David Heibeler Dr. Aaron Weiskittel. Outline. Introduction SIRS model basics Model framework Framework to R code Sample output/results Conclusion.
E N D
Spatial Modeling Using R: A case study with three landscape types and a disturbance agent NEMO 2010 Joseph Pekol Dr. David Heibeler Dr. Aaron Weiskittel
Outline • Introduction • SIRS model basics • Model framework • Framework to R code • Sample output/results • Conclusion
Using R as a spatial modeling tool • Event driven models • Accessible matrix manipulation • Quickly write, debug, and test code • Easy to produce graphical output • Built in statistical theory • No existing R packages available
SIRS basics • Susceptible-Infected-Recovered-Susceptible • Simulates a birth-death model of continuous-time spatial populations • Long-distance and local interactions between pathogen and target • Assumptions based on Poisson distribution and random occurrence on exponential dist. • Poisson used for simplicity – memoryloss property
Linking Poisson processes and simulations At time t, assume N(t) infected sites. Each infected site: • Produces offspring at rate φ • Dies off at rate μ E.g. each site acts as Poisson with rate φ+μ, so the total population acts as a Poisson process with rate λT = (φ+μ)N
Simulating an SIRS model in 3 steps • Choose when the next event will occur from an exponential distribution with mean 1/ λT. • Choose the origin of the event from all occupied sites with an equal probability. • Simulate an attempted infection of another site (pbirth) or else the death of the current site. If a birth occurs, choose a target site (specified as a percentage of local/long-distance infection)
Case study outline Landscape generator SIRS model Input 3 different sets of infection parameters Output visual landscape, residual proportions of stand types • Three landscape types • Hydric, Mesic, Xeric
However, we first need a landscape to work with…. Using a modified landscape from R script by David Hiebeler: • Three site types on a 100x100 lattice • Input site type densities and clustering level • Output includes matrix of stand types and visual representation • Total lines of code: 273
What R needs to know… [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [1,] 1 1 1 0 0 0 2 2 1 1 1 0 0 0 2 [2,] 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 [3,] 1 1 0 2 2 2 2 1 1 1 1 1 1 1 2 [4,] 1 1 1 1 2 2 2 1 1 1 1 1 1 1 2 [5,] 1 1 1 1 2 2 2 2 0 0 0 1 2 2 2 [6,] 1 1 1 1 1 0 0 0 0 0 2 2 2 2 2 [7,] 1 1 1 1 1 0 0 0 0 0 2 2 2 2 2 [8,] 1 1 1 0 0 1 1 1 1 1 2 2 2 2 2 [9,] 2 1 1 2 0 1 1 1 1 1 2 2 2 2 2 [10,] 2 1 1 2 1 1 2 1 1 1 2 2 2 2 2 [11,] 0 2 1 1 1 1 2 2 0 0 2 2 0 2 1 [12,] 0 0 0 0 1 1 2 2 0 0 0 0 0 2 1 [13,] 1 2 0 1 1 0 2 2 1 1 1 0 2 2 2 [14,] 1 1 1 1 0 0 0 1 1 1 1 1 1 2 2 [15,] 2 1 1 1 1 2 2 1 1 1 1 1 1 2 2 [16,] 2 1 1 1 1 2 2 1 1 1 1 1 1 2 2 [17,] 2 2 2 2 2 2 0 1 1 1 1 1 1 2 2 [18,] 1 1 2 2 1 1 2 0 2 0 0 1 1 0 2 [19,] 0 2 2 2 1 1 1 0 0 0 0 1 1 1 1 [20,] 0 0 2 2 1 1 1 0 0 0 0 1 1 1 1 =
SIRS model • Accepts landscape matrix • Input infection, dist. Site infection, recovery, and death rates • Randomly infects initial individuals • Outputs visual landscape; data frame of residual stand types proportions • Total lines of code: 335
SIRS model framework No – Record final site types, calculate proportion site type left, output data Update pop. vectors & Landscape matrix Infections remain? Yes – Calculate rates based on infected population. Choose where event originates and susceptible target. Long distance interaction? - Choose random S coord. Local interaction? - Choose neighbor using offset vector If target Susc. (=0,1,2) Infect target based on P(susceptible site) If target recovered. (=4) Infect target based on P(disturbed site) Infected indiv. loses infection with rate γ Site dies if runif(1) < P(Site recov.) Site recovers if runif(1) > P(Site recov.)
The key is good code for efficient performance • Keeping track of population totals • Choosing sites • Keeping track of infected locations • Updating population indices
Keeping track of everything Initialize population vectors chucksize = 500 # keeps track of how long are state vectors currently are vectorlengths = chunksize # Create vectors S = rep(0,chunksize) I = rep(0,chunksize) R = rep(0,chunksize) D = rep(0,chunksize) et = rep(0,chunksize) …….. if (i == vectorlengths) { vectorlengths = vectorlengths + chunksize length(S) = vectorlengths length(I) = vectorlengths ….. } Sets vector length, improves efficiancy so R doesn’t have to ‘grow as it goes’. Increase vector lengths as needed…
Choosing event locations Choose ‘origin’ site Vector containing neighborhood offsets: nhoodOffsets=matrix(c(0,-1,0,1,-1,0,1,0),nrow=4,byrow=TRUE) xoffsets = nhoodOffsets[,1] yoffsets = nhoodOffsets[,2] coordInd = floor(runif(1)*currentI)+1 x = Ixv[coordInd] y = Iyv[coordInd] X and Y vectors hold coordinates for infected individuals. Random choice of index determines event locations Or just… xoffsets = c(0,-1,0,1) yoffsets = c(-1,0,1,0) Choose ‘interaction’ site as local or long distance if (runif(1) < alpha) { # long-distance contact otherx = floor(runif(1)*L)+1 othery = floor(runif(1)*L)+1 } else { # local contact randInd = floor(runif(1)*4) + 1 otherx = x + xoffsets[randInd] othery = y + yoffsets[randInd] otherx = ((otherx - 1 + L) %% L) + 1 othery = ((othery - 1 + L) %% L) + 1 } Wrap around code allows movement off one edge of matrix and onto the opposite edge.
Tracking event occurrences Check whether an infection attempt is successful if (stateArray[otherx,othery] == 0 || stateArray[otherx,othery] == 1 || stateArray[otherx,othery] == 2 && (runif(1)< pColSuccess[hab[x,y]+1])) { # target susceptible currentS = currentS - 1 currentI = currentI + 1 stateArray[otherx,othery] = 3 Ixv[currentI] = otherx Iyv[currentI] = othery siteTypeInf[currentI] = hab[otherx,othery] } Update current population totals Update matrix with new stand classification Update location of infection for indexing in next iteration
Results – Proportions by trial Here, only one attempt per trial… However, a simple script can run the model multiple times, producing means, standard error, etc.
Results – Output of populations over time • deSolve package • Provides functions to solve first order, ordinary differential equations (ODE) among others • Modeled populations vs estimated ODE population curves
Conclusions – Pros and Cons of event driven models in R Pros • Quick and easy data manipulation • No need to compile: writing, testing, and debugging code is simplified • Opportunity for very robust results • Excellent packages make life easier • Graphical and statistical output very accessible Cons • Much slower than standard programming languages • Increased model complexity leads to trickier programming • Efficient coding a must for quick run-times But…the biggest Pro of all… R is free!!!
Thank you! Questions/Comments?