analysis caap attrib pct.Rmd

---
title: "pneumo israel"
author: "Dan Weinberger"
date: "2/11/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(readxl)
library(reshape2)
library(ggplot2)
library(lubridate)
library(dplyr)
library(INLA)
library(tidyr)
library(htmlTable)
#install.packages("waffle", repos = "https://cinc.rud.is")  
library(waffle)
library(abind)
library(pbapply)
library(rjags)
library(HDInterval)
library(patchwork)

source('./R/format_carr.R')
source('./R/format_caap.R')
source('./R/format_virus.R')
source('./R/format_ipd_func.R')
source('./R/inla_model3.R')
source('./R/Plot_Obs_exp_counterfact.R')
source('./R/attrib_pct_func.R')
source('./R/attrib_pct_fmt_func.R')
source('./R/attrib_pct_fmt_func_age.R')
source('./R/attrib_pct_waffle_age.R')
source('./R/format_post_func1.R')
source('./R/fitted_ci_mcmc.R')
source('./R/attrib_pct_age_mcmc.R')
source('./R/jags_model.R')
source('./R/plot_Obs_exp_counterfact_mcmc.R')
source('./R/gen_pred_interval_inla_ridge_ar1.R')

```

-CAAP and non-CAAP LRI 2016-2021 (southern Israel, children <5y)

-IPD 2016-2021 (Nationwide children <5y)

-All viruses (all nasal washes obtained, southern Israel, children <5y)

-CAAP virus coinfections 2016-2021 (Children <5y with CAAP and virus positive, southern Israel)

-Non-CAAP LRI virus coinfections 2016-2021 (Children <5y with non-CAAP LRI and virus positive, southern Israel)

-Carriage plus density 2016-2021 (Children <3y, no respiratory diagnosis, southern Israel)
 

```{r}
carr.formatted <- format_carr()

carr2.c <- carr.formatted$carr2.c

carr2.c <- carr2.c %>%
          group_by(ethnicity, agec) %>%
          mutate(N_swab=POS+NEG,
                  lag1_POS=lag(POS),
                 lag2_pOS=lag(POS,2),
                 lag1_N_swab=lag(N_swab),
                 lag2_N_swab=lag(N_swab,2),
                 prev_movave= (POS+lag1_POS+lag2_pOS)/(N_swab+lag1_N_swab+lag2_N_swab))

p1 <- ggplot(carr2.c, aes(x=date, y=prev_movave, group=ethnicity, col=ethnicity))+
  geom_line() +
  theme_classic()+
  facet_wrap(~agec) +
  ylim(0, 1)

p1
```

Prevalence over time, subsetted by density. The change in prevalence seen in 2020 is due to drop in low-density carriage 
```{r}
carr.formatted$plot.density
```

Viruses

```{r, fig.width=6, fig.height=4}
virus_formatted <- format_virus()

vir2.c <- virus_formatted$vir2.c

virus_formatted$virus.plot
```

## CAAP

```{r, fig.width=6, fig.height=2}
c3 <- format_caap()

p3 <- ggplot(c3, aes(x=date, y=CAAP, group=ethnicity, col=ethnicity)) +
  geom_line()+
  theme_classic() +
  facet_wrap(~agec, scales='free') +
  ggtitle('CAAP by age group and ethnicity')
p3


p4 <- ggplot(c3, aes(x=date, y=nonCAAP, group=ethnicity, col=ethnicity)) +
  geom_line()+
  theme_classic() +
  facet_wrap(~agec, scales='free') +
  ggtitle('non-CAAP by age group and ethnicity')
p4

```

IPD
This script just brings in data for kids <1 year
```{r}
ipd1 <- format_ipd_func()

ipd1$p.ipd_other

ipd1$p.ipd_pneu
```


Combine the virus, carriage and X-rays datasets
```{r}
all1 <- merge(c3, carr2.c, by=c('ethnicity', 'agec', 'date'), all=T)

all1 <- all1[!is.na(all1$CAAP) & all1$date>='2016-01-01',]

all1$prev[is.na(all1$prev)] <-0

all1 <- merge(all1, ipd1$ds[,c('agec','ethnicity','date','bacteremic.pneumonia','other.IPD')], by=c('ethnicity', 'agec', 'date'), all.x=T)


all1 <- merge(all1, vir2.c, by=c('ethnicity', 'date'), all=T)


all1$ethnicity <- as.factor(all1$ethnicity)

all1$agec <- as.factor(all1$agec)

all1$caap.pre <- all1$CAAP

all1$caap.pre[all1$date>=as.Date('2020-03-01')] <- NA

all1$non.caap.pre <- all1$nonCAAP
all1$non.caap.pre[all1$date>=as.Date('2020-03-01')] <- NA

all1$t <- interval(min(all1$date), all1$date) %/% months(1) 
all1$sin12 <- sin(2*pi*all1$t/12)
all1$cos12 <- cos(2*pi*all1$t/12)

all1$sin6 <- sin(2*pi*all1$t/6)
all1$cos6 <- cos(2*pi*all1$t/6)

all1$sin3 <- sin(2*pi*all1$t/3)
all1$cos3 <- cos(2*pi*all1$t/3)

all1$year <- year(all1$date)
all1$month <- month(all1$date)

all1$epiyr <- all1$year
all1$epiyr[all1$month<=6] <- all1$year[all1$month<=6]  - 1
all1$epiyr[all1$epiyr >2019] <- 2019
all1$epiyr <- as.factor(all1$epiyr)

all1$month <- as.factor(all1$month)

all1 <- all1[!is.na(all1$agec),]

```


```{r}

```



## INLA MODELS

## CAAP, fit to pre-COVID data only

Fit model to pre data only or to any 

```{r, eval=F}
res1.caap.pre <- jags_func(all1,  outcome.var='caap.pre', plot.var='CAAP')
res1.caap.all <- jags_func(all1,  outcome.var='CAAP', plot.var='CAAP')
res1.noncaap.pre <- jags_func(all1,  outcome.var='non.caap.pre', plot.var='nonCAAP')
res1.noncaap.all <- jags_func(all1,  outcome.var='nonCAAP', plot.var='nonCAAP')
res1.ipd.pneu.all <- jags_func(all1,  outcome.var='bacteremic.pneumonia', plot.var='bacteremic.pneumonia')
res1.ipd.other.all <- jags_func(all1,  outcome.var='other.IPD', plot.var='other.IPD')

mod.results <- list( 'res1.caap.pre'=res1.caap.pre,'res1.caap.all'=res1.caap.all,'res1.noncaap.pre'=res1.noncaap.pre,'res1.noncaap.all'=res1.noncaap.all,'res1.ipd.pneu.all'=res1.ipd.pneu.all,'res1.ipd.other.all'=res1.ipd.other.all)
saveRDS(mod.results,'./Results/model_fits.rds')
```

```{r}
mod.results <- readRDS('./Results/model_fits.rds')
res1.caap.pre <- mod.results$res1.caap.pre
res1.caap.all <- mod.results$res1.caap.all
res1.noncaap.pre <- mod.results$res1.noncaap.pre
res1.noncaap.all <- mod.results$res1.noncaap.all
res1.ipd.pneu.all <- mod.results$res1.ipd.pneu.all
res1.ipd.other.all <- mod.results$res1.ipd.other.all
```

Fit to pre-pandemic CAAP only

```{r, fig.width=9, fig.height=3}

plot_Obs_exp_counterfact_mcmc(ds.in=res1.caap.pre, plot.var='CAAP', post.ci=T)
```





## CAAP, fit to ALL data

```{r, fig.width=9, fig.height=3}
p1 <- plot_Obs_exp_counterfact_mcmc(ds.in=res1.caap.all, plot.var='CAAP') +
  ggtitle('CAAP')
p1
```



## non-CAAP, fit to all data


```{r, fig.width=9, fig.height=3}
p2 <- plot_Obs_exp_counterfact_mcmc(ds.in=res1.noncaap.all, plot.var='nonCAAP') +
    ggtitle('non-CAAP')
p2
```


## IPD pneumonia 


```{r, fig.width=9, fig.height=3}
p3 <- plot_Obs_exp_counterfact_mcmc(ds.in=res1.ipd.pneu.all, plot.var='bacteremic.pneumonia') +
      ggtitle('IPD (pneumonia)')
p3
```

## IPD non-pneumonia 

```{r, fig.width=9, fig.height=3}
p4 <- plot_Obs_exp_counterfact_mcmc(ds.in=res1.ipd.other.all, plot.var='other.IPD') +
  ggtitle('IPD (other)')
p4
```


```{r, fig.height=6, fig.width=6}
p1/ p2/ p3 /p4

```


## Compare attributable percents for different outcomes, etc
```{r}

ap1 <- attrib_pct_fmt_func_mcmc('res1.caap.all','CAAP_all')

ap3 <- attrib_pct_fmt_func_mcmc('res1.noncaap.all','NonCAAP_all')


ap5 <- attrib_pct_fmt_func_mcmc('res1.ipd.pneu.all','IPD_pneumonia_all')

ap6 <- attrib_pct_fmt_func_mcmc('res1.ipd.other.all','IPD_other_all')


#htmlTable(cbind(ap2, ap4[,2], ap1[,2], ap3[,2]))


tab1 <- htmlTable(cbind.data.frame(ap1, ap3[,-1, drop=F], ap5[,2, drop=F], ap6[,2, drop=F]))
tab1

out.df <- as.data.frame(cbind(ap1, ap3[,2], ap5[,2], ap6[,2]))

write.csv(out.df, './Results/overall_attrib.csv')

```


```{r}
ap2 <- attrib_pct_fmt_func_mcmc('res1.caap.pre','CAAP_pre')

ap2

write.csv(ap2, './Results/attrib_pct_caap_pre.csv')

```

Age
```{r}

age_ap1 <- attrib_pct_fmt_func_age_mcmc('res1.caap.all','CAAP_all')

age_ap3 <- attrib_pct_fmt_func_age_mcmc('res1.noncaap.all','NonCAAP_all')


age_ap5 <- attrib_pct_fmt_func_age_mcmc('res1.ipd.pneu.all','IPD_pneumonia_all')

age_ap6 <- attrib_pct_fmt_func_age_mcmc('res1.ipd.other.all','IPD_other_all')


#htmlTable(cbind(ap2, ap4[,2], ap1[,2], ap3[,2]))

combined.df <- cbind(age_ap1, age_ap3[,3, drop=F], age_ap5[,3, drop=F], age_ap6[,3, drop=F])

combined.df <- combined.df[order(combined.df$Virus,combined.df$agec),]
row.names(combined.df) <- NULL

combined.df$agec <- factor(combined.df$agec)
#htmlTable(combined.df[,-c(2)], rgroup=c('Any Virus','RSV','hMPV','Flu','Adeno', 'paraflu'), n.rgroup=rep(3,6) ,rnames=F)

#combined.df$virus <- factor(combined.df$virus, levels=c('any','rsv','hmpv','flu','adeno','paraflu'))


combined.df2 <- combined.df[order(combined.df$agec, combined.df$Virus),]

htmlTable(combined.df2[,-1], rgroup=c('<12m','12-23m','2-4y'), n.rgroup=rep(6,3), rnames=F )


combined.df3 <- as.data.frame(combined.df2)
write.csv(combined.df3, './Results/age_stratified_attrib.csv')

```

```{r, fig.width=6, fig.height=6}
attrib_pct_waffle_age_mcmc('res1.caap.all')

```



Rhino--only available in recent times
```{r}

all1$out.scale <- scale(all1$CAAP/all1$pop*1000)
modr1 <- lm(CAAP ~ agec*ethnicity *(RSV + influenza.a + influenza.b + hMPV + parainfluenza + rhinovirus ) , data=all1)

modr2 <- lm(CAAP ~ agec + ethnicity *(RSV + influenza.a + influenza.b + hMPV + parainfluenza + rhinovirus ) , data=all1)

modr3 <- lm(CAAP ~ ethnicity +agec*(RSV + influenza.a + influenza.b + hMPV + parainfluenza + rhinovirus ) , data=all1)

modr4 <- lm(CAAP ~ ethnicity*(RSV + influenza.a + influenza.b + hMPV + parainfluenza + rhinovirus ) +
              agec*(RSV + influenza.a + influenza.b + hMPV + parainfluenza + rhinovirus ) , data=all1)


AIC(modr1, modr2, modr3,modr4)
summary(modr)
```