if(!require("pacman")) install.packages("pacman")

pacman::p_load(broom,epiR,knitr,tidyverse)


## Primray Efficacy ##
rate_rr_1 <- as.table(matrix(c(85,85/11.74,167,167/23.64), nrow = 2, byrow = TRUE))
rr_relativeraterisk <- epi.2by2(dat = rate_rr_1, method = "cohort.time", conf.level = 0.95, units = 100, outcome = "as.columns")
rr_relativeraterisk
##              Outcome +    Time at risk        Inc rate *
## Exposed +           85            7.24              1174
## Exposed -          167            7.06              2364
## Total              252           14.30              1762
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc rate ratio                               0.50 (0.38, 0.65)
## Attrib rate *                                -1190.00 (-1626.85, -753.15)
## Attrib rate in population *                  -602.32 (-1021.67, -182.96)
## Attrib fraction in exposed (%)               -101.36 (-164.64, -54.17)
## Attrib fraction in population (%)            -34.19 (-38.02, -30.19)
## -------------------------------------------------------------------
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 units of population time at risk
tidy(rr_relativeraterisk, parameters = "moa") %>% kable
term estimate conf.low conf.high
IRR.strata.wald 0.4966159 0.3778774 0.6486139
ARate.strata.wald -1190.0000000 -1626.8524487 -753.1475513
PARate.strata.wald -602.3168882 -1021.6742229 -182.9595536
AFRate.strata.wald -1.0136286 -1.6463606 -0.5417493
PAFRate.strata.wald -0.3418985 -0.3802131 -0.3018822
tidy(rr_relativeraterisk, parameters = "stat") %>% kable
term statistic df p.value
chisq.strata 1.704546 1 0.1916945
ve_relativeraterisk <- 1 -  rr_relativeraterisk[["massoc"]][["IRR.strata.wald"]]
ve_relativeraterisk[c(1,3,2)]%>% kable
est upper lower
0.5033841 0.3513861 0.6221226
## Secondary Efficacy ##
rate_rr_2 <- as.table(matrix(c(7,7/0.97,31,31/4.39), nrow = 2, byrow = TRUE))
rr_relativeraterisk_2 <- epi.2by2(dat = rate_rr_2, method = "cohort.time", conf.level = 0.95, units = 100, outcome = "as.columns")
rr_relativeraterisk_2
##              Outcome +    Time at risk        Inc rate *
## Exposed +            7            7.22                97
## Exposed -           31            7.06               439
## Total               38           14.28               266
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc rate ratio                               0.22 (0.08, 0.51)
## Attrib rate *                                -342.00 (-512.43, -171.57)
## Attrib rate in population *                  -172.86 (-349.04, 3.33)
## Attrib fraction in exposed (%)               -352.58 (-1117.69, -95.52)
## Attrib fraction in population (%)            -64.95 (-70.58, -58.37)
## -------------------------------------------------------------------
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 units of population time at risk
tidy(rr_relativeraterisk_2, parameters = "moa") %>% kable
term estimate conf.low conf.high
IRR.strata.wald 0.2209567 0.0821224 0.5114515
ARate.strata.wald -342.0000000 -512.4261582 -171.5738418
PARate.strata.wald -172.8562500 -349.0440343 3.3315343
AFRate.strata.wald -3.5257732 -11.1769449 -0.9552197
PAFRate.strata.wald -0.6494845 -0.7057766 -0.5837344
tidy(rr_relativeraterisk_2, parameters = "stat") %>% kable
term statistic df p.value
chisq.strata 5.408648 1 0.0200372
ve_relativeraterisk_2 <- 1 - rr_relativeraterisk_2[["massoc"]][["IRR.strata.wald"]]
ve_relativeraterisk_2[c(1,3,2)] %>% kable
est upper lower
0.7790433 0.4885485 0.9178776

Source: Otavio Ranzani - otavio_ranzani