Analysis.Rmd

---
title: "NetworkConstruction"
output: html_document
date: "2023-12-07"
---

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

```{r Load data and packages}
#load packages
#install_infomap(target_folder = NULL, source = "binary")
#devtools::install_github('Ecological-Complexity-Lab/infomap_ecology_package', force=T)
library(tidyverse); library(scales);library(muxViz);library(infomapecology);library(bipartite);library(ggpubr);library(igraph);library(metaforest);library(ggnet);library(viridis);library(MuMIn)


```

```{r Edgelist Construction}
#Assign village
village <- metadata %>% select(social_netid, village) %>% rename("id1" = "social_netid") %>% distinct()

#Environmental network
enviro_adj2 <- enviro_adj %>% filter(network!="all", network!="village") %>% group_by(id1, id2, village) %>% summarise(weight=sum(weight))
enviro_adj2 <- enviro_adj2[!grepl("Z", enviro_adj2$id1),]
enviro_adj2 <- enviro_adj2[!grepl("Z", enviro_adj2$id2),]
enviro_vaillage_m_edge <- enviro_adj2 %>% filter(village =="M_Village") %>% select(id1, id2, weight)
enviro_vaillage_s_edge <- enviro_adj2 %>% filter(village =="S_Village") %>% select(id1, id2, weight)
enviro_vaillage_a_edge <- enviro_adj2 %>% filter(village =="A_Village") %>% select(id1, id2, weight) 

#Close Contact network
contact_adj <- left_join(contact_adj, village, by="id1")
contact_adj2 <- contact_adj[!grepl("Z", contact_adj$id1),]
contact_adj2 <- contact_adj2[!grepl("Z", contact_adj2$id2),]
contact_vaillage_m_edge <- contact_adj %>% filter(village =="M_Village") %>% select(id1, id2, close_contact_edge) %>% rename("weight" = "close_contact_edge") %>% filter(weight!=0)
contact_vaillage_s_edge <- contact_adj %>% filter(village =="S_Village") %>% select(id1, id2, close_contact_edge) %>% rename("weight" = "close_contact_edge") %>% filter(weight!=0)
contact_vaillage_a_edge <- contact_adj %>% filter(village =="A_Village") %>% select(id1, id2, close_contact_edge) %>% rename("weight" = "close_contact_edge") %>% filter(weight!=0)

#Naming network

name_undirected_man2 <-name_undirected_man$n.named
name_undirected_vaillage_s2 <-name_undirected_vaillage_s$n.named
name_undirected_vaillage_a2 <-name_undirected_vaillage_a$n.named

name_vaillage_m_edge <- igraph::as_data_frame(name_undirected_man2) %>% rename(id1 = from, id2 = to)
name_vaillage_s_edge <- igraph::as_data_frame(name_undirected_vaillage_s2)%>% rename(id1 = from, id2 = to)
name_vaillage_a_edge <- igraph::as_data_frame(name_undirected_vaillage_a2)%>% rename(id1 = from, id2 = to)

#Household distance network
house_dist <- house_dist[house_dist$dist!=0,]
house_dist$dist <- 1/house_dist$dist

distance_edges <- left_join(house_dist, village, by = "id1")%>% filter(dist !=0)
distance_vaillage_m_edge <- distance_edges %>% filter(village=="M_Village")
distance_vaillage_a_edge <- distance_edges %>% filter(village=="A_Village")
distance_vaillage_s_edge <- distance_edges %>% filter(village=="S_Village")


```

```{r Build Networks: M_Village}
#Setup intra-layer edge lists
enviro_vaillage_m_edge <- enviro_vaillage_m_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =1, layer.to =1) %>% select(node.from, layer.from, node.to, layer.to, weight) %>% mutate(weight = rescale(weight,to = c(0.01,1)))

contact_vaillage_m_edge <- contact_vaillage_m_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =2, layer.to =2)%>% select(node.from, layer.from, node.to, layer.to, weight)%>% mutate(weight = rescale(weight,to = c(0.01,1)))

name_vaillage_m_edge <- name_vaillage_m_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =3, layer.to =3)%>% select(node.from, layer.from, node.to, layer.to, weight)%>% mutate(weight = rescale(weight,to = c(0.01,1)))


ggplot(data= enviro_vaillage_m_edge, aes(x = weight))+ geom_histogram()
ggplot(data= contact_vaillage_m_edge, aes(x = weight))+ geom_histogram()
ggplot(data= name_vaillage_m_edge, aes(x = weight))+ geom_histogram()

mean(enviro_vaillage_m_edge$weight)
mean(contact_vaillage_m_edge$weight)
mean(name_vaillage_m_edge$weight)
mean(c(enviro_vaillage_m_edge$weight, contact_vaillage_m_edge$weight, name_vaillage_m_edge$weight))

name_vaillage_m_edge <- name_vaillage_m_edge %>% filter(node.to %in% node.from)#filter out people named in different villages 

```

```{r Build Networks: S_Village}
#Setup intra-layer edge lists
enviro_vaillage_s_edge <- enviro_vaillage_s_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =1, layer.to =1) %>% select(node.from, layer.from, node.to, layer.to, weight) %>% mutate(weight = rescale(weight,to = c(0.01,1)))
contact_vaillage_s_edge <- contact_vaillage_s_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =2, layer.to =2)%>% select(node.from, layer.from, node.to, layer.to, weight)%>% mutate(weight = rescale(weight,to = c(0.01,1)))
name_vaillage_s_edge <- name_vaillage_s_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =3, layer.to =3)%>% select(node.from, layer.from, node.to, layer.to, weight)%>% mutate(weight = rescale(weight,to = c(0.01,1)))

name_vaillage_s_edge <- name_vaillage_s_edge %>% filter(node.to %in% node.from)#filter out people named in different villages 



```

```{r Build Networks: A_Village}
#Setup intra-layer edge lists
enviro_vaillage_a_edge <- enviro_vaillage_a_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =1, layer.to =1) %>% select(node.from, layer.from, node.to, layer.to, weight) %>% mutate(weight = rescale(weight,to = c(0.01,1)))
contact_vaillage_a_edge <- contact_vaillage_a_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =2, layer.to =2)%>% select(node.from, layer.from, node.to, layer.to, weight)%>% mutate(weight = rescale(weight,to = c(0.01,1)))
name_vaillage_a_edge <- name_vaillage_a_edge %>% rename(node.from =id1, node.to = id2)%>% mutate(layer.from =3, layer.to =3)%>% select(node.from, layer.from, node.to, layer.to, weight)%>% mutate(weight = rescale(weight,to = c(0.01,1)))

name_vaillage_a_edge <- name_vaillage_a_edge %>% filter(node.to %in% node.from)#filter out people named in different villages 

```

```{r Single-layer Centrality}
#GET SINGLE-LAYER GRAPHS
#M_Village
enviro_vaillage_m_g <- graph_from_edgelist(as.matrix(enviro_vaillage_m_edge[,c(1,3)]), directed = FALSE)
name_vaillage_m_g <- graph_from_edgelist(as.matrix(name_vaillage_m_edge[,c(1,3)]), directed = FALSE)
contact_vaillage_m_g <- graph_from_edgelist(as.matrix(contact_vaillage_m_edge[,c(1,3)]), directed = FALSE)
distance_vaillage_m_g <- graph_from_edgelist(as.matrix(distance_vaillage_m_edge[,c(1,2)]), directed = FALSE)

enviro_vaillage_m_g<-simplify(enviro_vaillage_m_g, remove.multiple = FALSE, remove.loops = TRUE)
name_vaillage_m_g<-simplify(name_vaillage_m_g, remove.multiple = FALSE, remove.loops = TRUE)
contact_vaillage_m_g<-simplify(contact_vaillage_m_g, remove.multiple = FALSE, remove.loops = TRUE)
distance_vaillage_m_g<-simplify(distance_vaillage_m_g, remove.multiple = FALSE, remove.loops = TRUE)

E(enviro_vaillage_m_g)$weight <- enviro_vaillage_m_edge$weight
E(name_vaillage_m_g)$weight <- name_vaillage_m_edge$weight
E(contact_vaillage_m_g)$weight <- contact_vaillage_m_edge$weight
E(distance_vaillage_m_g)$weight <- distance_vaillage_m_edge$dist

#S_Village
enviro_vaillage_s_g <- graph_from_edgelist(as.matrix(enviro_vaillage_s_edge[,c(1,3)]), directed = FALSE)
name_vaillage_s_g <-graph_from_edgelist(as.matrix(name_vaillage_s_edge[,c(1,3)]), directed = FALSE)
contact_vaillage_s_g <- graph_from_edgelist(as.matrix(contact_vaillage_s_edge[,c(1,3)]), directed = FALSE)
distance_vaillage_s_g <- graph_from_edgelist(as.matrix(distance_vaillage_s_edge[,c(1,2)]), directed = FALSE)

enviro_vaillage_s_g<-simplify(enviro_vaillage_s_g, remove.multiple = FALSE, remove.loops = TRUE)
name_vaillage_s_g<-simplify(name_vaillage_s_g, remove.multiple = FALSE, remove.loops = TRUE)
contact_vaillage_s_g<-simplify(contact_vaillage_s_g, remove.multiple = FALSE, remove.loops = TRUE)
distance_vaillage_s_g<-simplify(distance_vaillage_s_g, remove.multiple = FALSE, remove.loops = TRUE)

E(enviro_vaillage_s_g)$weight <- enviro_vaillage_s_edge$weight
E(name_vaillage_s_g)$weight <- name_vaillage_s_edge$weight
E(contact_vaillage_s_g)$weight <- contact_vaillage_s_edge$weight
E(distance_vaillage_s_g)$weight <- distance_vaillage_s_edge$dist

#A_Village
enviro_vaillage_a_g <- graph_from_edgelist(as.matrix(enviro_vaillage_a_edge[,c(1,3)]), directed = FALSE)
name_vaillage_a_g <- graph_from_edgelist(as.matrix(name_vaillage_a_edge[,c(1,3)]), directed = FALSE)
contact_vaillage_a_g <- graph_from_edgelist(as.matrix(contact_vaillage_a_edge[,c(1,3)]), directed = FALSE)
distance_vaillage_a_g <- graph_from_edgelist(as.matrix(distance_vaillage_a_edge[,c(1,2)]), directed = FALSE)

E(enviro_vaillage_a_g)$weight <- enviro_vaillage_a_edge$weight
E(name_vaillage_a_g)$weight <- name_vaillage_a_edge$weight
E(contact_vaillage_a_g)$weight <- contact_vaillage_a_edge$weight
E(distance_vaillage_a_g)$weight <- distance_vaillage_a_edge$dist

enviro_vaillage_a_g<-simplify(enviro_vaillage_a_g, remove.multiple = FALSE, remove.loops = TRUE)
name_vaillage_a_g<-simplify(name_vaillage_a_g, remove.multiple = FALSE, remove.loops = TRUE)
contact_vaillage_a_g<-simplify(contact_vaillage_a_g, remove.multiple = FALSE, remove.loops = TRUE)
distance_vaillage_a_g<-simplify(distance_vaillage_a_g, remove.multiple = FALSE, remove.loops = TRUE)



#CALCULATE CENTRALITY
#Environmental networks
centrality_enviro_vaillage_m <- as.data.frame(page_rank(enviro_vaillage_m_g)$vector) 
centrality_enviro_vaillage_m$Node <- rownames(centrality_enviro_vaillage_m)
colnames(centrality_enviro_vaillage_m)[1] <- "PageRank"

centrality_enviro_vaillage_s <- as.data.frame(page_rank(enviro_vaillage_s_g)$vector) 
centrality_enviro_vaillage_s$Node <- rownames(centrality_enviro_vaillage_s)
colnames(centrality_enviro_vaillage_s)[1] <- "PageRank"

centrality_enviro_vaillage_a <- as.data.frame(page_rank(enviro_vaillage_a_g)$vector) 
centrality_enviro_vaillage_a$Node <- rownames(centrality_enviro_vaillage_a)
colnames(centrality_enviro_vaillage_a)[1] <- "PageRank"

#Close contact networks
centrality_contact_vaillage_m <- as.data.frame(page_rank(contact_vaillage_m_g)$vector) 
centrality_contact_vaillage_m$Node <- rownames(centrality_contact_vaillage_m)
colnames(centrality_contact_vaillage_m)[1] <- "PageRank"

centrality_contact_vaillage_s <- as.data.frame(page_rank(contact_vaillage_s_g)$vector) 
centrality_contact_vaillage_s$Node <- rownames(centrality_contact_vaillage_s)
colnames(centrality_contact_vaillage_s)[1] <- "PageRank"

centrality_contact_vaillage_a <- as.data.frame(page_rank(contact_vaillage_a_g)$vector) 
centrality_contact_vaillage_a$Node <- rownames(centrality_contact_vaillage_a)
colnames(centrality_contact_vaillage_a)[1] <- "PageRank"

#Naming networks
centrality_name_vaillage_m <- as.data.frame(page_rank(name_vaillage_m_g)$vector) 
centrality_name_vaillage_m$Node <- rownames(centrality_name_vaillage_m)
colnames(centrality_name_vaillage_m)[1] <- "PageRank"

centrality_name_vaillage_s <- as.data.frame(page_rank(name_vaillage_s_g)$vector) 
centrality_name_vaillage_s$Node <- rownames(centrality_name_vaillage_s)
colnames(centrality_name_vaillage_s)[1] <- "PageRank"

centrality_name_vaillage_a <- as.data.frame(page_rank(name_vaillage_a_g)$vector) 
centrality_name_vaillage_a$Node <- rownames(centrality_name_vaillage_a)
colnames(centrality_name_vaillage_a)[1] <- "PageRank"

#Distance networks
centrality_distance_vaillage_m <- as.data.frame(page_rank(distance_vaillage_m_g)$vector) 
centrality_distance_vaillage_m$Node <- rownames(centrality_distance_vaillage_m)
colnames(centrality_distance_vaillage_m)[1] <- "PageRank"

centrality_distance_vaillage_s <- as.data.frame(page_rank(distance_vaillage_s_g)$vector) 
centrality_distance_vaillage_s$Node <- rownames(centrality_distance_vaillage_s)
colnames(centrality_distance_vaillage_s)[1] <- "PageRank"

centrality_distance_vaillage_a <- as.data.frame(page_rank(distance_vaillage_a_g)$vector) 
centrality_distance_vaillage_a$Node <- rownames(centrality_distance_vaillage_a)
colnames(centrality_distance_vaillage_a)[1] <- "PageRank"


centrality_name_vaillage_m$Strength <- igraph::strength(name_vaillage_m_g)
centrality_name_vaillage_m$Betweenness <- igraph::betweenness(name_vaillage_m_g, weights = 1/E(name_vaillage_m_g)$weight, normalized = T)
centrality_name_vaillage_m$Eigen <- log(igraph::eigen_centrality(name_vaillage_m_g)$vector)

centrality_name_vaillage_s$Strength <- igraph::strength(name_vaillage_s_g)
centrality_name_vaillage_s$Betweenness <- igraph::betweenness(name_vaillage_s_g, weights =1/ E(name_vaillage_s_g)$weight, normalized = T)
centrality_name_vaillage_s$Eigen <- log(igraph::eigen_centrality(name_vaillage_s_g)$vector)

centrality_name_vaillage_a$Strength <- igraph::strength(name_vaillage_a_g)
centrality_name_vaillage_a$Betweenness <- igraph::betweenness(name_vaillage_a_g, weights = 1/E(name_vaillage_a_g)$weight, normalized = T)
centrality_name_vaillage_a$Eigen <- log(igraph::eigen_centrality(name_vaillage_a_g)$vector)


centrality_contact_vaillage_m$Strength <- igraph::strength(contact_vaillage_m_g)
centrality_contact_vaillage_m$Betweenness <- igraph::betweenness(contact_vaillage_m_g, weights = 1/E(contact_vaillage_m_g)$weight, normalized = T)
centrality_contact_vaillage_m$Eigen <- igraph::eigen_centrality(contact_vaillage_m_g)$vector

centrality_contact_vaillage_s$Strength <- igraph::strength(contact_vaillage_s_g)
centrality_contact_vaillage_s$Betweenness <- igraph::betweenness(contact_vaillage_s_g, weights = 1/E(contact_vaillage_s_g)$weight, normalized = T)
centrality_contact_vaillage_s$Eigen <- igraph::eigen_centrality(contact_vaillage_s_g)$vector

centrality_contact_vaillage_a$Strength <- igraph::strength(contact_vaillage_a_g)
centrality_contact_vaillage_a$Betweenness <- igraph::betweenness(contact_vaillage_a_g, weights = 1/E(contact_vaillage_a_g)$weight, normalized = T)
centrality_contact_vaillage_a$Eigen <- igraph::eigen_centrality(contact_vaillage_a_g)$vector


centrality_enviro_vaillage_m$Strength <- igraph::strength(enviro_vaillage_m_g)
centrality_enviro_vaillage_m$Betweenness <- igraph::betweenness(enviro_vaillage_m_g, weights = 1/E(enviro_vaillage_m_g)$weight, normalized = T)
centrality_enviro_vaillage_m$Eigen <- igraph::eigen_centrality(enviro_vaillage_m_g)$vector

centrality_enviro_vaillage_s$Strength <- igraph::strength(enviro_vaillage_s_g)
centrality_enviro_vaillage_s$Betweenness <- igraph::betweenness(enviro_vaillage_s_g, weights = 1/E(enviro_vaillage_s_g)$weight, normalized = T)
centrality_enviro_vaillage_s$Eigen <- igraph::eigen_centrality(enviro_vaillage_s_g)$vector

centrality_enviro_vaillage_a$Strength <- igraph::strength(enviro_vaillage_a_g)
centrality_enviro_vaillage_a$Betweenness <- igraph::betweenness(enviro_vaillage_a_g, weights = 1/E(enviro_vaillage_a_g)$weight, normalized = T)
centrality_enviro_vaillage_a$Eigen <- igraph::eigen_centrality(enviro_vaillage_a_g)$vector


centrality_distance_vaillage_m$Strength <- igraph::strength(distance_vaillage_m_g)
centrality_distance_vaillage_m$Betweenness <- igraph::betweenness(distance_vaillage_m_g, weights = 1/E(distance_vaillage_m_g)$weight, normalized = T)
centrality_distance_vaillage_m$Eigen <- igraph::eigen_centrality(distance_vaillage_m_g)$vector

centrality_distance_vaillage_s$Strength <- igraph::strength(distance_vaillage_s_g)
centrality_distance_vaillage_s$Betweenness <- igraph::betweenness(distance_vaillage_s_g, weights = 1/E(distance_vaillage_s_g)$weight, normalized = T)
centrality_distance_vaillage_s$Eigen <- igraph::eigen_centrality(distance_vaillage_s_g)$vector

centrality_distance_vaillage_a$Strength <- igraph::strength(distance_vaillage_a_g)
centrality_distance_vaillage_a$Betweenness <- igraph::betweenness(distance_vaillage_a_g, weights = 1/E(distance_vaillage_a_g)$weight, normalized = T)
centrality_distance_vaillage_a$Eigen <- igraph::eigen_centrality(distance_vaillage_a_g)$vector


```

```{r}

centrality_contact_vaillage_m$Type <- "Contact"
centrality_name_vaillage_m$Type <- "Name"
centrality_enviro_vaillage_m$Type <- "Enviro"
centrality_distance_vaillage_m$Type <- "House"

contact_df_vaillage_m <- rbind(centrality_contact_vaillage_m, centrality_name_vaillage_m,centrality_enviro_vaillage_m,centrality_distance_vaillage_m)
contact_df_vaillage_m$Village <- "M"

centrality_contact_vaillage_a$Type <- "Contact"
centrality_name_vaillage_a$Type <- "Name"
centrality_enviro_vaillage_a$Type <- "Enviro"
centrality_distance_vaillage_a$Type <- "House"

contact_df_vaillage_a <- rbind(centrality_contact_vaillage_a, centrality_name_vaillage_a,centrality_enviro_vaillage_a,centrality_distance_vaillage_a)
contact_df_vaillage_a$Village <- "A"

centrality_contact_vaillage_s$Type <- "Contact"
centrality_name_vaillage_s$Type <- "Name"
centrality_enviro_vaillage_s$Type <- "Enviro"
centrality_distance_vaillage_s$Type <- "House"

contact_df_vaillage_s <- rbind(centrality_contact_vaillage_s, centrality_name_vaillage_s,centrality_enviro_vaillage_s,centrality_distance_vaillage_s)
contact_df_vaillage_s$Village <- "S"

contact_df <- rbind(contact_df_vaillage_m, contact_df_vaillage_a, contact_df_vaillage_s)

colnames(contact_df) <- c( "PR", "Node","Stg", "Btwn", "Egn", "Type", "Village")

contact_df2 <-contact_df %>% pivot_wider(names_from = Type, values_from = c(PR, Stg, Btwn, Egn))

 
centrality_vars <- contact_df2[ , -c(1,2)]

write.csv(contact_df2, "contact_df2.csv")


```

```{r Explore Centrality Correlations}
library(corrplot)

corrplot(cor((centrality_name_vaillage_m[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_name_vaillage_s[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_name_vaillage_a[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)

corrplot(cor((centrality_enviro_vaillage_m[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_enviro_vaillage_s[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_enviro_vaillage_a[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)

corrplot(cor((centrality_contact_vaillage_m[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_contact_vaillage_s[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_contact_vaillage_a[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)

corrplot(cor((centrality_distance_vaillage_m[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_distance_vaillage_s[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)
corrplot(cor((centrality_distance_vaillage_a[,c(1,3:5)])),type = "lower", tl.col = "black",tl.srt = 0)

```

```{r Explore Centrality Correlations2}
cor_name_vaillage_m <- as.data.frame(cor((centrality_name_vaillage_m[,c(1,3:5)])))
cor_name_vaillage_s <-as.data.frame(cor((centrality_name_vaillage_s[,c(1,3:5)])))
cor_name_vaillage_a <-as.data.frame(cor((centrality_name_vaillage_a[,c(1,3:5)])))

cor_enviro_vaillage_m <-as.data.frame(cor((centrality_enviro_vaillage_m[,c(1,3:5)])))
cor_enviro_vaillage_s <-as.data.frame(cor((centrality_enviro_vaillage_s[,c(1,3:5)])))
cor_enviro_vaillage_a <-as.data.frame(cor((centrality_enviro_vaillage_a[,c(1,3:5)])))

cor_contact_vaillage_m <-as.data.frame(cor((centrality_contact_vaillage_m[,c(1,3:5)])))
cor_contact_vaillage_s <-as.data.frame(cor((centrality_contact_vaillage_s[,c(1,3:5)])))
cor_contact_vaillage_a <-as.data.frame(cor((centrality_contact_vaillage_a[,c(1,3:5)])))

cor_distance_vaillage_m <-as.data.frame(cor((centrality_distance_vaillage_m[,c(1,3:5)])))
cor_distance_vaillage_s <-as.data.frame(cor((centrality_distance_vaillage_s[,c(1,3:5)])))
cor_distance_vaillage_a <-as.data.frame(cor((centrality_distance_vaillage_a[,c(1,3:5)])))


cor_name_vaillage_m$Village <- "M_Village"
cor_name_vaillage_s$Village <- "S_Village"
cor_name_vaillage_a$Village <- "A_Village"

cor_enviro_vaillage_m$Village <- "M_Village"
cor_enviro_vaillage_s$Village <- "S_Village"
cor_enviro_vaillage_a$Village <- "A_Village"

cor_contact_vaillage_m$Village <- "M_Village"
cor_contact_vaillage_s$Village <- "S_Village"
cor_contact_vaillage_a$Village <- "A_Village"

cor_distance_vaillage_m$Village <- "M_Village"
cor_distance_vaillage_s$Village <- "S_Village"
cor_distance_vaillage_a$Village <- "A_Village"

cor_name_vaillage_m[2,1] <- NA
cor_name_vaillage_m[1,3] <- NA
cor_name_vaillage_m[1,4] <- NA
cor_name_vaillage_m[2,3] <- NA
cor_name_vaillage_m[2,4] <- NA
cor_name_vaillage_m[4,3] <- NA

cor_name_vaillage_s[2,1] <- NA
cor_name_vaillage_s[1,3] <- NA
cor_name_vaillage_s[1,4] <- NA
cor_name_vaillage_s[2,3] <- NA
cor_name_vaillage_s[2,4] <- NA
cor_name_vaillage_s[4,3] <- NA

cor_name_vaillage_a[2,1] <- NA
cor_name_vaillage_a[1,3] <- NA
cor_name_vaillage_a[1,4] <- NA
cor_name_vaillage_a[2,3] <- NA
cor_name_vaillage_a[2,4] <- NA
cor_name_vaillage_a[4,3] <- NA

cor_name <- rbind(cor_name_vaillage_m, cor_name_vaillage_a, cor_name_vaillage_s)
cor_name$name2 <- rep(c("PageRank", "Strength", "Betweenness", "Eigenvector"),3)
cor_name2 <- cor_name %>% pivot_longer(cols = PageRank:Eigen)
cor_name2$name[cor_name2$name=="Eigen"] <- "Eigenvector"

cor_name_vaillage_m_plot <- ggplot(data = cor_name2[cor_name2$Village=="M_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    name = "Correlation",
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) 

cor_name_vaillage_a_plot <- ggplot(data = cor_name2[cor_name2$Village=="A_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    name = "Correlation",
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) 

cor_name_vaillage_s_plot <- ggplot(data = cor_name2[cor_name2$Village=="S_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2,)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    name = "Correlation",
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) 
  

cor_name_plot <- ggarrange(cor_name_vaillage_m_plot, cor_name_vaillage_s_plot, cor_name_vaillage_a_plot, common.legend = TRUE, ncol = 3)




cor_enviro_vaillage_m[2,1] <- NA
cor_enviro_vaillage_m[1,3] <- NA
cor_enviro_vaillage_m[1,4] <- NA
cor_enviro_vaillage_m[2,3] <- NA
cor_enviro_vaillage_m[2,4] <- NA
cor_enviro_vaillage_m[4,3] <- NA

cor_enviro_vaillage_s[2,1] <- NA
cor_enviro_vaillage_s[1,3] <- NA
cor_enviro_vaillage_s[1,4] <- NA
cor_enviro_vaillage_s[2,3] <- NA
cor_enviro_vaillage_s[2,4] <- NA
cor_enviro_vaillage_s[4,3] <- NA

cor_enviro_vaillage_a[2,1] <- NA
cor_enviro_vaillage_a[1,3] <- NA
cor_enviro_vaillage_a[1,4] <- NA
cor_enviro_vaillage_a[2,3] <- NA
cor_enviro_vaillage_a[2,4] <- NA
cor_enviro_vaillage_a[4,3] <- NA

cor_enviro <- rbind(cor_enviro_vaillage_m, cor_enviro_vaillage_a, cor_enviro_vaillage_s)
cor_enviro$name2 <- rep(c("PageRank", "Strength", "Betweenness", "Eigenvector"),3)
cor_enviro2 <- cor_enviro %>% pivot_longer(cols = PageRank:Eigen)
cor_enviro2$name[cor_enviro2$name=="Eigen"] <- "Eigenvector"

cor_enviro_vaillage_m_plot <- ggplot(data = cor_enviro2[cor_enviro2$Village=="M_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )

cor_enviro_vaillage_a_plot <- ggplot(data = cor_enviro2[cor_enviro2$Village=="A_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )

cor_enviro_vaillage_s_plot <- ggplot(data = cor_enviro2[cor_enviro2$Village=="S_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2,)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )
  

cor_enviro_plot <-ggarrange(cor_enviro_vaillage_m_plot, cor_enviro_vaillage_s_plot, cor_enviro_vaillage_a_plot,ncol = 3)







cor_distance_vaillage_m[2,1] <- NA
cor_distance_vaillage_m[1,3] <- NA
cor_distance_vaillage_m[1,4] <- NA
cor_distance_vaillage_m[2,3] <- NA
cor_distance_vaillage_m[2,4] <- NA
cor_distance_vaillage_m[4,3] <- NA

cor_distance_vaillage_s[2,1] <- NA
cor_distance_vaillage_s[1,3] <- NA
cor_distance_vaillage_s[1,4] <- NA
cor_distance_vaillage_s[2,3] <- NA
cor_distance_vaillage_s[2,4] <- NA
cor_distance_vaillage_s[4,3] <- NA

cor_distance_vaillage_a[2,1] <- NA
cor_distance_vaillage_a[1,3] <- NA
cor_distance_vaillage_a[1,4] <- NA
cor_distance_vaillage_a[2,3] <- NA
cor_distance_vaillage_a[2,4] <- NA
cor_distance_vaillage_a[4,3] <- NA

cor_distance <- rbind(cor_distance_vaillage_m, cor_distance_vaillage_a, cor_distance_vaillage_s)
cor_distance$name2 <- rep(c("PageRank", "Strength", "Betweenness", "Eigenvector"),3)
cor_distance2 <- cor_distance %>% pivot_longer(cols = PageRank:Eigen)
cor_distance2$name[cor_distance2$name=="Eigen"] <- "Eigenvector"

cor_distance_vaillage_m_plot <- ggplot(data = cor_distance2[cor_distance2$Village=="M_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )

cor_distance_vaillage_a_plot <- ggplot(data = cor_distance2[cor_distance2$Village=="A_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )

cor_distance_vaillage_s_plot <- ggplot(data = cor_distance2[cor_distance2$Village=="S_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2,)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )
  

cor_distance_plot <-ggarrange(cor_distance_vaillage_m_plot, cor_distance_vaillage_s_plot, cor_distance_vaillage_a_plot,ncol = 3)



cor_contact_vaillage_m[2,1] <- NA
cor_contact_vaillage_m[1,3] <- NA
cor_contact_vaillage_m[1,4] <- NA
cor_contact_vaillage_m[2,3] <- NA
cor_contact_vaillage_m[2,4] <- NA
cor_contact_vaillage_m[4,3] <- NA

cor_contact_vaillage_s[2,1] <- NA
cor_contact_vaillage_s[1,3] <- NA
cor_contact_vaillage_s[1,4] <- NA
cor_contact_vaillage_s[2,3] <- NA
cor_contact_vaillage_s[2,4] <- NA
cor_contact_vaillage_s[4,3] <- NA

cor_contact_vaillage_a[2,1] <- NA
cor_contact_vaillage_a[1,3] <- NA
cor_contact_vaillage_a[1,4] <- NA
cor_contact_vaillage_a[2,3] <- NA
cor_contact_vaillage_a[2,4] <- NA
cor_contact_vaillage_a[4,3] <- NA

cor_contact <- rbind(cor_contact_vaillage_m, cor_contact_vaillage_a, cor_contact_vaillage_s)
cor_contact$name2 <- rep(c("PageRank", "Strength", "Betweenness", "Eigenvector"),3)
cor_contact2 <- cor_contact %>% pivot_longer(cols = PageRank:Eigen)
cor_contact2$name[cor_contact2$name=="Eigen"] <- "Eigenvector"

cor_contact_vaillage_m_plot <- ggplot(data = cor_contact2[cor_contact2$Village=="M_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )

cor_contact_vaillage_a_plot <- ggplot(data = cor_contact2[cor_contact2$Village=="A_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )

cor_contact_vaillage_s_plot <- ggplot(data = cor_contact2[cor_contact2$Village=="S_Village",], aes(x = name, y=name2,  fill=value))+
  geom_tile(color="white", size=2,)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )
  

cor_contact_plot <-ggarrange(cor_contact_vaillage_m_plot, cor_contact_vaillage_s_plot, cor_contact_vaillage_a_plot,ncol = 3)

cor_plot_total <-ggarrange(cor_name_plot, cor_contact_plot, cor_distance_plot,cor_enviro_plot, nrow = 4)

cor_contact3 <- cor_contact2 %>% drop_na(value) %>% group_by(name2, name) %>% summarise(value= mean(value))

cor_contact3$color_group <- cor_contact3$value
cor_contact3$color_group[cor_contact3$color_group==1] <- "one"
cor_contact3$color_group[cor_contact3$color_group!="one"] <- "other"

cor_contact_plot2 <- ggplot(data = cor_contact3, aes(x = name, y=name2,  fill=value, alpha = color_group))+
  geom_tile(color="white", size=2,)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )+
  ggtitle("Close Contact")
  


cor_name_plot <-ggarrange(cor_name_vaillage_m_plot, cor_name_vaillage_s_plot, cor_name_vaillage_a_plot,ncol = 3)

cor_plot_total <-ggarrange(cor_name_plot, cor_name_plot, cor_distance_plot,cor_enviro_plot, nrow = 4)

cor_name3 <- cor_name2 %>% drop_na(value) %>% group_by(name2, name) %>% summarise(value= mean(value))

cor_name3$color_group <- cor_name3$value
cor_name3$color_group[cor_name3$color_group==1] <- "one"
cor_name3$color_group[cor_name3$color_group!="one"] <- "other"

cor_name_plot2 <- ggplot(data = cor_name3, aes(x = name, y=name2,  fill=value, alpha = color_group))+
  geom_tile(color="white", size=2,)+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 2)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )+
  ggtitle("Social")


cor_enviro_plot <-ggarrange(cor_enviro_vaillage_m_plot, cor_enviro_vaillage_s_plot, cor_enviro_vaillage_a_plot,ncol = 3)

cor_plot_total <-ggarrange(cor_name_plot, cor_enviro_plot, cor_distance_plot,cor_enviro_plot, nrow = 4)

cor_enviro3 <- cor_enviro2 %>% drop_na(value) %>% group_by(name2, name) %>% summarise(value= mean(value))

cor_enviro3$color_group <- cor_enviro3$value
cor_enviro3$color_group[cor_enviro3$color_group==1] <- "one"
cor_enviro3$color_group[cor_enviro3$color_group!="one"] <- "other"

cor_enviro_plot2 <- ggplot(data = cor_enviro3, aes(x = name, y=name2,  fill=value, alpha= color_group))+
  geom_tile(size=2, color="white")+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 3)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )+
  ggtitle("Environmental")

cor_distance_plot <-ggarrange(cor_distance_vaillage_m_plot, cor_distance_vaillage_s_plot, cor_distance_vaillage_a_plot,ncol = 3)

cor_plot_total <-ggarrange(cor_name_plot, cor_distance_plot, cor_distance_plot,cor_enviro_plot, nrow = 4)

cor_distance3 <- cor_distance2 %>% drop_na(value) %>% group_by(name2, name) %>% summarise(value= mean(value))

cor_distance3$Correlation <- cor_distance3$value

cor_distance3$color_group <- cor_distance3$value
cor_distance3$color_group[cor_distance3$color_group==1] <- "one"
cor_distance3$color_group[cor_distance3$color_group!="one"] <- "other"

cor_distance_plot2 <- ggplot(data = cor_distance3, aes(x = name, y=name2,  fill=value, alpha= color_group))+
  geom_tile(size=2, color="white")+
  theme_classic()+
  xlab("Centrality")+
  ylab("Centrality")+
  scale_fill_viridis_c(
    direction = -1,
    na.value = "white",
    limits = c(-0.06, 1) 
  )+ geom_text(
    aes(label = round(value, 3)), 
    color = "white", 
    na.rm = TRUE, 
    size = 5 
  ) +
  theme(
    legend.position = "none" # Remove the legend
  )+
  ggtitle("Household")


correlation_plot_total <- ggarrange(cor_name_plot2, cor_contact_plot2, cor_distance_plot2, cor_enviro_plot2, labels = c("a", "b", "c", "d"))

```

```{r Explore Centrality Correlations- Scatterplots}

library(GGally)

centrality_name_vaillage_m$Village <- "M"
centrality_name_vaillage_a$Village <- "A"
centrality_name_vaillage_s$Village <- "S"

centrality_name <- rbind(centrality_name_vaillage_m, centrality_name_vaillage_a, centrality_name_vaillage_s)
colnames(centrality_name)[5] <- "Eigenvector"

scatmat_name <- ggscatmat(centrality_name, columns = c(1,3,4,5), color = "Village", alpha=0.05)+
  theme_classic()+
  ylab("Centrality Value")+
  xlab("Centrality Value")+
  guides(
    fill = "none", # Ensure fill legend is removed
    color = guide_legend(
      title = "Village",
      override.aes = list(
        alpha = 1, 
        size = 4  
      )
    )
  ) +
  scale_color_manual(
    values = c(A = "darkorchid", M = "darkgreen", S = "darkblue"),
    breaks = c("A", "M", "S"), 
    labels = c("A", "M", "S") 
  )+
  theme(legend.position = "top")
  

centrality_enviro_vaillage_m$Village <- "M"
centrality_enviro_vaillage_a$Village <- "A"
centrality_enviro_vaillage_s$Village <- "S"

centrality_enviro <- rbind(centrality_enviro_vaillage_m, centrality_enviro_vaillage_a, centrality_enviro_vaillage_s)
colnames(centrality_enviro)[5] <- "Eigenvector"

scatmat_enviro <- ggscatmat(centrality_enviro, columns = c(1,3,4,5), color = "Village", alpha=0.05)+
  theme_classic()+
  ylab("Centrality Value")+
  xlab("Centrality Value")+
  guides(
    fill = "none", # Ensure fill legend is removed
    color = guide_legend(
      title = "Village",
      override.aes = list(
        alpha = 1, 
        size = 4  
      )
    )
  ) +
  scale_color_manual(
    values = c(A = "darkorchid", M = "darkgreen", S = "darkblue"),
    breaks = c("A", "M", "S"), 
    labels = c("A", "M", "S") 
  )+
  theme(legend.position = "top")




centrality_contact_vaillage_m$Village <- "M"
centrality_contact_vaillage_a$Village <- "A"
centrality_contact_vaillage_s$Village <- "S"

centrality_contact <- rbind(centrality_contact_vaillage_m, centrality_contact_vaillage_a, centrality_contact_vaillage_s)
colnames(centrality_contact)[5] <- "Eigenvector"

scatmat_contact <- ggscatmat(centrality_contact, columns = c(1,3,4,5), color = "Village", alpha=0.05)+
  theme_classic()+
  ylab("Centrality Value")+
  xlab("Centrality Value")+
  guides(
    fill = "none", # Ensure fill legend is removed
    color = guide_legend(
      title = "Village",
      override.aes = list(
        alpha = 1, 
        size = 4  
      )
    )
  ) +
  scale_color_manual(
    values = c(A = "darkorchid", M = "darkgreen", S = "darkblue"),
    breaks = c("A", "M", "S"), 
    labels = c("A", "M", "S") 
  )+
  theme(legend.position = "top")




centrality_distance_vaillage_m$Village <- "M"
centrality_distance_vaillage_a$Village <- "A"
centrality_distance_vaillage_s$Village <- "S"

centrality_distance <- rbind(centrality_distance_vaillage_m, centrality_distance_vaillage_a, centrality_distance_vaillage_s)
colnames(centrality_distance)[5] <- "Eigenvector"

scatmat_distance <- ggscatmat(centrality_distance, columns = c(1,3,4,5), color = "Village", alpha=0.05)+
  theme_classic()+
  ylab("Centrality Value")+
  xlab("Centrality Value")+
  guides(
    fill = "none", # Ensure fill legend is removed
    color = guide_legend(
      title = "Village",
      override.aes = list(
        alpha = 1, 
        size = 4  
      )
    )
  ) +
  scale_color_manual(
    values = c(A = "darkorchid", M = "darkgreen", S = "darkblue"),
    breaks = c("A", "M", "S"), 
    labels = c("A", "M", "S") 
  )+
  theme(legend.position = "top")


```

```{r Explore Centrality Mean Values}
mean(centrality_name_vaillage_m$Strength)
mean(centrality_name_vaillage_a$Strength)
mean(centrality_name_vaillage_s$Strength)

mean(centrality_name_vaillage_m$PageRank)
mean(centrality_name_vaillage_a$PageRank)
mean(centrality_name_vaillage_s$PageRank)

mean(centrality_name_vaillage_m$Betweenness)
mean(centrality_name_vaillage_a$Betweenness)
mean(centrality_name_vaillage_s$Betweenness)

mean(centrality_name_vaillage_m$Eigen)
mean(centrality_name_vaillage_a$Eigen)
mean(centrality_name_vaillage_s$Eigen)

mean(centrality_enviro_vaillage_m$Strength)
mean(centrality_enviro_vaillage_a$Strength)
mean(centrality_enviro_vaillage_s$Strength)

mean(centrality_enviro_vaillage_m$PageRank)
mean(centrality_enviro_vaillage_a$PageRank)
mean(centrality_enviro_vaillage_s$PageRank)

mean(centrality_enviro_vaillage_m$Betweenness)
mean(centrality_enviro_vaillage_a$Betweenness)
mean(centrality_enviro_vaillage_s$Betweenness)

mean(centrality_enviro_vaillage_m$Eigen)
mean(centrality_enviro_vaillage_a$Eigen)
mean(centrality_enviro_vaillage_s$Eigen)

mean(centrality_contact_vaillage_m$Strength)
mean(centrality_contact_vaillage_a$Strength)
mean(centrality_contact_vaillage_s$Strength)

mean(centrality_contact_vaillage_m$PageRank)
mean(centrality_contact_vaillage_a$PageRank)
mean(centrality_contact_vaillage_s$PageRank)

mean(centrality_contact_vaillage_m$Betweenness)
mean(centrality_contact_vaillage_a$Betweenness)
mean(centrality_contact_vaillage_s$Betweenness)

mean(centrality_contact_vaillage_m$Eigen)
mean(centrality_contact_vaillage_a$Eigen)
mean(centrality_contact_vaillage_s$Eigen)

mean(centrality_distance_vaillage_m$Strength)
mean(centrality_distance_vaillage_a$Strength)
mean(centrality_distance_vaillage_s$Strength)

mean(centrality_distance_vaillage_m$PageRank)
mean(centrality_distance_vaillage_a$PageRank)
mean(centrality_distance_vaillage_s$PageRank)

mean(centrality_distance_vaillage_m$Betweenness)
mean(centrality_distance_vaillage_a$Betweenness)
mean(centrality_distance_vaillage_s$Betweenness)

mean(centrality_distance_vaillage_m$Eigen)
mean(centrality_distance_vaillage_a$Eigen)
mean(centrality_distance_vaillage_s$Eigen)

name_strength_density <- ggplot(data = corr_df_name, aes(x = Strength, color=Village))+
  geom_density()+
  ylab("Density")+
  theme_classic()

name_pagerank_density <- ggplot(data = corr_df_name, aes(x = PageRank, color=Village))+
  geom_density()+
  ylab("Density")+
  theme_classic()

name_eigen_density <- ggplot(data = corr_df_name, aes(x = Eigen, color=Village))+
  geom_density()+
  ylab("Density")+
  theme_classic()

name_betweenness_density <- ggplot(data = corr_df_name, aes(x = Betweenness, color=Village))+
  geom_density()+
  ylab("Density")+
  theme_classic()


name_betweenness_box <- ggplot(data = corr_df_name, aes(y = Betweenness, x=Village, color=Village))+
  geom_boxplot()+
  theme_classic()

enviro_betweenness_box <- ggplot(data = corr_df_enviro, aes(y = Betweenness, x=Village, color=Village))+
  geom_boxplot()+
  theme_classic()

contact_betweenness_box <- ggplot(data = corr_df_contact, aes(y = Betweenness, x=Village, color=Village))+
  geom_boxplot()+
  theme_classic()

distance_betweenness_box <- ggplot(data = corr_df_distance, aes(y = Betweenness, x=Village, color=Village))+
  geom_boxplot()+
  theme_classic()

```

```{r Scatterplots}

name_strength_pagerank <- ggplot(data = centrality_name_vaillage_m, aes(x = Strength, y= PageRank))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

name_strength_eigen<- ggplot(data = centrality_name_vaillage_m, aes(x = Strength, y= log(Eigen)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

name_strength_betweenness<- ggplot(data = centrality_name_vaillage_m, aes(x = Strength, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

name_eigen_pagerank<- ggplot(data = centrality_name_vaillage_m, aes(x = PageRank, y= log(Eigen)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

name_betweenness_pagerank<- ggplot(data = centrality_name_vaillage_m, aes(x = PageRank, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

name_eigen_betweenness <- ggplot(data = centrality_name_vaillage_m, aes(x = log(Betweenness), y= log(Eigen)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

ggarrange(name_strength_pagerank, name_strength_eigen, name_strength_betweenness, name_eigen_pagerank, name_eigen_pagerank, name_betweenness_pagerank, name_eigen_betweenness)




enviro_strength_pagerank <- ggplot(data = centrality_enviro_vaillage_m, aes(x = Strength, y= PageRank))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

enviro_strength_eigen<- ggplot(data =centrality_enviro_vaillage_m, aes(x = Strength, y= Eigen))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

enviro_strength_betweenness<- ggplot(data =centrality_enviro_vaillage_m, aes(x = Strength, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

enviro_eigen_pagerank<- ggplot(data = centrality_enviro_vaillage_m, aes(x = PageRank, y= Eigen))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

enviro_betweenness_pagerank<- ggplot(data = centrality_enviro_vaillage_m, aes(x = PageRank, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

enviro_eigen_betweenness <- ggplot(data = centrality_enviro_vaillage_m, aes(x = log(Betweenness), y= Eigen))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

ggarrange(enviro_strength_pagerank, enviro_strength_eigen, enviro_strength_betweenness, enviro_eigen_pagerank, enviro_eigen_pagerank, enviro_betweenness_pagerank, enviro_eigen_betweenness)




contact_strength_pagerank <- ggplot(data = centrality_contact_vaillage_m, aes(x = Strength, y= PageRank))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

contact_strength_eigen<- ggplot(data = centrality_contact_vaillage_m, aes(x = Strength, y= Eigen))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

contact_strength_betweenness<- ggplot(data = centrality_contact_vaillage_m, aes(x = Strength, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

contact_eigen_pagerank<- ggplot(data = centrality_contact_vaillage_m, aes(x = PageRank, y= Eigen))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

contact_betweenness_pagerank<- ggplot(data = centrality_contact_vaillage_m, aes(x = PageRank, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

contact_eigen_betweenness <- ggplot(data = centrality_contact_vaillage_m, aes(y = log(Betweenness), x= Eigen))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

ggarrange(contact_strength_pagerank, contact_strength_eigen, contact_strength_betweenness, contact_eigen_pagerank, contact_eigen_pagerank, contact_betweenness_pagerank, contact_eigen_betweenness)



distance_strength_pagerank <- ggplot(data = centrality_distance_vaillage_m, aes(x = Strength, y= PageRank))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

distance_strength_eigen<- ggplot(data = centrality_distance_vaillage_m, aes(x = Strength, y= log(Eigen)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

distance_strength_betweenness<- ggplot(data = centrality_distance_vaillage_m, aes(x = Strength, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

distance_eigen_pagerank<- ggplot(data = centrality_distance_vaillage_m, aes(x = PageRank, y= log(Eigen)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

distance_betweenness_pagerank<- ggplot(data = centrality_distance_vaillage_m, aes(x = PageRank, y= log(Betweenness)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

distance_eigen_betweenness <- ggplot(data = centrality_distance_vaillage_m, aes(y = log(Betweenness), x= log(Eigen)))+
  geom_point(alpha=0.1)+
  geom_smooth(method="lm")+
  theme_classic()

ggarrange(distance_strength_pagerank, distance_strength_eigen, distance_strength_betweenness, distance_eigen_pagerank, distance_eigen_pagerank, distance_betweenness_pagerank, distance_eigen_betweenness)


```
 
```{r Generate market integration indices}

# wall construction
metadata$material_wall_index <- 0
metadata$material_wall_index[metadata$material_wall =="bamboo"] <- 0
metadata$material_wall_index[metadata$material_wall =="rafia"] <- 0
metadata$material_wall_index[metadata$material_wall =="ravenala"] <- 0
metadata$material_wall_index[metadata$material_wall =="mud"] <- 0
metadata$material_wall_index[metadata$material_wall =="compacted_earth"] <- 0
metadata$material_wall_index[metadata$material_wall =="wood_planks"] <- 1
metadata$material_wall_index[metadata$material_wall =="brick_unfired"] <- 2
metadata$material_wall_index[metadata$material_wall =="brick_fired"] <- 2
metadata$material_wall_index[metadata$material_wall =="metal_sheets"] <- 3
metadata$material_wall_index[metadata$material_wall =="cement"] <- 4

# roof construction
metadata$material_roof_index <- 0
metadata$material_roof_index[metadata$material_roof=="bamboo"] <- 0
metadata$material_roof_index[metadata$material_roof=="thatch"] <- 0
metadata$material_roof_index[metadata$material_roof=="metal_sheets"] <- 1
metadata$material_roof_index[metadata$material_roof=="cement"] <- 2

# floor construction
metadata$material_floor_index <- 0
metadata$material_floor_index[metadata$material_floor=="dirt"] <- 0
metadata$material_floor_index[metadata$material_floor=="bamboo"] <- 0
metadata$material_floor_index[metadata$material_floor=="rafia"] <- 0
metadata$material_floor_index[metadata$material_floor=="ravinala"] <- 0
metadata$material_floor_index[metadata$material_floor=="wood_planks"] <- 1
metadata$material_floor_index[metadata$material_floor=="cement"] <- 02

# household SOL index
metadata$house_sol <- metadata$material_wall_index + metadata$material_roof_index + metadata$material_floor_index

#################################
## Durable Commercial Goods Index
#################################

# change "Yes" responses to 1s
metadata$own_cellphone[metadata$own_cellphone=="Yes"] <- 1
metadata$own_cellphone[metadata$own_cellphone!=1] <- 0
metadata$own_cellphone <- as.numeric(metadata$own_cellphone)

metadata$own_tv[metadata$own_tv=="Yes"] <- 1
metadata$own_tv[metadata$own_tv!=1] <- 0
metadata$own_tv <- as.numeric(metadata$own_tv)

metadata$own_bicycle[metadata$own_bicycle=="Yes"] <- 1
metadata$own_bicycle[metadata$own_bicycle!=1] <- 0
metadata$own_bicycle <- as.numeric(metadata$own_bicycle)

metadata$own_refrigerator[metadata$own_refrigerator=="Yes"] <- 1
metadata$own_refrigerator[metadata$own_refrigerator!=1] <- 0
metadata$own_refrigerator <- as.numeric(metadata$own_refrigerator)

metadata$own_motorcycle[metadata$own_motorcycle=="Yes"] <- 1
metadata$own_motorcycle[metadata$own_motorcycle!=1] <- 0
metadata$own_motorcycle <- as.numeric(metadata$own_motorcycle)

metadata$own_computer[metadata$own_computer=="Yes"] <- 1
metadata$own_computer[metadata$own_computer!=1] <- 0
metadata$own_computer <- as.numeric(metadata$own_computer)

metadata$own_generator[metadata$own_generator=="Yes"] <- 1
metadata$own_generator[metadata$own_generator!=1] <- 0
metadata$own_generator <- as.numeric(metadata$own_generator)
  
metadata$commercial_goods <- metadata$own_cellphone + metadata$own_tv + metadata$own_bicycle + metadata$own_refrigerator + metadata$own_motorcycle + metadata$own_computer + metadata$own_generator


write.csv(metadata,"~/Documents/GitHub/human_animal_networks/gps_networks/metadata2.csv")

```

```{r Centrality and Demography}
#metadata$age <- scale(metadata$age)
#metadata$house_sol <- scale(metadata$house_sol)
#metadata$commercial_goods <- scale(metadata$commercial_goods)
#metadata$household_size <- scale(metadata$household_size)

metadata$village[metadata$village=="Ampandrana"] <- "A_Village"

metadata$bmi <- (metadata$weight)/(metadata$height^2)
#metadata$bmi[metadata$bmi<18.5] <- 1
#metadata$bmi[metadata$bmi>1 & metadata$bmi<25] <- 2
#metadata$bmi[metadata$bmi>2 & metadata$bmi<30] <- 3
#metadata$bmi[metadata$bmi>30] <- 4

metadata$bmi <- scale(metadata$bmi)

metadata$meanDiastolic <- (metadata$diastolic1 +metadata$diastolic2 + metadata$diastolic3)/3
metadata$meanSystolic <- (metadata$systolic1 +metadata$systolic2 + metadata$systolic3)/3

metadata$bp <- metadata$meanSystolic / metadata$meanDiastolic
#metadata$bp[metadata$meanSystolic<90 | metadata$meanDiastolic <60] <- 0
metadata$bp[metadata$meanSystolic<120 & metadata$meanDiastolic <80 ] <- 1
metadata$bp[metadata$meanSystolic<=139 & metadata$meanDiastolic <=89 & metadata$bp!=1] <- 2
metadata$bp[metadata$meanSystolic>139 | metadata$meanDiastolic>89] <- 3
# Look into this https://www.cdc.gov/bloodpressure/about.htm#:~:text=Blood%20pressure%20is%20measured%20using,your%20heart%20rests%20between%20beats.
# low is just below 120
# make it ordinail 

#metadata$bp <- scale(metadata$bp)

# make education ordinal and add that 

# Look at networks in the different season 

#metadata$landsize_in_daba <- scale(metadata$landsize_in_daba)

metadata$school_level[metadata$school_level=="None"] <- 1
metadata$school_level[metadata$school_level=="Primary"] <- 2
metadata$school_level[metadata$school_level=="Secondary"] <- 3
metadata$school_level[metadata$school_level=="Higher"] <- 4

metadata$school_level <- as.numeric(metadata$school_level)

#metadata$school_level <- scale(metadata$school_level)

metadata <-metadata %>% mutate(livestock = farm_animals_cows +farm_animals_goats+farm_animals_pigs+farm_animals_sheep)
metadata$livestock <- scale(metadata$livestock)

cor.test(metadata$house_sol, metadata$commercial_goods)

cor.test(metadata$age, metadata$bp)

```

```{r Single-Layer Centrality Models Data Preparation}

metadata$Node <- metadata$social_netid

centrality_name_vaillage_m2 <- left_join(centrality_name_vaillage_m, metadata, by ="Node")
centrality_name_vaillage_s2 <- left_join(centrality_name_vaillage_s, metadata, by ="Node")
centrality_name_vaillage_a2 <- left_join(centrality_name_vaillage_a, metadata, by ="Node")

centrality_contact_vaillage_m2 <- left_join(centrality_contact_vaillage_m, metadata, by ="Node")
centrality_contact_vaillage_s2 <- left_join(centrality_contact_vaillage_s, metadata, by ="Node")
centrality_contact_vaillage_a2 <- left_join(centrality_contact_vaillage_a, metadata, by ="Node")

centrality_enviro_vaillage_m2 <- left_join(centrality_enviro_vaillage_m, metadata, by ="Node")
centrality_enviro_vaillage_s2 <- left_join(centrality_enviro_vaillage_s, metadata, by ="Node")
centrality_enviro_vaillage_a2 <- left_join(centrality_enviro_vaillage_a, metadata, by ="Node")

centrality_distance_vaillage_m2 <- left_join(centrality_distance_vaillage_m, metadata, by ="Node")
centrality_distance_vaillage_s2 <- left_join(centrality_distance_vaillage_s, metadata, by ="Node")
centrality_distance_vaillage_a2 <- left_join(centrality_distance_vaillage_a, metadata, by ="Node")


centrality_name_total2 <- rbind(centrality_name_vaillage_m2, centrality_name_vaillage_s2,centrality_name_vaillage_a2) %>% select(Node, PageRank, gender,  age, house_sol, commercial_goods, village, school_level, bmi, landsize_in_daba,  livestock, season, household_size, Strength, Betweenness, Eigen) %>% drop_na() #getting rid of all NA values removed 77 individuals (~6% of all nodes)

centrality_name_total2$season <- as.factor(centrality_name_total2$season)

centrality_name_vaillage_m_df <- centrality_name_total2[centrality_name_total2$village=="M_Village",]
centrality_name_vaillage_s_df <- centrality_name_total2[centrality_name_total2$village=="S_Village",]
centrality_name_vaillage_a_df <- centrality_name_total2[centrality_name_total2$village=="A_Village",]

centrality_contact_total2 <- rbind(centrality_contact_vaillage_m2, centrality_contact_vaillage_s2,centrality_contact_vaillage_a2) %>% select(Node, PageRank, gender,  age, house_sol, commercial_goods, village, school_level, bmi, landsize_in_daba, livestock, season, household_size, Strength, Betweenness, Eigen) %>%drop_na()#getting rid of all NA values removed 36 individuals (~4% of all nodes)

centrality_contact_total2$season <- as.factor(centrality_contact_total2$season)


centrality_contact_vaillage_m_df <- centrality_contact_total2[centrality_contact_total2$village=="M_Village",]
centrality_contact_vaillage_s_df <- centrality_contact_total2[centrality_contact_total2$village=="S_Village",]
centrality_contact_vaillage_a_df <- centrality_contact_total2[centrality_contact_total2$village=="A_Village",]

centrality_enviro_total2 <- rbind(centrality_enviro_vaillage_m2, centrality_enviro_vaillage_s2,centrality_enviro_vaillage_a2) %>% select(Node, PageRank, gender,  age, house_sol, commercial_goods, village, school_level, bmi, landsize_in_daba, livestock, season, household_size, Strength, Betweenness, Eigen) %>%drop_na()#getting rid of all NA values removed 42 individuals (~5% of all nodes)



centrality_enviro_vaillage_m_df <- centrality_enviro_total2[centrality_enviro_total2$village=="M_Village",]
centrality_enviro_vaillage_s_df <- centrality_enviro_total2[centrality_enviro_total2$village=="S_Village",]
centrality_enviro_vaillage_a_df <- centrality_enviro_total2[centrality_enviro_total2$village=="A_Village",]


centrality_enviro_total2$season <- as.factor(centrality_enviro_total2$season)

centrality_distance_total2 <- rbind(centrality_distance_vaillage_m2, centrality_distance_vaillage_s2,centrality_distance_vaillage_a2) %>% select(Node, PageRank, gender,  age, house_sol, commercial_goods, village, school_level, bmi, landsize_in_daba, livestock, season, household_size, Strength, Betweenness, Eigen) %>%drop_na()#getting rid of all NA values removed 67 individuals (~5% of all nodes)

centrality_distance_total2$season <- as.factor(centrality_distance_total2$season)

centrality_distance_vaillage_m_df <- centrality_distance_total2[centrality_distance_total2$village=="M_Village",]
centrality_distance_vaillage_s_df <- centrality_distance_total2[centrality_distance_total2$village=="S_Village",]
centrality_distance_vaillage_a_df <- centrality_distance_total2[centrality_distance_total2$village=="A_Village",]

```

```{r Centrality Model Function}

perform_model_averaging <- function(data, season) {
  
  data$Strength <- scale(data$Strength,center = TRUE, scale = TRUE)
  data$PageRank <- scale(data$PageRank,center = TRUE, scale = TRUE)
  data$Betweenness <- scale(data$Betweenness,center = TRUE, scale = TRUE)
  data$Eigen <- scale(data$Eigen,center = TRUE, scale = TRUE)
  
  data <- data[data$season==season,]
  
  data$age <- scale(data$age)
  data$house_sol <- scale(data$house_sol)
  data$commercial_goods <- scale(data$commercial_goods)
  data$school_level <- scale(data$school_level)
  data$bmi <- scale(data$bmi)
  data$landsize_in_daba <- scale(data$landsize_in_daba)
  data$livestock <- scale(data$livestock)
  data$household_size <- scale(data$household_size)
  
  #PageRank
  pagerank_model <- lm(PageRank ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size, data = data )

  
  options(na.action = "na.fail")
  pagerank_model_selection <- dredge(pagerank_model)
  options(na.action = "na.omit")
  
  pagerank_model_average <- summary(model.avg(pagerank_model_selection))
  
  pagerank_model_average_coef <- as.data.frame(pagerank_model_average$coefmat.full)
  pagerank_model_average_coef$model <- "PageRank"
  pagerank_model_average_coef$Variable <- rownames(pagerank_model_average_coef)
  
  #Strength
   strength_model <- lm(Strength ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size, data = data )
  
  options(na.action = "na.fail")
  strength_model_selection <- dredge(strength_model)
  options(na.action = "na.omit")
  
  strength_model_average <- summary(model.avg(strength_model_selection))
  
  strength_model_average_coef <- as.data.frame(strength_model_average$coefmat.full)
  strength_model_average_coef$model <- "Strength"
  strength_model_average_coef$Variable <- rownames(strength_model_average_coef)
  
    
  #Betweenness Centrality
   betweenness_model <- lm(Betweenness ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size, data = data )
  
  options(na.action = "na.fail")
  betweenness_model_selection <- dredge(betweenness_model)
  options(na.action = "na.omit")
  
  betweenness_model_average <- summary(model.avg(betweenness_model_selection))
  
  betweenness_model_average_coef <- as.data.frame(betweenness_model_average$coefmat.full)
  betweenness_model_average_coef$model <- "betweenness"
  betweenness_model_average_coef$Variable <- rownames(betweenness_model_average_coef)
  betweenness_model_average_coef$`Adjusted SE` <- NA
    betweenness_model_average_coef <-   betweenness_model_average_coef %>% select(colnames(  strength_model_average_coef))
  
   #Eigenvector Centrality
  eigen_model <- lm(Eigen ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size, data = data )
  
  options(na.action = "na.fail")
  eigen_model_selection <- dredge(eigen_model)
  options(na.action = "na.omit")
  
  eigen_model_average <- summary(model.avg(eigen_model_selection))
  
  eigen_model_average_coef <- as.data.frame(eigen_model_average$coefmat.full)
  eigen_model_average_coef$model <- "eigen"
  eigen_model_average_coef$Variable <- rownames(eigen_model_average_coef)
  

  top_model <- get.models(eigen_model_selection, 1)[[1]]
  plot(top_model)
  
  model_average_coef <- rbind(pagerank_model_average_coef, strength_model_average_coef, betweenness_model_average_coef,eigen_model_average_coef)
  
  model_average_coef$season <- season
  
  return(model_average_coef)
}

 
model_avg_name_vaillage_m_1 <- perform_model_averaging(centrality_name_vaillage_m_df, 1)
model_avg_name_vaillage_m_1$Village <- "M_Village"
model_avg_name_vaillage_m_1$Network <- "Naming"

model_avg_name_vaillage_s_1 <- perform_model_averaging(centrality_name_vaillage_s_df, 1)
model_avg_name_vaillage_s_1$Village <- "S_Village"
model_avg_name_vaillage_s_1$Network <- "Naming"

model_avg_name_vaillage_a_1 <- perform_model_averaging(centrality_name_vaillage_a_df, 1)
model_avg_name_vaillage_a_1$Village <- "A_Village"
model_avg_name_vaillage_a_1$Network <- "Naming"

model_avg_name_vaillage_m_2 <- perform_model_averaging(centrality_name_vaillage_m_df, 2)
model_avg_name_vaillage_m_2$Village <- "M_Village"
model_avg_name_vaillage_m_2$Network <- "Naming"


model_avg_name_vaillage_s_2 <- perform_model_averaging(centrality_name_vaillage_s_df, 2)
model_avg_name_vaillage_s_2$Village <- "S_Village"
model_avg_name_vaillage_s_2$Network <- "Naming"

model_avg_name_vaillage_a_2 <- perform_model_averaging(centrality_name_vaillage_a_df, 2)
model_avg_name_vaillage_a_2$Village <- "A_Village"
model_avg_name_vaillage_a_2$Network <- "Naming"

model_avg_name_vaillage_s_3 <- perform_model_averaging(centrality_name_vaillage_s_df, 3)
model_avg_name_vaillage_s_3$Village <- "S_Village"
model_avg_name_vaillage_s_3$Network <- "Naming"

model_avg_name_vaillage_a_3 <- perform_model_averaging(centrality_name_vaillage_a_df, 3)
model_avg_name_vaillage_a_3$Village <- "A_Village"
model_avg_name_vaillage_a_3$Network <- "Naming"

model_avg_contact_vaillage_m_1 <- perform_model_averaging(centrality_contact_vaillage_m_df, 1)
model_avg_contact_vaillage_m_1$Village <- "M_Village"
model_avg_contact_vaillage_m_1$Network <- "Close Contact"

model_avg_contact_vaillage_s_1 <- perform_model_averaging(centrality_contact_vaillage_s_df, 1)
model_avg_contact_vaillage_s_1$Village <- "S_Village"
model_avg_contact_vaillage_s_1$Network <- "Close Contact"

model_avg_contact_vaillage_a_1 <- perform_model_averaging(centrality_contact_vaillage_a_df, 1)
model_avg_contact_vaillage_a_1$Village <- "A_Village"
model_avg_contact_vaillage_a_1$Network <- "Close Contact"

model_avg_contact_vaillage_s_2 <- perform_model_averaging(centrality_contact_vaillage_s_df, 2)
model_avg_contact_vaillage_s_2$Village <- "S_Village"
model_avg_contact_vaillage_s_2$Network <- "Close Contact"

model_avg_contact_vaillage_a_2 <- perform_model_averaging(centrality_contact_vaillage_a_df, 2)
model_avg_contact_vaillage_a_2$Village <- "A_Village"
model_avg_contact_vaillage_a_2$Network <- "Close Contact"

model_avg_contact_vaillage_s_3 <- perform_model_averaging(centrality_contact_vaillage_s_df, 3)
model_avg_contact_vaillage_s_3$Village <- "S_Village"
model_avg_contact_vaillage_s_3$Network <- "Close Contact"

model_avg_contact_vaillage_a_3 <- perform_model_averaging(centrality_contact_vaillage_a_df, 3)
model_avg_contact_vaillage_a_3$Village <- "A_Village"
model_avg_contact_vaillage_a_3$Network <- "Close Contact"

model_avg_enviro_vaillage_m_1 <- perform_model_averaging(centrality_enviro_vaillage_m_df, 1)
model_avg_enviro_vaillage_m_1$Village <- "M_Village"
model_avg_enviro_vaillage_m_1$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_s_1 <- perform_model_averaging(centrality_enviro_vaillage_s_df, 1)
model_avg_enviro_vaillage_s_1$Village <- "S_Village"
model_avg_enviro_vaillage_s_1$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_a_1 <- perform_model_averaging(centrality_enviro_vaillage_a_df, 1)
model_avg_enviro_vaillage_a_1$Village <- "A_Village"
model_avg_enviro_vaillage_a_1$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_m_2 <- perform_model_averaging(centrality_enviro_vaillage_m_df, 2)
model_avg_enviro_vaillage_m_2$Village <- "M_Village"
model_avg_enviro_vaillage_m_2$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_s_2 <- perform_model_averaging(centrality_enviro_vaillage_s_df, 2)
model_avg_enviro_vaillage_s_2$Village <- "S_Village"
model_avg_enviro_vaillage_s_2$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_a_2 <- perform_model_averaging(centrality_enviro_vaillage_a_df, 2)
model_avg_enviro_vaillage_a_2$Village <- "A_Village"
model_avg_enviro_vaillage_a_2$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_s_3 <- perform_model_averaging(centrality_enviro_vaillage_s_df, 3)
model_avg_enviro_vaillage_s_3$Village <- "S_Village"
model_avg_enviro_vaillage_s_3$Network <- "Environmental Overlap"

model_avg_enviro_vaillage_a_3 <- perform_model_averaging(centrality_enviro_vaillage_a_df, 3)
model_avg_enviro_vaillage_a_3$Village <- "A_Village"
model_avg_enviro_vaillage_a_3$Network <- "Environmental Overlap"



model_avg_distance_vaillage_m_1 <- perform_model_averaging(centrality_distance_vaillage_m_df, 1)
model_avg_distance_vaillage_m_1$Village <- "M_Village"
model_avg_distance_vaillage_m_1$Network <- "Distance"

model_avg_distance_vaillage_s_1 <- perform_model_averaging(centrality_distance_vaillage_s_df, 1)
model_avg_distance_vaillage_s_1$Village <- "S_Village"
model_avg_distance_vaillage_s_1$Network <- "Distance"

model_avg_distance_vaillage_a_1 <- perform_model_averaging(centrality_distance_vaillage_a_df, 1)
model_avg_distance_vaillage_a_1$Village <- "A_Village"
model_avg_distance_vaillage_a_1$Network <- "Distance"

model_avg_distance_vaillage_m_2 <- perform_model_averaging(centrality_distance_vaillage_m_df, 2)
model_avg_distance_vaillage_m_2$Village <- "M_Village"
model_avg_distance_vaillage_m_2$Network <- "Distance"

model_avg_distance_vaillage_s_2 <- perform_model_averaging(centrality_distance_vaillage_s_df, 2)
model_avg_distance_vaillage_s_2$Village <- "S_Village"
model_avg_distance_vaillage_s_2$Network <- "Distance"

model_avg_distance_vaillage_a_2 <- perform_model_averaging(centrality_distance_vaillage_a_df, 2)
model_avg_distance_vaillage_a_2$Village <- "A_Village"
model_avg_distance_vaillage_a_2$Network <- "Distance"

model_avg_distance_vaillage_s_3 <- perform_model_averaging(centrality_distance_vaillage_s_df, 3)
model_avg_distance_vaillage_s_3$Village <- "S_Village"
model_avg_distance_vaillage_s_3$Network <- "Distance"

model_avg_distance_vaillage_a_3 <- perform_model_averaging(centrality_distance_vaillage_a_df, 3)
model_avg_distance_vaillage_a_3$Village <- "A_Village"
model_avg_distance_vaillage_a_3$Network <- "Distance"

```

```{r Pool Effect Sizes in a Meta-analysis Approach: Naming Network}

name_function <- function(model2, variable){
  
model_df <- rbind(model_avg_name_vaillage_m_1, model_avg_name_vaillage_s_1, model_avg_name_vaillage_a_1, model_avg_name_vaillage_m_2, model_avg_name_vaillage_s_2, model_avg_name_vaillage_a_2,
model_avg_name_vaillage_s_3, model_avg_name_vaillage_a_3) %>% filter(model == model2) %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable==variable)

metaanalysis_model <- rma(yi = model_df$Estimate,   
         sei = model_df$`Std. Error`, 
         method = "FE")

output <- c(metaanalysis_model$beta , metaanalysis_model$ci.lb, metaanalysis_model$ci.ub, metaanalysis_model$pval, metaanalysis_model$se, nrow(model_df))

return(output)

}

name_age_pr <- name_function(model2 = "PageRank", variable ="age")
name_sex_pr <- name_function("PageRank","genderMale")
name_bmi_pr <- name_function("PageRank", "bmi")
name_hs_pr <- name_function("PageRank", "household_size")
name_hsol_pr <- name_function("PageRank", "house_sol")
name_cg_pr <- name_function("PageRank", "commercial_goods")
name_l_pr <- name_function("PageRank", "livestock")
name_ls_pr <- name_function("PageRank", "landsize_in_daba")
names_sl_pr <- name_function("PageRank", "school_level")

name_age_str <- name_function("Strength", "age")
name_sex_str <- name_function("Strength","genderMale")
name_bmi_str <- name_function("Strength", "bmi")
name_hs_str <- name_function("Strength", "household_size")
name_hsol_str <- name_function("Strength", "house_sol")
name_cg_str <- name_function("Strength", "commercial_goods")
name_l_str <- name_function("Strength", "livestock")
name_ls_str <- name_function("Strength", "landsize_in_daba")
names_sl_str <- name_function("Strength", "school_level")

name_age_btwn <- name_function("betweenness", "age")
name_sex_btwn <- name_function("betweenness","genderMale")
name_bmi_btwn <- name_function("betweenness", "bmi")
name_hs_btwn <- name_function("betweenness", "household_size")
name_hsol_btwn <- name_function("betweenness", "house_sol")
name_cg_btwn <- name_function("betweenness", "commercial_goods")
name_l_btwn <- name_function("betweenness", "livestock")
name_ls_btwn <- name_function("betweenness", "landsize_in_daba")
names_sl_btwn <- name_function("betweenness", "school_level")

name_age_egn <- name_function("eigen", "age")
name_sex_egn <- name_function("eigen","genderMale")
name_bmi_egn <- name_function("eigen", "bmi")
name_hs_egn <- name_function("eigen", "household_size")
name_hsol_egn <- name_function("eigen", "house_sol")
name_cg_egn <- name_function("eigen", "commercial_goods")
name_l_egn <- name_function("eigen", "livestock")
name_ls_egn <- name_function("eigen", "landsize_in_daba")
name_sl_egn <- name_function("eigen", "school_level")

name_intercept_egn <- name_function("eigen", "(Intercept)")

```

```{r Pool Effect Sizes in a Meta-analysis Approach: Environmental Overlap Network}


enviro_function <- function(model2, variable){
  
model_df <- rbind(model_avg_enviro_vaillage_m_1, model_avg_enviro_vaillage_s_1, model_avg_enviro_vaillage_a_1, model_avg_enviro_vaillage_m_2, model_avg_enviro_vaillage_s_2, model_avg_enviro_vaillage_a_2,
model_avg_enviro_vaillage_s_3, model_avg_enviro_vaillage_a_3) %>% filter(model == model2) %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable==variable)

metaanalysis_model <- rma(yi = model_df$Estimate,   
         sei = model_df$`Std. Error`, 
         method = "FE")

output <- c(metaanalysis_model$beta , metaanalysis_model$ci.lb, metaanalysis_model$ci.ub, metaanalysis_model$pval, metaanalysis_model$se, nrow(model_df))


return(output)

}

enviro_age_pr <- enviro_function(model2 = "PageRank", variable ="age")
enviro_sex_pr <- enviro_function("PageRank","genderMale")
enviro_bmi_pr <- enviro_function("PageRank", "bmi")
enviro_hs_pr <- enviro_function("PageRank", "household_size")
enviro_hsol_pr <- enviro_function("PageRank", "house_sol")
enviro_cg_pr <- enviro_function("PageRank", "commercial_goods")
enviro_l_pr <- enviro_function("PageRank", "livestock")
enviro_ls_pr <- enviro_function("PageRank", "landsize_in_daba")
enviros_sl_pr <- enviro_function("PageRank", "school_level")

enviro_age_str <- enviro_function("Strength", "age")
enviro_sex_str <- enviro_function("Strength","genderMale")
enviro_bmi_str <- enviro_function("Strength", "bmi")
enviro_hs_str <- enviro_function("Strength", "household_size")
enviro_hsol_str <- enviro_function("Strength", "house_sol")
enviro_cg_str <- enviro_function("Strength", "commercial_goods")
enviro_l_str <- enviro_function("Strength", "livestock")
enviro_ls_str <- enviro_function("Strength", "landsize_in_daba")
enviros_sl_str <- enviro_function("Strength", "school_level")

enviro_age_btwn <- enviro_function("betweenness", "age")
enviro_sex_btwn <- enviro_function("betweenness","genderMale")
enviro_bmi_btwn <- enviro_function("betweenness", "bmi")
enviro_hs_btwn <- enviro_function("betweenness", "household_size")
enviro_hsol_btwn <- enviro_function("betweenness", "house_sol")
enviro_cg_btwn <- enviro_function("betweenness", "commercial_goods")
enviro_l_btwn <- enviro_function("betweenness", "livestock")
enviro_ls_btwn <- enviro_function("betweenness", "landsize_in_daba")
enviros_sl_btwn <- enviro_function("betweenness", "school_level")

enviro_age_egn <- enviro_function("eigen", "age")
enviro_sex_egn <- enviro_function("eigen","genderMale")
enviro_bmi_egn <- enviro_function("eigen", "bmi")
enviro_hs_egn <- enviro_function("eigen", "household_size")
enviro_hsol_egn <- enviro_function("eigen", "house_sol")
enviro_cg_egn <- enviro_function("eigen", "commercial_goods")
enviro_l_egn <- enviro_function("eigen", "livestock")
enviro_ls_egn <- enviro_function("eigen", "landsize_in_daba")
enviro_sl_egn <- enviro_function("eigen", "school_level")

```

```{r Pool Effect Sizes in a Meta-analysis Approach: Close Contact Network}


contact_function <- function(model2, variable){
  
model_df <- rbind(model_avg_contact_vaillage_m_1, model_avg_contact_vaillage_s_1, model_avg_contact_vaillage_a_1,  model_avg_contact_vaillage_s_2, model_avg_contact_vaillage_a_2,
model_avg_contact_vaillage_s_3, model_avg_contact_vaillage_a_3) %>% filter(model == model2) %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable==variable)

metaanalysis_model <- rma(yi = model_df$Estimate,   
         sei = model_df$`Std. Error`, 
         method = "FE")

output <- c(metaanalysis_model$beta , metaanalysis_model$ci.lb, metaanalysis_model$ci.ub, metaanalysis_model$pval, metaanalysis_model$se, nrow(model_df))


return(output)

}

contact_age_pr <- contact_function(model2 = "PageRank", variable ="age")
contact_sex_pr <- contact_function("PageRank","genderMale")
contact_bmi_pr <- contact_function("PageRank", "bmi")
contact_hs_pr <- contact_function("PageRank", "household_size")
contact_hsol_pr <- contact_function("PageRank", "house_sol")
contact_cg_pr <- contact_function("PageRank", "commercial_goods")
contact_l_pr <- contact_function("PageRank", "livestock")
contact_ls_pr <- contact_function("PageRank", "landsize_in_daba")
contacts_sl_pr <- contact_function("PageRank", "school_level")

contact_age_str <- contact_function("Strength", "age")
contact_sex_str <- contact_function("Strength","genderMale")
contact_bmi_str <- contact_function("Strength", "bmi")
contact_hs_str <- contact_function("Strength", "household_size")
contact_hsol_str <- contact_function("Strength", "house_sol")
contact_cg_str <- contact_function("Strength", "commercial_goods")
contact_l_str <- contact_function("Strength", "livestock")
contact_ls_str <- contact_function("Strength", "landsize_in_daba")
contacts_sl_str <- contact_function("Strength", "school_level")

contact_age_btwn <- contact_function("betweenness", "age")
contact_sex_btwn <- contact_function("betweenness","genderMale")
contact_bmi_btwn <- contact_function("betweenness", "bmi")
contact_hs_btwn <- contact_function("betweenness", "household_size")
contact_hsol_btwn <- contact_function("betweenness", "house_sol")
contact_cg_btwn <- contact_function("betweenness", "commercial_goods")
contact_l_btwn <- contact_function("betweenness", "livestock")
contact_ls_btwn <- contact_function("betweenness", "landsize_in_daba")
contacts_sl_btwn <- contact_function("betweenness", "school_level")

contact_age_egn <- contact_function("eigen", "age")
contact_sex_egn <- contact_function("eigen","genderMale")
contact_bmi_egn <- contact_function("eigen", "bmi")
contact_hs_egn <- contact_function("eigen", "household_size")
contact_hsol_egn <- contact_function("eigen", "house_sol")
contact_cg_egn <- contact_function("eigen", "commercial_goods")
contact_l_egn <- contact_function("eigen", "livestock")
contact_ls_egn <- contact_function("eigen", "landsize_in_daba")
contact_sl_egn <- contact_function("eigen", "school_level")

```

```{r Pool Effect Sizes in a Meta-analysis Approach: Distance Network}


distance_function <- function(model2, variable){
  
model_df <- rbind(model_avg_distance_vaillage_m_1,   model_avg_distance_vaillage_s_2, model_avg_distance_vaillage_a_2,
model_avg_distance_vaillage_s_3, model_avg_distance_vaillage_a_3) %>% filter(model == model2) %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable==variable)

metaanalysis_model <- rma(yi = model_df$Estimate,   
         sei = model_df$`Std. Error`, 
         method = "FE")

output <- c(metaanalysis_model$beta , metaanalysis_model$ci.lb, metaanalysis_model$ci.ub, metaanalysis_model$pval, metaanalysis_model$se, nrow(model_df))

return(output)

}

distance_age_pr <- distance_function(model2 = "PageRank", variable ="age")
distance_sex_pr <- distance_function("PageRank","genderMale")
distance_bmi_pr <- distance_function("PageRank", "bmi")
distance_hs_pr <- distance_function("PageRank", "household_size")
distance_hsol_pr <- distance_function("PageRank", "house_sol")
distance_cg_pr <- distance_function("PageRank", "commercial_goods")
distance_l_pr <- distance_function("PageRank", "livestock")
distance_ls_pr <- distance_function("PageRank", "landsize_in_daba")
distances_sl_pr <- distance_function("PageRank", "school_level")

distance_age_str <- distance_function("Strength", "age")
distance_sex_str <- distance_function("Strength","genderMale")
distance_bmi_str <- distance_function("Strength", "bmi")
distance_hs_str <- distance_function("Strength", "household_size")
distance_hsol_str <- distance_function("Strength", "house_sol")
distance_cg_str <- distance_function("Strength", "commercial_goods")
distance_l_str <- distance_function("Strength", "livestock")
distance_ls_str <- distance_function("Strength", "landsize_in_daba")
distances_sl_str <- distance_function("Strength", "school_level")

distance_age_btwn <- distance_function("betweenness", "age")
distance_sex_btwn <- distance_function("betweenness","genderMale")
distance_bmi_btwn <- distance_function("betweenness", "bmi")
distance_hs_btwn <- distance_function("betweenness", "household_size")
distance_hsol_btwn <- distance_function("betweenness", "house_sol")
distance_cg_btwn <- distance_function("betweenness", "commercial_goods")
distance_l_btwn <- distance_function("betweenness", "livestock")
distance_ls_btwn <- distance_function("betweenness", "landsize_in_daba")
distances_sl_btwn <- distance_function("betweenness", "school_level")

distance_age_egn <- distance_function("eigen", "age")
distance_sex_egn <- distance_function("eigen","genderMale")
distance_bmi_egn <- distance_function("eigen", "bmi")
distance_hs_egn <- distance_function("eigen", "household_size")
distance_hsol_egn <- distance_function("eigen", "house_sol")
distance_cg_egn <- distance_function("eigen", "commercial_goods")
distance_l_egn <- distance_function("eigen", "livestock")
distance_ls_egn <- distance_function("eigen", "landsize_in_daba")
distance_sl_egn <- distance_function("eigen", "school_level")

```

```{r Compile effect sizes}

meta_naming_df <- as.data.frame(rbind(name_age_pr ,name_sex_pr ,name_bmi_pr,name_hs_pr,name_hsol_pr,name_cg_pr, name_l_pr,name_ls_pr, names_sl_pr,
                                      name_age_str ,name_sex_str ,name_bmi_str,name_hs_str,name_hsol_str,name_cg_str, name_l_str,name_ls_str, names_sl_str,
                                      name_age_btwn ,name_sex_btwn ,name_bmi_btwn  ,name_hs_btwn,name_hsol_btwn,name_cg_btwn, name_l_btwn,name_ls_btwn, names_sl_btwn,
                                      name_age_egn ,name_sex_egn ,name_bmi_egn  ,name_hs_egn,name_hsol_egn,name_cg_egn, name_l_egn,name_ls_egn, name_sl_egn))

colnames(meta_naming_df) <- c("Estimate", "Lower", "Upper", "Pr(>|z|)","SE", "N")

meta_naming_df$Centrality <- c( rep("PageRank", 9), rep("Strength", 9), rep("Betweenness", 9), rep("Eigenvector", 9) )

meta_naming_df$Variable <- rep(c("Age", "Gender[Man]", "BMI", "Household Size", "House Materials", "Commercial Goods", "Livestock", "Land Size", "School Level"), 4)


meta_naming_df$Significance <- meta_naming_df$`Pr(>|z|)`
meta_naming_df$Significance[meta_naming_df$`Pr(>|z|)`<0.05] <- "p < 0.05"
meta_naming_df$Significance[meta_naming_df$`Pr(>|z|)`<0.1 & meta_naming_df$`Pr(>|z|)`> 0.05] <- "0.05 < p < 0.1"
meta_naming_df$Significance[meta_naming_df$`Pr(>|z|)`>=0.1] <- "p > 0.1"


meta_enviro_df <- as.data.frame(rbind(enviro_age_pr ,enviro_sex_pr ,enviro_bmi_pr  ,enviro_hs_pr,enviro_hsol_pr,enviro_cg_pr, enviro_l_pr,enviro_ls_pr, enviros_sl_pr,
                                      enviro_age_str ,enviro_sex_str ,enviro_bmi_str  ,enviro_hs_str,enviro_hsol_str,enviro_cg_str, enviro_l_str,enviro_ls_str, enviros_sl_str,
                                      enviro_age_btwn ,enviro_sex_btwn ,enviro_bmi_btwn ,enviro_hs_btwn,enviro_hsol_btwn,enviro_cg_btwn, enviro_l_btwn,enviro_ls_btwn, enviros_sl_btwn,
                                      enviro_age_egn ,enviro_sex_egn ,enviro_bmi_egn ,enviro_hs_egn,enviro_hsol_egn,enviro_cg_egn, enviro_l_egn,enviro_ls_egn, enviro_sl_egn))

colnames(meta_enviro_df) <- c("Estimate", "Lower", "Upper", "Pr(>|z|)","SE", "N")

meta_enviro_df$Centrality <- c( rep("PageRank", 9), rep("Strength", 9), rep("Betweenness", 9), rep("Eigenvector", 9) )

meta_enviro_df$Variable <- rep(c("Age", "Gender[Man]", "BMI",  "Household Size", "House Materials", "Commercial Goods", "Livestock", "Land Size", "School Level"), 4)


meta_enviro_df$Significance <- meta_enviro_df$`Pr(>|z|)`
meta_enviro_df$Significance[meta_enviro_df$`Pr(>|z|)`<0.05] <- "p < 0.05"
meta_enviro_df$Significance[meta_enviro_df$`Pr(>|z|)`<0.1 & meta_enviro_df$`Pr(>|z|)`> 0.05] <- "0.05 < p < 0.1"
meta_enviro_df$Significance[meta_enviro_df$`Pr(>|z|)`>=0.1] <- "p > 0.1"


meta_contact_df <- as.data.frame(rbind(contact_age_pr ,contact_sex_pr ,contact_bmi_pr  ,contact_hs_pr,contact_hsol_pr,contact_cg_pr, contact_l_pr,contact_ls_pr, contacts_sl_pr,
                                      contact_age_str ,contact_sex_str ,contact_bmi_str  ,contact_hs_str,contact_hsol_str,contact_cg_str, contact_l_str,contact_ls_str, contacts_sl_str,
                                      contact_age_btwn ,contact_sex_btwn ,contact_bmi_btwn ,contact_hs_btwn,contact_hsol_btwn,contact_cg_btwn, contact_l_btwn,contact_ls_btwn, contacts_sl_btwn,
                                      contact_age_egn ,contact_sex_egn ,contact_bmi_egn ,contact_hs_egn,contact_hsol_egn,contact_cg_egn, contact_l_egn,contact_ls_egn, contact_sl_egn))

colnames(meta_contact_df) <- c("Estimate", "Lower", "Upper", "Pr(>|z|)","SE", "N")

meta_contact_df$Centrality <- c( rep("PageRank", 9), rep("Strength", 9), rep("Betweenness", 9), rep("Eigenvector", 9) )

meta_contact_df$Variable <- rep(c("Age", "Gender[Man]", "BMI", "Household Size", "House Materials", "Commercial Goods", "Livestock", "Land Size", "School Level"), 4)


meta_contact_df$Significance <- meta_contact_df$`Pr(>|z|)`
meta_contact_df$Significance[meta_contact_df$`Pr(>|z|)`<0.05] <- "p < 0.05"
meta_contact_df$Significance[meta_contact_df$`Pr(>|z|)`<0.1 & meta_contact_df$`Pr(>|z|)`> 0.05] <- "0.05 < p < 0.1"
meta_contact_df$Significance[meta_contact_df$`Pr(>|z|)`>=0.1] <- "p > 0.1"


meta_distance_df <- as.data.frame(rbind(distance_age_pr ,distance_sex_pr ,distance_bmi_pr  ,distance_hs_pr,distance_hsol_pr,distance_cg_pr, distance_l_pr,distance_ls_pr, distances_sl_pr,
                                      distance_age_str ,distance_sex_str ,distance_bmi_str ,distance_hs_str,distance_hsol_str,distance_cg_str, distance_l_str,distance_ls_str, distances_sl_str,
                                      distance_age_btwn ,distance_sex_btwn ,distance_bmi_btwn  ,distance_hs_btwn,distance_hsol_btwn,distance_cg_btwn, distance_l_btwn,distance_ls_btwn, distances_sl_btwn,
                                      distance_age_egn ,distance_sex_egn ,distance_bmi_egn  ,distance_hs_egn,distance_hsol_egn,distance_cg_egn, distance_l_egn,distance_ls_egn, distance_sl_egn))

colnames(meta_distance_df) <- c("Estimate", "Lower", "Upper", "Pr(>|z|)","SE", "N")

meta_distance_df$Centrality <- c( rep("PageRank", 9), rep("Strength", 9), rep("Betweenness", 9), rep("Eigenvector", 9) )

meta_distance_df$Variable <- rep(c("Age", "Gender[Man]", "BMI",  "Household Size", "House Materials", "Commercial Goods", "Livestock", "Land Size", "School Level"), 4)


meta_distance_df$Significance <- meta_distance_df$`Pr(>|z|)`
meta_distance_df$Significance[meta_distance_df$`Pr(>|z|)`<0.05] <- "p < 0.05"
meta_distance_df$Significance[meta_distance_df$`Pr(>|z|)`<0.1 & meta_distance_df$`Pr(>|z|)`> 0.05] <- "0.05 < p < 0.1"
meta_distance_df$Significance[meta_distance_df$`Pr(>|z|)`>=0.1] <- "p > 0.1"


write.csv(meta_distance_df, "meta_distance_df.csv")
write.csv(meta_enviro_df, "meta_enviro_df.csv")
write.csv(meta_naming_df, "meta_naming_df.csv")
write.csv(meta_contact_df, "meta_contact_df.csv")

```

```{r Plotting Results}
meta_naming_df$Centrality <- factor(meta_naming_df$Centrality, levels=c("Betweenness" ,"Eigenvector", "PageRank", "Strength"))

meta_naming_df$lower_ci <- meta_naming_df$Estimate - 1.96 * meta_naming_df$SE
meta_naming_df$upper_ci <- meta_naming_df$Estimate + 1.96 * meta_naming_df$SE

meta_naming_df$lower_ci_90 <- meta_naming_df$Estimate - qnorm(0.95) * meta_naming_df$SE
meta_naming_df$upper_ci_90 <- meta_naming_df$Estimate + qnorm(0.95) * meta_naming_df$SE


meta_naming_df$Variable[meta_naming_df$Variable=="School Level"] <- "Education"
meta_naming_df$Variable[meta_naming_df$Variable=="Commercial Goods"] <- "Goods"
meta_naming_df$Variable[meta_naming_df$Variable=="Sex[Male]"] <- "Gender[Man]"
#meta_naming_df$Variable[meta_naming_df$Variable=="HSOL"] <- "House Materials"

meta_naming_df$Variable <- factor(meta_naming_df$Variable, levels = c("Livestock", "Land Size", "House Materials", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]"))

naming_plot2 <- ggplot(data=meta_naming_df, aes(x = Variable , y= Estimate))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.2) +
  geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.5) +
 geom_point(aes(x = Variable, y= Estimate, color = Variable, shape= Centrality, size=Centrality), position = position_dodge(1))+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Social")+
  coord_flip()+
  scale_shape_manual(values=c(19, 17, 15, 18))+
  scale_color_manual(values = c("darkorchid","springgreen3","orange2", "deeppink3" ,"cyan3", "darkgreen", "gold2" , "magenta3" , "blue3"), na.translate = F)+
  theme_classic()+
  scale_size_manual(values = c(1.5, 1.5, 1.5, 2.2))


meta_enviro_df$Centrality <- factor(meta_enviro_df$Centrality, levels=c("Betweenness" ,"Eigenvector", "PageRank", "Strength"))

meta_enviro_df$lower_ci <- meta_enviro_df$Estimate - 1.96 * meta_enviro_df$SE
meta_enviro_df$upper_ci <- meta_enviro_df$Estimate + 1.96 * meta_enviro_df$SE

meta_enviro_df$lower_ci_90 <- meta_enviro_df$Estimate - qnorm(0.95) * meta_enviro_df$SE
meta_enviro_df$upper_ci_90 <- meta_enviro_df$Estimate + qnorm(0.95) * meta_enviro_df$SE


meta_enviro_df$Variable[meta_enviro_df$Variable=="School Level"] <- "Education"
meta_enviro_df$Variable[meta_enviro_df$Variable=="Commercial Goods"] <- "Goods"
#meta_enviro_df$Variable[meta_enviro_df$Variable=="HSOL"] <- "House Materials"
meta_enviro_df$Variable[meta_enviro_df$Variable=="Sex[Male]"] <- "Gender[Man]"

meta_enviro_df$Variable <- factor(meta_enviro_df$Variable, levels = c("Livestock", "Land Size", "House Materials", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]"))

enviro_plot2 <- ggplot(data=meta_enviro_df, aes(x = Variable , y= Estimate))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.2) +
  geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.5) +
 geom_point(aes(x = Variable, y= Estimate, color = Variable, shape= Centrality, size=Centrality), position = position_dodge(1))+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Environmental")+
  coord_flip()+
  scale_shape_manual(values=c(19, 17, 15, 18))+
  scale_color_manual(values = c("darkorchid","springgreen3","orange2", "deeppink3" ,"cyan3", "darkgreen", "gold2" , "magenta3" , "blue3"), na.translate = F)+
  theme_classic()+
  scale_size_manual(values = c(1.5, 1.5, 1.5, 2.2))



meta_contact_df$Centrality <- factor(meta_contact_df$Centrality, levels=c("Betweenness" ,"Eigenvector", "PageRank", "Strength"))

meta_contact_df$lower_ci <- meta_contact_df$Estimate - 1.96 * meta_contact_df$SE
meta_contact_df$upper_ci <- meta_contact_df$Estimate + 1.96 * meta_contact_df$SE

meta_contact_df$lower_ci_90 <- meta_contact_df$Estimate - qnorm(0.95) * meta_contact_df$SE
meta_contact_df$upper_ci_90 <- meta_contact_df$Estimate + qnorm(0.95) * meta_contact_df$SE


meta_contact_df$Variable[meta_contact_df$Variable=="School Level"] <- "Education"
meta_contact_df$Variable[meta_contact_df$Variable=="Commercial Goods"] <- "Goods"
#meta_contact_df$Variable[meta_contact_df$Variable=="HSOL"] <- "House Materials"

meta_contact_df$Variable <- factor(meta_contact_df$Variable, levels = c("Livestock", "Land Size", "House Materials", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]"))

contact_plot2 <- ggplot(data=meta_contact_df, aes(x = Variable , y= Estimate))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.2) +
  geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.5) +
 geom_point(aes(x = Variable, y= Estimate, color = Variable, shape= Centrality, size=Centrality), position = position_dodge(1))+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Close Contact")+
  coord_flip()+
  scale_shape_manual(values=c(19, 17, 15, 18))+
  scale_color_manual(values = c("darkorchid","springgreen3","orange2", "deeppink3" ,"cyan3", "darkgreen", "gold2" , "magenta3" , "blue3"), na.translate = F)+
  theme_classic()+
  scale_size_manual(values = c(1.5, 1.5, 1.5, 2.2))



meta_distance_df$Centrality <- factor(meta_distance_df$Centrality, levels=c("Betweenness" ,"Eigenvector", "PageRank", "Strength"))

meta_distance_df$lower_ci <- meta_distance_df$Estimate - 1.96 * meta_distance_df$SE
meta_distance_df$upper_ci <- meta_distance_df$Estimate + 1.96 * meta_distance_df$SE

meta_distance_df$lower_ci_90 <- meta_distance_df$Estimate - qnorm(0.95) * meta_distance_df$SE
meta_distance_df$upper_ci_90 <- meta_distance_df$Estimate + qnorm(0.95) * meta_distance_df$SE


meta_distance_df$Variable[meta_distance_df$Variable=="School Level"] <- "Education"
meta_distance_df$Variable[meta_distance_df$Variable=="Commercial Goods"] <- "Goods"
#meta_distance_df$Variable[meta_distance_df$Variable=="HSOL"] <- "House Materials"

meta_distance_df$Variable <- factor(meta_distance_df$Variable, levels = c("Livestock", "Land Size", "House Materials", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]"))

distance_plot2 <- ggplot(data=meta_distance_df, aes(x = Variable , y= Estimate))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.2) +
  geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90, color = Variable, group=Centrality), position = position_dodge(1), linewidth=0.5) +
 geom_point(aes(x = Variable, y= Estimate, color = Variable, shape= Centrality, size=Centrality), position = position_dodge(1))+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Household")+
  coord_flip()+
  scale_shape_manual(values=c(19, 17, 15, 18))+
  scale_color_manual(values = c("darkorchid","springgreen3","orange2", "deeppink3" ,"cyan3", "darkgreen", "gold2" , "magenta3" , "blue3"), na.translate = F)+
  theme_classic()+
  scale_size_manual(values = c(1.5, 1.5, 1.5, 2.2))



ggarrange(naming_plot2,  contact_plot2, distance_plot2,enviro_plot2, common.legend = TRUE, labels = c("a", "b", "c", "d"))

write.csv(meta_naming_df, "meta_naming_df.csv")
write.csv(meta_contact_df, "meta_contact_df.csv")
write.csv(meta_enviro_df, "meta_enviro_df.csv")
write.csv(meta_distance_df, "meta_distance_df.csv")



```

```{r}

model_df1 <- rbind(model_avg_name_vaillage_m_1, model_avg_name_vaillage_s_1, model_avg_name_vaillage_a_1, model_avg_name_vaillage_m_2, model_avg_name_vaillage_s_2, model_avg_name_vaillage_a_2,
model_avg_name_vaillage_s_3, model_avg_name_vaillage_a_3) %>% filter(model == "betweenness") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")

model_df2 <- rbind(model_avg_name_vaillage_m_1, model_avg_name_vaillage_s_1, model_avg_name_vaillage_a_1, model_avg_name_vaillage_m_2, model_avg_name_vaillage_s_2, model_avg_name_vaillage_a_2,
model_avg_name_vaillage_s_3, model_avg_name_vaillage_a_3) %>% filter(model == "eigen") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")

model_df3 <- rbind(model_avg_name_vaillage_m_1, model_avg_name_vaillage_s_1, model_avg_name_vaillage_a_1, model_avg_name_vaillage_m_2, model_avg_name_vaillage_s_2, model_avg_name_vaillage_a_2,
model_avg_name_vaillage_s_3, model_avg_name_vaillage_a_3) %>% filter(model == "PageRank") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")

model_df4 <- rbind(model_avg_name_vaillage_m_1, model_avg_name_vaillage_s_1, model_avg_name_vaillage_a_1, model_avg_name_vaillage_m_2, model_avg_name_vaillage_s_2, model_avg_name_vaillage_a_2,
model_avg_name_vaillage_s_3, model_avg_name_vaillage_a_3) %>% filter(model == "Strength") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")



model_df1_enviro <- rbind(model_avg_enviro_vaillage_m_1, model_avg_enviro_vaillage_s_1, model_avg_enviro_vaillage_a_1, model_avg_enviro_vaillage_m_2, model_avg_enviro_vaillage_s_2, model_avg_enviro_vaillage_a_2,
model_avg_enviro_vaillage_s_3, model_avg_enviro_vaillage_a_3) %>% filter(model == "betweenness") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")

model_df2_enviro <- rbind(model_avg_enviro_vaillage_m_1, model_avg_enviro_vaillage_s_1, model_avg_enviro_vaillage_a_1, model_avg_enviro_vaillage_m_2, model_avg_enviro_vaillage_s_2, model_avg_enviro_vaillage_a_2,
model_avg_enviro_vaillage_s_3, model_avg_enviro_vaillage_a_3) %>% filter(model == "eigen") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")

model_df3_enviro <- rbind(model_avg_enviro_vaillage_m_1, model_avg_enviro_vaillage_s_1, model_avg_enviro_vaillage_a_1, model_avg_enviro_vaillage_m_2, model_avg_enviro_vaillage_s_2, model_avg_enviro_vaillage_a_2,
model_avg_enviro_vaillage_s_3, model_avg_enviro_vaillage_a_3) %>% filter(model == "PageRank") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")

model_df4 <- rbind(model_avg_enviro_vaillage_m_1, model_avg_enviro_vaillage_s_1, model_avg_enviro_vaillage_a_1, model_avg_enviro_vaillage_m_2, model_avg_enviro_vaillage_s_2, model_avg_enviro_vaillage_a_2,
model_avg_enviro_vaillage_s_3, model_avg_enviro_vaillage_a_3) %>% filter(model == "Strength") %>% select(Estimate, `Std. Error`, Variable, model, season, Village, Network) %>% filter(Variable== "genderMale")
```

```{r}

mand_eigen_name <- ggplot(data = centrality_name_vaillage_m, aes(x= log(Eigen)))+
  geom_histogram()+
  theme_classic()

sara_eigen_name <- ggplot(data = centrality_name_vaillage_s, aes(x= log(Eigen)))+
  geom_histogram()+
  theme_classic()

anda_eigen_name <- ggplot(data = centrality_name_vaillage_a, aes(x= log(Eigen)))+
  geom_histogram()+
  theme_classic()

ggarrange(mand_eigen_name, sara_eigen_name, anda_eigen_name, ncol=3, labels = c("M_Village", "S_Village", "A_Village"))


mand_eigen_name_scatter <- ggplot(data = centrality_name_vaillage_m, aes(x= Eigen, y= PageRank))+
  geom_point()+
  theme_classic()

sara_eigen_name_scatter <- ggplot(data = centrality_name_vaillage_s, aes(x= Eigen, y= PageRank))+
  geom_point()+
  theme_classic()

anda_eigen_name_scatter <- ggplot(data = centrality_name_vaillage_a, aes(x= Eigen, y= PageRank))+
  geom_point()+
  theme_classic()

ggarrange(mand_eigen_name_scatter, sara_eigen_name_scatter, anda_eigen_name_scatter, ncol=3, labels = c("M_Village", "S_Village", "A_Village"))


ggplot(data = centrality_distance_vaillage_m, aes(x= Eigen))+
  geom_histogram()+
  theme_classic()


mand_eigen_name_scatter <- ggplot(data = centrality_name_vaillage_m, aes(x= Eigen, y= PageRank))+
  geom_point()+
  theme_classic()

unique(igraph::components(name_vaillage_m_g)$membership)
unique(igraph::components(name_vaillage_a_g)$membership)
unique(igraph::components(name_vaillage_s_g)$membership)

```

```{r Regression(Not meta-analysis): Naming}

centrality_name_total2$Betweenness <- scale(centrality_name_total2$Betweenness)
centrality_name_total2$Eigen <- scale(centrality_name_total2$Eigen)
centrality_name_total2$PageRank <- scale(centrality_name_total2$PageRank)
centrality_name_total2$Strength <- scale(centrality_name_total2$Strength)

#centrality_name_total2$centrality <- (centrality_name_total2$Betweenness +centrality_name_total2$Eigen+centrality_name_total2$PageRank+centrality_name_total2$Strength)/4

#centrality_name_total2$centrality <- scale(centrality_name_total2$centrality)

centrality_name_total3 <- centrality_name_total2 %>% pivot_longer(names_to = "centrality", values_to = "value", cols = c(Betweenness, Eigen, PageRank, Strength))

  
  centrality_name_total3$age <- scale(centrality_name_total3$age)
  centrality_name_total3$house_sol <- scale(centrality_name_total3$house_sol)
  centrality_name_total3$commercial_goods <- scale(centrality_name_total3$commercial_goods)
  centrality_name_total3$school_level <- scale(centrality_name_total3$school_level)
  centrality_name_total3$bmi <- scale(centrality_name_total3$bmi)
  centrality_name_total3$landsize_in_daba <- scale(centrality_name_total3$landsize_in_daba)
  centrality_name_total3$livestock <- scale(centrality_name_total3$livestock)
  centrality_name_total3$household_size <- scale(centrality_name_total3$household_size)
  
  centrality_name_model <- lmer(value ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size + season + village +centrality + (1|Node), data= centrality_name_total3)

  
  options(na.action = "na.fail")
  name_model_selection <- dredge(centrality_name_model)
  options(na.action = "na.omit")
  
  name_model_average <- summary(model.avg(  name_model_selection))
  
 name_model_average_coef <- as.data.frame(name_model_average$coefmat.full)
 name_model_average_coef$Variable <- rownames(name_model_average_coef)
 
 
 importance <- as.data.frame(sw(model.avg(name_model_selection)))
  colnames(importance)[1] <- "Importance"
  importance$Variable <- rownames(importance)
  importance$Variable[importance$Variable=="gender"] <- "genderMale"
 
 name_model_average_coef <- left_join(name_model_average_coef, importance, by = "Variable")
 
 name_model_average_coef$Importance[name_model_average_coef$Variable=="season2"] <-  importance$Importance[importance$Variable=="season"]
 
  name_model_average_coef$Importance[name_model_average_coef$Variable=="season3"] <-  importance$Importance[importance$Variable=="season"]
  
  name_model_average_coef$Importance[name_model_average_coef$Variable=="villageM_Village"] <-  importance$Importance[importance$Variable=="village"]
   
 name_model_average_coef$Importance[name_model_average_coef$Variable=="villageS_Village"] <-  importance$Importance[importance$Variable=="village"]
     
  name_model_average_coef$Importance[name_model_average_coef$Variable=="centralityEigen"|name_model_average_coef$Variable=="centralityStrength"|name_model_average_coef$Variable=="centralityPageRank"] <-  importance$Importance[importance$Variable=="centrality"]
  
  #  name_model_average_coef$Importance[name_model_average_coef$Variable=="CentralityEigen:landsize_in_daba"|name_model_average_coef$Variable=="CentralityStrength:landsize_in_daba"|name_model_average_coef$Variable=="CentralityPageRank:landsize_in_daba"] <-  importance$Importance[importance$Variable=="Centrality:landsize_in_daba"]
 
 name_model_average_coef$lower_ci <-  name_model_average_coef$Estimate - 1.96 *  name_model_average_coef$`Std. Error`
 
 name_model_average_coef$upper_ci <-  name_model_average_coef$Estimate + 1.96 *  name_model_average_coef$`Std. Error`

 name_model_average_coef$lower_ci_90 <-  name_model_average_coef$Estimate - 1.644854*  name_model_average_coef$`Std. Error`
 
 name_model_average_coef$upper_ci_90 <-  name_model_average_coef$Estimate +1.644854*  name_model_average_coef$`Std. Error`

 
  name_model_average_coef$Variable <- c("Intercept","Season[2]", "Season[3]", "Village[M]", "Village[S]", "Gender[Man]","Livestock", "BMI","Education","Goods","Land Size","Age", "Household Size", "House Material" ,"Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]")

  name_model_average_coef$Variable <- factor(name_model_average_coef$Variable, levels = c("Livestock", "Land Size", "House Material", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]", "Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]", "Season[2]", "Season[3]", "Village[M]", "Village[S]", "Intercept"))
  
  name_model_average_coef <- name_model_average_coef%>% filter(Variable!="Centrality[Eigenvector]", 
                                                               Variable!="Centrality[PageRank]",
                                                               Variable!="Centrality[Strength]",
                                                               Variable!="Season[2]",
                                                               Variable!="Season[3]",
                                                               Variable!="Village[M]",
                                                               Variable!="Village[S]",
                                                               Variable!= "Intercept")
 
name_plot3 <- ggplot(data=name_model_average_coef, aes(x = Variable , y= Estimate, color = Importance))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci), position = position_dodge(1), linewidth=0.3) +
geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90), position = position_dodge(1), linewidth=1) +
 geom_point(aes(x = Variable, y= Estimate), position = position_dodge(1), size=2)+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Social")+
  coord_flip()+
  theme_classic()+
  scale_color_viridis_c(
    direction = -1, limits = c(0,1.01))


```

```{r Regression(Not meta-analysis): Distance}

centrality_distance_total2$Betweenness <- scale(centrality_distance_total2$Betweenness)
centrality_distance_total2$Eigen <- scale(centrality_distance_total2$Eigen)
centrality_distance_total2$PageRank <- scale(centrality_distance_total2$PageRank)
centrality_distance_total2$Strength <- scale(centrality_distance_total2$Strength)


#centrality_distance_total2$centrality <- (centrality_distance_total2$Betweenness +centrality_distance_total2$Eigen+centrality_distance_total2$PageRank+centrality_distance_total2$Strength)/4

#centrality_distance_total2$centrality <- scale(centrality_distance_total2$centrality)


centrality_distance_total3 <- centrality_distance_total2 %>% pivot_longer(names_to = "centrality", values_to = "value", cols = c(Betweenness, Eigen, PageRank, Strength))


  centrality_distance_total3$age <- scale(centrality_distance_total3$age)
  centrality_distance_total3$house_sol <- scale(centrality_distance_total3$house_sol)
  centrality_distance_total3$commercial_goods <- scale(centrality_distance_total3$commercial_goods)
  centrality_distance_total3$school_level <- scale(centrality_distance_total3$school_level)
  centrality_distance_total3$bmi <- scale(centrality_distance_total3$bmi)
  centrality_distance_total3$landsize_in_daba <- scale(centrality_distance_total3$landsize_in_daba)
  centrality_distance_total3$livestock <- scale(centrality_distance_total3$livestock)
  centrality_distance_total3$household_size <- scale(centrality_distance_total3$household_size)
  
    #centrality_distance_total3$value <- scale(centrality_distance_total3$value)
  
  centrality_distance_model <- lmer(value ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size + season + village+ centrality+ (1|Node)
                           , data= centrality_distance_total3)

  
  options(na.action = "na.fail")
  distance_model_selection <- dredge(centrality_distance_model)
  options(na.action = "na.omit")
  
  distance_model_average <- summary(model.avg(distance_model_selection))
  
 distance_model_average_coef <- as.data.frame(distance_model_average$coefmat.full)
 distance_model_average_coef$Variable <- rownames(distance_model_average_coef)
 
 
  
 importance <- as.data.frame(sw(model.avg(distance_model_selection)))
  colnames(importance)[1] <- "Importance"
  importance$Variable <- rownames(importance)
  importance$Variable[importance$Variable=="gender"] <- "genderMale"
 
 distance_model_average_coef <- left_join(distance_model_average_coef, importance, by = "Variable")
 
 distance_model_average_coef$Importance[distance_model_average_coef$Variable=="season2"] <-  importance$Importance[importance$Variable=="season"]
 
  distance_model_average_coef$Importance[distance_model_average_coef$Variable=="season3"] <-  importance$Importance[importance$Variable=="season"]
  
  distance_model_average_coef$Importance[distance_model_average_coef$Variable=="villageM_Village"] <-  importance$Importance[importance$Variable=="village"]
   
 distance_model_average_coef$Importance[distance_model_average_coef$Variable=="villageS_Village"] <-  importance$Importance[importance$Variable=="village"]
     
  distance_model_average_coef$Importance[distance_model_average_coef$Variable=="centralityEigen"|distance_model_average_coef$Variable=="centralityStrength"|distance_model_average_coef$Variable=="centralityPageRank"] <-  importance$Importance[importance$Variable=="centrality"]
  
  #distance_model_average_coef$Importance[distance_model_average_coef$Variable=="CentralityEigen:house_sol"|distance_model_average_coef$Variable=="CentralityStrength:house_sol"|distance_model_average_coef$Variable=="CentralityPageRank:house_sol"] <-  importance$Importance[importance$Variable=="Centrality:house_sol"]
     
 
 distance_model_average_coef$lower_ci <-  distance_model_average_coef$Estimate - 1.96 *  distance_model_average_coef$`Std. Error`
 
 distance_model_average_coef$upper_ci <-  distance_model_average_coef$Estimate + 1.96 *  distance_model_average_coef$`Std. Error`

 distance_model_average_coef$lower_ci_90 <-  distance_model_average_coef$Estimate - 1.644854*  distance_model_average_coef$`Std. Error`
 
 distance_model_average_coef$upper_ci_90 <-  distance_model_average_coef$Estimate +1.644854*  distance_model_average_coef$`Std. Error`

 
 
  distance_model_average_coef$Variable <- c("Intercept", "Season[2]", "Season[3]","Village[M]", "Village[S]","Goods", "Land Size", "Age", "House Material", "Household Size", "Gender[Man]", "Livestock","BMI", "Education","Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]" )

  distance_model_average_coef$Variable <- factor(distance_model_average_coef$Variable, levels = c("Livestock", "Land Size", "House Material", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]", "Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]", "Season[2]", "Season[3]", "Village[M]", "Village[S]", "Intercept"))
  
distance_model_average_coef <-distance_model_average_coef%>% filter(Variable!="Centrality[Eigenvector]", 
                                                               Variable!="Centrality[PageRank]",
                                                               Variable!="Centrality[Strength]",
                                                               Variable!="Season[2]",
                                                               Variable!="Season[3]",
                                                               Variable!="Village[M]",
                                                               Variable!="Village[S]",
                                                               Variable!= "Intercept")
 
distance_plot3 <- ggplot(data=distance_model_average_coef, aes(x = Variable , y= Estimate, color = Importance))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci), position = position_dodge(1), linewidth=0.3) +
geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90), position = position_dodge(1), linewidth=1) +
 geom_point(aes(x = Variable, y= Estimate), position = position_dodge(1), size=2)+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Household")+
  coord_flip()+
  theme_classic()+
  scale_color_viridis_c(
    direction = -1, limits = c(0,1.01))

```

```{r Regression(Not meta-analysis): Contact}

centrality_contact_total2$Betweenness <- scale(centrality_contact_total2$Betweenness)
centrality_contact_total2$Eigen <- scale(centrality_contact_total2$Eigen)
centrality_contact_total2$PageRank <- scale(centrality_contact_total2$PageRank)
centrality_contact_total2$Strength <- scale(centrality_contact_total2$Strength)


centrality_contact_total3 <- centrality_contact_total2 %>% pivot_longer(names_to = "centrality", values_to = "value", cols = c(Betweenness, Eigen, PageRank, Strength))

  
  centrality_contact_total3$age <- scale(centrality_contact_total3$age)
  centrality_contact_total3$house_sol <- scale(centrality_contact_total3$house_sol)
  centrality_contact_total3$commercial_goods <- scale(centrality_contact_total3$commercial_goods)
  centrality_contact_total3$school_level <- scale(centrality_contact_total3$school_level)
  centrality_contact_total3$bmi <- scale(centrality_contact_total3$bmi)
  centrality_contact_total3$landsize_in_daba <- scale(centrality_contact_total3$landsize_in_daba)
  centrality_contact_total3$livestock <- scale(centrality_contact_total3$livestock)
  centrality_contact_total3$household_size <- scale(centrality_contact_total3$household_size)
  
  centrality_contact_model <- lmer(value ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size + season + village +centrality +(1|Node) , data= centrality_contact_total3)

  
  options(na.action = "na.fail")
  contact_model_selection <- dredge(centrality_contact_model)
  options(na.action = "na.omit")
  
  contact_model_average <- summary(model.avg(contact_model_selection))
  

 contact_model_average_coef <- as.data.frame(contact_model_average$coefmat.full)
 contact_model_average_coef$Variable <- rownames(contact_model_average_coef)
 
  importance <- as.data.frame(sw(model.avg(contact_model_selection)))
  colnames(importance)[1] <- "Importance"
  importance$Variable <- rownames(importance)
  importance$Variable[importance$Variable=="gender"] <- "genderMale"
 
 contact_model_average_coef <- left_join(contact_model_average_coef, importance, by = "Variable")
 
 contact_model_average_coef$Importance[contact_model_average_coef$Variable=="season2"] <-  importance$Importance[importance$Variable=="season"]
 
  contact_model_average_coef$Importance[contact_model_average_coef$Variable=="season3"] <-  importance$Importance[importance$Variable=="season"]
  
   contact_model_average_coef$Importance[contact_model_average_coef$Variable=="villageM_Village"] <-  importance$Importance[importance$Variable=="village"]
   
     contact_model_average_coef$Importance[contact_model_average_coef$Variable=="villageS_Village"] <-  importance$Importance[importance$Variable=="village"]
     
  #contact_model_average_coef$Importance[contact_model_average_coef$Variable=="CentralityEigen:landsize_in_daba"|contact_model_average_coef$Variable=="CentralityStrength:landsize_in_daba"|contact_model_average_coef$Variable=="CentralityPageRank:landsize_in_daba"] <-  importance$Importance[importance$Variable=="Centrality:landsize_in_daba"]
  
    #contact_model_average_coef$Importance[contact_model_average_coef$Variable=="CentralityEigen:school_level"|contact_model_average_coef$Variable=="CentralityStrength:school_level"|contact_model_average_coef$Variable=="CentralityPageRank:school_level"] <-  importance$Importance[importance$Variable=="Centrality:school_level"]
 
     contact_model_average_coef$Importance[contact_model_average_coef$Variable=="centralityEigen"|contact_model_average_coef$Variable=="centralityStrength"|contact_model_average_coef$Variable=="centralityPageRank"] <-  importance$Importance[importance$Variable=="centrality"]
     
          #contact_model_average_coef$Importance[contact_model_average_coef$Variable=="CentralityEigen:household_size"|contact_model_average_coef$Variable=="CentralityStrength:household_size"|contact_model_average_coef$Variable=="CentralityPageRank:household_size"] <-  importance$Importance[importance$Variable=="Centrality:household_size"]
 
# contact_model_average_coef$Importance[is.na(contact_model_average_coef$Importance)] <- 0
 
 contact_model_average_coef$lower_ci <-  contact_model_average_coef$Estimate - 1.96 *  contact_model_average_coef$`Std. Error`
 
 contact_model_average_coef$upper_ci <-  contact_model_average_coef$Estimate + 1.96 *  contact_model_average_coef$`Std. Error`

 contact_model_average_coef$lower_ci_90 <-  contact_model_average_coef$Estimate - 1.644854*  contact_model_average_coef$`Std. Error`
 
 contact_model_average_coef$upper_ci_90 <-  contact_model_average_coef$Estimate +1.644854*  contact_model_average_coef$`Std. Error`

  contact_model_average_coef$Variable <- c("Intercept", "House Material","Season[2]", "Season[3]","Village[M]", "Village[S]",  "Age", "Gender[Man]", "Land Size", "Goods", "BMI", "Household Size", "Livestock", "Education","Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]")

  contact_model_average_coef$Variable <- factor(contact_model_average_coef$Variable, levels = c("Livestock", "Land Size", "House Material", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]", "Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]", "Season[2]", "Season[3]", "Village[M]", "Village[S]", "Intercept"))
  
contact_model_average_coef <- contact_model_average_coef%>% filter(Variable!="Centrality[Eigenvector]", 
                                                               Variable!="Centrality[PageRank]",
                                                               Variable!="Centrality[Strength]",
                                                               Variable!="Season[2]",
                                                               Variable!="Season[3]",
                                                               Variable!="Village[M]",
                                                               Variable!="Village[S]",
                                                               Variable!= "Intercept")
 
contact_plot3 <- ggplot(data=contact_model_average_coef, aes(x = Variable , y= Estimate, color=Importance))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci), position = position_dodge(1), linewidth=0.3) +
geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90), position = position_dodge(1), linewidth=1) +
 geom_point(aes(x = Variable, y= Estimate), position = position_dodge(1), size=2)+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Close Contact")+
  coord_flip()+
  theme_classic()+
  scale_color_viridis_c(
    direction = -1, limits = c(0,1.01))

```

```{r Regression(Not meta-analysis): Enviro}


centrality_enviro_total2$Betweenness <- scale(centrality_enviro_total2$Betweenness)
centrality_enviro_total2$Eigen <- scale(centrality_enviro_total2$Eigen)
centrality_enviro_total2$PageRank <- scale(centrality_enviro_total2$PageRank)
centrality_enviro_total2$Strength <- scale(centrality_enviro_total2$Strength)


#centrality_enviro_total2$centrality <- (centrality_enviro_total2$Betweenness +centrality_enviro_total2$Eigen+centrality_enviro_total2$PageRank+centrality_enviro_total2$Strength)/4
#centrality_enviro_total2$centrality <- scale(centrality_enviro_total2$centrality)


centrality_enviro_total3 <- centrality_enviro_total2 %>% pivot_longer(names_to = "centrality", values_to = "value", cols = c(Betweenness, Eigen, PageRank, Strength))
  
  centrality_enviro_total3$age <- scale(centrality_enviro_total3$age)
  centrality_enviro_total3$house_sol <- scale(centrality_enviro_total3$house_sol)
  centrality_enviro_total3$commercial_goods <- scale(centrality_enviro_total3$commercial_goods)
  centrality_enviro_total3$school_level <- scale(centrality_enviro_total3$school_level)
  centrality_enviro_total3$bmi <- scale(centrality_enviro_total3$bmi)
  centrality_enviro_total3$landsize_in_daba <- scale(centrality_enviro_total3$landsize_in_daba)
  centrality_enviro_total3$livestock <- scale(centrality_enviro_total3$livestock)
  centrality_enviro_total3$household_size <- scale(centrality_enviro_total3$household_size)
  
  centrality_enviro_model <- lmer(value ~gender + age + house_sol + commercial_goods+ school_level+bmi+landsize_in_daba +livestock + household_size + season + village +centrality + (1|Node)
                           , data= centrality_enviro_total3)
  
  summary_no_scale <- summary(centrality_enviro_model)
  summary_scale <- summary(centrality_enviro_model)
  
  options(na.action = "na.fail")
  enviro_model_selection <- dredge(centrality_enviro_model)
  options(na.action = "na.omit")
  
  enviro_model_average <- summary(model.avg(enviro_model_selection))
  
 enviro_model_average_coef <- as.data.frame(enviro_model_average$coefmat.full)
 enviro_model_average_coef$Variable <- rownames(enviro_model_average_coef)
 
 
 importance <- as.data.frame(sw(model.avg(enviro_model_selection)))
  colnames(importance)[1] <- "Importance"
  importance$Variable <- rownames(importance)
  importance$Variable[importance$Variable=="gender"] <- "genderMale"
 
 enviro_model_average_coef <- left_join(enviro_model_average_coef, importance, by = "Variable")
 
 enviro_model_average_coef$Importance[enviro_model_average_coef$Variable=="season2"] <-  importance$Importance[importance$Variable=="season"]
 
  enviro_model_average_coef$Importance[enviro_model_average_coef$Variable=="season3"] <-  importance$Importance[importance$Variable=="season"]
  
   enviro_model_average_coef$Importance[enviro_model_average_coef$Variable=="villageM_Village"] <-  importance$Importance[importance$Variable=="village"]
   
     enviro_model_average_coef$Importance[enviro_model_average_coef$Variable=="villageS_Village"] <-  importance$Importance[importance$Variable=="village"]
     
  #enviro_model_average_coef$Importance[enviro_model_average_coef$Variable=="age:CentralityEigen"|enviro_model_average_coef$Variable=="age:CentralityStrength"|enviro_model_average_coef$Variable=="age:CentralityPageRank"] <-  importance$Importance[importance$Variable=="age:Centrality"]
  
  enviro_model_average_coef$Importance[enviro_model_average_coef$Variable=="centralityEigen"|enviro_model_average_coef$Variable=="centralityStrength"|enviro_model_average_coef$Variable=="centralityPageRank"] <-  importance$Importance[importance$Variable=="centrality"]
     
 
 enviro_model_average_coef$lower_ci <-  enviro_model_average_coef$Estimate - 1.96 *  enviro_model_average_coef$`Std. Error`
 
 enviro_model_average_coef$upper_ci <-  enviro_model_average_coef$Estimate + 1.96 *  enviro_model_average_coef$`Std. Error`

 enviro_model_average_coef$lower_ci_90 <-  enviro_model_average_coef$Estimate - 1.644854*  enviro_model_average_coef$`Std. Error`
 
 enviro_model_average_coef$upper_ci_90 <-  enviro_model_average_coef$Estimate +1.644854*  enviro_model_average_coef$`Std. Error`

  enviro_model_average_coef$Variable <- c("Intercept", "Gender[Man]","Season[2]", "Season[3]","Village[M]", "Village[S]", "House Material","Goods", "Land Size", "Age", "Education",  "BMI","Household Size", "Livestock", "Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]")

  enviro_model_average_coef$Variable <- factor(enviro_model_average_coef$Variable, levels =  c("Livestock", "Land Size", "House Material", "Household Size", "Goods", "Education","BMI", "Age", "Gender[Man]", "Centrality[Eigenvector]", "Centrality[PageRank]", "Centrality[Strength]", "Season[2]", "Season[3]", "Village[M]", "Village[S]", "Intercept"))
  
enviro_model_average_coef <- enviro_model_average_coef%>% filter(Variable!="Centrality[Eigenvector]", 
                                                               Variable!="Centrality[PageRank]",
                                                               Variable!="Centrality[Strength]",
                                                               Variable!="Season[2]",
                                                               Variable!="Season[3]",
                                                               Variable!="Village[M]",
                                                               Variable!="Village[S]",
                                                               Variable!= "Intercept")
 
enviro_plot3 <- ggplot(data=enviro_model_average_coef, aes(x = Variable , y= Estimate, color =Importance))+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_linerange(aes(ymin = lower_ci, ymax = upper_ci), position = position_dodge(1), linewidth=0.3) +
geom_linerange(aes(ymin = lower_ci_90, ymax = upper_ci_90), position = position_dodge(1), linewidth=1) +
 geom_point(aes(x = Variable, y= Estimate), position = position_dodge(1), size=2)+
  xlab("")+
  ylab("Standardized Effect on Centrality")+
  ggtitle("Environmental")+
  coord_flip()+
  theme_classic()+
  scale_color_viridis_c(
    direction = -1, limits = c(0,1.01))

main_fig <- ggarrange(name_plot3, contact_plot3, distance_plot3, enviro_plot3, common.legend = TRUE, labels = c("a", "b", "c", "d"))

ggarrange(name_plot4, contact_plot4, distance_plot4, enviro_plot4, common.legend = TRUE, labels = c("a", "b", "c", "d"))


write.csv(enviro_model_average_coef, "enviro_model_average_coef.csv")
write.csv(contact_model_average_coef, "contact_model_average_coef.csv")
write.csv(name_model_average_coef, "naming_model_average_coef.csv")
write.csv(distance_model_average_coef, "distance_model_average_coef.csv")

ggsave("main_fig.png", plot = main_fig, width = 9, height = 7)



library(magick)

# Read in the PNG
img <- image_read("main_fig.png")

# Write out as JPEG
image_write(img, path = "main_fig.jpg", format = "jpeg")

image_write(img, path = "main_fig.tif", format = "tiff", compression = "lzw")

```

```{r}
# Extract the best model from the model selection process
best_model <- get.models(name_model_selection, subset = 1)[[1]]

summary(best_model )

car::vif(centrality_name_model)


resid_panel(centrality_name_model)
resid_panel(centrality_contact_model)
resid_panel(centrality_enviro_model)
resid_panel(centrality_distance_model)

summary(centrality_distance_model)

centrality_distance_model


hist(centrality_distance_model$residuals, breaks = 30, main = "Histogram of Model Residuals",
     xlab = "Residuals", col = "lightblue", border = "black")



```


```{r}

centrality_enviro_total3$Type <- "Env."
centrality_name_total3$Type <- "Name"
centrality_distance_total3$Type <- "House"
centrality_contact_total3$Type <- "Contact"

centrality_total <- rbind(centrality_enviro_total3, centrality_name_total3, centrality_distance_total3, centrality_contact_total3)%>%pivot_wider(values_from = value, names_from = centrality) %>% select(Node, PageRank, Strength, Betweenness, Eigen, Type)
colnames(centrality_total) <- c("Node", "PR", "Stg", "Btwn", "Egn", "Type")

centrality_total <- centrality_total %>% pivot_wider(names_from = Type, values_from = c(PR, Stg, Btwn, Egn))


#install.packages("psych")  # Run this only if you haven't already
library(psych)
 
centrality_vars <- centrality_total[ , -1]

centrality_vars2 <- centrality_vars %>% drop_na()

pca_result <- prcomp(centrality_vars2)

# View summary
summary(pca_result)

pca_result$x

# Scree plot of variance explained
plot(pca_result, type = "l")

# Biplot of PCA
biplot(pca_result)

library(psych)
fa.parallel(centrality_vars2, fa = "fa")  # Gives you a scree plot + eigenvalues

fa.parallel(centrality_vars2, fa = "fa", n.iter = 100, show.legend = TRUE, main = "Parallel Analysis Scree Plot")

scree(centrality_vars2) 

factor_analysis <- fa(centrality_vars2, 
                nfactors = 1,
                rotate = "varimax",
                fm = "ml",              # Maximum Likelihood method
                scores = "regression", missing = TRUE) 

factor_analysis$loadings
factor_analysis$communalities
factor_analysis$Vaccounted
print(factor_analysis)

fa.diagram(factor_analysis, cex = 1.2, e.size = 0.1, rsize = 0.2)
 


```