covid_data.R

library(ggplot2)
library(scales)
library(dplyr)

################## parameters  ##############
all <- 10000000 # population
sick <- 200000  # number of sick people

sensitivity <- 0.95   # how many SICK people will be detected correctly
specificity <- 0.95   # how many HEALTHY people will be detected correctly

#############################################

healthy <- all-sick

sick.positive <- sensitivity * sick
sick.negative <- sick - sick.positive

healthy.positive <- healthy * (1-specificity)
healthy.negative <- healthy - healthy.positive

data <- data.frame(
  list(
    reality = c("HEALTHY","HEALTHY","SICK","SICK"),
    detected = c("POSITIVE","NEGATIVE", "POSITIVE","NEGATIVE"),
    correct = c("NO","YES","YES","NO"),
    people = c(healthy.positive, healthy.negative, sick.positive, sick.negative)
  )
)

data

data %>% group_by(reality) %>% summarise( people = sum(people)) %>%
  ggplot(aes(x="REALITY",y=people,fill=reality))+
  geom_bar(stat="identity")+
  scale_fill_manual(values=c("lightgray","red"))+
  labs(x="",y="Number od people", fill="")+
  scale_y_continuous(labels=comma)+theme_bw()

ggplot(data, aes(x=reality,y=detected))+
  geom_point(aes(col=correct,fill=correct,size=people), shape=22)+
  geom_text(aes(label=format(people,format="d", big.mark = " ")), vjust=3)+
  scale_size_continuous(guide="none", range=c(5,60))+
  scale_fill_manual(values=c("pink","lightblue"), guide="none")+
  labs(x="Reality", y="Detected using PCR test", fill="correct detection")+
  theme_bw()


############  probability of at least one sick in group ###############

p <- sick / all
n <- seq(1,30)
p_fail <- 1-((1-p)^n)
data_fail <- data.frame(
  list(
    n = n,
    p_fail = p_fail
  )
)

data_fail$label=""
data_fail[data_fail$n<15,]$label=sprintf("%.1f%%",data_fail[data_fail$n<15,]$p_fail * 100)

ggplot(data_fail, aes(x=n,y=p_fail * 100))+
  geom_line(col="red")+
  geom_point(col="red")+  
  geom_text(aes(label=label), size=2.5, hjust=-.5)+
  scale_y_continuous(breaks=seq(0,100,5))+
  scale_x_continuous(breaks=seq(1,1000,1), limits=c(1,30))+
  labs(x="Number of people",y="Probability at least one is sick [%]")+
  theme_bw()