Modeling Corona Virus

Abstract

In this paper we will look at the SIR model for the mathematical modeling of Covid19 pendemic.We will be using R programming to solve the model

Introduction

The SIR Model is used in epidemiology to compute the amount of susceptible, infected, and recovered people in a population. It is also used to explain the change in the number of people needing medical attention during an epidemic. It is important to note that this model does not work with all diseases To know more of SIR model

Assumptions

The SIR Model is used in epidemiology to compute the amount of susceptible,infected, recovered people in a population. This model is an appropriate one to use under the following assumptions:
1. The population is fixed.
2. The only way a person can leave the susceptible group is to become infected.The only way a person can leave the infected group is to recover from the disease. Once a person has recovered, the person received immunity.
3. Age, sex, social status, and race do not affect the probability of being infected.
4. There is no inherited immunity.
5. The member of the population mix homogeneously (have the same interactions with one another to the same degree).

SIR Formulas

The model starts with some basic notation:

Woking Formula

Woking Formula

To Lean the mathematics behind them model see the link below Click here to see how SIR model works Process

Building model

Import data and clean data

library(dplyr)


library(reshape2)


library(googlesheets4)

confirmedCases <- read_sheet("https://docs.google.com/spreadsheets/d/1lirxEXoIeCTucFdrYmzwjG-nja7R94xWBPt6sIAT01M/edit#gid=0",1)
## Using an auto-discovered, cached token.
## To suppress this message, modify your code or options to clearly consent to the use of a cached token.
## See gargle's "Non-interactive auth" vignette for more details:
## https://gargle.r-lib.org/articles/non-interactive-auth.html
## The googlesheets4 package is using a cached token for krishnakumarshrestha00@gmail.com.
recoveredCases<- read_sheet("https://docs.google.com/spreadsheets/d/1lirxEXoIeCTucFdrYmzwjG-nja7R94xWBPt6sIAT01M/edit#gid=0",2)




deathCases <- read_sheet("https://docs.google.com/spreadsheets/d/1lirxEXoIeCTucFdrYmzwjG-nja7R94xWBPt6sIAT01M/edit#gid=0",3)

confirmedCases<-confirmedCases%>%select(-c(Lat,Long))%>%melt(id=c('Country/Region','Province/State'))
confirmedCases<-confirmedCases%>%group_by(`Country/Region`,variable)%>%summarise(Confirmed=sum(value))

deathCases<-deathCases%>%select(-c(Lat,Long))%>%melt(id=c('Country/Region','Province/State'))
deathCases<-deathCases%>%group_by(`Country/Region`,variable)%>%summarise(Deaths=sum(value))

recoveredCases<-recoveredCases%>%select(-c(Lat,Long))%>%melt(id=c('Country/Region','Province/State'))
recoveredCases<-recoveredCases%>%group_by(`Country/Region`,variable)%>%summarise(Recovered=sum(value))


# rename table columns
colnames(confirmedCases)<-c("Country","Date","Confirmed")
colnames(deathCases)<-c("Country","Date","Death")
colnames(recoveredCases)<-c("Country","Date","Recovered")


# merge all atbles together

mergedCases<-merge(confirmedCases,deathCases, by.y=c("Country","Date"))
mergedCases<-merge(mergedCases,recoveredCases, by.y=c("Country","Date"))

# convert factors to date format

mergedCases$Date<-as.Date(mergedCases$Date,"%m/%d/%y")
mergedCases<- mergedCases %>% filter(Country== "Nepal")

# summarize cases by date
x<-"03/23/2020"
x<-as.Date(x,"%m/%d/%y")

df1<-mergedCases %>% group_by(Date) %>% summarise_at(c("Confirmed","Recovered","Death"),sum)

The lockdown was started from 2020-03-23

Modeling

library(deSolve)
df2<-df1%>% filter(Date>x)



Infected <- as.integer(df2$Confirmed)
Day <- 1:(length(Infected))
N <- 10500000  # population 

#function for SIR model
SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * (I * S)
    dI <- beta/N *(I * S) - (gamma/25) * I
    dR <- (gamma/25) * I
    list(c(dS, dI, dR))
  })
}

init <- c(S = N*0.5-Infected[1], I = Infected[1], R = 1)

RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}


Opt <- optim(c(0.5, 0.2), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1)) # optimize with some sensible conditions

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
##      beta     gamma 
## 0.2280168 0.6358305

Predicting

library(DT)

t <- 1:80 # time in days(upto 80 days)
fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par))
Recovered=as.integer(fit$R)
Infecte=as.integer(fit$I)
Days=fit$time

 data<-data.frame(Days,Infecte,Recovered)

datatable(data = data)
col <- 1:2 # colour
ggplot(data=fit,aes(x=t))+
  geom_line(aes(y=I),col="red")+
  geom_line(aes(y=R),col="blue")+ggtitle("Covid 19 in Nepal")

a <- 1:180 # time in days
fit2 <- data.frame(ode(y = init, times = a, func = SIR, parms = Opt_par))
col1 <- 1:3 # colour

matplot(fit2$time, fit2[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col1);points(Day, Infected);
legend("topright", c("Susceptibles", "Infecteds", "Recovereds"), lty = 1, lwd = 2, col = col, inset = 0.05);title("SIR model COVID-19 Nepal", outer = TRUE, line = -2)

R0 <- setNames(Opt_par["beta"] / Opt_par["gamma"], "R0")



#fit[fit$I == max(fit$I), "I", drop = FALSE] # height of pandemic








d<-max(data$Infecte)* 0.033 # deaths with supposed 3.3% fatality rate i.e of India

Finding

So, according to this model,with in 30 days span number of cases will reach to 2186 and death will be 72.138 (considering death rate of india 3.3%). Basic reproduction number is 0.3586126 which is very low which means speed of transmission is slow in Nepal(World has 5.7).It will reach it’s peak around time span of 100 days if it still grows with this rate.

Conclusion

Do not panic! All of this is preliminary and hopefully (probably!) false.My model does not consider many factors which i have discussed below.But this is not a time to relax. Our government has not been able to increase testing and isolate people and open boader is making it much difficult.In addidtion to it contact tracing is not efficient and very flow. Hope fully in coming days Government functions properly with set plan and people obey them.

Problem with the model

All the finding that are produced by this model may not be accurate. The model is the general class of SIR differential equations used in epidemiology and is therefore not well suited for incorporating real world conditions at fine or large scale.
These include:

  1. Significant interactive local dynamics and travel restrictions that cannot be seen from aggregate quantities or averages across geographic locations

  2. Non-normal distributions of the number of infections per person (superspreader events) as well as the infection period,and

  3. Dynamic or stochastic values of parameters that arise from variations in sampling of distributions as well as the impact of changing social response efforts.

  4. this model doesnot consider some of the imortant measures take by Nepal government, our socio economic dynamice, population density , etc

  5. Another most important factor is open boader with india.

Suggestion

Despite not consedering many other factor , model can be used as a refrence and we can consider for worst case senerio.In order to Reduce the transmission rate and flat the curve as soon as possible ,here are few things needed to be done:**

  1. Number of testing should be increased rapidely(its not posssible to lockdown country forever) and isolate the infected.
  2. Contact Tracing should be done in fast track and Qurintine them.
  3. Social distancing should be strickly followed by general people.(As if now best way to reduce transmission when lockdown is over is social distancing and should be followed for at list 2-3 month prior the lockdown too).
  4. Personal sanitation should be kept at most priority (most important washing hand)
  5. It is better to bring all people in the border to Nepal.This way we can closely monitor them.

I am currently studying acturial science.This is not supposed to be perfect.Looking forward for the suggestion. Thank you.