fil0607
fil0607

Reputation: 101

Converting the outputs of a for loop from a list to a data frame

I have constructed a discrete time SIR model using a loop within a function (i have added my code below).

Currently the results of the iterations are coming out as a list which seems to show all the S values first followed by the I values and then the R values, which I have deduced myself from the nature of the values.

I need the output as a data frame with the column names: 'Iteration', 'S', 'I' and 'R' from left to right and the corresponding values underneath such that when a row is read it will tell you the iteration and values of S, I and R at that iteration.

I do not know how to construct a data frame that and returns the output values in this way, I have only started learning R a few weeks ago and so am not yet proficient so any help would be HUGELY appreciated.

Thank you in advance.

#INITIAL CONDITIONS 
S=999
I=1
R=0

#PARAMETERS 
beta  = 0.003 # infectious contact rate (/person/day)
gamma = 0.2 # recovery rate (/day)

#SIR MODEL WITH POISSON SAMPLING 
discrete_SIR_model <- function(){
for(i in 1:30){ #the number of iterations of loop indicates the 
  #duration of the model in days 
  # i.e. 'i in 1:30' constitutes 30 days

  deltaI<- rpois(1,beta * I * S) #rate at which individuals in the 
  #population are becoming infected
  
  deltaR<-rpois(1,gamma * I)#rate at which infected individuals are
  #recovering 
  
  S[i+1]<-S[i] -deltaI
  I[i+1] <-I[i] + deltaI -deltaR
  R[i+1]<-R[i]+deltaR
}
}
output <- list(c(S, I, R))
output

Upvotes: 0

Views: 53

Answers (1)

tpetzoldt
tpetzoldt

Reputation: 5813

If a foor loop is used, one can define vectors or a data frame beforehand where the results are stored:

beta  <- 0.001 # infectious contact rate (/person/day)
gamma <- 0.2   # recovery rate (/day)

S <- I <- R <- numeric(31)

S[1] <- 999
I[1] <- 1
R[1] <- 0

set.seed(123) # makes the example reproducible
for(i in 1:30){
  deltaI <- rpois(1, beta * I[i] * S[i])
  deltaR <- rpois(1, gamma * I[i])

  S[i+1] <- S[i] - deltaI
  I[i+1] <- I[i] + deltaI - deltaR
  R[i+1] <- R[i] + deltaR
}

output <- data.frame(S, I, R)
output
matplot(output)

As an alternative, it is also possible to employ a package for this. Package deSolve is intended for differential equations, but it can also solve the discrete case with method "euler":

library(deSolve)

discrete_SIR_model <- function(t, y, p) {
  with(as.list(c(y, p)), {
    deltaI <- rpois(1, beta * I * S)
    deltaR <- rpois(1, gamma * I)
    list(as.double(c(-deltaI, deltaI - deltaR,  deltaR)))
  })
}

y0 <- c(S = 999.0, I=1, R=0)

p <- c(
  beta  = 0.001, # infectious contact rate (/person/day)
  gamma = 0.2    # recovery rate (/day)
)

times <- 1:30

set.seed(576) # to make the example reproducible
output <- ode(y0, times, discrete_SIR_model, p, method="euler")

plot(output, mfrow=c(1,3))

Note: I reduced beta, otherwise the discrete model would become unstable.

discrete SIR

Upvotes: 1

Related Questions