CentralitySharing.Rmd

---
title: "Centrality Sharing"
author: "Jackson Pullman"
date: "2022-12-23"
output: html_notebook
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

```

```{r}
library(ggplot2)
library(ggpubr)
library(ggplotify)
```


Calculate Diversity
```{r}
#Calculate diversity indices
richness  <- rowSums(t(species_data) > 0)
shannon <- vegan::diversity(t(species_data), index = "shannon")
invsimpson <- vegan::diversity(t(species_data), index = "invsimpson")
diversity_df <- data.frame(respondent_master_id = names(richness),
                           richness = unname(richness),
                           shannon = unname(shannon),
                           invsimpson = unname(invsimpson))


```

Calculate Centrality and Similarity
```{r}
count <- 0
for(i in 1:length(village_names)){
  sn_vil <- SN %>% filter(village_code_w3 == village_names[i])
  sn_vil_all <- all_edges %>% filter(village_code_w3 == village_names[i])
  village_ids <- unique(c(sn_vil$ego, sn_vil$alter))
  count <- count + length(village_ids)
  sn_vil <- simplify(graph_from_data_frame(sn_vil, directed = FALSE))
  sn_vil_all <- simplify(graph_from_data_frame(sn_vil_all, directed = FALSE))

  eigen_vil <- eigen_centrality(sn_vil)$vector
  between_vil <- igraph::betweenness(sn_vil, normalized = TRUE)
  degree_vil <- igraph::degree(sn_vil)
  clustering_coef <- igraph::transitivity(sn_vil, type = "local", isolates = 'zero')
  names(clustering_coef) <- names(V(sn_vil))

  eigen_vil_all <- eigen_centrality(sn_vil_all)$vector
  eigen_vil_all <- eigen_vil_all[match(names(eigen_vil), names(eigen_vil_all))]
  between_vil_all <- igraph::betweenness(sn_vil_all, normalized = TRUE)
  between_vil_all <- between_vil_all[match(names(between_vil), names(between_vil_all))]
  degree_vil_all <- igraph::degree(sn_vil_all)
  degree_vil_all <- degree_vil_all[match(names(degree_vil), names(degree_vil_all))]
  clustering_coef_all <- igraph::transitivity(sn_vil_all, type = "local", isolates = 'zero')
  names(clustering_coef_all) <- names(V(sn_vil_all))
  clustering_coef_all <- clustering_coef_all[match(names(clustering_coef), names(clustering_coef_all))]

  
  strain_vil <- strain_rate[rownames(strain_rate) %in% village_ids,
                            colnames(strain_rate) %in% village_ids]
  diag(strain_vil) <- NA
  strain_sim <- colMeans(strain_vil,na.rm = TRUE)
  
  strain_sim <- strain_sim[match(names(eigen_vil), names(strain_sim))]
  
  if(i == 1){
    centrality_df <- data.frame(respondent_master_id = names(V(sn_vil)),
                                eigen_c = unname(eigen_vil),
                                between_c = unname(between_vil),
                                degree_c = unname(degree_vil),
                                cluster_coef = unname(clustering_coef),
                                eigen_c_all = unname(eigen_vil_all),
                                between_c_all = unname(between_vil_all),
                                degree_c_all = unname(degree_vil_all),
                                cluster_coef_all = unname(clustering_coef_all),
                                village = rep(village_names[i], length(eigen_vil)),
                                strain_sim_all = unname(strain_sim)
                                )
  }
  else{
    centrality_df <- rbind(centrality_df,
                                data.frame(respondent_master_id = names(V(sn_vil)),
                                eigen_c = unname(eigen_vil),
                                between_c = unname(between_vil),
                                degree_c = unname(degree_vil),
                                cluster_coef = unname(clustering_coef),
                                eigen_c_all = unname(eigen_vil_all),
                                between_c_all = unname(between_vil_all),
                                degree_c_all = unname(degree_vil_all),
                                cluster_coef_all = unname(clustering_coef_all),
                                village = rep(village_names[i], length(eigen_vil)),
                                strain_sim_all = unname(strain_sim)
                                ))
  }
}

```


Merge data to SN
```{r}

diversity_df_ego <- diversity_df
names(diversity_df_ego) <- paste0(names(diversity_df_ego),"_ego")
names(diversity_df_ego)[1] <- "ego"
diversity_df_alter <- diversity_df
names(diversity_df_alter) <- paste0(names(diversity_df_alter),"_alter")
names(diversity_df_alter)[1] <- "alter"

SN_richness <- merge(SN, diversity_df_ego) %>% merge(diversity_df_alter)

SN_Graph_all <- graph_from_data_frame(SN_richness %>% distinct(pair_key,.keep_all = TRUE),
                                      directed = FALSE)
sim_deg_one <- c()
name_vec <- c()

for( i in 1:length(V(SN_Graph_all))){
  test_incident <- incident(SN_Graph_all, V(SN_Graph_all)[i])
  name_vec <-c(name_vec, names(V(SN_Graph_all)[i]))
  sim_deg_one <- c(sim_deg_one, mean(edge_attr(SN_Graph_all, "strain_sharing_rate", test_incident)))
  vert_neighbors <- names(neighbors(SN_Graph_all, V(SN_Graph_all)[i]))
}


sim_deg_one_df <- data.frame(sim_deg_one = sim_deg_one,
                            respondent_master_id = name_vec)

centrality_sim_df <- merge(centrality_df,sim_deg_one_df)
centrality_sim_df <- merge(centrality_sim_df, diversity_df)
centrality_sim_df_covars <- merge(centrality_sim_df, covariates)


centrality_sim_df_covars_all <- inner_join(extraction_batch, sequencing_depth, by=join_by(SampleID)) |> 
  inner_join(shipping_batch, by=join_by(SampleID)) |> 
  inner_join(centrality_sim_df_covars, by=join_by(SampleID==mb_m05))
```


Calculate average similarity to first degree by clustering coefficient
```{r}
#Want to filter to people who have more than one connection, that might be driving it
median(centrality_sim_df$sim_deg_one[centrality_sim_df$cluster_coef_all <= .25], na.rm = TRUE)

median(centrality_sim_df$sim_deg_one[centrality_sim_df$cluster_coef_all >= .25 &
                                       centrality_sim_df$cluster_coef_all <= .5], na.rm = TRUE)

median(centrality_sim_df$sim_deg_one[centrality_sim_df$cluster_coef_all >= .5 &
                                       centrality_sim_df$cluster_coef_all <= .75], na.rm = TRUE)
median(centrality_sim_df$sim_deg_one[centrality_sim_df$cluster_coef_all >= .75], na.rm = TRUE)

centrality_sim_df$cluster_bucket <- case_when(
  centrality_sim_df$cluster_coef_all <= .25 ~ "0 - .25",
  centrality_sim_df$cluster_coef_all <= .5 ~ ".25 - .5",
  centrality_sim_df$cluster_coef_all <= .75 ~ ".5 - .75",
  centrality_sim_df$cluster_coef_all <= 1 ~ ".75 - 1"
  )

cluster_pvals <- compare_means(sim_deg_one ~ cluster_bucket,
                               centrality_sim_df,
                               p.adjust.method = "BH")

cluster_pvals <- cluster_pvals %>%
  filter(p.signif != "ns") %>%
  mutate(y.position = c(25, 27.5,30,32.5,35))

cluster_pvals$p.adj <- paste0("p = ",cluster_pvals$p.adj)

cluster_bucket_fig <- ggboxplot(
  centrality_sim_df,
  x = "cluster_bucket",
  y = "sim_deg_one",
  outlier.shape = NA,
  color = "cluster_bucket",
  show.legend = FALSE
) +
  # geom_jitter(
  #   aes(color = cluster_bucket),
  #   alpha = .5,
  #   width = .25,
  #   show.legend = FALSE
  # ) +
  stat_pvalue_manual(cluster_pvals, label = "p.adj", tip.length = 0) +
  xlab("Clustering Coefficient") +
  ylab("Average Strain-Sharing Rate\nWith First Degree Connections (%)") +
  theme_pubr() +
  labs_pubr() +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0,40))
```

Make some plots
```{r}


richness_plot <-
  ggplot(centrality_sim_df, aes(x = richness_deg_one, y = richness)) +
  geom_point() +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  ylab("Host Species Richness") +
  xlab("Average Richness of Connections 1st Degree Connections") +
  theme_minimal()

deg_c_sharing_all <-
  ggplot(centrality_sim_df, aes(x = degree_c, y = strain_sim_all)) +
  geom_point() +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  ylab("Average Strain-Sharing Rate\nWith All Village Members (%)") +
  xlab("Degree Centrality") +
  theme_pubr() +
  labs_pubr()


richness_deg_c <- ggplot(centrality_sim_df, aes(x = degree_c, y = richness)) +
  geom_point() +
  ylab("Host Species Richness") +
  xlab("Degree Centrality") +
  theme_minimal()

shannon_deg_c <- ggplot(centrality_sim_df, aes(x = degree_c, y = shannon)) +
  geom_point() +
  ylab("Host Shannon Diversity") +
  xlab("Degree Centrality") +
  theme_minimal()

invsimpson_deg_c <- ggplot(centrality_sim_df, aes(x = degree_c, y = invsimpson)) +
  geom_point() +
  ylab("Host Inverse Simpson Diversity") +
  xlab("Degree Centrality") +
  theme_minimal()



deg_c_sharing_first <-
  ggplot(centrality_sim_df, aes(x = degree_c, y = sim_deg_one)) +
  geom_point() +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  ylab("Average Strain-Sharing Rate\nWith First Degree Connections (%)") +
  xlab("Degree Centrality") +
  theme_pubr() +
  labs_pubr()


clustering_coef_fig <-
  ggplot(centrality_sim_df, aes(x = cluster_coef_all, y = sim_deg_one)) +
  geom_point() +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  ylab("Average Strain-Sharing Rate with First Degree Connections") +
  xlab("Degree Centrality") +
  theme_minimal()


species_pop_deg_c <-
  ggplot(centrality_sim_df, aes(x = degree_c, y = avg_species_prev)) +
  geom_point() +
  stat_smooth(method = "lm",
              formula = y ~ x,
              geom = "smooth") +
  ylab("Average Prevalence of Host Species") +
  xlab("Degree Centrality") +
  theme_minimal()

```
Test for correlation, Figure 4B,D
```{r}
broom::tidy(cor.test(centrality_sim_df$degree_c, centrality_sim_df$strain_sim_all)) |> 
  select(cor=estimate, p.value) |> 
  mutate(test='degree_c~strain_sim_all') |> 
  bind_rows(
broom::tidy(cor.test(centrality_sim_df$degree_c, centrality_sim_df$sim_deg_one)) |> 
  select(cor=estimate, p.value) |> 
  mutate(test='degree_c~strain_sim_first_degree'))

```


No diversity reelationship
```{r}
svglite("../Figures/SFigure10/sfigure10_full.svg")

ggarrange(richness_deg_c, shannon_deg_c, invsimpson_deg_c,
          labels = c("A", "B", "C"), nrow = 1)
dev.off()


```

```{r}
set.seed(22)
test <- sn_vil_graph_15


vertex_attr(test)$similarity <- centrality_sim_df$strain_sim_all[match(names(V(test)),
                                           centrality_sim_df$respondent_master_id)]


palette <- heat.colors(3)

palette <- rev(palette)


#Trichotomize sharing
vertex_attr(test)$color <- case_when(
  vertex_attr(test)$similarity <= quantile(vertex_attr(test)$similarity, .33) ~ palette[1],
  vertex_attr(test)$similarity <= quantile(vertex_attr(test)$similarity, .66) ~ palette[2],
  TRUE ~ palette[3],
)

vertex_attr(test)$size <- (centrality_sim_df$degree_c[match(names(V(test)),
                                           centrality_sim_df$respondent_master_id)] +1)/3

par(mar=c(0,0,0,0))
set.seed(22)
plot(test,
     vertex.label = NA,
     layout = layout_with_fr)
legend('topleft',
       legend = c(">= 4.6% (Top Tertile)",
                  ">= 3.9% (Middle Tertile)",
                  ">= 0% (Bottom Tertile)"),
       #fill = rev(palette),
       pt.bg  = rev(palette),
       pch    = 21,
       cex    = 1,
       border = NA,
       bty    = "n",
       title = "Average Strain-Sharing Rate With\nAll Village Members")

set.seed(22)
strain_sim_cent_network <- as.ggplot(expression(
  plot(
    test,
    vertex.label = NA,
    layout = layout_with_fr(test)
  ),
  legend(x = -.5,
         y = 1.4,
    #'topleft',
    legend = c(">= 4.6% (Top Tertile)",
                  ">= 3.9% (Middle Tertile)",
                  ">= 0% (Bottom Tertile)"),
    pt.bg  = rev(palette),
    pch    = 21,
    cex    = .75,
    bty    = "n",
    y.intersp = 1.25,
    title = "Average Strain-Sharing Rate With\nAll Village Members"
  )
))

set.seed(22)

vertex_attr(test)$similarity <- centrality_sim_df$sim_deg_one[match(names(V(test)),
                                           centrality_sim_df$respondent_master_id)]
palette <- heat.colors(3)
palette <- rev(palette)

vertex_attr(test)$color <- palette[order(vertex_attr(test)$similarity)]

#Trichotomize sharing
vertex_attr(test)$color <- case_when(
  vertex_attr(test)$similarity <= quantile(vertex_attr(test)$similarity, .33) ~ palette[1],
  vertex_attr(test)$similarity <= quantile(vertex_attr(test)$similarity, .66) ~ palette[2],
  TRUE ~ palette[3],
)

vertex_attr(test)$size <- (centrality_sim_df$degree_c[match(names(V(test)),
                                           centrality_sim_df$respondent_master_id)] +1)/3

set.seed(22)
plot(test,
     vertex.label = NA,
#     vertex.size = 5,
     layout = layout_with_fr(test))
legend('topleft',
       legend = c(">= 9.7% (Top Tertile)",
                  ">= 7.3% (Middle Tertile)",
                  ">= 0% (Bottom Tertile)"),
       pt.bg  = rev(palette),
       pch    = 21,
       cex    = 1,
       bty    = "n",
       title = "Average Strain-Sharing Rate With\nFirst Degree Connections")

set.seed(22)
sim_deg_one_network <- as.ggplot(expression(
  plot(
    test,
    vertex.label = NA,
    layout = layout_with_fr
  ),
  legend(x = -.5,
         y = 1.4,
    #'topleft',
    legend = c(">= 9.7% (Top Tertile)",
                  ">= 7.3% (Middle Tertile)",
                  ">= 0% (Bottom Tertile)"),
    pt.bg  = rev(palette),
    pch    = 21,
    cex    = .75,
    bty    = "n",
    y.intersp = 1.25,
    title = "Average Strain-Sharing Rate With\nFirst Degree Connections"
  )
))

```

```{r}


sim_scatter_fig <- ggarrange(deg_c_sharing_all,
          deg_c_sharing_first,
          nrow = 2,
          labels = c("B", "D"))

sim_network_fig <- ggarrange(strain_sim_cent_network,
          sim_deg_one_network,
          nrow = 2,
          labels = c("C", "E"))

#This comes from Sharing_by_degree
ssr_by_degree_plot <- ggarrange(ssr_by_degree_plot, labels = "A")

svglite("../FiguresNew/Figure4/full_fig4_test.svg",
        width = 14,
        height = 8)
ggarrange(ssr_by_degree_plot,sim_scatter_fig, sim_network_fig,nrow = 1)
dev.off()

ggarrange(ssr_by_degree_plot,sim_scatter_fig, sim_network_fig,nrow = 1)
```




```{r}
sub_covar_centrality <- centrality_sim_df_covars_all |> select(BMI, age_at_survey, household_wealth_index_w3, gender, DDS, starts_with('med_'), l0400, mb_m0300)
sub_covar_centrality$gender <- ifelse(sub_covar_centrality$gender == 'man', 1, 0)
sub_covar_centrality$mb_m0300 <- as.integer(sub_covar_centrality$mb_m0300)
sub_covar_centrality$l0400 <- as.integer(sub_covar_centrality$l0400)

saveRDS(centrality_sim_df_covars_all, 'centrality_sim_df_covars_all.RDS')
saveRDS(sub_covar_centrality, 'sub_covar_centrality.RDS')
```

# Social Network Position and Strain Network Position

## Simple
### Average strain-sharing rate to the rest of the villages, betweeness centrality
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ between_c_all + age_at_survey +household_wealth_index_w3 +
                         gender + DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

### Average strain-sharing rate to the rest of the village, eigenvector centrality
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ eigen_c_all + age_at_survey +household_wealth_index_w3 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

### Average strain-sharing rate to the rest of the village, degree
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ degree_c_all + age_at_survey +household_wealth_index_w3 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads + (1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```


### First degree connections, clustering coefficient
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ cluster_coef_all + age_at_survey +household_wealth_index_w3 +
                         gender +  DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village_name) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

### First degree connections, eigenvector centrality
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ eigen_c_all + age_at_survey +household_wealth_index_w3 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

### First degree connections, betweeness centrality
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ between_c_all + age_at_survey +household_wealth_index_w3 +
                         gender+  DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

### First degree connections, degree
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ degree_c_all + age_at_survey +household_wealth_index_w3 +
                         gender+  DNA_conc + BMI + l0400 + mb_m0300 + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

## Separate covariates

## Average strain-sharing rate to the rest of the villages, betweeness centrality, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ between_c_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender + DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

## Average strain-sharing rate to the rest of the village, eigenvector centrality, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ eigen_c_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

## Average strain-sharing rate to the rest of the village, degree, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ degree_c_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads + (1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```


## First degree connections, clustering coefficient, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ cluster_coef_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender +  DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village_name) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

## First degree connections, eigenvector centrality, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ eigen_c_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

## First degree connections, betweeness centrality, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ between_c_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender+  DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```

## First degree connections, degree, use separate covariates
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ degree_c_all + age_at_survey +household_wealth_index_w3 +DDS +med_painkillers +med_antibiotics+ med_antidiarrheal+ med_antiparasite+ med_vitamins+ med_zinc+ med_antifungal+ med_antihypertensives+ med_antidiabetics+ med_antiacids+ med_laxatives+ l0400+ mb_m0300 +
                         gender+  DNA_conc + BMI + l0400 + mb_m0300 + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```


## Mahalanobis
## Average strain-sharing rate to the rest of the villages, betweeness centrality, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ between_c_all + maha_social + DNA_conc + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
)
summary(lm.1)
```

## Average strain-sharing rate to the rest of the village, eigenvector centrality, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ eigen_c_all + maha_social + DNA_conc + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
)
summary(lm.1)
```

## Average strain-sharing rate to the rest of the village, degree, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(strain_sim_all ~ degree_c_all + maha_social + DNA_conc + sampling_date + shipment_batch + Plate + WIKI + Reads + (1|village) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
)
summary(lm.1)
```


## First degree connections, clustering coefficient, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ cluster_coef_all + maha_social + DNA_conc + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village_name) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
)
summary(lm.1)
```

## First degree connections, eigenvector centrality, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ eigen_c_all + maha_social + DNA_conc + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
)
summary(lm.1)
```

## First degree connections, betweeness centrality, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ between_c_all + maha_social + DNA_conc  + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
)
summary(lm.1)
```

## First degree connections, degree, use Mahalanobis distance
```{r}
lm.1 <- lmerTest::lmer(sim_deg_one ~ degree_c_all + maha_social + DNA_conc + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = bind_cols(centrality_sim_df_covars_all, maha_social = mahalanobis(sub_covar_centrality, FALSE, var(sub_covar_centrality, na.rm=TRUE)))
             )

summary(lm.1)
```



## Inverse Simpson diversity, betweeness centralitys
```{r}
lm.1 <- lmerTest::lmer(invsimpson ~ between_c_all + age_at_survey +household_wealth_index_w3 +
                         gender+ DNA_conc + BMI + sampling_date + shipment_batch + Plate + WIKI + Reads +(1|village) + (1|building_id),
             data = centrality_sim_df_covars_all)
summary(lm.1)
```




See if women share more to their direct contacts than men
```{r}
library(ggplot2)


SN_Graph_all <- graph_from_data_frame(SN %>% distinct(pair_key, .keep_all = TRUE),
                                               directed = FALSE)
woman_list_strain <- c()
man_list_strain <- c()

woman_list_bray <- c()
man_list_bray <- c()

woman_list_jaccard <- c()
man_list_jaccard <- c()

sum(hmb_census$respondent_master_id == names(V(SN_Graph_all)[i]))

for( i in 1:length(V(SN_Graph_all))){
  test_incident <- incident(SN_Graph_all, V(SN_Graph_all)[i])
  
  if(hmb_census$gender[hmb_census$respondent_master_id == names(V(SN_Graph_all)[i])] == "woman") {
    woman_list_strain <- c(woman_list_strain,
                    mean(edge_attr(SN_Graph_all, "strain_sharing_rate", test_incident)))
    woman_list_bray <- c(woman_list_bray,
                    mean(edge_attr(SN_Graph_all, "bray_curtis_sim", test_incident)))
    woman_list_jaccard <- c(woman_list_jaccard,
                          mean(edge_attr(SN_Graph_all, "jaccard_sim", test_incident)))
  }
  
  else{
     man_list_strain <- c(man_list_strain,
                   mean(edge_attr(SN_Graph_all, "strain_sharing_rate", test_incident)))
     man_list_bray <- c(man_list_bray,
                    mean(edge_attr(SN_Graph_all, "bray_curtis_sim", test_incident)))
     man_list_jaccard <- c(man_list_jaccard,
                    mean(edge_attr(SN_Graph_all, "jaccard_sim", test_incident)))
  }
}


man_woman_df <- data.frame("Strain_Sharing_Rate" = c(man_list_strain, woman_list_strain),
                           "Bray_Curtis_Similarity" = c(man_list_bray, woman_list_bray),
                           "Jaccard_Similarity" = c(man_list_jaccard, woman_list_jaccard),
                           "Gender" = c(rep("Male", length(man_list_strain)),
                                        rep("Female", length(woman_list_strain))))

man_woman_df.m <- man_woman_df |>
  pivot_longer(-Gender, names_to = 'variable', values_to = 'value')

medians_mw <- aggregate(value ~ variable + Gender, data = man_woman_df.m, median, na.rm = TRUE)

#Want all pairwise comparisons
man_woman_df.m$all_group <- paste0(man_woman_df.m$Gender, man_woman_df.m$variable)
man_woman_df.m$variable <- factor(man_woman_df.m$variable, levels = c('Strain_Sharing_Rate','Bray_Curtis_Similarity','Jaccard_Similarity'))
mw_share_pvals  <- compare_means(value ~ all_group, data = man_woman_df.m)

mw_share_pvals <- mw_share_pvals[c(3, 8, 12), ]

# If you want the two levels of event plotted side by side
mw_share  <- ggplot(man_woman_df.m, aes(Gender, value, colour=variable)) +
  facet_wrap(~variable, scales = "free") +
  geom_boxplot(width=0.7, show.legend = FALSE) +
  theme_minimal() +
  stat_compare_means(comparisons = list(c("Female", "Male")),
    method = "wilcox.test",
    size = 3,
    tip.length = 0) +
  geom_text(data = medians_mw,
            aes(label = sprintf("%0.2f", round(value, digits = 2)),
                y = value + 0.02),
            show.legend = FALSE)+
  theme(axis.title.y =element_blank(),
        axis.title.x =element_blank(),
        plot.title = element_text(hjust = .5))+
  ggtitle("Average Microbiome Similarity to Direct Connections by Gender")

# mw_share <- ggplot(man_woman_df, aes(x = Gender, fill = Gender)) +
#   geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
#   facet_grid(rows = vars(Strain_Sharing_Rate, Bray_Curtis_Similarity, Jaccard_Similarity)) +
#   coord_cartesian(ylim = c(0, 1)) +
#   #stat_compare_means(method = "wilcox.test",label.y = .3) +
#   theme_minimal() +
#   ylab("Average Strain-Sharing Rate to Direct Links") +
#   ggtitle("Average Strain-Sharing Rate to Direct Links") +
#   theme(plot.title = element_text(hjust = .5, size = 10))
```




```{r}
#Remove spouse_pairs_inherently gendered
spouse_pairs <- SN %>% filter(relationship == "Partner") %>% pull(pair_key) %>% unique()
SN_no_spouse <- SN %>% filter(!(pair_key %in% spouse_pairs)) %>% distinct(pair_key, .keep_all = TRUE)


SN_no_spouse$pair_gender <- paste0(SN_no_spouse$gender_ego, SN_no_spouse$gender_alter)

SN_no_spouse$pair_gender <- case_when(
  SN_no_spouse$pair_gender == "womanman" | SN_no_spouse$pair_gender == "manwoman" ~ "Female - Male",
  SN_no_spouse$pair_gender == "manman" ~ "Male - Male",
  SN_no_spouse$pair_gender == "womanwoman" ~ "Female - Female"
)

SN$pair_gender <- paste0(SN$gender_ego, SN$gender_alter)

SN$pair_gender <- case_when(
  SN$pair_gender == "womanman" | SN$pair_gender == "manwoman" ~ "Female - Male",
  SN$pair_gender == "manman" ~ "Male - Male",
  SN$pair_gender == "womanwoman" ~ "Female - Female"
)

my_comparisons <- list( c("Female - Male", "Male - Male"),
                        c("Female - Male", "Female - Female"),
                        c("Male - Male", "Female - Female") )

means_gender <- aggregate(strain_sharing_rate ~  pair_gender, SN_no_spouse, median, na.rm = TRUE)



```




```{r}
SN_tri_gender <- SN_no_spouse %>%
  filter(!(relationship %in% c("Partner", "Mother", "Father", "Child"))) %>%
  distinct(pair_key, relationship, .keep_all = TRUE)


ann_text <- aggregate(strain_sharing_rate ~  relationship + pair_gender,
                     SN_tri_gender, median, na.rm = TRUE)
ann_text_2 <- aggregate(strain_sharing_rate ~  relationship + pair_gender,
                     SN_tri_gender, length)
ann_text_3 <- aggregate(strain_sharing_rate ~  relationship + pair_gender,
                     SN_tri_gender, median, na.rm = TRUE)

names(ann_text_2) <- c("relationship", "pair_gender", "count")

ann_text <- merge(ann_text, ann_text_2)

ann_text$lab <- paste0(as.character(round(ann_text$strain_sharing_rate,3)),
                       "\n(N = ", ann_text$count, ")")


gender_split_share_ssr <- ggboxplot(SN_tri_gender, color = "pair_gender",
                                x = "pair_gender", y = "strain_sharing_rate",
                                facet.by = "relationship",outlier.shape = NA) +
  coord_cartesian(ylim = c(0,35)) +
  theme(#axis.text.x = element_text(angle = 45, hjust=1),
    axis.text.x = element_blank(),
    legend.title = element_blank(),
    legend.position= "bottom",
    legend.margin=margin(0,0,0,0, unit='cm')) +
  stat_compare_means(label.y = c(20,25,30), tip.length = 0,
                     comparisons = my_comparisons) +
  xlab("") +
  ylab("Strain-Sharing Rate") +
  geom_text(data = ann_text, label = ann_text$lab,
            size = 3)

```


```{r}
ann_text <- aggregate(bray_curtis_sim ~  relationship + pair_gender,
                     SN_tri_gender, median)
ann_text_2 <- aggregate(bray_curtis_sim ~  relationship + pair_gender,
                     SN_tri_gender, length)


names(ann_text_2) <- c("relationship", "pair_gender", "count")

ann_text <- merge(ann_text, ann_text_2)

ann_text$lab <- paste0(as.character(round(ann_text$bray_curtis_sim,3)),
                       "\n(N = ", ann_text$count, ")")


gender_split_share_bc <- ggboxplot(SN_tri_gender, color = "pair_gender",
                                x = "pair_gender", y = "bray_curtis_sim",
                                facet.by = "relationship",outlier.shape = NA) +
  theme(#axis.text.x = element_text(angle = 45, hjust=1),
    axis.text.x = element_blank(),
    legend.title = element_blank(),
    legend.position= "bottom",
    legend.margin=margin(0,0,0,0, unit='cm')) +
  stat_compare_means(label.y = c(.85,.9,.95),
                     comparisons = my_comparisons, tip.length = 0,
                     method.args = list(method.args = list(p.adjust.method = "BH"))) +
  xlab("") +
  ylab("Bray-Curtis Similarity") +
  geom_text(data = ann_text, label = ann_text$lab,
            size = 3)
```


```{r}
ann_text <- aggregate(jaccard_sim ~  relationship + pair_gender,
                     SN_tri_gender, median)
ann_text_2 <- aggregate(jaccard_sim ~  relationship + pair_gender,
                     SN_tri_gender, length)


names(ann_text_2) <- c("relationship", "pair_gender", "count")

ann_text <- merge(ann_text, ann_text_2)

ann_text$lab <- paste0(as.character(round(ann_text$jaccard_sim,3)),
                       "\n(N = ", ann_text$count, ")")

compare_means(jaccard_sim ~  pair_gender,
                     SN_tri_gender, p.adjust.method = "BH", )


gender_split_share_jac <- ggboxplot(SN_tri_gender, color = "pair_gender",
                                x = "pair_gender", y = "jaccard_sim",
                                facet.by = "relationship",outlier.shape = NA) +
  theme(#axis.text.x = element_text(angle = 45, hjust=1),
    axis.text.x = element_blank(),
    legend.title = element_blank(),
    legend.position= "bottom",
    legend.margin=margin(0,0,0,0, unit='cm')) +
  stat_compare_means(label.y = c(.78,.8,.85),
                     tip.length = 0,
                     comparisons = my_comparisons,
                     method = 'wilcox.test',
                     method.args = list(method.args = list(p.adjust.method = "BH"))) +
  xlab("") +
  ylab("Jaccard Similarity") +
  geom_text(data = ann_text, label = ann_text$lab,
            size = 3)

```


Make plot
```{r}

gender_fig <- ggarrange(ggarrange(gender_split_share_jac, gender_split_share_bc, labels = c('A','B')),
ggarrange(gender_split_share_ssr, mw_share, labels = c('C','D')), nrow = 2)
          
svglite('Figures/supp_fig9.svg', fix_text_size = FALSE, height = 16, width = 12)
gender_fig
dev.off()

```

Mother-Child sharing by father sharing and age differential

```{r}
SN_parent <- SN %>%
  filter(relationship %in% c("Child","Mother", "Father"))

SN_parent <- SN_parent %>% distinct(pair_key, .keep_all = TRUE)

for(i in 1:nrow(SN_parent)){
  if(SN_parent$relationship[i] == "Mother"){
    SN_parent$relationship[i] <- "Child"
    SN_parent$gender_ego[i] <- "Mother"
    SN_parent$age_at_survey_alter[i] <- SN_parent$age_at_survey_ego[i]
  }
  else if(SN_parent$relationship[i] == "Father"){
    SN_parent$relationship[i] <- "Child"
    SN_parent$gender_ego[i] <- "Father"
    SN_parent$age_at_survey_alter[i] <- SN_parent$age_at_survey_ego[i]
  }
}
SN_parent$gender_ego[SN_parent$gender_ego == "woman"] <- "Mother"
SN_parent$gender_ego[SN_parent$gender_ego == "man"] <- "Father"


mc_scatter_ssr <- ggscatter(SN_parent %>% filter(gender_ego == "Mother"),
          x = "age_at_survey_alter", y = "strain_sharing_rate",
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE) +
  coord_cartesian(ylim = c(0,.75)) +
  stat_cor(method = "pearson", label.x = 20, label.y = .75,
           aes(label = paste(after_stat(rr.label), ..p.label.., sep = "~`,`~"))) +
  xlab("Child Age") +
  ylab("Strain-Sharing Rate") +
  ggtitle("Mother-Child Strain-Sharing Rate by Child's Age") +
  theme(plot.title = element_text(size = 10))

mc_scatter_bc <- ggscatter(SN_parent %>% filter(gender_ego == "Mother"),
          x = "age_at_survey_alter", y = "bray_curtis_sim",
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE) + 
  coord_cartesian(ylim = c(0,.75)) +
  stat_cor(method = "pearson", label.x = 20, label.y = .75,
           aes(label = paste(after_stat(rr.label), ..p.label.., sep = "~`,`~"))) +
  xlab("Child Age") +
  ylab("Bray-Curtis Similarity") +
  ggtitle("Mother-Child Bray-Curtis Similarity by Child's Age")+
  theme(plot.title = element_text(size = 10))
  


mc_scatter_j <- ggscatter(SN_parent %>% filter(gender_ego == "Mother"),
          x = "age_at_survey_alter", y = "jaccard_sim",
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE) + 
  coord_cartesian(ylim = c(0,.75)) +
  stat_cor(method = "pearson", label.x = 20, label.y = .75,
           aes(label = paste(after_stat(rr.label), ..p.label.., sep = "~`,`~"))) +
  xlab("Child Age") +
  ylab("Jaccard Similarity") +
  ggtitle("Mother-Child Jaccard Similarity by Child's Age")+
  theme(plot.title = element_text(size = 10))



fc_scatter_ssr <- ggscatter(SN_parent %>% filter(gender_ego == "Father"),
          x = "age_at_survey_alter", y = "strain_sharing_rate",
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE) +
  coord_cartesian(ylim = c(0,.75)) +
  stat_cor(method = "pearson", label.x = 20, label.y = .75,
           aes(label = paste(after_stat(rr.label), ..p.label.., sep = "~`,`~"))) +
  xlab("Child Age") +
  ylab("Strain-Sharing Rate") +
  ggtitle("Father-Child Strain-Sharing Rate by Child's Age") +
  theme(plot.title = element_text(size = 10))

fc_scatter_bc <- ggscatter(SN_parent %>% filter(gender_ego == "Father"),
          x = "age_at_survey_alter", y = "bray_curtis_sim",
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE) + 
  coord_cartesian(ylim = c(0,.75)) +
  stat_cor(method = "pearson", label.x = 20, label.y = .75,
           aes(label = paste(after_stat(rr.label), ..p.label.., sep = "~`,`~"))) +
  xlab("Child Age") +
  ylab("Bray-Curtis Similarity") +
  ggtitle("Father-Child Bray-Curtis Similarity by Child's Age")+
  theme(plot.title = element_text(size = 10))
  


fc_scatter_j <- ggscatter(SN_parent %>% filter(gender_ego == "Father"),
          x = "age_at_survey_alter", y = "jaccard_sim",
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE) + 
  coord_cartesian(ylim = c(0,.75)) +
  stat_cor(method = "pearson", label.x = 20, label.y = .75,
           aes(label = paste(after_stat(rr.label), ..p.label.., sep = "~`,`~"))) +
  xlab("Child Age") +
  ylab("Jaccard Similarity") +
  ggtitle("Father-Child Jaccard Similarity by Child's Age")+
  theme(plot.title = element_text(size = 10))


parent_share_ssr <-ggplot(SN_parent, aes(x = gender_ego,
                                     y = strain_sharing_rate,
                                     fill = gender_ego)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  stat_compare_means(label.y = .40) +
  coord_cartesian(ylim = c(0,.50)) +
  ylab("Strain-Sharing Rate (%)") +
  xlab("Parent Gender") +
  ggtitle("Child to Parent Strain-Sharing Rate") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = .5, size = 10)) +
  theme_pubr()

parent_share_bc <-ggplot(SN_parent, aes(x = gender_ego,
                                     y = bray_curtis_sim,
                                     fill = gender_ego)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  stat_compare_means(label.y = .75) +
  coord_cartesian(ylim = c(0,1)) +
  ylab("Bray-Curtis Similarity") +
  xlab("Parent Gender") +
  ggtitle("Child to Parent Bray-Curtis Similarity") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = .5, size = 10)) +
  theme_pubr()

parent_share_j <-ggplot(SN_parent, aes(x = gender_ego,
                                     y = jaccard_sim,
                                     fill = gender_ego)) +
  geom_boxplot(outlier.shape = NA, show.legend = FALSE) +
  stat_compare_means(label.y = .75) +
  coord_cartesian(ylim = c(0,1)) +
  ylab("Jaccard Similarity") +
  xlab("Parent Gender") +
  ggtitle("Child to Parent Jaccard Similarity") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = .5, size = 10)) +
  theme_pubr()


```


```{r}
sfig5_1 <- ggarrange(mc_scatter_ssr, mc_scatter_bc, mc_scatter_j, ncol = 3, labels = c("A", "B", "C"))

sfig5_2 <- ggarrange(fc_scatter_ssr, fc_scatter_bc, fc_scatter_j, ncol = 3, labels = c("D", "E", "F"))

sfig5_3 <- ggarrange(parent_share_ssr, parent_share_bc, parent_share_j, ncol = 3, labels = c("G", "H", "I"))


svglite("../FiguresNew/SFigure6/sfigure6_full.svg",
        width = 12, height = 10)
ggarrange(sfig5_1, sfig5_2,sfig5_3, ncol = 1)
dev.off()

```