interaction_analysis_clean.rmd

---
title: "Interaction_analysis_clean"
output: pdf_document
date: "2024-01-30"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(ggplot2)
library(rigr)
library(readr)
library(lubridate)
library(ggthemes)
library(Polychrome)
library(ggbeeswarm)
library(stringr)
library(tidyverse)
library(smplot2)

```



# Repository
```{r}

# set directory
repo_root <- "/Users/eabousam" 
repo      <- sprintf('%s/SFS-coinfection-interactions/',     repo_root)
setwd(repo)

```


# Read Data

```{R}

#full encounter sample data
full_encounter_sample = readRDS('model_data.RDS')

full_encounter_sample = full_encounter_sample %>%
  filter(organism %in% c("Rhinovirus", "SARS_CoV_2")) 



#Clean Vaccination data
vacc_data = read.csv("SCAN_vacc_data.csv")

#Assign vaccination status
vacc_data$vaccine_status <- ifelse(vacc_data$covidshot1 == '["yes"]' & vacc_data$covidshot2 == '["yes"]', "fully_vaccinated", vacc_data$vaccine_status)
vacc_data$vaccine_status <- ifelse(vacc_data$covidshot1 == '["no"]' & vacc_data$covidshot2 == '["no"]', "not_vaccinated", vacc_data$vaccine_status)
vacc_data$vaccine_status <- ifelse(vacc_data$covidshot1 == '["yes"]' & vacc_data$covidshot2 == '["yes"]' & vacc_data$covidshot3 == '["yes"]', "boosted", vacc_data$vaccine_status)
#Assign vaccination status for unknown classification
vacc_data <- vacc_data %>%
  mutate(vaccine_status = case_when(
    vaccine_status == "unknown" & covidshot1date != "" & covidshot2date == "" ~ "partially_vaccinated",
    vaccine_status == "unknown" & covidshot1date != "" & covidshot2date != "" ~ "fully_vaccinated",
    TRUE ~ vaccine_status
  ))
#Assign vaccination status for empty classification
vacc_data <- vacc_data %>%
  mutate(vaccine_status = case_when(
    vaccine_status == "" & covidshot1date != "" & covidshot2date == "" ~ "partially_vaccinated",
    vaccine_status == "" & covidshot1date != "" & covidshot2date != "" ~ "fully_vaccinated",
    TRUE ~ vaccine_status
  ))
#Remove incomplete/invalid/null data
vacc_data = vacc_data %>%
  filter(vaccine_status != "invalid" &
         vaccine_status != "unknown" &
         vaccine_status != "na" &
         vaccine_status != "")
  


# Join datasets by sfs sample identification
full_encounter_sample_vacc = full_encounter_sample %>%
  left_join(vacc_data, by = "sample")




# Adding vaccine status at encounter variable


full_encounter_sample_vacc <- full_encounter_sample_vacc %>%
  mutate(vacc_status_at_enc = case_when(
    covidshot1date > date ~ "not_vaccinated",
    covidshot1date < date & covidshot2date > date ~ "partially_vaccinated",
    covidshot1date < date & covidshot2date < date ~ "fully_vaccinated",
    covidshot1date < date & covidshot2date < date & covidshot3date < date ~ "boosted",
    vaccine_status == "not_vaccinated" ~ "not_vaccinated",
    TRUE ~ NA_character_  # Handles any other cases
  ))


```




#Hypothesis: RV decrease/increase covid infectivity odds

```{r}
range(full_encounter_sample_vacc$date, na.rm = T)

encounter_with_multiple_data_rv = full_encounter_sample_vacc %>%
  filter(organism == "Rhinovirus") %>%
  dplyr::select(individual, present, date ) %>%
  rename(present_rv = present,
         encountered_date_rv = date) %>%
  na.omit()

encounter_with_multiple_data_cv = full_encounter_sample_vacc %>%
  filter(organism %in% c("SARS_CoV_2")) %>%
  dplyr::select(individual, present, date ) %>%
  rename(present_cv = present,
         encountered_date_cv = date) %>%
  na.omit()


encounter_rv_cv_pairs = encounter_with_multiple_data_rv %>%
  inner_join(encounter_with_multiple_data_cv)

encounter_rv_cv_pairs = encounter_rv_cv_pairs %>%
  #filter(encountered_week_cv > encountered_week_rv) %>%
  mutate(present_cv_n = as.integer(as.logical(present_cv))) %>%
  mutate(present_rv_n = as.integer(as.logical(present_rv))) %>%
  #mutate(rv_cv_diff = encountered_week_cv- encountered_week_rv)
  mutate(rv_cv_diff = encountered_date_cv - encountered_date_rv)



encounter_rv_cv_pairs$rv_cv_diff <- as.numeric(encounter_rv_cv_pairs$rv_cv_diff, units = "days")


#Cleaning vaccination dataset
full_encounter_sample_vacc_analysis = full_encounter_sample_vacc %>%
  select(individual, sex, age,covidshot1date, covidshot2date, covidshot3date, date ) %>%
  na.omit("individual") %>%
  rename(encountered_date_cv = date )


#Merge covariates and vaccination data
encounter_rv_cv_pairs_conf = encounter_rv_cv_pairs %>%
  left_join(full_encounter_sample_vacc_analysis, by = c("individual","encountered_date_cv")) %>%
  distinct()

#Replace missing dates with non-vaccination
encounter_rv_cv_pairs_conf$covidshot1date[encounter_rv_cv_pairs_conf$covidshot1date==""]<-"Notvaccinated"
encounter_rv_cv_pairs_conf$covidshot2date[encounter_rv_cv_pairs_conf$covidshot2date==""]<-"Notvaccinated"
encounter_rv_cv_pairs_conf$covidshot3date[encounter_rv_cv_pairs_conf$covidshot3date==""]<-"Notvaccinated"

encounter_rv_cv_pairs_conf = encounter_rv_cv_pairs_conf %>%
  mutate(across(starts_with("covid"), ~ gsub("\\[|\\]|\\(|\\)|\"", "", .))) %>%
  na.omit()

encounter_rv_cv_pairs_conf = encounter_rv_cv_pairs_conf %>%
  mutate(across(c(covidshot1date, covidshot2date, covidshot3date), na_if, "Notvaccinated"))


#Mutate vaccination status variable
encounter_rv_cv_pairs_conf <- encounter_rv_cv_pairs_conf %>%
 mutate(
    covidshot1date = as.Date(covidshot1date),
    covidshot2date = as.Date(covidshot2date),
    covidshot3date = as.Date(covidshot3date)
  ) %>%
  mutate(
    vacc_status_at_enc_cv = case_when(
      covidshot1date > encountered_date_cv ~ "not_vaccinated",
      covidshot1date < encountered_date_cv & covidshot2date > encountered_date_cv ~ "partially_vaccinated",
      covidshot1date < encountered_date_cv & covidshot2date < encountered_date_cv ~ "fully_vaccinated",
      covidshot1date < encountered_date_cv & covidshot2date < encountered_date_cv & covidshot3date < encountered_date_cv ~ "boosted",
      TRUE ~ NA_character_  # Handles any other cases
    )
  ) 


encounter_rv_cv_pairs_conf = encounter_rv_cv_pairs_conf %>%
    mutate(vacc_status_at_enc_cv = replace_na(vacc_status_at_enc_cv, "not_vaccinated")) %>%
  select(-covidshot1date, -covidshot2date, -covidshot3date)

table(encounter_rv_cv_pairs_conf$vacc_status_at_enc_cv)



#Defining periods

encounter_rv_cv_pair_2weeks = encounter_rv_cv_pairs_conf %>%
  filter(rv_cv_diff < 14 &
           rv_cv_diff > 0) %>%
  na.omit()

encounter_rv_cv_pair_4weeks = encounter_rv_cv_pairs_conf %>%
  filter(rv_cv_diff < 30 &
           rv_cv_diff > 0) %>%
  na.omit()

encounter_rv_cv_pair_8weeks = encounter_rv_cv_pairs_conf %>%
  filter(rv_cv_diff < 60 &
           rv_cv_diff > 0) %>%
  na.omit()

encounter_rv_cv_pair_8weeks %>%
  na.omit()
#Running models

#within 14 days
regress("odds", formula = present_cv_n~present_rv_n*rv_cv_diff  + sex +age+vacc_status_at_enc_cv, data = encounter_rv_cv_pair_2weeks)
#within 30 days
regress("odds", formula = present_cv_n~present_rv_n*rv_cv_diff + sex +age +vacc_status_at_enc_cv , data = encounter_rv_cv_pair_4weeks)

names(encounter_rv_cv_pair_4weeks)
#within 60 days
regress("odds", formula = present_cv_n~present_rv_n*rv_cv_diff + sex +age+ vacc_status_at_enc_cv, data = encounter_rv_cv_pair_8weeks)

#x axis rv_Cs_diff faceting by presence rv, y = presence for cv

  
table(encounter_rv_cv_pair_2weeks$rv_cv_diff)


encounter_rv_cv_pair_8weeks %>%
  filter(rv_cv_diff=="0")


# Add infection window category variable

```


# Sample Size Calculations RV>CV
```{r}
size_sum(encounter_rv_cv_pair_2weeks) 
size_sum(encounter_rv_cv_pair_4weeks)
size_sum(encounter_rv_cv_pair_8weeks)

table(encounter_rv_cv_pair_8weeks$present_rv_n)

```


```{r}
#Visualizing 

encounter_rv_cv_pair_2weeks %>%
  na.omit() %>%
  ggplot(aes(x = rv_cv_diff, y = present_cv_n, group = present_rv_n, color = as.factor(present_rv_n))) +
  #geom_point()+
  geom_smooth()+
  #eom_smooth(method = "glm", method.args = list(family = "binomial"), se = T) +  # Use logistic regression line
  theme_minimal() +  # Use a minimal theme
  labs(title = "Association Between RV and COVID Infections within 8 weeks",
       x = "Time Since RV Infection (Days)",
       y = "COVID Infection",
       color = "RV Presence") +
  scale_color_manual(values = c("0" = "blue", "1" = "red")) +  # Customize color scale
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12)) 




encounter_rv_cv_pair_8weeks_band <- transform(encounter_rv_cv_pair_8weeks,
                                            time_period = cut(rv_cv_diff,
                                                             breaks = c(0, 5, 10, 15, 20, 30, 40, 50, 60),
                                                             labels = c("0-4 days", "5-9 days", "10-14 days", "15-20 days", "20-30 days", "30-40 days", "40-50 days", "50-60 days")))


encounter_rv_cv_pair_8weeks_band <- transform(encounter_rv_cv_pair_8weeks,
                                            time_period = cut(rv_cv_diff,
                                                             breaks = c(-1, 0, 5, 10, 15, 20, 30, 40, 50, 60),
                                                             labels = c("0", "1-4 days", "5-9 days", "10-14 days", "15-20 days", "20-30 days", "30-40 days", "40-50 days", "50-60 days")))
table(encounter_rv_cv_pair_8weeks_band$time_period)


encounter_rv_cv_pair_8weeks_band %>%
  na.omit() %>%
  ggplot(aes(x = as.factor(time_period), y = present_cv_n, color = as.factor(present_rv_n),
             group = present_rv_n)) +
  #geom_violin(alpha = 0.6) +
  #stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +  
  sm_pointplot(legends = T) +  
  scale_color_manual(values = sm_palette(2))+
  labs(title = "Proportion comparison between RV and COVID Infections within 8 weeks",
       x = "Time Since RV Infection (Days)",
    #   y = "Probability of COVID Infection",
       color = "RV Presence") +
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12))

encounter_rv_cv_pair_8weeks_band %>%
  na.omit() %>%
  ggplot(aes(x = rv_cv_diff)) +
  #geom_violin(alpha = 0.6) +
  #stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +  
  sm_hist(legends = T) +  
  scale_color_manual(values = sm_palette(2))+
  labs(title = "distribution of available events",
       x = "Time Since RV Infection (Days)",
    #   y = "Probability of COVID Infection",
       color = "RV Presence") +
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12))



```








# Hypothesis: Covid decrease/increase RV infectivity odds


```{r}


encounter_with_multiple_data_rv = full_encounter_sample_vacc %>%
  filter(organism == "Rhinovirus") %>%
  select(individual, present, date ) %>%
  rename(present_rv = present,
         encountered_date_rv = date) %>%
  na.omit()

encounter_with_multiple_data_cv = full_encounter_sample_vacc %>%
  filter(organism == "SARS_CoV_2") %>%
  select(individual, present, date ) %>%
  rename(present_cv = present,
         encountered_date_cv = date) %>%
  na.omit()


encounter_cv_rv_pairs = encounter_with_multiple_data_rv %>%
  inner_join(encounter_with_multiple_data_cv, by="individual")


encounter_cv_rv_pairs = encounter_cv_rv_pairs %>%
  #filter(encountered_week_cv > encountered_week_rv) %>%
  mutate(present_cv_n = as.integer(as.logical(present_cv))) %>%
  mutate(present_rv_n = as.integer(as.logical(present_rv))) %>%
  #mutate(rv_cv_diff = encountered_week_cv- encountered_week_rv)
  mutate(cv_rv_diff = encountered_date_rv - encountered_date_cv)


encounter_cv_rv_pairs$cv_rv_diff <- as.numeric(encounter_cv_rv_pairs$cv_rv_diff, units = "days")


#Cleaning vaccination dataset
full_encounter_sample_vacc_analysis_rv = full_encounter_sample_vacc %>%
  select(individual, sex, age,covidshot1date, covidshot2date, covidshot3date, date ) %>%
  na.omit("individual") %>%
  rename(encountered_date_rv = date )


#Merge covariates and vaccination data
encounter_cv_rv_pairs_conf = encounter_cv_rv_pairs %>%
  left_join(full_encounter_sample_vacc_analysis_rv, by = c("individual","encountered_date_rv")) %>%
  distinct()

#Replace missing dates with non-vaccination
encounter_cv_rv_pairs_conf$covidshot1date[encounter_cv_rv_pairs_conf$covidshot1date==""]<-"Notvaccinated"
encounter_cv_rv_pairs_conf$covidshot2date[encounter_cv_rv_pairs_conf$covidshot2date==""]<-"Notvaccinated"
encounter_cv_rv_pairs_conf$covidshot3date[encounter_cv_rv_pairs_conf$covidshot3date==""]<-"Notvaccinated"

encounter_cv_rv_pairs_conf = encounter_cv_rv_pairs_conf %>%
  mutate(across(starts_with("covid"), ~ gsub("\\[|\\]|\\(|\\)|\"", "", .))) %>%
  na.omit()

encounter_cv_rv_pairs_conf = encounter_cv_rv_pairs_conf %>%
  mutate(across(c(covidshot1date, covidshot2date, covidshot3date), na_if, "Notvaccinated"))


#Mutate vaccination status variable
encounter_cv_rv_pairs_conf <- encounter_cv_rv_pairs_conf %>%
 mutate(
    covidshot1date = as.Date(covidshot1date),
    covidshot2date = as.Date(covidshot2date),
    covidshot3date = as.Date(covidshot3date)
  ) %>%
  mutate(
    vacc_status_at_enc_rv = case_when(
      covidshot1date > encountered_date_rv ~ "not_vaccinated",
      covidshot1date < encountered_date_rv & covidshot2date > encountered_date_rv ~ "partially_vaccinated",
      covidshot1date < encountered_date_rv & covidshot2date < encountered_date_rv ~ "fully_vaccinated",
      covidshot1date < encountered_date_rv & covidshot2date < encountered_date_rv & covidshot3date < encountered_date_rv ~ "boosted",
      TRUE ~ NA_character_  # Handles any other cases
    )
  ) 


encounter_cv_rv_pairs_conf = encounter_cv_rv_pairs_conf %>%
    mutate(vacc_status_at_enc_rv = replace_na(vacc_status_at_enc_rv, "not_vaccinated")) %>%
  select(-covidshot1date, -covidshot2date, -covidshot3date)


#Defining periods

encounter_cv_rv_pair_2weeks = encounter_cv_rv_pairs_conf %>%
  filter(cv_rv_diff < 14 &
           cv_rv_diff > 0) %>%
  na.omit()

table(encounter_rv_cv_pair_2weeks$rv_cv_diff)

encounter_cv_rv_pair_4weeks = encounter_cv_rv_pairs_conf %>%
  filter(cv_rv_diff < 30 &
           cv_rv_diff > 0) %>%
  na.omit()


encounter_cv_rv_pair_8weeks = encounter_cv_rv_pairs_conf %>%
  filter(cv_rv_diff < 60 &
          cv_rv_diff > -1) %>%
  na.omit()

#Running models

#within 14 days
regress("odds", formula = present_rv_n~present_cv_n*cv_rv_diff + sex  +age +vacc_status_at_enc_rv, data = encounter_cv_rv_pair_2weeks)

#within 30 days
regress("odds", formula = present_rv_n~present_cv_n*cv_rv_diff + sex  +age +vacc_status_at_enc_rv, data = encounter_cv_rv_pair_4weeks)

#within 60 days
regress("odds", formula = present_rv_n~present_cv_n*cv_rv_diff + sex +age + vacc_status_at_enc_rv, data = encounter_cv_rv_pair_8weeks)

#x axis rv_Cs_diff faceting by presence rv, y = presence for cv

encounter_cv_rv_pair_4weeks %>%
  na.omit() %>%
  ggplot(aes( x = cv_rv_diff, y =  present_rv_n, group = present_cv_n , color = as.factor(present_cv_n))) +
  #geom_point()+
  geom_smooth() +
  theme_classic()# +
  #acet_wrap(~Race)

encounter_cv_rv_pair_8weeks %>%
  na.omit() %>%
  ggplot(aes(x = cv_rv_diff, y = present_rv_n, group = present_cv_n, color = as.factor(present_cv_n))) +
  #geom_point()+
  geom_smooth()+
  #eom_smooth(method = "glm", method.args = list(family = "binomial"), se = T) +  # Use logistic regression line
  theme_minimal() +  # Use a minimal theme
  labs(title = "Association Between RV and COVID Infections within 8 weeks",
       x = "Time Since CV Infection (Days)",
       y = "RV Infection Probability",
       color = "RV Presence") +
  scale_color_manual(values = c("0" = "blue", "1" = "red")) +  # Customize color scale
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12)) 

#Idea
#Hospitalization data? admittance?? ask about


encounter_cv_rv_pair_8weeks %>%
  na.omit() %>%
  ggplot(aes(x = cv_rv_diff, y = present_rv_n, color = as.factor(present_cv_n))) +
  #geom_point(alpha = 0.6) +
  stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +  
  theme_minimal() +  
  labs(title = "Association Between RV and COVID Infections within 8 weeks",
       x = "Time Since COVID Infection (Days)",
       y = "Probability of RV Infection",
       color = "RV Presence") +
  scale_color_manual(values = c("0" = "blue", "1" = "red")) +
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12))





```



# Sample Size Calculations RV>CV
```{r}
size_sum(encounter_cv_rv_pair_2weeks) 
size_sum(encounter_cv_rv_pair_4weeks)
size_sum(encounter_cv_rv_pair_8weeks)

table(encounter_rv_cv_pair_8weeks$present_rv_n)
range(encounter_rv_cv_pair_8weeks$encountered_date_cv)
```








```{r}


encounter_cv_rv_pair_8weeks_band <- transform(encounter_cv_rv_pair_8weeks,
                                            time_period = cut(cv_rv_diff,
                                                             breaks = c(-1, 0, 5, 10, 15, 20, 30, 40, 50, 60),
                                                             labels = c("0", "1-4 days", "5-9 days", "10-14 days", "15-20 days", "20-30 days", "30-40 days", "40-50 days", "50-60 days")))




encounter_cv_rv_pair_8weeks_band %>%
  na.omit() %>%
  ggplot(aes(x = as.factor(time_period), y = present_rv_n, color = as.factor(present_cv_n),
             group = present_cv_n)) +
  #geom_violin(alpha = 0.6) +
  #stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +  
  sm_pointplot(legends = T) +  
  scale_color_manual(values = sm_palette(2))+
  labs(title = "Proportion comparison between RV and COVID Infections within 8 weeks",
       x = "Time Since CV Infection (Days)",
    #   y = "Probability of COVID Infection",
       color = "CV Presence") +
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12))

```

# Hypothesis: Vaccination effectWithin 3-weeks

```{r}
#within 3 weeks, count of present (RV)

vacc_encouter_data = full_data %>%
  select(individual, present, date, covidshot1date,covidshot2date, covidshot3date, organism, sex, age, 
         Race, Income) %>%
  na.omit()

vacc_encouter_data = vacc_encouter_data %>%
  mutate(Race = fct_recode(
    Race, "White Non-hispanic" = "{white}", 
    "Black Non-hispanic" = "{blackOrAfricanAmerican}",
    "Black Non-hispanic" = "{blackOrAfricanAmerican,white}",
    "Black Non-hispanic"= "{blackOrAfricanAmerican,other}",
    "Black Non-hispanic"= "{blackOrAfricanAmerican,nativeHawaiian}",
    "Black Non-hispanic"= "{blackOrAfricanAmerican,nativeHawaiian,white}",
    "Asian" = "{asian}",
    "Asian" = "{asian,white}",
   "Asian" = "{asian,blackOrAfricanAmerican}",
   "Asian"= "{asian,nativeHawaiian}",
   "Asian"= "{asian,blackOrAfricanAmerican,white}",
   "Asian" = "{asian,other}",
   "Asian"= "{asian,other,white}",
   "Asian"= "{asian,nativeHawaiian,white}",
    "Other" = "{other}", 
    "dont_say" = "{}",
    "Other"= "{nativeHawaiian}",
    "Other"= "{americanIndianOrAlaskaNative}",
    "Other"= "{americanIndianOrAlaskaNative}",
   "Other" = "{americanIndianOrAlaskaNative,asian,white}",
   "Other"= "{asian,blackOrAfricanAmerican,other,white}",
   "Other"= "{other,white}",
   "Other"= "{americanIndianOrAlaskaNative,white}",
   "Other"="{americanIndianOrAlaskaNative,blackOrAfricanAmerican,white}",
   "Other" = "{nativeHawaiian,white}"
    )) 


vacc_encouter_data = vacc_encouter_data %>%
  filter(Race %in% c("White Non-hispanic", "Black Non-hispanic", "Asian", "dont_say", "Other")) %>%
  filter(Income %in% c("100k_125k", "125k_150k", "25k_50k", "50k_75k", "75k_100k", "dont_know", "dont_say",
                       "less_25k", "more_150k")) %>%
  na.omit()





vacc_encouter_data <- vacc_encouter_data %>%
  #filter(vacc_status_at_enc == "fully_vaccinated") %>%
  filter(organism %in% c("Rhinovirus")) %>%
  mutate(present = as.integer(as.logical(present))) %>%
  mutate(days_since_shot1 = date - as.Date(covidshot1date),
         days_since_shot2 = date - as.Date(covidshot2date)) %>%
  na.omit() %>%
  filter(days_since_shot1 < 180 &
           days_since_shot1 > 0) %>%
  mutate(within_30_days_sh1 = ifelse(days_since_shot1 < 30, TRUE, FALSE),
         within_30_days_sh2 = ifelse(days_since_shot2 < 30, TRUE, FALSE),
         within_30_days = within_30_days_sh1 | within_30_days_sh2)
  

#Need to find most recent vaccination followed by a rhinovirus infection  
#Treat vaccination as infection, time since last vaccination


regress("odds", formula = present~within_30_days+age+Income+Race+sex, data =  vacc_encouter_data)

regress("odds", formula = present~within_30_days_sh2+age+sex+Race+Income , data =  vacc_encouter_data)





#vizualizing distribution

vacc_encouter_data %>%
  ggplot(aes(x = present, group = within_30_days, fill = within_30_days)) +
  geom_histogram(aes(y = ..density..), bins = 10, position = "dodge", alpha = 0.7) +
  facet_wrap(~within_30_days) +
  scale_fill_viridis_d(option = "magma", begin = 0.2, end = 0.8) +  # Use viridis color palette
  theme_classic() +
  labs(title = "Distribution of 'RV present' Within 30 Days of CV vaccination",
       x = "Present",
       y = "Density") +
  theme(legend.position = "top")+
  xlab("RV Presence")


vacc_encouter_data %>%
  ggplot(aes(x = present, group = within_30_days_sh1, fill = within_30_days_sh1)) +
  geom_histogram(aes(y = ..density..), bins = 10, position = "dodge", alpha = 0.7) +
  facet_wrap(~within_30_days_sh1) +
  scale_fill_viridis_d(option = "magma", begin = 0.2, end = 0.8) +  # Use viridis color palette
  theme_classic() +
  labs(title = "Distribution of 'present' Within 30 Days",
       x = "Present",
       y = "Density") +
  theme(legend.position = "top")+
  xlab("RV Presence")
 



```



# Hypothesis: Ct value on y axis for covid and RV ct value on y axis for rhinovirus


##Data for CT
```{r}

#selecting variables of interest for the full data
full_data_with_ct = full_data1 %>%
  #filter(!is.na(week)) %>%
  dplyr::select(individual, organism, present, sample , site, date,ct)  %>%
  mutate(present = as.character(present))


#Merging multiple encounters with original tracking data


```

```{r}

encounter_with_multiple_data_rv_ct = full_data_with_ct %>%
  filter(organism == "Rhinovirus") %>%
  dplyr::select(individual, present, date, ct ) %>%
  rename(present_rv = present,
         encountered_date_rv = date,
         rv_ct = ct) %>%
  na.omit("individual")

encounter_with_multiple_data_cv_ct = full_data_with_ct %>%
  filter(organism == "SARS_CoV_2") %>%
  dplyr::select(individual, present, date, ct ) %>%
  rename(present_cv = present,
         encountered_date_cv = date,
         covid_ct = ct) %>%
  na.omit("individual")


encounter_rv_cv_pairs_ct = encounter_with_multiple_data_rv_ct %>%
  inner_join(encounter_with_multiple_data_cv_ct)

encounter_rv_cv_pairs_ct = encounter_rv_cv_pairs_ct %>%
  #filter(encountered_week_cv > encountered_week_rv) %>%
  mutate(present_cv_n = as.integer(as.logical(present_cv))) %>%
  mutate(present_rv_n = as.integer(as.logical(present_rv))) %>%
  #mutate(rv_cv_diff = encountered_week_cv- encountered_week_rv)
  mutate(rv_cv_diff = encountered_date_rv- encountered_date_cv)


encounter_rv_cv_pairs_ct$rv_cv_diff <- as.numeric(encounter_rv_cv_pairs_ct$rv_cv_diff, units = "days")


#adding confounders

encounter_rv_cv_pairs_conf_ct = encounter_rv_cv_pairs_ct %>%
  inner_join(full_data_conf_clean, by = "individual")

#Defining periods

encounter_rv_cv_pair_2weeks = encounter_rv_cv_pairs_conf_ct %>%
  filter(rv_cv_diff < 14 &
           rv_cv_diff > 0)


encounter_rv_cv_pair_4weeks = encounter_rv_cv_pairs_conf %>%
  filter(rv_cv_diff < 30 &
           rv_cv_diff > 0)

encounter_rv_cv_pair_8weeks_ct = encounter_rv_cv_pairs_ct %>%
  filter(rv_cv_diff < 60 &
           rv_cv_diff > -30)


rv_cv_tidy = encounter_rv_cv_pair_8weeks_ct %>%
  select(present_cv_n, rv_cv_diff, rv_ct) %>%
  na.omit()

write_tsv(rv_cv_tidy, file = "rv_cv_tidy.txt")

#Running models

#within 14 days
regress("odds", formula = present_cv_n~present_rv_n*rv_cv_diff + sex + Income +Race+age , data = encounter_rv_cv_pair_2weeks)
#within 30 days
regress("odds", formula = present_cv_n~present_rv_n*rv_cv_diff + sex +age+ Race + Income, data = encounter_rv_cv_pair_4weeks)

#within 60 days
regress("odds", formula = present_cv_n~present_rv_n*rv_cv_diff + sex +age+ Race + Income, data = encounter_rv_cv_pair_8weeks_ct)

#x axis rv_Cs_diff faceting by presence rv, y = presence for cv

  
table(encounter_rv_cv_pair_2weeks$present_rv_n)

table(encounter_rv_cv_pair_8weeks_ct$rv_cv_diff)
encounter_rv_cv_pair_8weeks_ct %>%
  ggplot(aes(x = rv_cv_diff, y = rv_ct, color = as.factor(present_cv_n))) +
  geom_point()+
  facet_wrap(~present_cv_n) +
  #geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +  # Use logistic regression line
  theme_minimal() +  # Use a minimal theme
  labs(title = "Association Between RV and COVID Infections within 2 weeks",
       x = "Time Since RV Infection (Days)",
       y = "COVID Infection Probability",
       color = "RV Presence") +
  scale_color_manual(values = c("0" = "blue", "1" = "red")) +  # Customize color scale
  theme(legend.title = element_text(face = "bold", size = 10),
        legend.text = element_text(size = 8),
        axis.title = element_text(face = "bold", size = 10),
        axis.text = element_text(size = 8),
        plot.title = element_text(face = "bold", size = 12)) 

encounter_rv_cv_pair_8weeks_ct %>%
  filter(present_rv_n == "1" & present_cv_n == "1")

```



```{r}
# Create a new variable for time periods
library(ggplot2)

# Create a new variable for time periods
encounter_rv_cv_pair_8weeks_ct <- transform(encounter_rv_cv_pair_8weeks_ct,
                                            time_period = cut(rv_cv_diff,
                                                             breaks = c(0, 5, 10, 15, 20, 30, 40, 50, 60),
                                                             labels = c("0-4 days", "5-9 days", "10-14 days", "15-20 days", "20-30 days", "30-40 days", "40-50 days", "50-60 days")))

# Create the violin plot with mean points
ggplot(encounter_rv_cv_pair_8weeks_ct, aes(x = time_period, y = rv_ct, fill = as.factor(present_cv_n))) +
  geom_violin(trim = FALSE, scale = "width", width = 0.8, show.legend = FALSE) +
  stat_summary(fun.data = "mean_cl_boot", geom = "point", shape = 18, size = 3, color = "black") +
  labs(title = "Violin Plot of covid_ct for Different Time Periods",
       x = "Time Period",
       y = "rv_ct") +
  facet_wrap(~present_cv_n) +
  theme_minimal() +
  ylim(0,60)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

```