---
title: "figure_17_12_24"
format: html
editor: visual
---

```{r}
library(tidyverse)
library(Epi) 
library(popEpi) 
library(survival) 
library(mgcv) 
library(modelr)
library(mice)
library(marginaleffects)
library(patchwork)
library(tidycmprsk)
library(ggsurvfit)
```

read in mace data, with all-cause death as competing event.



```{r}
d1 <- read_csv("d1_scrambled.csv")
d1$vac <- if_else(d1$vaccinated == "yes", 1, 0)
d1$sex <- as.factor(d1$sex)
#recode categorical variables as factors (for the mgcv::gam function)
d1$vac <- if_else(d1$vaccinated == "yes", 1, 0)
d1$education <- as.factor(d1$education)
d1$heart <- as.factor(d1$heart)
d1$lung <- as.factor(d1$lung)
d1$diabetes <- as.factor(d1$diabetes)
d1$vaccinated <- as.factor(d1$vaccinated)
```


# Figure 1 was not made in R.

# Figure 2.

MACE rates for the vaccinated and the unvaccinated over the study follow-up period. The vaccinated rates are in green and the unvaccinated are in red. Shaded areas represent 95% confidence intervals. A. MACE race for the full study population. B.-E. MACE rates for populations stratified by sex and age. B. \< 70-year-old males; C. \< 70-year-old females; D. 70+ males; E. 70+ females. The X-axis shows time from baseline (T0) in days, and the Y-axis plots "instantaneous" mace rates per day, which were transformed for clarity to 100-year periods (each day's rates were multiplied by 36500).

```{r}
#make lexis object
lex <- Lexis(exit = list(tfe = time), 
            exit.status = factor(status, labels = c("Alive", "MACE", "Dead")),
             data = d1)
#split the lexis object into unequal slices of time 
sL <- splitMulti(lex, tfe = c(0, 10, 20, 30, 40, 60, 80, 100, 120, 150, 200, 250, 300, 350, 360, 365))

sL$vaccinated <- as.factor(sL$vaccinated)

#calculate IP weights from a logistic regression fit. data encompasses full study population 
m_logistic2 <- gam(vac ~ education + s(age, by=sex) + sex + tongue * county + diabetes * heart*sex + lung + cal_time1, data = d1 , family = binomial)
dat2 <- predictions(m_logistic2, newdata = d1) 
#predicts probability of being vaccinated
dat2$wts <- ifelse(dat2$vac == 1, 1 / dat2$estimate, 1 / (1 - dat2$estimate))
#calculates IP weights
d1$wts <- dat2$wts

wghts <- d1 %>% select(id, wts)
sL <- sL %>% left_join(wghts)
#adds IP weights to the long-form table used for survival modelling (the marginal structural model)

k5_mace_ipw <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = sL, weights = wts)
#the MSM

#predict rates over time
pr_ipw <- predictions(
    k5_mace_ipw,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")

#transform rates from "per day" to "per 100 years"
pr_ipw <- as.data.frame(pr_ipw) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

#plot transformed rates
fig_mace_all <- ggplot(pr_ipw, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype = vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab("wIR")+
  xlab("Days")+
  #ylim(0, 8e-05)+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

#separate male's and female's data
sL_male <- sL %>% filter(sex =="M")
sL_female <- sL %>% filter(sex =="N")
d1_male <- d1 %>% filter(sex =="M")
d1_female <- d1 %>% filter(sex =="N")
#younger males
d_male_young <- sL_male %>% filter(age <70)
d1_male_young <- d1_male %>% filter(age <70)
#again IP weights
m_log_male_young <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data =  d1_male_young, family = binomial)
dat_m <- predictions(m_log_male_young, newdata = d1_male_young)
dat_m$wts <- ifelse(dat_m$vac == 1, 1 / dat_m$estimate, 1 / (1 - dat_m$estimate))
d1_male_young$wts <- dat_m$wts

wghts_male_young <- d1_male_young %>% select(id, wts)
d_male_young <- d_male_young %>% left_join(wghts_male_young)

k5_mace_M <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = d_male_young, weights = wts)

pr_ipw_M <- predictions(
    k5_mace_M,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_M <- as.data.frame(pr_ipw_M) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)
mace_y_M <- ggplot(pr_ipw_M, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype = vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  xlab(NULL)+
  ggtitle("<70 males")+
 # xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

#younger females
d_female_young <- sL_female %>% filter(age <70)
d1_female_young <- d1_female %>% filter(age <70)
m_log_female_young <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data =  d1_female_young, family = binomial)
dat_n <- predictions(m_log_female_young, newdata = d1_female_young)
dat_n$wts <- ifelse(dat_n$vac == 1, 1 / dat_n$estimate, 1 / (1 - dat_n$estimate))
d1_female_young$wts <- dat_n$wts

wghts_female_young <- d1_female_young %>% select(id, wts)
d_female_young <- d_female_young %>% left_join(wghts_female_young)

mace_N <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated), 
             family = poisreg, data = d_female_young, weights = wts)

pr_ipw_N <- predictions(
    mace_N,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_N <- as.data.frame(pr_ipw_N) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

mace_y_F <- ggplot(pr_ipw_N, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype = vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  xlab(NULL)+
  ggtitle("<70 females")+
  #xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

#older females
d_female_old <- sL_female %>% filter(age >=70)
d1_female_old <- d1_female %>% filter(age >=70)

m_log_female_old <- gam(vac ~ education + s(age) + tongue + county  + diabetes * heart + lung + cal_time1, data =  d1_female_old, family = binomial)
dat_n2 <- predictions(m_log_female_old, newdata = d1_female_old)
dat_n2$wts <- ifelse(dat_n2$vac == 1, 1 / dat_n2$estimate, 1 / (1 - dat_n2$estimate))
d1_female_old$wts1 <- dat_n2$wts

wghts_female_old <- d1_female_old %>% select(id, wts1)
d_female_old <- d_female_old %>% left_join(wghts_female_old)

k5_mace_N3 <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = d_female_old, weights = wts1)

pr_ipw_N3 <- predictions(
    k5_mace_N3,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts1")
pr_ipw_N3 <- as.data.frame(pr_ipw_N3) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

mace_old_N <- ggplot(pr_ipw_N3, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype = vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  xlab(NULL)+
  ggtitle("70+ females")+
  #xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

#older males
d_male_old <- sL_male %>% filter(age >=70)
d1_male_old <- d1_male %>% filter(age >=70)
m_log_male_old <- gam(vac ~ education + s(age) + tongue + county  + diabetes * heart + lung + cal_time1, data =  d1_male_old, family = binomial)
dat_m2 <- predictions(m_log_male_old, newdata = d1_male_old)
dat_m2$wts <- ifelse(dat_m2$vac == 1, 1 / dat_m2$estimate, 1 / (1 - dat_m2$estimate))
d1_male_old$wts <- dat_m2$wts

wghts_male_old <- d1_male_old %>% select(id, wts)
d_male_old <- d_male_old %>% left_join(wghts_male_old)

k5_mace_M2 <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = d_male_old, weights = wts)

pr_ipw_M2 <- predictions(
    k5_mace_M2,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_M2 <- as.data.frame(pr_ipw_M2) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)
mace_old_M <- ggplot(pr_ipw_M2, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype = vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  xlab(NULL)+
  ggtitle("70+ males")+
  #xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

#making the composite figure with the patchwork packege
fig_mace_all /(mace_y_M| mace_y_F)/ (mace_old_M| mace_old_N) + plot_annotation(tag_levels = 'A')
ggsave("Fig2y.pdf", units = "cm", width = 12, height = 14) 
```

# Fig 3.

All-cause death rates rates for the vaccinated and the unvaccinated over the study follow-up period. The vaccinated rates are in green and the unvaccinated are in red. Shaded areas represent 95% confidence intervals. A. MACE race for the full study population. B.-E. MACE rates for populations stratified by sex and age. B. \< 70-year-old males; C. \< 70-year-old females; D. 70+ males; E. 70+ females. The X-axis shows time from baseline (T0) in days, and the Y-axis plots "instantaneous" mace rates per day, which were transformed for clarity to 100-year periods (each day's rates were multiplied by 36500).

```{r}
#read in death data (mace is not a competing risk)
dd <- read_csv("deaths_scrambled.csv")
#see Fig 2 for code comments
lex <- Lexis(exit = list(tfe = time), 
            exit.status = factor(status, labels = c("Alive", "Dead")),
             data = dd)
sL <- splitMulti(lex, tfe = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 145, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 365))
sL$vac <- if_else(sL$vaccinated == "yes", 1, 0)
sL$vaccinated <- as.factor(sL$vaccinated)

m_logistic <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data = dd, family = binomial)
dat2 <- predictions(m_logistic, newdata = dd)
dat2$wts <- ifelse(dat2$vac == 1, 1 / dat2$estimate, 1 / (1 - dat2$estimate))
dd$wts <- dat2$wts

wts_all <- dd %>% select(id, wts)
sL <- sL %>% left_join(wts_all)

k5_death <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
                family = poisreg, data = sL, weights = wts)

pred_ipw_d <- predictions(
    k5_death,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pred_ipw_d <- as.data.frame(pred_ipw_d) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

fig_death_all <- ggplot(pred_ipw_d, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab("wIR")+
  xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen")) 

sL_male <- sL %>% filter(sex =="M")
sL_female <- sL %>% filter(sex =="N")
dd_male <- dd %>% filter(sex =="M")
dd_female <- dd %>% filter(sex =="N")
d_male_young <- sL_male %>% filter(age <70)
d1_male_young <- dd_male %>% filter(age <70)

m_log_male_young <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data =  d1_male_young, family = binomial)
dat_m <- predictions(m_log_male_young, newdata = d1_male_young)
dat_m$wts <- ifelse(dat_m$vac == 1, 1 / dat_m$estimate, 1 / (1 - dat_m$estimate))
d1_male_young$wts <- dat_m$wts
wts_m_y <- d1_male_young %>% select(id, wts)
d_male_young <- d_male_young %>% left_join(wts_m_y)

k5_mace_M <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = d_male_young, weights = wts)

pr_ipw_M <- predictions(
    k5_mace_M,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_M <- as.data.frame(pr_ipw_M) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)
fig_death_Y_M <- ggplot(pr_ipw_M, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  ggtitle("<70 males")+
  xlab(NULL)+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

d_female_young <- sL_female %>% filter(age <70)
d1_female_young <- dd_female %>% filter(age <70)

m_log_female_young <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data =  d1_female_young, family = binomial)
dat_n <- predictions(m_log_female_young, newdata = d1_female_young)
dat_n$wts <- ifelse(dat_n$vac == 1, 1 / dat_n$estimate, 1 / (1 - dat_n$estimate))
d1_female_young$wts <- dat_n$wts

wts_f_y <- d1_female_young %>% select(id, wts)
d_female_young <- d_female_young %>% left_join(wts_f_y)

mace_N <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=6), 
             family = poisreg, data = d_female_young, weights = wts)

pr_ipw_N <- predictions(
    mace_N,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_N <- as.data.frame(pr_ipw_N) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)
fig_death_Y_N <- ggplot(pr_ipw_N, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  ggtitle("<70 females")+
  xlab(NULL)+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

d_female_old <- sL_female %>% filter(age >=70)
d1_female_old <- dd_female %>% filter(age >=70)

m_log_female_old <- gam(vac ~ education + s(age) + tongue + county  + diabetes * heart + lung + cal_time1, data =  d1_female_old, family = binomial)
dat_n2 <- predictions(m_log_female_old, newdata = d1_female_old)
dat_n2$wts <- ifelse(dat_n2$vac == 1, 1 / dat_n2$estimate, 1 / (1 - dat_n2$estimate))
d1_female_old$wts <- dat_n2$wts

wts_f_o <- d1_female_old %>% select(id, wts)
d_female_old <- d_female_old %>% left_join(wts_f_o)

k5_mace_N3 <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = d_female_old, weights = wts)

pr_ipw_N3 <- predictions(
    k5_mace_N3,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_N3 <- as.data.frame(pr_ipw_N3) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)
fig_death_old_N <- ggplot(pr_ipw_N3, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  ggtitle("70+ females")+
  xlab(NULL)+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

d_male_old <- sL_male %>% filter(age >=70)
d1_male_old <- dd_male %>% filter(age >=70)

m_log_male_old <- gam(vac ~ education + s(age) + tongue + county  + diabetes * heart + lung + cal_time1, data =  d1_male_old, family = binomial)
dat_m2 <- predictions(m_log_male_old, newdata = d1_male_old)
dat_m2$wts <- ifelse(dat_m2$vac == 1, 1 / dat_m2$estimate, 1 / (1 - dat_m2$estimate))
d1_male_old$wts <- dat_m2$wts

wts_m_o <- d1_male_old %>% select(id, wts)
d_male_old <- d_male_old %>% left_join(wts_m_o)

k5_mace_M2 <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = d_male_old, weights = wts)

pr_ipw_M2 <- predictions(
    k5_mace_M2,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_ipw_M2 <- as.data.frame(pr_ipw_M2) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)
fig_death_old_M <- ggplot(pr_ipw_M2, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab(NULL)+
  ggtitle("70+ males")+
  xlab(NULL)+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

fig_death_all /(fig_death_Y_M| fig_death_Y_N)/ (fig_death_old_M| fig_death_old_N) + 
  plot_annotation(tag_levels = 'A')
ggsave("Fig3y.pdf", units = "cm", width = 12, height = 14) 
```

# Supplementary figures

# Figure S1 DAG, was not done in R.

# Fig S2.

Standardized mean differences of the IPW model counfounders calculated on the actual dataset (red dots) and on the virtual dataset obtained by implementing IPW weights (blue dots). The range typeically seen as not significally different from zero effect, from -0.1 to 0.1, is denoted by dotted lines. Upper panel -- the weights are calculated from full data; middle panel -- data is restricted to males; lower panel -- data is restricted to females. The tidysmd R package was used to calculate and plot the smd-s.

```{r}
library(tidysmd)

m_logistic2 <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung, data = d1, family = binomial)
dat2 <- predictions(m_logistic2, newdata = d1)
wts <- ifelse(dat2$vac == 1, 1 / dat2$estimate, 1 / (1 - dat2$estimate))
d1$wts <- wts

fig_smd_full <- tidy_smd(d1, c(age, sex, tongue, education, heart, lung, diabetes, county), .group = vaccinated, .wts = wts, na.rm = TRUE) %>% ggplot(aes(smd, variable, color=method))+ 
  geom_point(size=3)+
  geom_vline(xintercept = c(-0.1, 0.1), linetype=2)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  ggtitle("full data")+
  ylab(NULL)+
  xlab(NULL)+
  xlim(-0.2, 0.43)

fig_smd_full <- tidy_smd(d1, c(age, sex, tongue, education, heart, lung, diabetes, county), .group = vaccinated, .wts = wts, na.rm = TRUE) %>% 
  ggplot(aes(smd, variable, shape = method))+ 
  scale_shape_discrete(solid=F) +
   geom_point(size=3)+
  geom_vline(xintercept = c(-0.1, 0.1), linetype=2)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  ggtitle("full data")+
  ylab(NULL)+
  xlab(NULL)+
  theme(legend.position = "none")+
  xlim(-0.2, 0.43)
d_mehed <- d1 %>% filter(sex=="M")
m_mehed <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung, data = d_mehed, family = binomial)
dat_M <- predictions(m_mehed, newdata = d_mehed)
wts <- ifelse(dat_M$vac == 1, 1 / dat_M$estimate, 1 / (1 - dat_M$estimate))
d_mehed$wts <- wts

fig_smd_mehed <- tidy_smd(d_mehed, c(age, tongue, education, heart, lung, diabetes, county), .group = vaccinated, .wts = wts, na.rm = TRUE) %>% 
  ggplot(aes(smd, variable, color=method))+ 
  geom_point(size=3)+
  geom_vline(xintercept = c(-0.1, 0.1), linetype=2)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  ggtitle("males")+
  ylab(NULL)+
  xlab(NULL)+
  xlim(-0.2, 0.43)

fig_smd_mehed <- tidy_smd(d_mehed, c(age, tongue, education, heart, lung, diabetes, county), .group = vaccinated, .wts = wts, na.rm = TRUE) %>% 
  ggplot(aes(smd, variable, shape = method))+ 
  scale_shape_discrete(solid=F) +
  geom_point(size=3)+
  geom_vline(xintercept = c(-0.1, 0.1), linetype=2)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  ggtitle("males")+
  ylab(NULL)+
  xlab(NULL)+
  theme(legend.title=element_blank())+
  xlim(-0.2, 0.43)

d_naised <- d1 %>% filter(sex=="N")
m_naised <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung, data = d_naised, family = binomial)
dat_N <- predictions(m_naised, newdata = d_naised)
wts <- ifelse(dat_N$vac == 1, 1 / dat_N$estimate, 1 / (1 - dat_N$estimate))
d_naised$wts <- wts

fig_smd_naised <- tidy_smd(d_naised, c(age, tongue, education, heart, lung, diabetes, county), .group = vaccinated, .wts = wts, na.rm = TRUE) %>% 
  ggplot(aes(smd, variable, shape = method))+ 
  scale_shape_discrete(solid=F) +
  geom_point(size=3)+
  geom_vline(xintercept = c(-0.1, 0.1), linetype=2)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  ggtitle("females")+
  ylab(NULL)+
  theme(legend.position = "none")+
  xlim(-0.2, 0.43)
library(patchwork)
fig_smd_full / fig_smd_mehed + fig_smd_naised
ggsave("FigS2y.pdf", units = "cm", width = 12, height = 12)
```

# Fig S3.

Unadjusted cumulative incidence estimates for all-cause Death, stratified by sex and age group. Left panels show females and right panels show males. The risks for the vaccinated are shown in green and for the unvaccinated in red. Shaded areas denote 95% confidence intervals.

```{r fig.height=5, fig.width=5}

d1 <- d1 %>% mutate(age_binned = case_when(age <= 70 ~ "<70",
                                        age >70  ~ "71+"))
a <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="M", age_binned =="71+")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "2"))+
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "males 71+")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

b <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="N", age_binned =="71+")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "2")) +
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "females 71+")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

c <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="M", age_binned =="<70")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "2"))+
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "males <70")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

d <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="N", age_binned =="<70")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "2")) +
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "females <70")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

a+b+c+d + plot_layout(nrow = 2) 
ggsave("FigS3y.pdf", units = "cm", width = 12, height = 12)
```

# Fig S4.

Unadjusted cumulative incidence estimates, given as percentages, for MACE, stratified by sex and age group. Left panels show females and right panels show males. The risks for the vaccinated are shown in green and for the unvaccinated in red. Shaded areas denote 95% confidence intervals.

```{r fig.height=5, fig.width=5}

a <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="M", age_binned =="71+")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "1"))+
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "males 71+")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

b <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="N", age_binned =="71+")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "1")) +
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "females 71+")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

c <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="M", age_binned =="<70")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "1"))+
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "males <70")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

d <- cuminc(Surv(time, status, type="mstate") ~ vaccinated + age_binned, 
       data = d1 %>% filter(sex=="N", age_binned =="<70")) %>% 
  ggsurvfit::ggcuminc(outcome = c( "1")) +
  add_confidence_interval() +
  scale_ggsurvfit()+
  labs(subtitle = "females <70")+
  xlab(NULL)+
  ylab(NULL)+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))+
  scale_x_continuous(breaks = seq(0, 365, by = 100))

a+b+c+d + plot_layout(nrow = 2) 
ggsave("FigS4y.pdf", units = "cm", width = 12, height = 12)
```

# Fig S5

For data restricted to non-hard covid-19 infections, females' (panels A and C) and males' (B and D) mace rates (A and B) and all-cause death rates (C and D) over the study follow-up period. The rates were calculated using IP weights. The vaccinated are shown in green and the unvaccinated are in red. Shaded areas represent 95% confidence inrtervals.

firstly mace:

```{r}
d_lightC19 <- read_csv("d_lightC19_scrambled.csv")

lex_l <- Lexis(exit = list(tfe = time), 
            exit.status = factor(status, labels = c("Alive", "MACE", "Dead")),
             data = d_lightC19)
sL_l <- splitMulti(lex_l, tfe = c(0, 10, 20, 30, 40, 60, 80, 100, 150, 200, 300, 310, 320, 330, 340, 350, 360, 365))
sL$vaccinated <- as.factor(sL$vaccinated)

sL_male3 <- sL_l %>% filter(sex =="M")
sL_female3 <- sL_l %>% filter(sex =="N")
d_lightC19_F <- d_lightC19 %>% filter(sex=="N")
d_lightC19_M <- d_lightC19 %>% filter(sex=="M")

# 1.  females
m_logistic3_F_young <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, 
                           data = d_lightC19_F, 
                           family = binomial)
dat3_F_y <- predictions(m_logistic3_F_young, newdata = d_lightC19_F)
dat3_F_y$wts <- ifelse(dat3_F_y$vac == 1, 1 / dat3_F_y$estimate, 1 / (1 - dat3_F_y$estimate))

wts_F <- dat3_F_y %>% select(id, wts)
sL_female3 <- sL_female3 %>% left_join(wts_F)

k5_mace_F_light_y <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5),
                         family = poisreg, data = sL_female3, weights = wts)


pr_F_y <- predictions(
    k5_mace_F_light_y,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_F_y <- as.data.frame(pr_F_y) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)


A <- ggplot(pr_F_y, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
   ylab("MACE rates")+
  xlab("Days")+
  ggtitle("females")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

# 2. males
m_logistic3_M <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, 
                           data = d_lightC19_M, 
                           family = binomial)
dat3_M <- predictions(m_logistic3_M, newdata = d_lightC19_M)
dat3_M$wts <- ifelse(dat3_M$vac == 1, 1 / dat3_M$estimate, 1 / (1 - dat3_M$estimate))


wts_M <- dat3_M %>% select(id, wts)
sL_male3 <- sL_male3 %>% left_join(wts_M)

k5_mace_F_light_o <- gam(cbind(lex.Xst == "MACE", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5),
                         family = poisreg, data = sL_male3, weights = wts)

pr_M <- predictions(
    k5_mace_F_light_o,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pr_M <- as.data.frame(pr_M) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

B <- ggplot(pr_M, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
 ylab("MACE rates")+
  xlab("Days")+
  ggtitle("males")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))
```

and deaths...

```{r}
d_light <- read_csv("d_light_death_scrambled.csv")
lex_ld <- Lexis(exit = list(tfe = time), 
            exit.status = factor(status, labels = c("Alive", "Dead")),
             data = d_light)
sL_ld <- splitMulti(lex_ld, tfe = c(0, 10, 20, 30, 40, 60, 80, 100, 150, 200, 300, 310, 320, 330, 340, 350, 360, 365))
sL_ld$vac <- if_else(sL_ld$vaccinated == "yes", 1, 0)
sL_ld$vac <- if_else(sL_ld$vaccinated == "yes", 1, 0)
sL_ld$vaccinated <- as.factor(sL_ld$vaccinated)

sL_male <- sL_ld %>% filter(sex =="M")
sL_female <- sL_ld %>% filter(sex =="N")
d_male <- d_light %>% filter(sex =="M")
d_female <- d_light %>% filter(sex =="N")

# females and males

m_logistic_female <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data = d_female, family = binomial)
dat2_N <- predictions(m_logistic_female, newdata = d_female)
dat2_N$wts <- ifelse(dat2_N$vac == 1, 1 / dat2_N$estimate, 1 / (1 - dat2_N$estimate))

d_female$wts <- dat2_N$wts
weights <- d_female %>% select(id, wts)
sL_female <- sL_female %>% left_join(weights)

m_logistic_male <- gam(vac ~ education + s(age) + tongue * county + diabetes * heart + lung + cal_time1, data = d_male, family = binomial)
dat2_M <- predictions(m_logistic_male, newdata = d_male)
dat2_M$wts <- ifelse(dat2_M$vac == 1, 1 / dat2_M$estimate, 1 / (1 - dat2_M$estimate))

d_male$wts <- dat2_M$wts
weights <- d_male %>% select(id, wts)
sL_male <- sL_male %>% left_join(weights)

m_female_death <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = sL_female, weights = wts)
m_male_death <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ vaccinated + s(tfe, by= vaccinated, k=5), 
             family = poisreg, data = sL_male, weights = wts)

pred_ipw_female_d <- predictions(
    m_female_death,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pred_ipw_female_d <- as.data.frame(pred_ipw_female_d) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

pred_ipw_male_d <- predictions(
    m_male_death,
    type = "response",
    by = c("vaccinated", "tfe"),
    wts = "wts")
pred_ipw_male_d <- as.data.frame(pred_ipw_male_d) %>% 
  mutate(estimate = estimate*36500, 
         conf.low = conf.low*36500,
         conf.high = conf.high*36500)

C <- ggplot(pred_ipw_female_d, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab("death rates")+
  xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))
D <- ggplot(pred_ipw_male_d, aes(tfe, estimate)) + 
  geom_line(aes(color=vaccinated, linetype=vaccinated)) + 
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill= vaccinated), alpha=0.3)+
  ylab("death rates")+
  xlab("Days")+
  theme_bw()+
  theme(legend.position="none")+
   scale_fill_manual(values=c("firebrick", "darkgreen"))

A+B+C+D + plot_annotation(tag_levels = 'A')
ggsave("figS5y.pdf", width = 8, height = 4)
```
