Vaccine Infection Prevention Efficacy via SIMULATION.

vipe.R

#!/usr/bin/Rscript --no-init-file

# Vaccine Infection Prevention Efficacy
#
# Eficácia da vacina na prevenção de infecção no contato com vacinados.
#
# Extraído do artigo "Estimating the probability that a vaccinated person still
# infects others with Covid-19" de 19/01/2021, disponível no R-blogger.
#

library(dtplyr)
library(dplyr)
library(magrittr)

line.separator <- ifelse(checkmate::testOS("windows"), "\r\n", "\n")

PRINT <- function (...) {
  for (arg in list(...)) if (is.character(arg)) cat(arg) else print(arg)
}

##  ROTINA DE SIMULAÇÃO

simulate.meeting.data <- function(
  init.sick.share = 0.3,  # share of initially sick people
  infect.prob = 0.1,      # probability of infection
  vipe = 0.7,             # vaccine infection prevention efficacy
  vacc.share = 0.3,       # share of vaccinated people
  prob.unobserved = 0,    # probability that a meeting is unobserved
  n = 100000,             # number of persons
  m = 100000              # number of meetings
) {
  # Draw persons that meet
  m1 = sample.int(n, 2*m, replace=TRUE)
  m2 = sample.int(n, 2*m, replace=TRUE)

  # Only keep first m meetings
  # where different persons meet
  keep = m1 != m2
  m1 = m1[keep][1:m]
  m2 = m2[keep][1:m]

  # Which persons are vaccinated
  vacc = (runif(n) <= vacc.share)

  # Which persons are initially infected
  # for vaccinated persons that may be the hypothetical
  # infection if they would not be vaccinated
  start.sick = (runif(n) <= init.sick.share)

  # Simulate which meetings would transfer an infection
  infect.dice = runif(m)
  infect1.threshold = infect.prob*ifelse(vacc[m2], 1-vipe,1)
  infect2.threshold = infect.prob*ifelse(vacc[m1], 1-vipe,1)

  infected1 = (infect.dice <= infect1.threshold) & start.sick[m2]
  infected2 = (infect.dice <= infect2.threshold) & start.sick[m1]

  # Which subjects are sick after all infections
  post.sick1 = m1[infected1]
  post.sick2 = m2[infected2]
  post.sick = start.sick
  post.sick[c(post.sick1,post.sick2)] = TRUE

  # Is the meeting recorded / observed
  record = runif(m) > prob.unobserved

  # in every meeting let s be sender
  # and r be the receiver of a potential infection
  dat = bind_rows(
    tibble(s = m1, r = m2, record=record),
    tibble(s = m2, r = m1, record=record)
  ) %>%
  filter(record) %>%
  mutate(
    start.sick.s = start.sick[s],
    start.sick.r = start.sick[r],
    post.sick.s = post.sick[s],
    post.sick.r = post.sick[r],
    infected = !start.sick.r & post.sick.r,
    vacc.s = vacc[s],
    vacc.r = vacc[r]
  ) %>%
  # only keep people that are initially sick
  filter(!start.sick.r) %>%
  # aggregate meeting data on person level
  lazy_dt() %>%
  group_by(r) %>%
  summarize(
    infected = 1L*first(infected),
    m = n(),
    ms = sum(start.sick.s & !vacc.s),
    mv = sum(vacc.s)
  ) %>%
  rename(subjectid = r) %>%
  as_tibble()

  dat   # return value
}

##  PARÂMETROS ARBITRÁRIOS DA SIMULAÇÃO

ISH <- 0.1                  # init.sick.share

NUMBER.PERSONS <- 500000    # número de pessoas a simular
MEETINGS.PER.PERSON <- 2    # número de encontros por pessoa

##  SIMULAÇÃO DOS DADOS

set.seed(1)
dat <- simulate.meeting.data(
  init.sick.share = ISH,                  # percentual de infectados no início
  infect.prob = 0.1,                      # Pr(infecção entre pessoas)
  vipe = 0.7,                             # eficácia da vacina na prevenção …
  vacc.share = 0.5,                       # taxa de vacinação
  prob.unobserved = 0.1,                  # Pr(não observar encontro)
  n = NUMBER.PERSONS,                     # tamanho da população simulada
  m = NUMBER.PERSONS*MEETINGS.PER.PERSON  # número de encontros da população
)
PRINT("=== Dados Simulados", line.separator, line.separator, head(dat))

##  TRANSFORMAÇÃO DOS DADOS SIMULADOS

agg <- dat %>%
group_by(mv, ms) %>%
summarize(
  n = n(),
  prob.infected = mean(infected),
  y = log(1-prob.infected)
) %>%
ungroup()
PRINT(line.separator, "=== Dados Adequados", line.separator, line.separator,
  head(agg))

##  MODELO LINEAR

agg <- filter(agg, is.finite(y))
reg <- lm(y~mv+ms, data=agg, weights=agg$n)
beta <- coef(reg)
PRINT(line.separator, "=== Modelo Linear", line.separator, line.separator,
  formula(reg), line.separator, beta, line.separator, anova(reg))

##  ESTIMATIVA DAS PROBABILIDADES

names(beta) <- NULL
pi.v <- 1-exp(beta[2])          # Pr(infecção no contato com vacinado)
pi.s <- 1-exp(beta[3])          # Pr(infecção no contato com não vacinado)
vipe <- 1-(pi.v / ISH / pi.s)   # Vaccine Infection Prevention Efficacy

PRINT(line.separator, "=== Estimativa das Probabilidades e Eficácia",
  line.separator, line.separator,
  c('pi.s(%)'=pi.s, 'pi.v(%)'=pi.v, 'vipe(%)'=vipe)*100)

Contains the code from Estimating the probability that a vaccinated person still infects others with Covid-19 in 2021/01/19, available at R-blogger.

Warning: THIS IS A SIMULATION!