Load library and the data

The data is scrapped from the website of the Nigerian Centre for Disease Control (NCDC) i.e (NCDC 2020). The scrapping was done with some bits of tricks. Please see my post on that. The BREAKS were established from the visual inspection of the data (see (Nmadu, Yisa, and Mohammed 2009))

library(tidyverse)
library(splines)
library(Metrics)
library(scales)
library(readxl)
library(patchwork)
library(Dyn4cast)
BREAKS <- c(70, 131, 173, 228, 274, 326)
z.s <- readRDS("data/zs.RDS")
cumper <- readRDS("data/cumper.RDS")
niz <- readRDS("data/niz.RDS")
States_county <- read_excel("data/States.xlsx")

Like mentioned earlier, some tricks were used to obtain the daily data from the website of NCDC. The NCDC uploads the cumulative data on daily basis, so each day, the data were scrapped and then the daily incidence were recorded. For example, below is the data scrapped from the website on March 15, 2021.

Deaths <- readr::read_csv("data/today.csv")
names(Deaths) <- c("Region", "Confirmed", "Active", "Recovered", "Deaths")
DDmap <- Deaths
Deaths
## # A tibble: 37 × 5
## Region Confirmed Active Recovered Deaths
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Abia 1629 44 1564 21
## 2 Adamawa 942 641 270 31
## 3 Akwa Ibom 1707 483 1210 14
## 4 Anambra 1909 64 1826 19
## 5 Bauchi 1482 198 1267 17
## 6 Bayelsa 828 37 765 26
## 7 Benue 1188 575 591 22
## 8 Borno 1321 83 1200 38
## 9 Cross River 344 14 313 17
## 10 Delta 2593 780 1744 69
## # … with 27 more rows

Creating date items

The dates are created starting with the day COVID-19 was first recorded in Nigeria, which is February 29, 2020. The sequence of dates were created fron the day to the last day the dataset was uploaded, which is March 15, 2021. In addition, sequence of dates were created from March 16, 2021 to march the number of days the virus have been recorded in Nigeria.

ss <- seq(1:length(cumper$DailyCase))
Dss <- seq(as.Date("2020-02-29"), by = "day", length.out = length(cumper$DailyCase))
Dsf <- seq(as.Date(Dss[length(Dss)] + lubridate::days(1)),
by = "day", length.out = length(cumper$DailyCase))
day01 <- format(Dss[1], format = "%B %d, %Y")
daylast <- format(Dss[length(Dss)], format = "%B %d, %Y")
casesstarts <- paste("Starting from", day01, "to", daylast,
collapse=" ")
lastdayfo <- format(Dsf[length(Dsf)], format = "%B %d, %Y")

Modelling and forecasting of COVID-19 caes

The incidence of the virus was model using five machine learning algorithms which are:

Splines model without Knots
Smooth spline
Splines model with Knots (or breaks)
Quadratic Polynomial
ARIMA

Then using the ensemble technology, another set of four machine learning models were estimated as follows:

Essembled based on summed weight
Essembled based on weight of fit
Essembled based on weight
Essembled with equal weight

The esemble technology is an attempt to combine models into one especially were the individual models are weak (Mahoney 2021).

fit0 <- lm(Case ~ bs(Day, knots = NULL), data = niz)
fit <- lm(Case ~ bs(Day, knots = BREAKS), data = niz)
fit1 <- smooth.spline(niz$Day, niz$Case)
fita <- forecast::auto.arima(niz$Case)
fitpi <- lm(Case ~ Day + I(Day^2), data = niz)
kk <- forecast::forecast(fitted.values(fit),
h = length(Dsf))
kk0 <- forecast::forecast(fitted.values(fit0),
h = length(Dsf))
kk1 <- forecast::forecast(fitted.values(fit1),
h = length(Dsf))
kk10 <- forecast::forecast(fitted.values(fitpi),
h = length(Dsf))
kk2 <- forecast::forecast(forecast::auto.arima(z.s$Case), h = length(Dsf))
kk30 <- (fitted.values(fit) + fitted.values(fit0) +
fitted.values(fit1) + fitted.values(fitpi) +
fita[["fitted"]])/5
kk31 <- forecast::forecast(kk30, h = length(Dsf))
kk40 <- lm(niz[,2]~fitted.values(fit)*fitted.values(fit0)*
fitted.values(fit1)*fitted.values(fitpi)*
fita[["fitted"]])
kk41 <- forecast::forecast(fitted.values(kk40),
h = length(Dsf))
kk60 <- lm(niz[,2]~fitted.values(fit)+fitted.values(fit0)+
fitted.values(fit1)+fitted.values(fitpi)+
fita[["fitted"]])
kk61 <- forecast::forecast(fitted.values(kk60),
h = length(Dsf))
KK <- as.data.frame(cbind("Date" = Dsf,"Day" = ss, "Without Knots" = kk[["mean"]], "Smooth spline" = kk0[["mean"]], "With Knots" = kk1[["mean"]], "Polynomial" = kk10[["mean"]], "Lower ARIMA" = kk2[["lower"]], "Upper ARIMA" =
kk2[["upper"]]))
KK <- KK[,-c(7,9)]
names(KK) <- c("Date", "Day", "Without Knots", "Smooth spline", "With Knots", "Polynomial", "Lower ARIMA", "Upper ARIMA")
WK <- sum(KK$`Without Knots`)
WKs <- sum(KK$`With Knots`)
SMsp <- sum(KK$`Smooth spline`)
LA <- sum(KK$`Lower ARIMA`)
UA <- sum(KK$`Upper ARIMA`)
POL <- sum(KK$Polynomial)
RMSE <- c("Without knots" = rmse(niz[,2],
fitted.values(fit)),
"Smooth Spline" = rmse(niz[,2], fitted.values(fit0)),
"With knots" = rmse(niz[,2], fitted.values(fit1)),
"Polynomial" = rmse(niz[,2], fitted.values(fitpi)),
"Lower ARIMA" = rmse(niz[,2], fita[["fitted"]]),
"Upper ARIMA" = rmse(niz[,2], fita[["fitted"]]))
#RMSE <- 1/RMSE
RMSE_weight <- as.list(RMSE / sum(RMSE))
KK$Date <- as.Date(KK$Date, origin = as.Date("1970-01-01"))
DDf <- c("Without knots", "Smooth Spline",
"With knots", "Quadratic Polynomial",
"Lower ARIMA", "Upper ARIMA",
"Essembled with equal weight",
"Essembled based on weight",
"Essembled based on summed weight",
"Essembled based on weight of fit" )
KK$`Essembled with equal weight` <- kk31[["mean"]]
KK$`Essembled based on weight` <- kk41[["mean"]]
KK$`Essembled based on summed weight` <- kk61[["mean"]]
ESS <- sum(KK$`Essembled with equal weight`)
ESM <- sum(KK$`Essembled based on weight`)
ESMS <- sum(KK$`Essembled based on summed weight`)
P_weight <-
(fitted.values(fit) * RMSE_weight$`Without knots`) + (fitted.values(fit0) *
RMSE_weight$`Smooth Spline`) +
(fitted.values(fit1) * RMSE_weight$`With knots`) + (fitted.values(fitpi) *
RMSE_weight$Polynomial) + (fita[["fitted"]] * RMSE_weight$`Lower ARIMA`)
kk51 <- forecast::forecast(P_weight, h = length(Dsf))
KK$`Essembled based on weight of fit` <-
kk51[["mean"]]
ESF <- sum(KK$`Essembled based on weight of fit`)
RMSE$`Essembled with equal weight` <- rmse(niz[,2], kk30)
RMSE$`Essembled based on weight` <- rmse(niz[,2], fitted.values(kk40))
RMSE$`Essembled based on summed weight` <- rmse(niz[,2], fitted.values(kk60))
RMSE$`Essembled based on weight of fit` <- rmse(niz[,2], P_weight)
Forcasts <- colSums(KK[,-c(1,2)])
RMSE_f <- as.data.frame(cbind("Model" = DDf,
"Confirmed cases" =
comma(round(Forcasts, 0))))
Deaths$County <- States_county$County
Deaths <- Deaths[order(-Deaths$Deaths),]
Deaths$Region <- factor(Deaths$Region, levels = Deaths$Region)
Deaths$Ratio <- Deaths$Deaths/Deaths$Confirmed*100
Deaths <- Deaths[order(-Deaths$Ratio),]
Deaths$Region <- factor(Deaths$Region, levels = Deaths$Region)
Deaths$RecRatio <- Deaths$Recovered/Deaths$Confirmed*100
Deaths <- Deaths[order(-Deaths$RecRatio),]
Deaths$Region <- factor(Deaths$Region, levels = Deaths$Region)
PerDeath <- round(sum((Deaths$Deaths)/
sum(Deaths$Confirmed)*100), 2)
PerRec <- round(sum((Deaths$Recovered)/
sum(Deaths$Confirmed)*100), 2)
RMSE_f$Recoveries <- comma(round(Forcasts * PerRec/100, 0))
RMSE_f$Fatalities <- round(Forcasts * PerDeath/100, 0)
RMSE_f$RMSE <- RMSE
RMSE_f <- RMSE_f[order(-RMSE_f$Fatalities),]
RMSE_f$Model <- factor(RMSE_f$Model, levels = RMSE_f$Model)
RMSE_f$Fatalities <- comma(RMSE_f$Fatalities)
RMSE_fk <- knitr::kable(RMSE_f, row.names = FALSE, "html")

The report

kableExtra::kable_styling(RMSE_fk, "striped")

Take a further look at some of the visuals

plot(z.s[,1], z.s[,2], main = " ", col = "red", type = "l", lwd = 2, xlab = casesstarts, ylab = "COVID-19 confirmed cases")
grid(10,20,lty=1,lwd=1)
abline(v = BREAKS, lwd = 2, lty = 2, col = "black")
#lines(density(z.s$Case), lwd = 2, type = "l", col = "blue", bw = 69.22)
lines(fitted.values(smooth.spline(z.s[,1],z.s[,2])), type = "l",
col = "blue", lwd = 2)
DDf <- c("Observed", "Cut points", "Smooth spline") #, "Density")
legend("topleft",
inset = c(0, 0), bty = "n", x.intersp = 0.5,
xjust = 0, yjust = 0,
legend = DDf, col = c("red", "black", "blue"),
lwd = 2, cex = 0.75, xpd = TRUE, horiz = FALSE)

The above plot shows that as at March 15, 2021; Nigeria has gone through two waves and there were more cases of COVID-19 in the second wave.

plot(cumper$Per, main = " ",
xlab = casesstarts, ylab = "Per Cent",
ylim = c(-50,100),
type = "l", lwd = 2, col="forestgreen")
grid(10,20,lty=1,lwd=1)
lines(cumper$Per, lty=1,lwd=2, col= "forestgreen")
lines(cumper$PerCum, lty=1,lwd=2, col= "orangered")
lines(log(cumper$DailyCase), lty=1,lwd=2, col= "magenta")
lines(cumper$LPer, lty=1,lwd=2, col= "brown")
lines(cumper$LPerCum, lty=1,lwd=2, col= "blue")
legend("topright",
inset = c(0, 0, 0, 0), bty = "n", x.intersp = 0.5,
xjust = 0, yjust = 0,
"topright", col = c("forestgreen", "orangered", "magenta", "brown", "blue"),
c("Per cent of total cases","Percent of cummulative cases",
"Daily cases, log transformed",
"Difference from previsous rate",
"Difference from previous cummulative rate"),
lwd = 2, cex = 0.75, xpd = TRUE, horiz = FALSE)

The above two visuals show the log transformation of the distribution and the difference from previous day. The plots were mage with plot and ggplot functions respectively.

The plot on the right is the cumulative cases of COVID-19 recorded in Nigeria while on the left is the forecast using spline and polynomial models. The forecast is based on equal duation of the lenght of the data.