Analysis of proportion of variants 501Y.V1 and 501Y.V2 in Belgium

VariantAnalysis.R

# remotes::install_github("rvlenth/emmeans")
# use the latest version of emmeans, as the current cran version 1.5.4 containsa bug
library(dplyr)
library(tidyr)
library(ggplot2)
library(lme4)
library(emmeans)

# create fits function:
create_fit <- function(x, variant = "nv1",n = 100){
  fit <- glmer(cbind(npos, n-npos) ~ weeknum + (1|weeknum),
               data = x, family = binomial)
  sigma <- as.data.frame(VarCorr(fit))$sdcor
  grid <- list(weeknum = seq(min(x$weeknum),
                             max(x$weeknum) + 0.1,
                             length.out = n))
  ci <- emmeans(fit, ~ weeknum,
                at = grid,
                type = "response")
  out <- as.data.frame(ci, bias.adjust = TRUE, sigma = sigma)
  out$variant <- variant
  return(out)
}

# Make date sequences for fitting and plots
weekseq <- seq.POSIXt(from = as.POSIXct("2020-12-6"),
                    to = as.POSIXct("2021-2-14"),
                    by = "1 week")
weeknum <- scale(as.numeric(weekseq))

# Construct data from weekly epidemiological report Sciensano
# dated february 19, 2021
thedata <- 
  tibble(
    n = c(33,109,92,22,80,98,104,470,460,501,103),
    nv1 = c(0,0,0,3,0,7,8,63,101,192,34),
    nv2 = c(0,0,0,0,0,0,2,22,32,32,4),
    week = weekseq,
    weeknum = as.vector(weeknum)
  ) 
processed <- pivot_longer(thedata,
                         cols = c("nv1","nv2"),
                         names_to = "variant",
                         values_to = "npos") 

# Construct models and CI for every variant

fitnv1 <- create_fit(
  filter(processed, variant == "nv1"),
  variant = "501Y.V1"
)
fitnv2 <- create_fit(
  filter(processed, variant == "nv2"),
  variant = "501Y.V2"
)
fits <- rbind(fitnv1, fitnv2)
fits$date <- as.POSIXct(fits$weeknum* attr(weeknum, "scaled:scale") + attr(weeknum, "scaled:center"), origin = "1970-01-01")

# Plot construction

pdata <- data.frame(
  variant = processed$variant,
  date = processed$week,
  Hmisc::binconf(processed$npos, processed$n)) %>%
  mutate(variant = factor(variant,
                          labels= c("501Y.V1", "501Y.V2")))



ggplot(pdata, aes(x = date, color = variant)) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper),
                width = 3e5,
                position = position_dodge2(width = 3e5)) +
  geom_point(aes(y = PointEst, color = variant),
             position = position_dodge2(width = 3e5))+
  geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL, fill = variant),
              data = fits,
              alpha = 0.2, linetype = 0) +
  geom_line(aes(y=prob), data = fits) + 
  labs(x = "Datum",y = "fractie",
       title = "Fractie per variant: data + logistisch model",
       subtitle = "CI: Wilson score interval",
       caption = "Data van Sciensano epidemiologisch bulletin 19/2/2021 - @JorisMeys")