###################################################
## R syntax for superspreading potential        ###
## ONEHEALTH LAB, 5 June 2025, onehealth.or.kr  ###
## Catholic University of Korea                 ###
###################################################

## Instrall MASS packages ## 
install.packages("MASS", dependencies = FALSE)

# offspring distribution
offspring <- c(rep(0,199),rep(1,57),rep(2,18),rep(3,6),rep(4,5),5,rep(6,2),10,11)
offspring
# frequency table
table(offspring)

# plot histogram
hist(offspring)
hist(offspring, right = F)

# fitting negative binomial distribution
library(MASS)
fit <- fitdistr(offspring, "negative binomial")
fit

# plot PDF
hist(offspring, prob=T, col="gray",right = F)
xvalue <- 0:max(offspring)
yvalue <- dnbinom(xvalue, 
                  size = fit$estimate["size"], 
                  mu = fit$estimate["mu"])
lines(xvalue,yvalue)

# 95% CI for R0
# Estimate 95th percentile of size and mu   
conf_interval <- confint(fit, level = 0.95)
print(conf_interval)