### Analysis script 1

#Table 1, Table S1, S2 and S3

library(tidyverse)
library(rstudioapi)
library(MatchIt)
library(lubridate)
library(scales)
library(epitools)

rm(list = ls()) 
setwd(dirname(getActiveDocumentContext()$path))

##### all COVID records between 01-MAR-2020 and 31-DEC-2020
covid_all_cleaned = read_csv("covid_all_cleaned.csv")
included_tenants = read_csv("included_tenants.csv")

#age category

covid_all_cleaned %>% group_by(agecat) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 5 x 3
# agecat      n     freq
# <chr>   <int>    <dbl>
#   1 <25     59293 0.159   
# 2 >=65    92620 0.248   
# 3 25-49  130914 0.351   
# 4 50-64   90116 0.242   
# 5 NA         89 0.000239


#gender
covid_all_cleaned %>% group_by(gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 3 x 3
# gender      n    freq
# <chr>   <int>   <dbl>
#   1 F      201325 0.540  
# 2 M      170771 0.458  
# 3 NA        936 0.00251

#deceased
covid_all_cleaned %>% group_by(deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# deceased      n   freq
# <lgl>     <int>  <dbl>
#   1 FALSE    355582 0.953 
# 2 TRUE      17450 0.0468

#by region
covid_all_cleaned %>% group_by(zip_code)%>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 11 x 3
# zip_code     n   freq
# <chr>    <int>  <dbl>
#   1 0        12742 0.0342
# 2 1        14898 0.0399
# 3 2        63631 0.171 
# 4 3        68043 0.182 
# 5 4        10704 0.0287
# 6 5        12056 0.0323
# 7 6        28363 0.0760
# 8 7        57538 0.154 
# 9 8        63053 0.169 
# 10 9        35407 0.0949
# 11 NA        6597 0.0177 

#### 50+ in included_tenants for recombinant adjuvanted varicella zoster vaccine analysis
#effect of vaccine reducing deaths
covid_zoster_data_cleaned = read_csv("covid_zoster_data_cleaned.csv")

####Table 1
covid_zoster_data_cleaned %>% group_by(gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 3 x 3
# gender     n    freq
# <chr>  <int>   <dbl>
#   1 F      87428 0.517  
# 2 M      81155 0.480  
# 3 NA       423 0.00250

covid_zoster_data_cleaned %>% group_by(recombzostervax) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# recombzostervax      n   freq
# <dbl>  <int>  <dbl>
#   1               0 161739 0.957 
# 2               1   7267 0.0430

covid_zoster_data_cleaned %>% group_by(recombzostervax, gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 6 x 4
# # Groups:   recombzostervax [2]
# recombzostervax gender     n    freq
# <dbl> <chr>  <int>   <dbl>
#   1               0 F      83181 0.514  
# 2               0 M      78144 0.483  
# 3               0 NA       414 0.00256
# 4               1 F       4247 0.584  
# 5               1 M       3011 0.414  
# 6               1 NA         9 0.00124

covid_zoster_data_cleaned %>% group_by(deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# deceased      n   freq
# <chr>     <int>  <dbl>
#   1 0        153701 0.909 
# 2 1         15305 0.0906

covid_zoster_data_cleaned %>% group_by(emergency_inpatient) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# emergency_inpatient      n  freq
# <chr>                <int> <dbl>
#   1 0                    30894 0.183
# 2 1                   138112 0.817

covid_zoster_data_cleaned %>% group_by(fluvx) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# fluvx      n  freq
# <dbl>  <int> <dbl>
#   1     0 123370 0.730
# 2     1  45636 0.270

#odds ratio: all records
result = covid_zoster_data_cleaned %>% group_by(recombzostervax, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$recombzostervax==0,]$n,
                            result[result$deceased==0 & result$recombzostervax==1,]$n,
                            result[result$deceased==1 & result$recombzostervax==0,]$n,
                            result[result$deceased==1 & result$recombzostervax==1,]$n), dim=c(2,2))
oddsratio(contingency_table)
# $data
# Outcome
# Predictor  Disease1 Disease2  Total
# Exposed1   146711    15028 161739
# Exposed2     6990      277   7267
# Total      153701    15305 169006
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower     upper
# Exposed1 1.0000000        NA        NA
# Exposed2 0.3870899 0.3421413 0.4360548
# 
# $p.value
# two-sided
# Predictor  midp.exact fisher.exact
# Exposed1         NA           NA
# Exposed2          0 3.360276e-70
# two-sided
# Predictor   chi.square
# Exposed1          NA
# Exposed2 4.34284e-57
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: history of influenza vx
result = covid_zoster_data_cleaned %>% filter(fluvx==1) %>% group_by(recombzostervax, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$recombzostervax==0,]$n,
                            result[result$deceased==0 & result$recombzostervax==1,]$n,
                            result[result$deceased==1 & result$recombzostervax==0,]$n,
                            result[result$deceased==1 & result$recombzostervax==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2 Total
# Exposed1    35657     4025 39682
# Exposed2     5725      229  5954
# Total       41382     4254 45636
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower     upper
# Exposed1 1.0000000        NA        NA
# Exposed2 0.3545937 0.3086642 0.4052784
# 
# $p.value
# two-sided
# Predictor  midp.exact fisher.exact   chi.square
# Exposed1         NA           NA           NA
# Exposed2          0  3.37521e-66 9.329746e-55
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: emergency department and inpatient only
result = covid_zoster_data_cleaned %>% filter(emergency_inpatient==1) %>% group_by(recombzostervax, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$recombzostervax==0,]$n,
                            result[result$deceased==0 & result$recombzostervax==1,]$n,
                            result[result$deceased==1 & result$recombzostervax==0,]$n,
                            result[result$deceased==1 & result$recombzostervax==1,]$n), dim=c(2,2))
oddsratio(contingency_table)
# $data
# Outcome
# Predictor  Disease1 Disease2  Total
# Exposed1   117885    14940 132825
# Exposed2     5011      276   5287
# Total      122896    15216 138112
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower    upper
# Exposed1 1.0000000        NA       NA
# Exposed2 0.4348523 0.3839198 0.490431
# 
# $p.value
# two-sided
# Predictor  midp.exact fisher.exact   chi.square
# Exposed1         NA           NA           NA
# Exposed2          0 3.384331e-51 6.982035e-43
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: Female gender
result = covid_zoster_data_cleaned %>% filter(gender=='F') %>% group_by(recombzostervax, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$recombzostervax==0,]$n,
                            result[result$deceased==0 & result$recombzostervax==1,]$n,
                            result[result$deceased==1 & result$recombzostervax==0,]$n,
                            result[result$deceased==1 & result$recombzostervax==1,]$n), dim=c(2,2))
oddsratio(contingency_table)



#### 25-64y in included_tenants for MMR vaccine analysis
#effect of vaccine reducing deaths
covid_mmr_data_cleaned = read_csv("covid_mmr_data_cleaned.csv")

####Table 1
covid_mmr_data_cleaned %>% group_by(gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 3 x 3
# gender      n    freq
# <chr>   <int>   <dbl>
#   1 F      113399 0.547  
# 2 M       93415 0.451  
# 3 NA        517 0.00249

covid_mmr_data_cleaned %>% group_by(mmr) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# mmr        n   freq
# <chr>  <int>  <dbl>
#   1 0     203683 0.982 
# 2 1       3648 0.0176

covid_mmr_data_cleaned %>% group_by(deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# deceased      n   freq
# <lgl>     <int>  <dbl>
#   1 FALSE    203106 0.980 
# 2 TRUE       4225 0.0204

covid_mmr_data_cleaned %>% group_by(emergency_inpatient) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# emergency_inpatient      n  freq
# <chr>                <int> <dbl>
#   1 0                    55915 0.270
# 2 1                   151416 0.730

covid_mmr_data_cleaned %>% group_by(fluvx) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 2 x 3
# fluvxpositive      n  freq
# <chr>          <int> <dbl>
#   1 0             167801 0.809
# 2 1              39530 0.191

#gender, vaccinated
covid_mmr_data_cleaned %>% group_by(mmr, gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
# # A tibble: 6 x 4
# # Groups:   mmr [2]
# mmr   gender      n    freq
# <chr> <chr>   <int>   <dbl>
#   1 0     F      109958 0.540  
# 2 0     M       93217 0.458  
# 3 0     NA        508 0.00249
# 4 1     F        3441 0.943  
# 5 1     M         198 0.0543 
# 6 1     NA          9 0.00247

#odds ratio: all records
result = covid_mmr_data_cleaned %>% group_by(mmr, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$mmr==0,]$n,
                            result[result$deceased==0 & result$mmr==1,]$n,
                            result[result$deceased==1 & result$mmr==0,]$n,
                            result[result$deceased==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2
# Exposed1   199470     4213
# Exposed2     3636       12
# Total      203106     4225
# Outcome
# Predictor   Total
# Exposed1 203683
# Exposed2   3648
# Total    207331
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate
# Exposed1 1.0000000
# Exposed2 0.1583893
# odds ratio with 95% C.I.
# Predictor       lower
# Exposed1         NA
# Exposed2 0.08453063
# odds ratio with 95% C.I.
# Predictor      upper
# Exposed1        NA
# Exposed2 0.2662038
# 
# $p.value
# two-sided
# Predictor  midp.exact
# Exposed1         NA
# Exposed2          0
# two-sided
# Predictor  fisher.exact
# Exposed1           NA
# Exposed2 2.714088e-19
# two-sided
# Predictor    chi.square
# Exposed1           NA
# Exposed2 1.704139e-13
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: with influenza vx history
result = covid_mmr_data_cleaned %>% filter(fluvx == 1) %>% group_by(mmr, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$mmr==0,]$n,
                            result[result$deceased==0 & result$mmr==1,]$n,
                            result[result$deceased==1 & result$mmr==0,]$n,
                            result[result$deceased==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2
# Exposed1    37280      956
# Exposed2     1287        7
# Total       38567      963
# Outcome
# Predictor  Total
# Exposed1 38236
# Exposed2  1294
# Total    39530
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate      lower
# Exposed1 1.0000000         NA
# Exposed2 0.2170211 0.09239537
# odds ratio with 95% C.I.
# Predictor      upper
# Exposed1        NA
# Exposed2 0.4220955
# 
# $p.value
# two-sided
# Predictor   midp.exact
# Exposed1          NA
# Exposed2 1.08413e-07
# two-sided
# Predictor  fisher.exact
# Exposed1           NA
# Exposed2 1.824564e-07
# two-sided
# Predictor    chi.square
# Exposed1           NA
# Exposed2 6.917281e-06
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: emergency department and inpatient only
result = covid_mmr_data_cleaned %>% filter(emergency_inpatient == 1) %>% group_by(mmr, deceased) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$deceased==0 & result$mmr==0,]$n,
                            result[result$deceased==0 & result$mmr==1,]$n,
                            result[result$deceased==1 & result$mmr==0,]$n,
                            result[result$deceased==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2
# Exposed1   143833     4198
# Exposed2     3373       12
# Total      147206     4210
# Outcome
# Predictor   Total
# Exposed1 148031
# Exposed2   3385
# Total    151416
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate      lower
# Exposed1 1.0000000         NA
# Exposed2 0.1235681 0.06592792
# odds ratio with 95% C.I.
# Predictor      upper
# Exposed1        NA
# Exposed2 0.2076793
# 
# $p.value
# two-sided
# Predictor  midp.exact
# Exposed1         NA
# Exposed2          0
# two-sided
# Predictor  fisher.exact
# Exposed1           NA
# Exposed2 5.126205e-27
# two-sided
# Predictor    chi.square
# Exposed1           NA
# Exposed2 3.880138e-18
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"


#### general population 50+ in included_tenants for recombinant adjuvanted varicella zoster vaccine analysis
#preventative effect of vaccine against COVID infection

general_over50_data_cleaned = read_csv("general_over50_data_cleaned.csv")

#basic demographic statistics
general_over50_data_cleaned %>% group_by(gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 3 x 3
# gender       n    freq
# <chr>    <int>   <dbl>
#   1 F      3170002 0.565  
# 2 M      2420788 0.431  
# 3 NA       19622 0.00350

general_over50_data_cleaned %>% group_by(recombzostervx) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# zoster       n   freq
# <dbl>   <int>  <dbl>
#   1        0 5312513 0.947 
# 2        1  297899 0.0531

general_over50_data_cleaned %>% group_by(gender, recombzostervx) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 6 x 4
# # Groups:   gender [3]
# gender zoster       n   freq
# <chr>     <dbl>   <int>  <dbl>
#   1 F             0 2992382 0.944 
# 2 F             1  177620 0.0560
# 3 M             0 2301058 0.951 
# 4 M             1  119730 0.0495
# 5 NA            0   19073 0.972 
# 6 NA            1     549 0.0280

general_over50_data_cleaned %>% group_by(fluvax) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# fluvax       n  freq
# <dbl>   <int> <dbl>
#   1      0 4257920 0.759
# 2      1 1352492 0.241

general_over50_data_cleaned %>% group_by(emergency_inpatient) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# emergency_inpatient       n  freq
# <dbl>   <int> <dbl>
#   1                   0 4186121 0.746
# 2                   1 1424291 0.254

#odds ratio: all individuals
result = general_over50_data_cleaned %>% group_by(recombzostervx, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$recombzostervx==0,]$n,
                            result[result$covid==0 & result$recombzostervx==1,]$n,
                            result[result$covid==1 & result$recombzostervx==0,]$n,
                            result[result$covid==1 & result$recombzostervx==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  5154217   158296 5312513
# Exposed2   290678     7221  297899
# Total     5444895   165517 5610412
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower     upper
# Exposed1 1.0000000        NA        NA
# Exposed2 0.8088859 0.7897165 0.8283479
# 
# $p.value
# two-sided
# Predictor  midp.exact fisher.exact   chi.square
# Exposed1         NA           NA           NA
# Exposed2          0 4.927176e-72 3.918704e-68
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: individuals with influenza vx history
result = general_over50_data_cleaned %>% filter(fluvax==1) %>% group_by(recombzostervx, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$recombzostervx==0,]$n,
                            result[result$covid==0 & result$recombzostervx==1,]$n,
                            result[result$covid==1 & result$recombzostervx==0,]$n,
                            result[result$covid==1 & result$recombzostervx==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  1077145    38970 1116115
# Exposed2   230463     5914  236377
# Total     1307608    44884 1352492
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower     upper
# Exposed1 1.0000000        NA        NA
# Exposed2 0.7093085 0.6898552 0.7291892
# 
# $p.value
# two-sided
# Predictor  midp.exact  fisher.exact    chi.square
# Exposed1         NA            NA            NA
# Exposed2          0 9.837945e-141 1.647067e-131
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: emergency department and inpatient only
result = general_over50_data_cleaned %>% filter(emergency_inpatient==1) %>% group_by(zoster, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$zoster==0,]$n,
                            result[result$covid==0 & result$zoster==1,]$n,
                            result[result$covid==1 & result$zoster==0,]$n,
                            result[result$covid==1 & result$zoster==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  1271553   111387 1382940
# Exposed2    38057     3294   41351
# Total     1309610   114681 1424291
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower    upper
# Exposed1 1.0000000        NA       NA
# Exposed2 0.9881173 0.9529026 1.024258
# 
# $p.value
# two-sided
# Predictor  midp.exact fisher.exact chi.square
# Exposed1         NA           NA         NA
# Exposed2  0.5159057    0.5209033  0.5149889
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#### general population 25-64y in included_tenants for MMR vaccine analysis
#preventative effect of vaccine against COVID infection

general_25_64_data_cleaned = read_csv("general_25_64_data_cleaned.csv")

#basic demographic statistics
general_25_64_data_cleaned %>% group_by(gender) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 3 x 3
# gender       n    freq
# <chr>    <int>   <dbl>
#   1 F      3776648 0.587  
# 2 M      2635151 0.410  
# 3 NA       21422 0.00333

general_25_64_data_cleaned %>% group_by(mmr) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# mmr       n   freq
# <dbl>   <int>  <dbl>
#   1     0 6306799 0.980 
# 2     1  126422 0.0197


general_25_64_data_cleaned %>% group_by(gender, mmr) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 6 x 4
# # Groups:   gender [3]
# gender   mmr       n    freq
# <chr>  <dbl>   <int>   <dbl>
#   1 F          0 3657677 0.968  
# 2 F          1  118971 0.0315 
# 3 M          0 2627954 0.997  
# 4 M          1    7197 0.00273
# 5 NA         0   21168 0.988  
# 6 NA         1     254 0.0119 

general_25_64_data_cleaned %>% group_by(fluvax) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# fluvax       n  freq
# <dbl>   <int> <dbl>
#   1      0 5337852 0.830
# 2      1 1095369 0.170

general_25_64_data_cleaned %>% group_by(emergency_inpatient) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))

# # A tibble: 2 x 3
# emergency_inpatient       n  freq
# <dbl>   <int> <dbl>
#   1                   0 4489263 0.698
# 2                   1 1943958 0.302

#odds ratio: all individuals
result = general_25_64_data_cleaned %>% group_by(mmr, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$mmr==0,]$n,
                            result[result$covid==0 & result$mmr==1,]$n,
                            result[result$covid==1 & result$mmr==0,]$n,
                            result[result$covid==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  6106969   199830 6306799
# Exposed2   122786     3636  126422
# Total     6229755   203466 6433221
# 
# $measure
# odds ratio with 95% C.I.
# Predictor  estimate     lower    upper
# Exposed1 1.000000        NA       NA
# Exposed2 0.905017 0.8752432 0.935481
# 
# $p.value
# two-sided
# Predictor    midp.exact fisher.exact  chi.square
# Exposed1           NA           NA          NA
# Exposed2 2.416225e-09 2.558455e-09 4.05127e-09
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: individuals with influenza vaccination history
result = general_25_64_data_cleaned %>% filter(fluvax == 1) %>% group_by(mmr, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$mmr==0,]$n,
                            result[result$covid==0 & result$mmr==1,]$n,
                            result[result$covid==1 & result$mmr==0,]$n,
                            result[result$covid==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  1017351    37819 1055170
# Exposed2    38913     1286   40199
# Total     1056264    39105 1095369
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower    upper
# Exposed1 1.0000000        NA       NA
# Exposed2 0.8891184 0.8398887 0.940359
# 
# $p.value
# two-sided
# Predictor    midp.exact fisher.exact   chi.square
# Exposed1           NA           NA           NA
# Exposed2 3.350457e-05 3.526927e-05 4.425126e-05
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: emergency department and inpatient only
result = general_25_64_data_cleaned %>% filter(emergency_inpatient == 1) %>% group_by(mmr, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$mmr==0,]$n,
                            result[result$covid==0 & result$mmr==1,]$n,
                            result[result$covid==1 & result$mmr==0,]$n,
                            result[result$covid==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  1760012   114353 1874365
# Exposed2    66929     2664   69593
# Total     1826941   117017 1943958
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower     upper
# Exposed1 1.0000000        NA        NA
# Exposed2 0.6126501 0.5889795 0.6369899
# 
# $p.value
# two-sided
# Predictor  midp.exact  fisher.exact    chi.square
# Exposed1         NA            NA            NA
# Exposed2          0 3.701622e-153 2.785317e-135
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"

#odds ratio: female only
result = general_25_64_data_cleaned %>% filter(gender == 'F') %>% group_by(mmr, covid) %>% summarise(n=n()) %>% mutate(freq = n / sum(n))
contingency_table = array(c(result[result$covid==0 & result$mmr==0,]$n,
                            result[result$covid==0 & result$mmr==1,]$n,
                            result[result$covid==1 & result$mmr==0,]$n,
                            result[result$covid==1 & result$mmr==1,]$n), dim=c(2,2))
oddsratio(contingency_table)

# $data
# Outcome
# Predictor  Disease1 Disease2   Total
# Exposed1  3549572   108105 3657677
# Exposed2   115541     3430  118971
# Total     3665113   111535 3776648
# 
# $measure
# odds ratio with 95% C.I.
# Predictor   estimate     lower    upper
# Exposed1 1.0000000        NA       NA
# Exposed2 0.9747743 0.9415811 1.008818
# 
# $p.value
# two-sided
# Predictor  midp.exact fisher.exact chi.square
# Exposed1         NA           NA         NA
# Exposed2  0.1452043    0.1486363   0.145989
# 
# $correction
# [1] FALSE
# 
# attr(,"method")
# [1] "median-unbiased estimate & mid-p exact CI"