How to apply a system of ordinary differential equations to a cell grid?

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

How to apply a system of ordinary differential equations to a cell grid?

Marine
Hello,

I am developing an agent-based model to simulate the spread of infectious diseases in heterogeneous landscapes composed of habitat polygons (or clumps of connected cells). To simplify the model, I consider a habitat grid (or raster) containing the polygon ID of each cell. In addition, I have epidemiological parameters associated with each polygon ID. At each time step, the parameter values change in the polygon. Thus, the data frame �landscape� (see below) is updated at each time step. Here is an example at t = 0:

landscape <- data.frame(polygon_ID = seq(1, 10, by = 1),
                        beta = sample(c(100, 200, 400, 600), 10, replace = TRUE),
                        gamma = sample(c(25, 26, 27, 28), 10, replace = TRUE))

To study the disease dynamics, I also am developing a compartmental model based on a system of ordinary differential equations (ODEs). Here is an example to represent the system of ODEs:

solve_sir_model <- function (times, parameters) {

  sir_model <- function (times, states, parameters) {

    with(as.list(c(states, parameters)), {

      dSdt <- -beta*S*I
      dIdt <- beta*S*I-gamma*I
      dRdt <- gamma*I
      dNdt <- dSdt + dIdt + dRdt

      return(list(c(dSdt, dIdt, dRdt, dNdt)))

    })
  }

  states <- c(S = 99, I = 1, R = 0, N = 100)
  return(ode(y = states, times = times, func = sir_model, parms = parameters))
}

require(deSolve)
output <- as.data.frame(solve_sir_model(times = seq(0, 5, by = 1), parameters = c(beta = 400, gamma = 28)))

Here is my question: at each time step, is it possible to apply the system of ODEs to each habitat polygon (thus each row) in the data frame �landscape�? I am using lsoda as an ODE solver. Do I need to use another solver to apply the ODEs at each time step?

Thank you very much for your advice.
Have a nice day
Marine


        [[alternative HTML version deleted]]


______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: How to apply a system of ordinary differential equations to a cell grid?

David Winsemius

> On Jun 21, 2017, at 12:48 PM, Marine Regis <[hidden email]> wrote:
>
> Hello,
>
> I am developing an agent-based model to simulate the spread of infectious diseases in heterogeneous landscapes composed of habitat polygons (or clumps of connected cells). To simplify the model, I consider a habitat grid (or raster) containing the polygon ID of each cell. In addition, I have epidemiological parameters associated with each polygon ID. At each time step, the parameter values change in the polygon. Thus, the data frame �landscape� (see below) is updated at each time step. Here is an example at t = 0:
>
> landscape <- data.frame(polygon_ID = seq(1, 10, by = 1),
>                        beta = sample(c(100, 200, 400, 600), 10, replace = TRUE),
>                        gamma = sample(c(25, 26, 27, 28), 10, replace = TRUE))
>
> To study the disease dynamics, I also am developing a compartmental model based on a system of ordinary differential equations (ODEs). Here is an example to represent the system of ODEs:
>
> solve_sir_model <- function (times, parameters) {
>
>  sir_model <- function (times, states, parameters) {
>
>    with(as.list(c(states, parameters)), {
>
>      dSdt <- -beta*S*I
>      dIdt <- beta*S*I-gamma*I
>      dRdt <- gamma*I
>      dNdt <- dSdt + dIdt + dRdt
>
>      return(list(c(dSdt, dIdt, dRdt, dNdt)))
>
>    })
>  }
>
>  states <- c(S = 99, I = 1, R = 0, N = 100)
>  return(ode(y = states, times = times, func = sir_model, parms = parameters))
> }
>
> require(deSolve)
> output <- as.data.frame(solve_sir_model(times = seq(0, 5, by = 1), parameters = c(beta = 400, gamma = 28)))
>
> Here is my question: at each time step, is it possible to apply the system of ODEs to each habitat polygon (thus each row) in the data frame �landscape�? I am using lsoda as an ODE solver. Do I need to use another solver to apply the ODEs at each time step?
>
> Thank you very much for your advice.
> Have a nice day
> Marine
>

There's also ode.2D in the same package {deSolve} and it's help page has a 2-d diffusion example that might be cognate.

>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: How to apply a system of ordinary differential equations to a cell grid?

Marine
Thank you very much for your answer. I'm not sure if I can use the function ode.2 because it's a solver for 2-D partial differential equation problems. My equations don't contain diffusion parameters.


Thank you for your help

Marine


________________________________
De : David Winsemius <[hidden email]>
Envoyé : mercredi 21 juin 2017 22:30
À : Marine Regis
Cc : [hidden email]
Objet : Re: [R] How to apply a system of ordinary differential equations to a cell grid?


> On Jun 21, 2017, at 12:48 PM, Marine Regis <[hidden email]> wrote:
>
> Hello,
>
> I am developing an agent-based model to simulate the spread of infectious diseases in heterogeneous landscapes composed of habitat polygons (or clumps of connected cells). To simplify the model, I consider a habitat grid (or raster) containing the polygon ID of each cell. In addition, I have epidemiological parameters associated with each polygon ID. At each time step, the parameter values change in the polygon. Thus, the data frame �landscape� (see below) is updated at each time step. Here is an example at t = 0:
>
> landscape <- data.frame(polygon_ID = seq(1, 10, by = 1),
>                        beta = sample(c(100, 200, 400, 600), 10, replace = TRUE),
>                        gamma = sample(c(25, 26, 27, 28), 10, replace = TRUE))
>
> To study the disease dynamics, I also am developing a compartmental model based on a system of ordinary differential equations (ODEs). Here is an example to represent the system of ODEs:
>
> solve_sir_model <- function (times, parameters) {
>
>  sir_model <- function (times, states, parameters) {
>
>    with(as.list(c(states, parameters)), {
>
>      dSdt <- -beta*S*I
>      dIdt <- beta*S*I-gamma*I
>      dRdt <- gamma*I
>      dNdt <- dSdt + dIdt + dRdt
>
>      return(list(c(dSdt, dIdt, dRdt, dNdt)))
>
>    })
>  }
>
>  states <- c(S = 99, I = 1, R = 0, N = 100)
>  return(ode(y = states, times = times, func = sir_model, parms = parameters))
> }
>
> require(deSolve)
> output <- as.data.frame(solve_sir_model(times = seq(0, 5, by = 1), parameters = c(beta = 400, gamma = 28)))
>
> Here is my question: at each time step, is it possible to apply the system of ODEs to each habitat polygon (thus each row) in the data frame �landscape�? I am using lsoda as an ODE solver. Do I need to use another solver to apply the ODEs at each time step?
>
> Thank you very much for your advice.
> Have a nice day
> Marine
>

There's also ode.2D in the same package {deSolve} and it's help page has a 2-d diffusion example that might be cognate.

>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
R-help Info Page - Homepage - SfS – Seminar for Statistics<https://stat.ethz.ch/mailman/listinfo/r-help>
stat.ethz.ch
The main R mailing list, for announcements about the development of R and the availability of new code, questions and answers about problems and solutions using R ...



> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA


        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Hunting a histogram variant

S Ellison-2
In reply to this post by David Winsemius

I'm looking for a histogram variant in which data points are plotted as labelled rectangles 'piled up' to form a histogram. I've posted an example at https://www.dropbox.com/s/ozi8bhdn5kqaufm/labelled_histogram.png?dl=0

It seems to have a long pedigree, as I see it (as in this example) in documents going back beyond the '80s. But I've not seen it in recent textbooks. So it may be one of those older graphical displays that's just fallen out of use.

General questions:
i) Does this thing have a name?
ii) Can anyone point me to a literature source for it?

... and the R question:
ii) Is it already hiding somewhere in an R package?*

S Ellison

*If it's not, I'll be adding it to one, hence the hunt for due credit/sources



*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Hunting a histogram variant

Bert Gunter-2
?stem

for something close and built in.

Cheers,
Bert


Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Wed, Jun 21, 2017 at 9:01 PM, S Ellison <[hidden email]> wrote:

>
> I'm looking for a histogram variant in which data points are plotted as labelled rectangles 'piled up' to form a histogram. I've posted an example at https://www.dropbox.com/s/ozi8bhdn5kqaufm/labelled_histogram.png?dl=0
>
> It seems to have a long pedigree, as I see it (as in this example) in documents going back beyond the '80s. But I've not seen it in recent textbooks. So it may be one of those older graphical displays that's just fallen out of use.
>
> General questions:
> i) Does this thing have a name?
> ii) Can anyone point me to a literature source for it?
>
> ... and the R question:
> ii) Is it already hiding somewhere in an R package?*
>
> S Ellison
>
> *If it's not, I'll be adding it to one, hence the hunt for due credit/sources
>
>
>
> *******************************************************************
> This email and any attachments are confidential. Any u...{{dropped:8}}

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Hunting a histogram variant

Richard M. Heiberger
In reply to this post by S Ellison-2
Two things to look at are

?monthplot
which shows seasonal subseries plotted on an overall plot.
The parallel to your example is small histograms plotted on an overall plot.
That goes back (through the reference in the Blue Book (Becker,
Chambers, and Wilks)) to
Bill Cleveland and Irma Terpenning in JASA 77, 377, 52-62 (1982).

Starting from scratch I would be inclined to build your design
directly in grid, to place the
subgraphs in the main graph, and then I would use lattice to design
the subgraphs.
See Paul Murrell's book Figure 7.18 (First Edition) showing plots of
temperature for
cities in Australia.

Rich

On Thu, Jun 22, 2017 at 12:01 AM, S Ellison <[hidden email]> wrote:

>
> I'm looking for a histogram variant in which data points are plotted as labelled rectangles 'piled up' to form a histogram. I've posted an example at https://www.dropbox.com/s/ozi8bhdn5kqaufm/labelled_histogram.png?dl=0
>
> It seems to have a long pedigree, as I see it (as in this example) in documents going back beyond the '80s. But I've not seen it in recent textbooks. So it may be one of those older graphical displays that's just fallen out of use.
>
> General questions:
> i) Does this thing have a name?
> ii) Can anyone point me to a literature source for it?
>
> ... and the R question:
> ii) Is it already hiding somewhere in an R package?*
>
> S Ellison
>
> *If it's not, I'll be adding it to one, hence the hunt for due credit/sources
>
>
>
> *******************************************************************
> This email and any attachments are confidential. Any u...{{dropped:8}}

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Hunting a histogram variant

Jim Lemon-4
In reply to this post by S Ellison-2
Hi,
This looks very much like a waffle plot, which can be performed with
color2D.matplot in the plotrix package. At the moment, the function
fill.corner is used to produce the matrix, but this could be
generalized to produce any shape of "waffle". The value display is
already part of color2D.matplot.

Jim

On Thu, Jun 22, 2017 at 2:01 PM, S Ellison <[hidden email]> wrote:

>
> I'm looking for a histogram variant in which data points are plotted as labelled rectangles 'piled up' to form a histogram. I've posted an example at https://www.dropbox.com/s/ozi8bhdn5kqaufm/labelled_histogram.png?dl=0
>
> It seems to have a long pedigree, as I see it (as in this example) in documents going back beyond the '80s. But I've not seen it in recent textbooks. So it may be one of those older graphical displays that's just fallen out of use.
>
> General questions:
> i) Does this thing have a name?
> ii) Can anyone point me to a literature source for it?
>
> ... and the R question:
> ii) Is it already hiding somewhere in an R package?*
>
> S Ellison
>
> *If it's not, I'll be adding it to one, hence the hunt for due credit/sources
>
>
>
> *******************************************************************
> This email and any attachments are confidential. Any u...{{dropped:8}}

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Hunting a histogram variant

Jim Lemon-4
In reply to this post by S Ellison-2
and I almost forgot Derek Ogle's histStack again in the plotrix package.

Jim

On Thu, Jun 22, 2017 at 2:01 PM, S Ellison <[hidden email]> wrote:

>
> I'm looking for a histogram variant in which data points are plotted as labelled rectangles 'piled up' to form a histogram. I've posted an example at https://www.dropbox.com/s/ozi8bhdn5kqaufm/labelled_histogram.png?dl=0
>
> It seems to have a long pedigree, as I see it (as in this example) in documents going back beyond the '80s. But I've not seen it in recent textbooks. So it may be one of those older graphical displays that's just fallen out of use.
>
> General questions:
> i) Does this thing have a name?
> ii) Can anyone point me to a literature source for it?
>
> ... and the R question:
> ii) Is it already hiding somewhere in an R package?*
>
> S Ellison
>
> *If it's not, I'll be adding it to one, hence the hunt for due credit/sources
>
>
>
> *******************************************************************
> This email and any attachments are confidential. Any u...{{dropped:8}}

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Hunting a histogram variant

Peter Dalgaard-2
In reply to this post by Bert Gunter-2
Hmmno... The labels on a stem-and-leaf plots are the values. This is just the measurement number: Observations #2,5,6,9 from level 1 had a temperature between 89 and 90, making up the penultimate column of that histogram.

I would conjecture that, like stem-and-leaf, this has fallen out of favour because it doesn't scale well to larger samples. It is fine with 64 observations like this, but with (say) four times as many boxes, you'd lose all legibility.

-pd

> On 22 Jun 2017, at 06:16 , Bert Gunter <[hidden email]> wrote:
>
> ?stem
>
> for something close and built in.
>
> Cheers,
> Bert
>
>
> Bert Gunter
>
> "The trouble with having an open mind is that people keep coming along
> and sticking things into it."
> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>
>
> On Wed, Jun 21, 2017 at 9:01 PM, S Ellison <[hidden email]> wrote:
>>
>> I'm looking for a histogram variant in which data points are plotted as labelled rectangles 'piled up' to form a histogram. I've posted an example at https://www.dropbox.com/s/ozi8bhdn5kqaufm/labelled_histogram.png?dl=0
>>
>> It seems to have a long pedigree, as I see it (as in this example) in documents going back beyond the '80s. But I've not seen it in recent textbooks. So it may be one of those older graphical displays that's just fallen out of use.
>>
>> General questions:
>> i) Does this thing have a name?
>> ii) Can anyone point me to a literature source for it?
>>
>> ... and the R question:
>> ii) Is it already hiding somewhere in an R package?*
>>
>> S Ellison
>>
>> *If it's not, I'll be adding it to one, hence the hunt for due credit/sources
>>
>>
>>
>> *******************************************************************
>> This email and any attachments are confidential. Any u...{{dropped:8}}
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: How to apply a system of ordinary differential equations to a cell grid?

Marine
In reply to this post by David Winsemius
It seems that the method `iteration` in the function `ode` can be useful in my case:

> Method "iteration" is special in that here the function func should
> return the new value of the state variables rather than the rate of
> change. This can be used for individual based models, for difference
> equations, or in those cases where the integration is performed within
> func).

I have tested the method but I don't understand why it doesn't work with one time step:

    solve_sir_model <- function (times, parameters) {

      sir_model <- function (times, states, parameters) {

        with(as.list(c(states, parameters)), {

          dSdt <- -beta*S*I
          dIdt <- beta*S*I-gamma*I
          dRdt <- gamma*I
          dNdt <- dSdt + dIdt + dRdt

          return(list(c(dSdt, dIdt, dRdt, dNdt)))

        })
      }

      states <- c(S = 99, I = 1, R = 0, N = 100)
      return(ode(y = states, times = times, func = sir_model, parms = parameters, method = "iteration"))
    }

    require(deSolve)
    output <- as.data.frame(solve_sir_model(times = 0, parameters = c(beta = 400, gamma = 28)))

 Error in iteration(y, times, func, parms, ...) :
   times should be equally spaced In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

Thank you very much for your help.
Marine



________________________________
De : David Winsemius <[hidden email]>
Envoyé : mercredi 21 juin 2017 22:30
À : Marine Regis
Cc : [hidden email]
Objet : Re: [R] How to apply a system of ordinary differential equations to a cell grid?


> On Jun 21, 2017, at 12:48 PM, Marine Regis <[hidden email]> wrote:
>
> Hello,
>
> I am developing an agent-based model to simulate the spread of infectious diseases in heterogeneous landscapes composed of habitat polygons (or clumps of connected cells). To simplify the model, I consider a habitat grid (or raster) containing the polygon ID of each cell. In addition, I have epidemiological parameters associated with each polygon ID. At each time step, the parameter values change in the polygon. Thus, the data frame �landscape� (see below) is updated at each time step. Here is an example at t = 0:
>
> landscape <- data.frame(polygon_ID = seq(1, 10, by = 1),
>                        beta = sample(c(100, 200, 400, 600), 10, replace = TRUE),
>                        gamma = sample(c(25, 26, 27, 28), 10, replace = TRUE))
>
> To study the disease dynamics, I also am developing a compartmental model based on a system of ordinary differential equations (ODEs). Here is an example to represent the system of ODEs:
>
> solve_sir_model <- function (times, parameters) {
>
>  sir_model <- function (times, states, parameters) {
>
>    with(as.list(c(states, parameters)), {
>
>      dSdt <- -beta*S*I
>      dIdt <- beta*S*I-gamma*I
>      dRdt <- gamma*I
>      dNdt <- dSdt + dIdt + dRdt
>
>      return(list(c(dSdt, dIdt, dRdt, dNdt)))
>
>    })
>  }
>
>  states <- c(S = 99, I = 1, R = 0, N = 100)
>  return(ode(y = states, times = times, func = sir_model, parms = parameters))
> }
>
> require(deSolve)
> output <- as.data.frame(solve_sir_model(times = seq(0, 5, by = 1), parameters = c(beta = 400, gamma = 28)))
>
> Here is my question: at each time step, is it possible to apply the system of ODEs to each habitat polygon (thus each row) in the data frame �landscape�? I am using lsoda as an ODE solver. Do I need to use another solver to apply the ODEs at each time step?
>
> Thank you very much for your advice.
> Have a nice day
> Marine
>

There's also ode.2D in the same package {deSolve} and it's help page has a 2-d diffusion example that might be cognate.

>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
R-help -- Main R Mailing List: Primary help - Homepage - SfS<https://stat.ethz.ch/mailman/listinfo/r-help>
stat.ethz.ch
The main R mailing list, for announcements about the development of R and the availability of new code, questions and answers about problems and solutions using R ...



> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA


        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.