# Statistical method for modelling of disease incidence for NGO

I am trying to provide help to an NGO in the non-profit sector that is running a disease screening program:

The program visits thousands of villages a year. A village has a population (on average around 800), with an unknown number of infected individuals. Each village is visited by the NGO for mass screening.

For a given population in a given year:

• An infected person will attend the NGO screening and get cured with a probability of S (e.g. 50%)
• An infected person might independently self-diagnose himself, seek treatment and get cured in a health clinic outside the screening program with probability P (e.g. 10%)
• An infected person will infect another person per year with probability R (e.g. 50%)
• An infected person can live for 4 years.
• There might a fixed probability that a new infection will jump over from another village (J)

These probabilities (S, P, R) will differ per village (within a estimated prior distribution), but are constant per village over time. The NGO has data going back 7 years for each village.

The NGO would like to predict the number of infected people in the population for each village in year 8.

I am trying to help the NGO, but have struggled to apply statistical learning to this problem, due to the stochastic nature and path-dependency. Most of my approaches try to fit a max. likelihood estimate to a forward projection of a starting population variable in year 0 with sampling of parameters from the prior distributions. I understand we are limited by the many hidden parameters and very limited data per village, but I feel there must be at least some predictive value in the data, but struggle to find the right approach.

Hence – I was wondering whether someone can give their perspectives on what is the best approach / are the best tools to perform predictive modelling on such a problem, and help me understand which might be possible in terms of predictive accuracy?

• One new thought - i wonder whether you could model this as different priors (for each parameter and for the original level of incidence in year 1), and subsequently update each prior at the same time in some way by using the 7 likelihood functions of observed data Commented May 13, 2014 at 13:46
• If you have trouble writing it down on one line, I would suggests simulations. Simulate the next 8 years thousands of times, and record what the distribution of the infected population is Commented May 15, 2014 at 23:47
• Are you familiar with compartmental models used in edpidemiology (en.wikipedia.org/wiki/Compartmental_models_in_epidemiology)? Your description sounds similar to SIR models (Susceptible, Infected, Recovered) albeit your data sounds more complex. Commented May 16, 2014 at 15:01

Your problem sounds like a pretty classic example of a dynamical system. There is actually quite a lot of literature, mainly aimed at engineering and the physical sciences, about techniques for dynamical systems modeling, but unfortunately considerably less of it (still more than zero though!) aimed at those working in epidemiology. If you use search terms like "dynamical systems modeling epidemiology" on either google or amazon, you'll definitely come up with some useful hits, e.g., items such as B. F. Finkenstadt and B. T. Grenfell, "Time Series Modeling of Childhood Diseases: A Dynamical Systems Approach", Journal of the Royal Statistical Society: Series C, Volume 49, Issue 2, pg. 187-205 (2000), or Frank C. Hoppensteadt, "Mathematical Methods for Analysis of a Complex Disease, American Mathematical Society, 2011.

I haven't read either of those works, nor have I ever worked in epidemiology, but I have modeled dynamical systems in a science & engineering context, and here's how I would approach constructing a mathematical model to solve your problem.

First, I'd start by borrowing three high level concepts that are frequently applied in other dynamical systems modeling-related contexts: 1.) model parameters, 2.) state space, and 3.) observation states. Let's talk about the latter two first, and then we can talk about obtaining actual data-driven estimates for the model parameters at the end.

There is a pretty good diagram here which roughly shows the relationship between state space and observation space in a dynamical systems problem. In a classic type of problem, we would typically use the concept of "state space" to mean "the true state of the underlying physical system, which in general isn't directly observable", while the "observation state" refers to whatever we actually are able to see using some type of measurement or test process. To make this more concrete, let's suppose that I am tracking an airplane, using some kind of automated video camera system. In that problem, the "system state" (a vector within the "state space") would be "the position and velocity of the airplane at a given moment in time", while the "observation state" would be something like, "the pixel on the video camera focal plane array where the airplane was detected at a given moment in time". The reason we make such a careful distinction between "state space" vs. "observation space" is that often what we are able to measure (i.e., the "observation"; the crucial real-world data that keeps our model well-grounded) is quite distinct and different from the underlying dynamical, time-evolving, hidden variable or variables that we would actually like to know. Furthermore, the "state space" is often the place where it is easiest to model the "physics" (in the case of the airplane, processes such as gravity, thrust, lift, drag, etc.) of the system correctly.

So what would the "state space" be for an epidemiological problem like this? In the simplest case, it would just be a vector with two components: the number of infected individuals in a given population (for example, the population of one of the villages, let's say) vs. the number of uninfected individuals. If you had a more complicated case, such as a disease which affects men and women in different proportions, or adults vs. children in different proportions, you might need a different state space model, e.g., a four component vector: (number of infected men, number of healthy men, number of infected women, number of healthy women), etc.

A crucial part of your dynamical system model will be figuring out the correct equations to use for your state transition function. In a typical dynamical systems problem, this is a function which tells you, given the state of the system at time $t_{n}$, what will be the expected state at time $t_{n+1}$? In the airplane tracking example that I talked about previously, this would literally involve a bunch of physics: take the current position and velocity, include the effects of gravity, thrust, lift and drag, and use Newtonian kinematics to predict the new position of the plane a second or two later.

In your problem, instead of using physics to update your system state (i.e., the relative populations of infected vs. uninfected people) what you will instead want to do is update those estimates using your model parameters; i.e., the variables which you have already labeled as $S$, $P$, $R$ and $J$, etc. Don't worry just yet that you don't actually know the proper values yet for $S$, $P$, $R$ and $J$, instead, just write a function (for example, in R or Python or whatever data analysis language you may be using) which would allow you to predict the infected vs. uninfected populations in year $(n+1)$, provided you knew them for year $n$, given any hypothetical combination of values for $S$, $P$, $R$ and $J$. In the diagram that I mentioned previously, the state transition or update function would correspond to the horizontal, rightward-pointing arrows moving from one time step to the next.

Now let's talk about the observation state. Just as with the airplane tracking problem, where what we really "observed" were not the actual position and velocity of the airplane itself, but rather, successive frames in a video camera, you also technically don't really observe your state space directly either. Instead, what you observe, literally, are the outcomes of a screening test. At first glance, these may seem like the same things to you, but in fact, the screening test will have some non-zero rate of both false positives and false negatives, which are closely related to the concept of type I and type II errors in statistics, or alternatively, to sensitivity vs. specificity in medicine and biology. The point is, you are going to need to define a second type of transition function (technically, it's usually called an "observation matrix" or "observation function"), similar to the one that we already talked about above for transitions between successive moments in time, in order to go back and forth between the state space (i.e., the true number of infected/uninfected individuals) vs. the observation space (the number who show up as positives and negatives on the screening test). If your screening test doesn't have a super high sensitivity/specificity rate (likely to be the case for a test administered in a village in the developing world, where cost is a larger factor than perfect accuracy) then this part of the model is potentially very important and should not be ignored.

A couple of more points should be mentioned about the observation function. First, the sensitivity/specificity parameters may be numbers that the test manufacturer has already estimated for you, in which case you should be able to just use those numbers directly in your model. However, if the manufacturer doesn't know them very precisely either, then you may be forced instead to add them to your larger model as an additional pair of parameters which you yourself will have to estimate, directly from the data. Another key point to understand is that it is the numbers which will be coming out of the observation function which will allow you to anchor your model to the seven years of existing archival data. Literally, the observation function is responsible for producing predictions, from the underlying hidden state variables (the "true" number of infected/uninfected individuals) which can be compared directly to the "CF_S" and "CF_D" values in your data table. Lastly, in terms of the dynamical systems diagram that I previously shared above, the observation function will correspond to the vertically oriented, downward pointing arrows; i.e., pointing from x(t) (system state) to y(t) (observation state) in the diagram.

So, now we come to the last part: estimating the model parameters. You model parameters will include the following items: 1.) initial numbers of infected/uninfected people, referenced either to year 0, the year before the NGO started taking data, or to year 1, the first year in which observations were kept, 2.) all of the state update-related parameters; i.e., $S$, $P$, $R$, $J$, basically all of the things that would have defined the "physics" of the system if we had been talking about a physical system model, and lastly 3.) any state-to-observation translational parameters (e.g., sensitivity and specificity), if they haven't already been provided to you by the screening test manufacturer.

To solve this type of problem and obtain an "answer" (i.e., estimates of your model parameters, which you can then use to predict infection rates for unknown future years), you can use a non-linear regression solver. An example of this type of package would be something like nls in R, or curve_fit in Python, or nlinfit in MATLAB.

All of these packages work pretty similarly. You give the fit package a top-level "fit function" which you will also have to write, which takes your model parameters as input, and makes predictions (using the state transition function and the observation function, which we've already defined above) for what the expected observations CF_S and CF_D would be over the course of 7 consecutive years, for an arbitrary combination of settings of the model parameters. The fit package, in effect, tries lots and lots of different combinations of input parameters (albeit in a somewhat structured and intelligent way, using algorithms which have been designed to produce an optimally fast search of the potentially available parameter space) until it comes up with an "optimal" combination of parameter settings, which result in a model that most closely reproduces the observed CF_S and CF_D data during the 7 year observation period.

Finally, once you have a data-driven solution for those "optimal" fitted parameter settings for $S$, $P$, $R$, $J$, etc., it is then a relatively simple matter to run the state transition and observation functions through one last final iteration, in order to predict your best guess of the observation values in future years for which you haven't yet acquired any data.

Also, one last thing: this probably almost goes without saying, but if you have a lousy model that doesn't adequately represent the true underlying dynamics of disease transmission among your study population, then not only will you obtain a poor fit to your data, but also your model won't have much predictive power either. So, don't be afraid to play around with different types of dynamical models, perhaps even ones based upon other parameters besides or in addition to $S$, $P$, $R$ and $J$ that you haven't even thought of yet. If you need some inspiration for how you might go about constructing alternative models, this is the point where a perusal of the existing journal and textbook literature (a few items of which I cited at the beginning) can really help you out.

Building off stachyra's answer, which is fairly extensive, as someone who does work in the modeling of infectious diseases. I'd contest that there's "considerably less of it" in any meaningful sense - mathematical models of disease are quite common, well developed, etc.

It sounds like what you're looking to do is construct a meta-population model, where each village has a separate infection process within it, and then terms governing the transmission of the disease between each village.

Depending on the disease in question (it would be helpful to know what disease you're thinking about) you can build a mechanistic model of the disease, and then your NGO's data provides observed values for some states within that model, and the rest can be obtained from the literature, or using a curve fitting technique like MCMC or maximum likelihood. It's a question that's somewhat far afield from "straight" statistical work, but isn't impossible to do, and there are certainly ways to get some predictions from those models, though there may be a fairly high degree of uncertainty to them,

Any major language you wish to use should be fine for this problem - I've seen these models done in C, R, Python, MATLAB, etc. In terms of content, a solid introduction to infectious disease modeling would be Keeling and Rohani's Modeling Infectious Disease in Humans and Animals, which has code in MATLAB and C++ I believe in an online companion, and for which someone has translated very nearly everything into Python.

Two notes, one self-serving, the other...not.

1. Like statistical analysis of many kinds, there's something to be said for subject matter expertise and experience. While what you're suggesting is reasonably straightforward, it's not trivially so, and you may want to consider consulting someone who does infectious disease modeling for a living.
2. While infectious diseases have all kinds of fun non-linear dynamics, it's possible that, if the incidence of disease is stable enough, you may be able to get perfectly worthwhile estimates from a poisson regression model.