title: SIR model error - can't find bug, need help in locating potential source of deviation?
body: <p>This question will be an interesting one. I was trying to replicate the results of a <a href="http://www.maths.qmul.ac.uk/%7Elatora/buscarino_jump_38002_epl08.pdf" rel="nofollow noreferrer">paper</a> which concerned disease transmission in a system of freely moving agents (sounds like the perfect job for NetLogo). I coded up a simple SIR model in NetLogo pretty easily according to the details given in the paper, made sure my model parameters matched those listed, and let the simulation run. Everything ran perfectly until I checked how the experimental results matched with the predicted values (according to the results of the paper). They were off, and by a pretty sizeable margin. Thinking there was an error somewhere in the code, I triple checked everything, only to find <em>nothing</em>. I then made sure the ordering of events was correct (as the order of movement, infection, and recovery matters), and these also matched the paper. I mulled over the problem for quite some time until finally I opened R, coded up the exact same program in RStudio, and let it run, only to find that the results matched the prediction perfectly! The R code does the same thing I <em>expect</em> the NetLogo code to be doing, so I think that something is going on behind the scenes in NetLogo or I've a misunderstanding somewhere that is the source of the deviation... Note that since the result in the paper is a mean-field approximation, you would have to run the program a few times in order for it to approach the theoretical result.</p>
  <p>I'm not sure where I'm going wrong, as my R code confirms the predicted values are correct, so I conclude that something <em>somewhere</em> in my NetLogo code is incorrect. I'm not too familiar with NetLogo, and I would really appreciate it if someone could help me find where in the following code the deviation may be occurring. The experimental average tends to be lower than the predicted one, suggesting that infection occurs faster than it should, but of all the changes I looked at, none of them solved this problem (e.g. infections do not occur one at a time per infectious turtle). Any suggestions/help would be very much appreciated.</p>
  <p>A slimmed-down version of my code is presented below. This should run in a regular interface with the standard setup/go buttons. Results are stored in lists that can be plotted, and anyone curious can see the deviation as the simulation progresses via the Plot object. Thank you in advance.</p>
  <pre><code>;; Simple SIR model
  globals [
    ;; variables for storing predictions
    predS
    predE
    predI
    predR
    oldPredS
    oldPredE
    oldPredI
    oldPredR
  
    ;; list to store experimental values
    Slist
    ;; list to store predicted values
    predSList 
    
    ;; model variables
    length-of-patch ;; length of habitat (a square of area length-of-patch^2)
    infection-radius ;; the distance from an infectious individual a susceptible agent has to be within
    ;; in order to risk getting infected
    total-pop ;; total population in the model
    force-of-infection ;; probability of infection if within infection-radius distance
    I0 ;; initial infected
    recovery-rate ;; probability of recovery
  ]
  
  turtles-own [
    infected-status ;; 0 susceptible, 1 infected, 2 recovered
  ]
  
  to setup
    ca ;; clear
    
    ;; define the variables
    set length-of-patch 31.62278 ;; the square root of 1000 (so the density is 1)
    set infection-radius 1
    set total-pop 1000
    set force-of-infection 0.1
    set I0 10
    set recovery-rate 0.05
    
    ;; setup simulation
    setup-patches
    setup-agents
    reset-ticks
    
    ;; initialize lists as empty
    set Slist []
    set predSList []
  end
  
  
  to go
    ;; update experimental values (density of susceptible individuals)
    set Slist lput ((count turtles with [infected-status = 0]) / (length-of-patch ^ 2)) Slist
    
    if (ticks = 0) ;; if ticks == 0, make sure initial value is the same as experimental
    [
      ;; update predicted values with densities of agents
      set predS ((count turtles with [infected-status = 0]) / (length-of-patch ^ 2)) 
      set predI ((count turtles with [infected-status = 1]) / (length-of-patch ^ 2))  
      set predR 0
      ;; placeholder variables for iterative process
      set oldPredS predS
      set oldPredI predI
      set oldPredR predR
  
      ;; store predicted S population in corresponding list
      set predSList lput (predS) predSList
    ]
    if (ticks &gt; 0) ;; if ticks &gt; 0, then update predicted values according to paper results
    [
      ;; update predicted values
      set predI (oldPredI + oldPredS * (1 -  (1 - force-of-infection * oldPredI) ^ (pi * (infection-radius ^ 2))) - recovery-rate * oldPredI)
      set predR (oldPredR + recovery-rate * oldPredI)
      set predS ((total-pop / (length-of-patch ^ 2)) - predI - predR)
      ;; placeholder variables
      set oldPredS predS
      set oldPredI predI
      set oldPredR predR
  
      ;; store values in corresponding list
      set predSList lput (oldPredS) predSList
    ]
  
  
    ;; perform movement, infection, and recovery, in that order
    move-agents
    infect-agents
    recover-agents
  
    if (count turtles with [infected-status = 1] = 0) [
      ;; if no one else is infected, stop
      stop
    ]
   
    tick
  end
  
  
  to setup-patches
    ;; resize the world to make it fit comfortably in the interface
    resize-world 0 length-of-patch 0 length-of-patch 
    set-patch-size 400 / (length-of-patch)
  end
  
  to setup-agents
    ;; create susceptible agents
    crt (total-pop - I0) [
      set infected-status 0
      setxy random-pxcor random-pycor
      set color 55 ;; green
      set size 2
    ]
    ;; create I0 infected agents
    crt I0 [
      set infected-status 1
      setxy random-pxcor random-pycor
      set color 15 ;; red
      set size 2
    ]
  end
  
  
  to move-agents ;; move all the agents
    ask turtles [
      setxy random-pxcor random-pycor
    ]
  end
  
  to infect-agents
    ;; iterate over infected turtles
    ask turtles with [infected-status = 1] [
      ;; check neighborhood around infected turtle for susceptible turtles...
      let numNeighbors count (turtles with [infected-status = 0] in-radius infection-radius)
      
      if (numNeighbors &gt; 0) [ ;; there are susceptibles around, so we perform infection
        ask (turtles with [infected-status = 0] in-radius infection-radius) [
          let %draw (random-float 1)
          if (%draw &lt;= force-of-infection) [ ;; probability of infection
            ;; infect one of the neighbors
            set infected-status 1
            set color 15 ;; red
          ]
        ] 
      ] ;; end of if numneighbors &gt; 0
    ]
  end
  
  to recover-agents
    ask turtles with [infected-status = 1] [
      let %draw (random-float 1)
      if (%draw &lt;= recovery-rate) [ ;; an agent recovered
        set infected-status 2
        set color 105
      ]
    ]
  end
  </code></pre>