COVID-19 Data of Germany

At first, we plot the size profile of the data and it clearly indicates that the logistic growth model seems to be an appropriate choice to start our analysis. But from the size profile, it is hard to conclude whether or not any variation is present in r. So, we plotted the ISRP profile of r. We obtained the empirical estimate of ISRP of r for logistic growth. To check whether there is any significant variation in r (with time), we fit different continuous functions given in Table 2 (Karim et al., 2022), and select the best model based on AIC, and RMSE values. As there is no sign of the presence of periodic variation in the data, we have not considered any kind of periodic variation in r. Based on the summary of AIC and RMSE values, we observe that r is subjected to hyperbolic variation. Hence, we infer that the actual data generation process is not logistic, but rather logistic growth with hyperbolically varying r. To support our claim, we fit several derived growth equations from the logistic model taken from Table 2 (Karim et al., 2022) to the cumulative number of cases. We observe that the logistic model with a hyperbolically varying growth coefficient is the best choice among all the models (both RMSE and AIC values support the conclusion). The analysis is carried out using the non-linear least-squares method using the nls2 function available in R.

# select the working directory where you stored the data 

if (getwd() != "D:/Covid_19_data") {

  setwd("D:/Covid_19_data")

}

Covid_data=read.csv(file="Covid_Germany_data.csv")

# Required packages

if(!require(pracma)){

  install.packages("pracma")

  library(pracma)

}

#taking 5 days moving average

data_1=movavg(Covid_data$total_cases,5,type="s")

plot(Covid_data$total_cases,col="grey",type="l",lwd=2,xlab="day",

     ylab="cumulative cases",main="Cumulative cases in Germany")

lines(data_1,col="red",lwd=2)

legend("topleft",legend=c("total cases","five days moving average"),

       lwd=c(2,2),col=c("black","red"))

#We can see that till day 63 there have been minor changes in the total number of #infected people so here we can neglect them. We took data from the 63rd day till last #from the moving average data

data=data_1[63:length(data_1)]

n=length(data)

t=0:(n-1)

h=numeric(length= n-1)

for(i in 1:length(h)){

  h[i]=t[i+1]-t[i]

}

#plot for the final data for Germany by taking five days moving average

plot(t,data,xlab = "Time (Day)", ylab = " Cumulative cases", col="red",lwd=1,cex=1.5)

#function to calculate Root Mean Square Error

rmse=function(predicted,observed){

  sqrt(sum((predicted-observed)^2)/length(predicted))

}

#fitting Logistic model

fun_logistic=function(x0,K,r,x){

  K/(1+(K/x0 - 1)*exp(-r*x))

}

# Required Packages

if(!require(nls2)){

  install.packages("nls2")

  library(nls2)

}

st1= expand.grid(x0=seq(0,100,len=10),K=seq(130000,200000,len=10),r=seq(0,3,len=10))

mod1 <- nls2(data~fun_logistic(x0,K,r,t), start = st1, algorithm = "brute-force")

# use nls object mod1 just calculated as starting value for nls optimization

logistic_fit=nls2(data~fun_logistic(x0,K,r,t), start = mod1)

summary(logistic_fit)

#plotting ISRP of parameter r of the logistic model fitted to data

r_hat_vals=numeric(length = n-2)

fun_r=function(x,y,z){

  (1/x-1/y)/(1/y-1/z)

}

for(i in 1:(n-2)){

  r_hat_vals[i]=log(fun_r(data[i],data[i+1],data[i+2]))/h[i]

}

par(mfrow=c(1,1), oma = c(1, 1, 2, 2), mar=c(4, 5, 2, 1))

plot(t[1:(length(t)-2)], r_hat_vals, col = "red", type = "p",

     ylab = expression(bold(widehat(r[j](Delta~t)))), xlab=expression(bold(Time~(Day))))

#fitting different variations of r on ISRP of r

#Linearly growing r

fun_1=function(r0,c,t){

  r0*(1+c*t)

}

param1=expand.grid(r0=seq(0,1,len=10),c=seq(0,2,len=10))

fit_1=nls2(r_hat_vals~fun_1(r0,c,t[1:(length(t)-2)]),start = param1, algorithm = "brute-force")

summary(fit_1)

#linearly decaying r

fun_2=function(r0,c,t){

  r0*(1-c*t)

}

param2=expand.grid(r0=seq(0,1,len=10),c=seq(0,2,len=10))

fit_2=nls2(r_hat_vals~fun_2(r0,c,t[1:(length(t)-2)]),start = param2, algorithm = "brute-force")

summary(fit_2)

#Extended logistic

fun_3=function(r0,c,t){

  r0*t^(c-1)

}

param3=expand.grid(r0=seq(0,1,len=10),c=seq(1,2,len=10))

fit_3=nls2(r_hat_vals~fun_3(r0,c,t[1:(length(t)-2)]),start = param3, algorithm = "brute-force")

summary(fit_3)

#Exponentially decaying r

fun_4=function(r0,c,t){

  r0*exp(-c*t)

}

param4=expand.grid(r0=seq(0,1,len=10),c=seq(-1,1,len=10))

fit_4=nls2(r_hat_vals~fun_4(r0,c,t[1:(length(t)-2)]),start = param4, algorithm = "brute-force")

summary(fit_4)

#Exponentially growing r

fun_5=function(r0,c,t){

  r0*exp(c*t)

}

param5=expand.grid(r0=seq(0,1,len=10),c=seq(-1,1,len=10))

fit_5=nls2(r_hat_vals~fun_5(r0,c,t[1:(length(t)-2)]),start = param5, algorithm = "brute-force")

summary(fit_5)

#Hyperbolically varying r

fun_6=function(r0,c,t){

  r0/(1+c*t)

}

param6=expand.grid(r0=seq(0,1,len=10),c=seq(0,2,len=10))

fit_6=nls2(r_hat_vals~fun_6(r0,c,t[1:(length(t)-2)]),start = param6, algorithm = "brute-force")

summary(fit_6)

#########################

#fitting respective variations of the logistic model on data

#linearly growing r

fun_logistic_1=function(x0,K,r0,c,t){

  K/(1+(K/x0 - 1)*exp(-r0*t*(1+c*t/2)))

}

st1= expand.grid(x0=seq(10,100,len=10),K=seq(100000,200000,len=20),

                 r0=seq(0,3,len=10),c=seq(0,3,len=10))

final1 <- nls2(data~fun_logistic_1(x0,K,r0,c,t), start = st1, algorithm = "brute-force")

summary(final1)

#linearly decaying r

fun_logistic_2=function(x0,K,r0,c,t){

  K/(1+(K/x0 - 1)*exp(-r0*t*(1-c*t/2)))

}

st2= expand.grid(x0=seq(10,100,len=10),K=seq(100000,200000,len=10),

                 r0=seq(0,3,len=10),c=seq(0,3,len=10))

final2 <- nls2(data~fun_logistic_2(x0,K,r0,c,t), start = st2, algorithm = "brute-force")

summary(final2)

#Extended logistic

fun_logistic_3=function(x0,K,r0,c,t){

  K/(1+(K/x0 - 1)*exp(-(r0/c)*t^c))

}

st3= expand.grid(x0=seq(10,100,len=10),K=seq(100000,200000,len=20),

                 r0=seq(0,3,len=10),c=seq(0.01,0.5,len=20))

final3 <- nls2(data~fun_logistic_3(x0,K,r0,c,t), start = st3, algorithm = "brute-force")

summary(final3)

#Exponentially decaying r

fun_logistic_4=function(x0,K,r0,c,t){

  K/(1+(K/x0 - 1)*exp(-(r0/c)*(1-exp(-c*t))))

}

st4= expand.grid(x0=seq(10,100,len=10),K=seq(100000,200000,len=10),

                 r0=seq(0,3,len=10),c=seq(-2,2,len=10))

final4 <- nls2(data~fun_logistic_4(x0,K,r0,c,t), start = st4, algorithm = "brute-force")

#exponentially growing r

fun_logistic_5=function(x0,K,r0,c,t){

  K/(1+(K/x0 - 1)*exp(-(r0/c)*(exp(c*t)-1)))

}

st6= expand.grid(x0=seq(10,100,len=10),K=seq(100000,200000,len=20),

                 r0=seq(0,3,len=10),c=seq(-3,3,len=10))

final5 <- nls2(data~fun_logistic_5(x0,K,r0,c,t), start = st6, algorithm = "brute-force")

summary(final5)

#Hyperbollically varying r

fun_logistic_6=function(x0,K,r0,c,t){

  K/(1+(K/x0 - 1)*(1+c*t)^(-r0/c))

}

st6= expand.grid(x0=seq(10,100,len=10),K=seq(100000,200000,len=20),

                 r0=seq(0,2,len=10),c=seq(0,2,len=10))

final6 <- nls2(data~fun_logistic_6(x0,K,r0,c,t), start = st6, algorithm = "brute-force")

summary(final6)

############################

#final fittings on real data

plot(t,data,xlab = "Time (Day)", ylab = " Cumulative cases", col="red",lwd=1,cex=1.5)

x=t

curve(fun_logistic(coef(logistic_fit)[1],coef(logistic_fit)[2],coef(logistic_fit)[3],x)

      ,col="magenta",lwd=3,add=TRUE)

curve(fun_logistic_1(coef(final1)[1],coef(final1)[2],coef(final1)[3],coef(final1)[4],x)

      ,col="green",lwd=3,add=TRUE)

curve(fun_logistic_2(coef(final2)[1],coef(final2)[2],coef(final2)[3],coef(final2)[4],x)

      ,col="yellow",lwd=3,add=TRUE)

curve(fun_logistic_3(coef(final3)[1],coef(final3)[2],coef(final3)[3],coef(final3)[4],x)

      ,col="chocolate",lwd=3,add=TRUE)

curve(fun_logistic_4(coef(final4)[1],coef(final4)[2],coef(final4)[3],coef(final4)[4],x)

      ,col="blue",lwd=3,add=TRUE)

curve(fun_logistic_5(coef(final5)[1],coef(final5)[2],coef(final5)[3],coef(final5)[4],x)

      ,col="cyan",lwd=3,add=TRUE)

curve(fun_logistic_6(coef(final6)[1],coef(final6)[2],coef(final6)[3],coef(final6)[4],x)

      ,col="black",lwd=3,add=TRUE)

legend("bottomright", legend = c("Logistic", "Linearly growing r",

        "Linearly decaying r","Extended logistic","Exponentially decaying r",

        "Exponentially growing r","Hyperbolically varying r")

       ,col=c("magenta","green","yellow","chocolate","blue","cyan","black"),lwd=3)

#final fittings on ISRP of r

plot(t[1:(length(t)-2)], r_hat_vals, col = "red",lwd=2, type = "p",

     xlim=c(-1,(length(t)-2)),ylab = expression(bold(widehat(r[j](Delta~t)))), 

     xlab = "Time (Day)")

x=t[1:(length(t)-2)]

curve(fun_1(coef(fit_1)[1],coef(fit_1)[2],x),col="green",lwd=3,add=TRUE)

curve(fun_3(coef(fit_3)[1],coef(fit_3)[2],x),col="chocolate",lwd=3,add=TRUE)

curve(fun_4(coef(fit_4)[1],coef(fit_4)[2],x),col="blue",lwd=3,add=TRUE)

curve(fun_6(coef(fit_6)[1],coef(fit_6)[2],x),col="black",lwd=3,add=TRUE)

legend("topright", legend = c( "Linearly varying r", 

      "Extended logistic","Exponentially varying r",

      "Hyperbolically varying r") ,col=c("green","yellow","cyan","black"), 

      lwd=2)

####################################

#AIC and RMSE values for fitted model onto data

final_result = matrix(data = NA, nrow = 7, ncol = 2)

colnames(final_result) = c("RMSE", "AIC")

rownames(final_result) = c("logistic", "Linearly growing r", "linearly decaying r", 

                            "Extended logistic","Exponentially decaying r",

                            "Exponentially growing r","Hyperbolically varying r")

final_result[,1] = c(rmse(data,fun_logistic(coef(logistic_fit)[1],coef(logistic_fit)[2],

                                            coef(logistic_fit)[3],t)),

                     rmse(data,fun_logistic_1(coef(final1)[1],coef(final1)[2],

                                              coef(final1)[3],coef(final1)[4],t)),

                     rmse(data,fun_logistic_2(coef(final2)[1],coef(final2)[2]

                                              ,coef(final2)[3],coef(final2)[4],t)),

                     rmse(data,fun_logistic_3(coef(final3)[1],coef(final3)[2],

                                              coef(final3)[3],coef(final3)[4],t)),

                     rmse(data,fun_logistic_4(coef(final4)[1],coef(final4)[2],

                                              coef(final4)[3],coef(final4)[4],t)),

                     rmse(data,fun_logistic_5(coef(final5)[1],coef(final5)[2],

                                              coef(final5)[3],coef(final5)[4],t)),

                     rmse(data,fun_logistic_6(coef(final6)[1],coef(final6)[2],

                                              coef(final6)[3],coef(final6)[4],t)))

final_result[,2] = c(AIC(logistic_fit),AIC(final1),AIC(final2),AIC(final3),

                     AIC(final4),AIC(final5),AIC(final6))

print("AIC and RMSE values for fitted model onto data")

print(final_result)

ISRP_fit_result= matrix(data = NA, nrow = 4, ncol = 2)

colnames(ISRP_fit_result) = c("RMSE", "AIC")

rownames(ISRP_fit_result) = c( "Linearly varying r", 

                               "Extended logistic","Exponentially varying r",

                               "Hyperbolically varying r")

ISRP_fit_result[,1] = c(rmse(r_hat_vals,fun_1(coef(fit_1)[1],coef(fit_1)[2],

                                              t[1:(length(t)-2)])),

                        rmse(r_hat_vals,fun_3(coef(fit_3)[1],coef(fit_3)[2],

                                              t[1:(length(t)-2)])),

                        rmse(r_hat_vals,fun_4(coef(fit_4)[1],coef(fit_4)[2],

                                              t[1:(length(t)-2)])),

                        rmse(r_hat_vals,fun_6(coef(fit_6)[1],coef(fit_6)[2],

                                              t[1:(length(t)-2)])))

ISRP_fit_result[,2] = c(AIC(fit_1),AIC(fit_3),AIC(fit_4),AIC(fit_6))

print("AIC and RMSE values for fitted model onto Interval specific estimates of r are")

print(ISRP_fit_result)