PoissonRegOveDis: Main Function on Over-dispersed Poisson Regression
In fblues/Test: Prediction region based on over-dispersed Poisson Regression

Description Usage Arguments Details Value Author(s) Examples

View source: R/PoissonRegOveDis.R

This is the main file which incorporates other functions.

PoissonRegOveDis(
  daycutpred = 137,
  maxday = 154,
  ordfit = 5,
  DayEffect = TRUE,
  Country = "USA",
  xival = NULL,
  dat = NULL,
  alpha = 0.05,
  yhighval = NULL,
  targetval = 1e+05,
  LegPos = c("topright", "left"),
  ToPlot = TRUE
)

`daycutpred`	threshold for prediction
`maxday`	maximum day for the prediction
`ordfit`	the order of polynomial of Poisson regression
`DayEffect`	day effect
`Country`	country name
`xival`	over-dispersion parameter, xi
`dat`	data
`alpha`	alpha
`yhighval`	graphical parameter
`targetval`	graphical parameter
`LegPos`	graphical parameter
`ToPlot`	plot or not

Okay

list

E. A. Pena

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (daycutpred = 137, maxday = 154, ordfit = 5, DayEffect = TRUE, 
    Country = "USA", xival = NULL, dat = NULL, alpha = 0.05, 
    yhighval = NULL, targetval = 1e+05, LegPos = c("topright", 
        "left"), ToPlot = TRUE) 
{
    if (is.null(dat)) {
        countrydata = UpdateToday()
    }
    else {
        countrydata = dat
    }
    dat = AnalysisByCountry2(Country, alldata = countrydata)
    daynum = dat$daynum
    deaths = dat$dailydeaths
    cumdeaths = dat$cumdeaths
    DateRep = dat$DateRep
    close.screen(all.screens = T)
    daynumc = daynum - daynum[1]
    X = NULL
    collab = NULL
    for (j in 0:ordfit) {
        X = cbind(X, daynumc^j)
        collab = c(collab, paste("DayNum", j))
    }
    colnames(X) = collab
    if (DayEffect) {
        dayclass = daynum%%7
        dayclass = factor(dayclass, levels = as.character(0:6), 
            labels = c("Monday", "Tuesday", "Wednesday", "Thursday", 
                "Friday", "Saturday", "Sunday"))
        ll = length(dayclass)
        desmat = model.matrix(~dayclass)[, -1]
        daylab = c("Tuesday", "Wednesday", "Thursday", "Friday", 
            "Saturday", "Sunday")
        X = cbind(X, desmat)
        colnames(X) = c(collab, daylab)
    }
    inpred = (daynum <= daycutpred)
    numdaysused = sum(inpred)
    Y = deaths[inpred]
    X = X[inpred, ]
    X0 = NULL
    dayspred = daynum[1]:maxday
    dayspredc = dayspred - daynum[1]
    for (j in 0:ordfit) {
        X0 = cbind(X0, dayspredc^j)
    }
    if (DayEffect) {
        dayspredclass = as.factor(dayspred%%7)
        dayspredclass = factor(dayspredclass, levels = 0:6, labels = c("Monday", 
            "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", 
            "Sunday"))
        ll = length(dayspredclass)
        desmat = model.matrix(~dayspredclass)[, -1]
        X0 = cbind(X0, desmat)
    }
    numdaystopred = maxday - daycutpred
    if (numdaystopred == 0) {
        alphaeff = alpha
    }
    else {
        alphaeff = 1 - (1 - alpha)^(1/numdaystopred)
    }
    out = OveDisPoiReg(Y, X, X0, alpha = alphaeff, xival = xival)
    outglm = out$outglm
    outglmsumm = out$outglmsumm
    thetahat = out$theta
    xihat = out$xi
    XiCov = out$XiCov
    Xmean = out$Xmean
    Xstd = out$Xstd
    FitPred = out$FitPred[, 1:3]
    if (is.null(yhighval)) {
        yhighval = max(c(deaths, FitPred[, 3]))
    }
    if (ToPlot) {
        plot(NULL, NULL, xlim = c(daynum[1], maxday), ylim = c(0, 
            yhighval), xlab = "Day Number Since 12/31/2019", 
            ylab = "Daily Deaths", main = paste("Daily Deaths due to Covid-19 in ", 
                Country))
        points(daynum[inpred], Y, col = "black", lwd = 3)
        matlines(dayspred, FitPred, col = c("blue", rep("red", 
            2)), lty = 1, pch = 1, lwd = 3)
        points(daynum[(!inpred) & (daynum <= maxday)], deaths[(!inpred) & 
            (daynum <= maxday)], col = "red", lwd = 3)
        abline(v = daycutpred, col = "black", lwd = 2)
        abline(v = maxday, col = "green", lwd = 2)
    }
    lendayspred = length(dayspred)
    cumdeathspred = rep(0, lendayspred)
    cumdeathspredlow = rep(0, lendayspred)
    cumdeathspredhig = rep(0, lendayspred)
    for (i in 1:lendayspred) {
        cumdeathspred[i] = sum(FitPred[1:i, 1])
        cumdeathspredlow[i] = sum(FitPred[1:i, 2])
        cumdeathspredhig[i] = sum(FitPred[1:i, 3])
    }
    if (daycutpred < maxday) {
        daysahead = (daycutpred + 1):maxday
        lendaysahead = length(daysahead)
        lennow = length(daynum[daynum <= daycutpred])
        currentval = cumdeaths[lennow]
        cumdeathspredahead = rep(currentval, lendaysahead)
        cumdeathspredaheadlow = rep(currentval, lendaysahead)
        cumdeathspredaheadhig = rep(currentval, lendaysahead)
        for (i in 1:lendaysahead) {
            cumdeathspredahead[i] = cumdeathspredahead[i] + sum(FitPred[(lennow + 
                1):(lennow + i), 1])
            cumdeathspredaheadlow[i] = cumdeathspredaheadlow[i] + 
                sum(FitPred[(lennow + 1):(lennow + i), 2])
            cumdeathspredaheadhig[i] = cumdeathspredaheadhig[i] + 
                sum(FitPred[(lennow + 1):(lennow + i), 3])
        }
        maxy = max(cumdeathspredaheadhig)
    }
    else {
        maxy = max(cumdeaths[inpred])
    }
    if (ToPlot) {
        scan()
        plot(NULL, NULL, xlim = range(dayspred), ylim = c(0, 
            maxy), xlab = "Day Number Since 12/31/2019", ylab = "Cumulative Deaths", 
            main = paste("Cumulative Deaths Due to Covid-19 in ", 
                Country))
        lines(dayspred, cumdeathspred, type = "l", col = c("blue"), 
            lty = 1, pch = 1, lwd = 3)
        points(daynum[inpred], cumdeaths[inpred], col = "black", 
            lwd = 3)
        points(daynum[(!inpred) & (daynum <= maxday)], cumdeaths[(!inpred) & 
            (daynum <= maxday)], col = "red", lwd = 3)
        if (daycutpred < maxday) {
            matlines(daysahead, cbind(cumdeathspredaheadlow, 
                cumdeathspredaheadhig), col = "magenta", lty = 1, 
                pch = 2, lwd = 3)
        }
        abline(v = daycutpred, col = "black", lwd = 2)
        abline(v = maxday, col = "green", lwd = 2)
    }
    summpreds = NULL
    if (daycutpred < maxday) {
        summpreds = data.frame(daynum = daysahead, Pred = cumdeathspredahead, 
            PILow = cumdeathspredaheadlow, PIHig = cumdeathspredaheadhig)
    }
    return(list(Country = Country, dat = dat, outglm = outglm, 
        outglmsumm = outglmsumm, daycutpred = daycutpred, numdaysused = numdaysused, 
        maxday = maxday, numdaystopred = numdaystopred, ordfit = ordfit, 
        alpha = alpha, alphaeff = alphaeff, thetahat = thetahat, 
        xihat = xihat, XiCov = XiCov, Xmean = Xmean, Xstd = Xstd, 
        summpreds = summpreds))
  }