Covid-19 Analysis

[1] -1

Load data

Iowa COVID-19 Summary 2020-08-08

Today we had 384 new positive tests which is 9.8% (9.3% in the last week) of the tests recorded in the last 24 hours. There are 10,337 people actively sick and 2% (229) are hospitalized. 25.3% (58) of those hospitalized are in the ICU.

View Data

So, it turns out to be very difficult to find past data in the State of Iowa since they re-publish everything daily. Here are the last 6 days of data if you are interested:

tail(clean)

Data Exploration

Known Active Cases

This is a graph of the known active cases in Iowa. These people have tested positive, but have not moved on to a stage where they are not infectious. I say it that way because we know that a lot of patients with mild and moderate COVID-19 are very sick even after they are no longer infectious.

null device 
          1 

Sigh. On June 30 IDPH decided to count anyone who was still not listed as recovered after 28 days as recovered. So, let’s go back into the historic data and try to add that expectation to past recoverd data. Fortunately it is possible to download the data from the website https://coronavirus.iowa.gov/#CurrentStatus

Percent Positive Active Cases

null device 
          1 

png 
  2 

Bremer County

I started recording daily data on Bremer county on May 29, 2020. There will be problems with the recovered data since I have yet to get tables of data for the county and the definition of recovered changed on June 30. Thus, the following graph omits that data.

null device 
          1 

png 
  2 

png 
  2 

There are 97 active cases and 2 hospitalized in Bremer County today. The 7 day positivity rate is 15.5%.

Positive Tests

Logistic Model

To try to fit a logistic model to the data, we want to fit a parabola to the rate of change, but the input variable is cases or deaths and not time. This is because the logistic model comes from the differential equation \(\displaystyle{\frac{dP}{dt}= kP(1-\frac{P}{M})}\). Note that this equation as a function of \(P\) is a parabola with zeroes at \(P=0\) and \(P=M\).

clean <- clean[-c(1,10),]
with(clean, plot(positive, New.Positive, pch=16, xlab = "Positive Counts", ylab = "Rate of Change of Positive", cex.lab = 1.3, col = "blue"))

Now, it doesn’t really look quadratic, but we are not going to let that stop us from mathematics. Next, we will fit a quadratic model to the data.

clean$P2 <- as.numeric(clean$positive)^2
quadratic <-lm(New.Positive ~ positive + P2-1, data = clean)
r <- quadratic$coefficients[2]*(-1)
M <- quadratic$coefficients[1]/r
Mprint <- format(M, scientific = FALSE)
predictedcounts <- with(clean,predict(quadratic,list(positive = positive, P2 = P2)))

Let’s plot the model:

So, we have a quadratic model that looks like \(\frac{dP}{dt}\) = 0.0275512\(P^2\) + -3.923460210^{-7}\(P\) = 3.923460210^{-7}\(P(1-P/\) 70221.79\()\). This is looking pretty good. In particular, because \(M=\) 70221.79, we can estimate that 2.23% of Iowa’s population will become infected. We would like to see what that looks like for the original data. When R performs logistic modeling, it usually has binomial data and thus expects numbers that are between zero and one. We will use our estimate of \(M\) to create positive counts that are a percentage of \(M\).

non-integer #successes in a binomial glm!

Okay, that looks quite good. What will it predict for the next 28 days?

This predicts that in 14 days, Iowa will have 55813 infected people.

SIR Model

This should be more interesting and probably won’t be done for a week. The differential equations for the SIR model are:

\(\frac{dS}{dt} = -r_1SI\), \(\frac{dR}{dt} = r_2I\), \(\frac{dI}{dt} = -\frac{dS}{dt} - \frac{dR}{dt}\)

So, we may have enough data to estimate \(r_1\) and \(r_2\). We would do this using a linear model on the change of recovered to esitmate \(r_2\) and then that model and the rate of change of infecteds to estimate \(r_1\). We may have to use \(M\) from the logistic growth model to estimate \(S\).

Step 1: Estimate \(r_2\)

SIR <- with(clean, data.frame(
    dS = (-1)*diff(positive),
    dR = diff(Recovered),
    dI = diff(Still.Sick),
    dD = diff(deaths)
))
SIR <- SIR %>% mutate(
    R = clean$Recovered[-1],
    I = clean$Still.Sick[-1],
    day = clean$date[-1],
    S = M - clean$positive[-1]
    )
SIR$day = as.POSIXct(SIR$day)
SIR <- SIR[complete.cases(SIR),]

Let’s first just graph our data

qplot(day, data = SIR) + 
    stat_smooth(aes( y = S),method = "loess") + 
    stat_smooth( aes( y = I), color = "red", method = "loess") +
    stat_smooth( aes( y = R), color = "green", method = "loess") +
    theme(legend.position = "right")

(Note 5/30/2020) As of about March 15 this no longer looks like the SIR model that you would expect. We have this unexpected linear growth in \(R\) where we would expect things to be more concave down at this point. Not good.

Next, let’s see about coefficients:

dRmodel <- lm(dR ~ 0+I, data = SIR)
r2 <- dRmodel$coefficients[1]
dImodel <- lm (dS ~ 0 + S*I, data = SIR)
r1 <- dImodel$coefficients[1]

Now we plug this into the DE solver.

require(deSolve)
state <- c(S = as.numeric(1-clean$positive[25]/M), 
           I = as.numeric(clean$Still.Sick[25]/M), 
           R = as.numeric(clean$Recovered[25]/M))
times <- seq(0, nrow(clean)-24, by = 1)
parameters <- c(r1=as.numeric(r1*M),
                r2=as.numeric(r2))

# R function to calculate the value of the derivatives at each time value
    # Use the names of the variables as defined in the vectors above
SIR <- function(time, state, parameters){
      with(as.list(c(state, parameters)), {
        dS = r1*S*I
        dR = r2*I
        dI = -r1*S*I -r2*I
        return(list(c(dS, dI, dR)))
      })
}

## Integration with 'ode'
out <- ode(y = state, times = times, func = SIR, parms = parameters, method = rk4)
    
    ## Ploting
    out.df = as.data.frame(out) # required by ggplot: data object must be a data frame
    require(reshape2)
    out.m = melt(out.df, id.vars='time') # this makes plotting easier by puting all variables in a single column
    
    p <- ggplot(out.m, aes(time, value, color = variable)) + geom_point()
    #print(p)

Exponential Model

Run a linear model on the logarithm of the positive cases.

qplot(date, log(positive), data =clean) + geom_point() + geom_smooth()+
  stat_smooth(method = "lm", col = "red") 

For the model \(positive = e^{k date}\), \(k=\) 0.0399945. Estimated doubling time for positive cases is 17.3310456 days. Here is what that looks like on the original data.

Modified Exponential Model

By now there is strong evidence that the curve is flattening. We can create a new linear model that uses the more recent data to improve our model. Let’s look to see what happens if we only use data from later when we were doing more testing.

For the model \(positive = e^{k date}\), \(k=\) 0.0182274. Estimated doubling time for positive cases is 23.7213056 days. Here is what that looks like on the original data.

Deaths

Exponential Model

We will start the data when we have non-zero deaths.

For the model \(deaths = e^{k date}\), \(k=\) 7.6213089. Estimated time for deaths to double is 21.3708173 days.

Later Data

As above, let’s look at the same models, but only after we got serious about testing.

For the model \(deaths = e^{k date}\), \(k=\) 0.0238919. Estimated time for deaths to double is 29.011853 days.

New Positives

At the beginning, we were modeling the new positive data as an exponential curve. As of about the middle of May we changed to use the logistic growth curve from the Logistic Model to predict new positive tests.

As of about the middle of May, it becomes clear that the new positive data is not exponential. Which is really good news. This is when this model gets downgraded to the Alternative Modeling section.

Accuracy

As with any test, there will be false positives (people who test positive, but aren’t infected) and false negatives (people who test negative, but are infected). On May 1, it was reported by the Governor’s office:

“Governor Reynolds says the State Hygienic Lab has validated the machines used as part of the Test Iowa program. “I’m pleased to announce that the State Hygienic Lab completed the Test Iowa validation process yesterday, achieving high ratings of 95% accuracy for determining positives and 99.7% accuracy for determining negatives.”"

This is a bit hard to interpret, but I am going to assume that it means that the false positive rate is 0.3% and the false negative rate is 5%. This agrees with earlier reports of testing where the tests used to show that you are currently infected have pretty low false positive rates and higher false negative rates because it comes back positive if it finds the virus. A false negative would be that you have the virus, but the test can’t find it which isn’t out of the range of possibility. So, let’s go back to our data and calculate the numbers:

temp <- 1/(.997*.95 - (1-.995)*(1-.95))
accuracy <- clean %>% 
    mutate(
    tp = .95*positive/temp, 
    fp = (1-.993)/temp*positive, 
    tn = 0.993/temp*negative, 
    fn = (1-.95)/temp*negative) %>%
    mutate( pos = round(tp+fn, digits = 0),
            neg = round(tn+fp, digits = 0))
lastrow = length(accuracy)

Today, the state reported 13675 total positive and 134 total negative. Due to test inaccuracy, these numbers could be as much as NA positive and NA negative.

Store Predictions

This section stores predictions.

Alternative Modeling

This didn’t turn out well, but it is still here for historical perspective:

Positives Linear Model

First run a linear model on positive test results.

qplot(date, positive, data =clean) + geom_point() + geom_smooth()+
  stat_smooth(method = "lm", col = "red")

Deaths Linear Model

qplot(date, deaths, data =cdeaths, rm.na=TRUE) + geom_point() + geom_smooth()+
  stat_smooth(method = "lm", col = "red")

Ignoring unknown parameters: rm.na

model <- lm(deaths~date, data=cdeaths)
#model$coefficients

New Positives

As of about the middle of May, it becomes clear that the new positive data is not exponential. Which is really good news. This is when this model gets downgraded to the Alternative Modeling section. For the model \(new positives = e^{k date}\), \(k=\) 0.0182274. Estimated time for new positives to double is 38.0277688 days.