Problems with delta(log(cumulative))

This is a follow-up again to my previous posts. I have since discovered that effects come in from the linear regression that make the use of cumulative cases even less valid.

Again, I am playing with simulated epidemics of 10 days free growth, followed by 10 days with \(NPI_1\) then 10 days with \(NPI_1\) and \(NPI_2\) What is key is that \(NPI_1\) lasts from day 11 to day 30 and \(NPI_2\) lasts from day 21 to day 30. This is exactly what was done in the original paper. The result is that the regression over-weights the \(NPI_1\) and so has to under-weight \(NPI_2\).

I could try to explain it in intuitive terms but I’m not sure that I’m right. Fundamentally it’s a result of apply linear regression to something that is just not linear.

Here are some examples of this going badly wrong. The highlight is where the NPI that causes growth to increase comes out better than the NPI that drops growth to 0!

library(knitr)
# This class represents and Non-pharmacological intervention.
# The NPI is in place from day |start| to day |end| inclusive.
setClass("npi",
         slots = list(
           start = "numeric",
           end = "numeric",
           factor = "numeric"
         ))
# Convenient constructor for NPIs.
NPI <-
  function(start, end, factor)
    new("npi",
        start = start,
        end = end,
        factor = factor)

# Returns true if a day is impacted by |npi|
in_npi <- function(npi, i) {
  i >= npi@start && i <= npi@end
}

# Returns the last day impacted by any NPI in |npis|
npis_end <- function(npis) {
  max(sapply(npis, function(x) {
    x@end
  }))
}

# A vector of 0/1 for each day simulated. A day gets a 1 if the NPI is active on that day, 0 otherwise.
# Used as a predictor in the linear regression.
npi_pred <- function(npis, i) {
  last <- npis_end(npis)
  npi = npis[[i]]
  pred <- numeric()
  value = 0
  for (d in 1:last) {
    pred <- append(pred, (if (in_npi(npi, d))
      1
      else
        0))
  }
  pred
}

# Scales a growth rate so that it becomes |rate| every |period| days.
scale <- function(rate, period) {
  exp(log(rate) / period)
}

# Returns a vector of growth rates by combining |npis| with a base growth rate.
# All growth is scaled over |period| days.
createGrowthRates <- function(r0, npis, period) {
  last <- npis_end(npis)
  rates <- numeric()
  for (i in 1:last) {
    rate <- r0
    for (npi in npis) {
      if (in_npi(npi, i)) {
        rate <- rate * npi@factor
      }
    }
    rates <- append(rates, scale(rate, period))
  }
  rates
}

# Create a timeseries of cases based on the supplied vector of daily growth rates.
createCases <- function(rates) {
  cases <- c(1)
  for (rate in rates) {
    cases <- append(cases, tail(cases, 1) * rate)
  }
  cases
}

# Returns the coefficients for the impact of |npis| by applying linear regression to diff(log(|cases|)).
regress <- function(npis, cases) {
  d_l = diff(log(cases))
  npi1 <- npi_pred(npis, 1)
  npi2 <- npi_pred(npis, 2)
  d = data.frame(y = d_l,
                 npi1 = npi1,
                 npi2 = npi2)
  r.model <- lm(y ~ npi1 + npi2, data = d)
  r.model$coefficients
}

# Dump out some info on the coefficients.
renderRates <- function(c, period, title) {
  rate_c = exp(c * period)
  d = data.frame(c, rate_c)
  print(kable(d, col.names=c("log(r)", paste("r ", period, " days")),
              caption=paste("Regression using ", title), valign='t'))
}

# Applying regression and print the results and comparisons.
regressAndRender <- function(npis, cases, period, title) {
  c <- regress(npis, cases)
  renderRates(c, period, title)
  c
}

# Apply the methodology to |npis| with a base growth rate of |r0|.
# Scale all rates to be per |period|.
# Returns a list of coefficient sets from the regressions it ran.
# It will skip daily-case regression if any NPI has a factor of 0.
doOne <- function(npis, r0, period) {
  rates <- createGrowthRates(r0, npis, period)
  daily = createCases(rates)
  cumulative = cumsum(daily)

  plot(daily, ylab="daily cases")
  plot(cumulative, ylab="cumulative cases")

  d_l_daily <- diff(log(daily))
  d_l_cumulative <- diff(log(cumulative))

  plot(d_l_daily,ylab="delta(log(daily cases))")
  plot(d_l_cumulative, ylab="delta(log(cumulative cases))")

  cs = list()
  if (min(sapply(npis, function(npi) {npi@factor})) > 0) {
    c1 = regressAndRender(npis, daily, period, "daily")
    cs = append(cs, c1)
  } else {
    print("Skipping daily due to 0-factor")
  }
  c2 = regressAndRender(npis, cumulative, period, "cumulative")
  cs = append(cs, c2)
}

Two equally effective NPIs

Consider an epidemic with growth rate of 2.5 per 7-days. Apply two NPIs, each of which cut the 7-day growth rate by 50%, one running from days 11-30 and the other from days 21-30. The final 10 days have both NPIs in place for a cut of 75%.

Using daily cases, linear regression assigns the expected coefficients from which we can recover the 7-day growth rate.

Using cumulative cases, linear regression assigns a much lower coefficient to the first NPI.

npis <- c(NPI(11, 30, 0.5), NPI(21, 30, 0.5))
cs = doOne(npis, 2.5, 7)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Regression using daily

	log(r)	r 7 days
(Intercept)	0.1308987	2.5
npi1	-0.0990210	0.5
npi2	-0.0990210	0.5

Regression using cumulative

	log(r)	r 7 days
(Intercept)	0.3136616	8.9856700
npi1	-0.2063621	0.2358559
npi2	-0.0646044	0.6362075

Two very unequal NPIs but the same outcome

Another epidemic where the first NPI reduces growth by 1% and the second NPI reduces it by 99% but still appears weaker in the cumulative regression

npis <- c(NPI(11, 30, .99), NPI(21, 30, .01))
cs = doOne(npis, 2.5, 7)

Regression using daily

	log(r)	r 7 days
(Intercept)	0.1308987	2.50
npi1	-0.0014358	0.99
npi2	-0.6578815	0.01

Regression using cumulative

	log(r)	r 7 days
(Intercept)	0.3136616	8.9856700
npi1	-0.1630651	0.3193539
npi2	-0.1335377	0.3926784

An NPI that makes things worse beats the NPI that ends the epidemic!

Finally the absurdity, where the first NPI increases growth by 5% and the second NPI stops growth entirely but still appears weaker from the regression!

npis <- c(NPI(11, 30, 1.05), NPI(21, 30, 0))
cs = doOne(npis, 2.5, 7)

[1] "Skipping daily due to 0-factor"

Regression using cumulative

	log(r)	r 7 days
(Intercept)	0.3136616	8.9856700
npi1	-0.1587407	0.3291686
npi2	-0.1549208	0.3380891