LinkPredictionFI_reciprocated.Rmd

---
title: "Link Prediction only on reciprocated ties"
author: "Francesco Beghini"
date: "2023-11-22"
output: html_notebook
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, dev = c('png','svglite'))
```

# Link Prediction Methods
```{r message=FALSE, warning=FALSE, include=FALSE}
Covars_ego <- covariates %>%
    unite('indigenous_status', mb_a0500a:mb_a0500i, na.rm = TRUE) |> 
  dplyr::rename(ego = respondent_master_id, ego_age = age_at_survey, ego_gender = gender,
          ego_household = building_id,
         ego_household_wealth = household_wealth_index_w3, ego_education = b0100, ego_indigenous = indigenous_status, 
         ego_religion = b0600, ego_length_in_village = b0800, ego_length_in_village_years = b0900,
         ego_bristol = mb_m0300,
         ego_DDS = DDS, ego_watersource = l0400,
         ego_med_painkillers = med_painkillers, 
         ego_med_antibiotics = med_antibiotics,  
         ego_med_antidiarrheal = med_antidiarrheal,  
         ego_med_antiparasite = med_antiparasite,  
         ego_med_vitamins = med_vitamins, 
         ego_med_zinc = med_zinc,  
         ego_med_antifungal = med_antifungal,  
         ego_med_antihypertensives = med_antihypertensives,  
         ego_med_antidiabetics = med_antidiabetics,  
         ego_med_antiacids = med_antiacids,  
         ego_med_laxatives = med_laxatives,
         ego_latitude = building_latitude, ego_longitude = building_longitude) %>%
  mutate(ego_education = recode(ego_education, "Have not completed any type of school" = 0,
                                "1st grade" = 1, "2nd grade" = 2, "3rd grade" = 3,
                                "4th grade" = 4, "5th grade" = 5, "6th grade" = 6,
                                "Some secondary" = 7, "Secondary" = 8, "More than secondary" = 9),
         ego_indigenous = recode(ego_indigenous, "Lenca" = 1, "Ch'orti'/Maya Chorti" = 1, .default = 0),
         ego_bristol = as.numeric(ego_bristol)) %>% 
  dplyr::select(-c("village_name"))

Covars_alter <- covariates %>% 
      unite('indigenous_status', mb_a0500a:mb_a0500i, na.rm = TRUE) |> 
  dplyr::rename(alter = respondent_master_id, alter_age = age_at_survey, alter_gender = gender,
         alter_household = building_id,
         alter_household_wealth = household_wealth_index_w3, alter_education = b0100, alter_indigenous = indigenous_status, 
         alter_religion = b0600, alter_length_in_village = b0800, alter_length_in_village_years = b0900,
         alter_bristol = mb_m0300,
         alter_DDS = DDS, alter_watersource = l0400, 
         alter_med_painkillers = med_painkillers,  
         alter_med_antibiotics = med_antibiotics,  
         alter_med_antidiarrheal = med_antidiarrheal,  
         alter_med_antiparasite = med_antiparasite,  
         alter_med_vitamins = med_vitamins, 
         alter_med_zinc = med_zinc,  
         alter_med_antifungal = med_antifungal,  
         alter_med_antihypertensives = med_antihypertensives,  
         alter_med_antidiabetics = med_antidiabetics,  
         alter_med_antiacids = med_antiacids,  
         alter_med_laxatives = med_laxatives,
         alter_latitude = building_latitude, alter_longitude = building_longitude) %>%
  mutate(alter_education = recode(alter_education, "Have not completed any type of school" = 0,
                                "1st grade" = 1, "2nd grade" = 2, "3rd grade" = 3,
                                "4th grade" = 4, "5th grade" = 5, "6th grade" = 6,
                                "Some secondary" = 7, "Secondary" = 8, "More than secondary" = 9),
         alter_indigenous = recode(alter_indigenous, "Lenca" = 1, "Ch'orti'/Maya Chorti" = 1, .default = 0),
         alter_bristol = as.numeric(alter_bristol)) %>% 
  dplyr::select(-c("village_name"))

SN_All_Covars_reciprocated_ties <- SN |>
                                    filter(pair_key %in% E(SN_mutual_conns)$pair_key) |>
                                    inner_join(Covars_ego, by = join_by("ego"=="ego")) |>
                                    inner_join(Covars_alter, by = join_by("alter"=="alter"))
```

```{r get_covariate_matrix_no_relationship, message=FALSE, warning=FALSE, include=FALSE}
#Get covariate matrix for people without a relationship
village_names <- unique(SN_All_Covars_reciprocated_ties$village_code_w3)
foreach(v = village_names, .combine = rbind) %do%
{
  #Split up matrix into villages
  SN_Village_w3 <- SN_All_Covars_reciprocated_ties %>% filter(village_code_w3 == v)
  #IDs of everyone in village
  Village_IDs_w3 <- unique(c(as.character(SN_Village_w3$ego), as.character(SN_Village_w3$alter)))
  data.frame(t(combn(Village_IDs_w3, 2, simplify = TRUE))) |> 
    dplyr::rename(ego = X1, alter = X2) |> 
    mutate(pair = vapply(stri_split_boundaries(str_c(ego, alter), type='character'), function(x) stri_c(x[stri_order(x)], collapse = ''), '')) |> 
    filter(!pair %in% SN_Village_w3$pair) |> 
    inner_join(Covars_ego, by = 'ego') |> 
    inner_join(Covars_alter, by = 'alter', suffix = c('_ego','_alter')) |> 
    mutate(
      relationship = "Un-Nominated Same Village",
      village_code_w3 = v,
      gender.mm = ego_gender == 'man' & alter_gender == 'man',
      gender.mf = ego_gender == "man" & alter_gender == "woman" | ego_gender == "woman" & alter_gender == "man",
      gender.ff = ego_gender == 'woman' & alter_gender == 'woman',
      indigenous.both = ego_indigenous == 1 & alter_indigenous == 1, 
      indigenous.one = ego_indigenous == 1 & alter_indigenous == 0 | ego_indigenous == 0 & alter_indigenous == 1,
      indigenous.none = ego_indigenous == 0 & alter_indigenous == 0,
      religion.same = ego_religion == alter_religion,
      age_difference_abs = abs(ego_age - alter_age),
      average_age = (ego_age + alter_age)/2,
      wealth_difference_abs = abs(ego_household_wealth - alter_household_wealth),
      average_wealth = (ego_household_wealth + alter_household_wealth)/2,
      education_difference_abs = abs(ego_education - alter_education),
      average_education = (ego_education + alter_education)/2,
      household_same = ego_household == alter_household,
      average_diet = (ego_DDS == alter_DDS) / 2,
      different_diet = abs(ego_DDS - alter_DDS),
      average_bristol = (ego_bristol == alter_bristol) / 2,
      different_bristol = abs(ego_bristol - alter_bristol),
      watersource_same = ego_watersource == alter_watersource,
      same_usage_painkillers = alter_med_painkillers == ego_med_painkillers,
      same_usage_antibiotics = alter_med_antibiotics == ego_med_antibiotics,
      same_usage_antidiarrheal = alter_med_antidiarrheal == ego_med_antidiarrheal,
      same_usage_antiparasite = alter_med_antiparasite == ego_med_antiparasite,
      same_usage_vitamins = alter_med_vitamins == ego_med_vitamins,
      same_usage_zinc = alter_med_zinc == ego_med_zinc,
      same_usage_antifungal = alter_med_antifungal == ego_med_antifungal,
      same_usage_antihypertensives = alter_med_antihypertensives == ego_med_antihypertensives,
      same_usage_antidiabetics = alter_med_antidiabetics == ego_med_antidiabetics,
      same_usage_antiacids = alter_med_antiacids == ego_med_antiacids,
      same_usage_laxatives = alter_med_laxatives == ego_med_laxatives
    ) |> 
    rowwise() |> 
    mutate(strain_sharing_rate = strain_rate[ego, alter],
           bray_curtis_sim = species_bray[ego, alter]) |>
    ungroup()
} -> All_NoRel_Covars_reciprocated


#Get dataframe of relationships into format for regression
SN_All_Covars_for_Regression_reciprocated_ties <- SN_All_Covars_reciprocated_ties %>% distinct(pair_key, .keep_all = TRUE) %>%
  mutate(gender.mm = ego_gender == "man" & alter_gender == "man",
         gender.mf = ego_gender == "woman" & alter_gender == "man" | ego_gender == "man" & alter_gender == "woman",
         gender.ff = ego_gender == "woman" & alter_gender == "woman",
         indigenous.both = ego_indigenous == 1 & alter_indigenous == 1,
         indigenous.one = ego_indigenous == 1 & alter_indigenous == 0 | ego_indigenous == 0 & alter_indigenous == 1,
         indigenous.none = ego_indigenous == 0 & alter_indigenous == 0,
         religion.same = ego_religion == alter_religion,
         age_difference_abs = abs(ego_age - alter_age),
         average_age = (ego_age + alter_age)/2,
         wealth_difference_abs = abs(ego_household_wealth - alter_household_wealth),
         average_wealth = (ego_household_wealth + alter_household_wealth)/2,
         education_difference_abs = abs(ego_education - alter_education),
         average_education = (ego_education + alter_education)/2,
         household_same = ego_household == alter_household,
         average_diet = (ego_DDS == alter_DDS) / 2,
         different_diet = abs(ego_DDS - alter_DDS),
         average_bristol = (ego_bristol == alter_bristol) / 2,
         different_bristol = abs(ego_bristol - alter_bristol),
         watersource_same = ego_watersource == alter_watersource, 
         same_usage_painkillers = alter_med_painkillers == ego_med_painkillers,
         same_usage_antibiotics = alter_med_antibiotics == ego_med_antibiotics,
         same_usage_antidiarrheal = alter_med_antidiarrheal == ego_med_antidiarrheal,
         same_usage_antiparasite = alter_med_antiparasite == ego_med_antiparasite,
         same_usage_vitamins = alter_med_vitamins == ego_med_vitamins,
         same_usage_zinc = alter_med_zinc == ego_med_zinc,
         same_usage_antifungal = alter_med_antifungal == ego_med_antifungal,
         same_usage_antihypertensives = alter_med_antihypertensives == ego_med_antihypertensives,
         same_usage_antidiabetics = alter_med_antidiabetics == ego_med_antidiabetics,
         same_usage_antiacids = alter_med_antiacids == ego_med_antiacids,
         same_usage_laxatives = alter_med_laxatives == ego_med_laxatives)
```

```{r message=FALSE, warning=FALSE, include=FALSE}
family_house_pairs_reciprocated <- SN_All_Covars_reciprocated_ties |>
                                filter(building_id_ego == building_id_alter | relationship %in% c("father", "mother", "sibling", "child_over12_other_house")) |>
                                pull(pair_key)

#Make dataframe with only non-kin in different house
SN_Non_Kin_Dif_House_Covars_for_Regression_reciprocated_ties<- SN_All_Covars_for_Regression_reciprocated_ties %>%
  filter(!(pair_key %in% family_house_pairs_reciprocated))


#select predictors
SN_All_Covars_for_Regression_reciprocated_ties <- SN_All_Covars_for_Regression_reciprocated_ties %>%
  mutate(related = 1) %>% 
  dplyr::select("ego","alter","gender.mm", "gender.mf", "indigenous.both", "indigenous.one", "religion.same",
         "age_difference_abs", "average_age", "wealth_difference_abs",
         "average_wealth", "education_difference_abs","average_education", "household_same",
         "average_diet","different_diet", "average_bristol","different_bristol", "watersource_same", starts_with('same_usage'),
         "related", "strain_sharing_rate", "village_code_w3", "bray_curtis_sim")

SN_Non_Kin_Dif_House_Covars_for_Regression_reciprocated_ties <- SN_Non_Kin_Dif_House_Covars_for_Regression_reciprocated_ties %>%
  mutate(related = 1) %>% 
  dplyr::select("ego","alter","gender.mm", "gender.mf", "indigenous.both", "indigenous.one", "religion.same",
         "age_difference_abs", "average_age", "wealth_difference_abs",
         "average_wealth", "education_difference_abs","average_education", "household_same",
         "average_diet","different_diet", "average_bristol","different_bristol", "watersource_same", starts_with('same_usage'),
         "related", "strain_sharing_rate", "village_code_w3", "bray_curtis_sim")



All_NoRel_Covars_for_Regression_reciprocated_ties<- All_NoRel_Covars_reciprocated %>%
  mutate(related = 0) %>%
  dplyr::select("ego","alter","gender.mm", "gender.mf", "indigenous.both", "indigenous.one", "religion.same",
                "age_difference_abs", "average_age", "wealth_difference_abs",
                "average_wealth", "education_difference_abs","average_education", "household_same",
                "average_diet","different_diet", "average_bristol","different_bristol", "watersource_same", starts_with('same_usage'),
                "related", "strain_sharing_rate", "village_code_w3", "bray_curtis_sim")


#Combine relationships and non relationship dataframe
SN_All_Covars_for_Regression_complete_reciprocated_ties <- SN_All_Covars_for_Regression_reciprocated_ties[complete.cases(SN_All_Covars_for_Regression_reciprocated_ties),]

SN_Non_Kin_Dif_House_Covars_for_Regression_complete_reciprocated_ties <- SN_Non_Kin_Dif_House_Covars_for_Regression_reciprocated_ties[complete.cases(SN_Non_Kin_Dif_House_Covars_for_Regression_reciprocated_ties),]

All_for_Regression_reciprocated_ties <- rbind(All_NoRel_Covars_for_Regression_reciprocated_ties, SN_All_Covars_for_Regression_reciprocated_ties)

Non_Kin_House_for_Regression_reciprocated_ties <- rbind(All_NoRel_Covars_for_Regression_reciprocated_ties,
                                      SN_Non_Kin_Dif_House_Covars_for_Regression_complete_reciprocated_ties)

All_for_Regression_complete_reciprocated_ties <- All_for_Regression_reciprocated_ties[complete.cases(All_for_Regression_reciprocated_ties), ]

All_NoRel_Covars_for_Regression_complete_reciprocated_ties <- All_NoRel_Covars_for_Regression_reciprocated_ties[complete.cases(All_NoRel_Covars_for_Regression_reciprocated_ties), ]

Non_Kin_House_for_Regression_complete_reciprocated_ties <-
  Non_Kin_House_for_Regression_reciprocated_ties[complete.cases(Non_Kin_House_for_Regression_reciprocated_ties), ]

edge_counts <- SN_All_Covars_for_Regression_complete_reciprocated_ties %>%
  group_by(village_code_w3) %>%
  summarize(edge_count = sum(related))

edge_counts <- edge_counts[match(village_names, edge_counts$village_code_w3), ] %>%
  filter(!is.na(village_code_w3)) |> 
  pull(edge_count, village_code_w3)

edge_counts_nkh <- SN_Non_Kin_Dif_House_Covars_for_Regression_complete_reciprocated_ties %>%
  group_by(village_code_w3) %>%
  summarize(edge_count = sum(related))

edge_counts_nkh <- edge_counts_nkh[match(village_names, edge_counts_nkh$village_code_w3), ] %>%
  filter(!is.na(village_code_w3)) |> 
  pull(edge_count, village_code_w3)



```


```{r Downsample, message=FALSE, warning=FALSE, include=FALSE}
#Downsample

foreach(v = names(edge_counts), .combine = rbind) %do% {
  All_NoRel_Covars_for_Regression_complete_reciprocated_ties |>
  filter(village_code_w3 == as.integer(v)) |>
  slice_sample(n = edge_counts[[v]])
} -> Downsamp_all_reciprocated_ties
Downsamp_all_reciprocated_ties <- as.data.frame(Downsamp_all_reciprocated_ties)

Downsamp_all_reciprocated_ties <- rbind(Downsamp_all_reciprocated_ties, SN_All_Covars_for_Regression_complete_reciprocated_ties)

Downsamp_all_reciprocated_ties <- Downsamp_all_reciprocated_ties[complete.cases(Downsamp_all_reciprocated_ties),]

#Non-kin different house

foreach(v = names(edge_counts_nkh), .combine = rbind) %do% {
  All_NoRel_Covars_for_Regression_complete_reciprocated_ties %>%
  filter(village_code_w3 == as.integer(v)) %>% slice_sample(n = edge_counts_nkh[[v]])
} -> Downsamp_all_reciprocated_ties_nkh

Downsamp_all_reciprocated_ties_nkh <- as.data.frame(Downsamp_all_reciprocated_ties_nkh)
Downsamp_all_reciprocated_ties_nkh <- rbind(Downsamp_all_reciprocated_ties_nkh, SN_Non_Kin_Dif_House_Covars_for_Regression_complete_reciprocated_ties)

Downsamp_all_reciprocated_ties_nkh <- Downsamp_all_reciprocated_ties_nkh[complete.cases(Downsamp_all_reciprocated_ties_nkh),]
```



```{r oos_5cv, message=FALSE, warning=FALSE, include=FALSE}
cl <- parallel::makeCluster(5)
doParallel::registerDoParallel(cl)
foreach(i = c(1:5)) %dopar% {
  print(i)
  inds <- c(1:nrow(Downsamp_all_reciprocated_ties))
  test1_inds <- sample(inds, size = floor(.33*nrow(Downsamp_all_reciprocated_ties)), replace = F)
  train1_inds <- inds[!(inds %in% test1_inds)]
  test2_inds <- sample(train1_inds, size = floor(.5*length(train1_inds)), replace = F)
  train2_inds <- inds[!(inds %in% test2_inds)]
  test3_inds <- train2_inds[!(train2_inds %in% test1_inds)]
  train3_inds <- inds[!(inds %in% test3_inds)]

  train1 <- Downsamp_all_reciprocated_ties[train1_inds,]
  test1 <- Downsamp_all_reciprocated_ties[test1_inds,]
  train2 <- Downsamp_all_reciprocated_ties[train2_inds,]
  test2 <- Downsamp_all_reciprocated_ties[test2_inds,]
  train3 <- Downsamp_all_reciprocated_ties[train3_inds,]
  test3 <- Downsamp_all_reciprocated_ties[test3_inds,]

  logit.test.1 <- lme4::glmer(related ~ strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2 <- lme4::glmer(related ~ strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3 <- lme4::glmer(related ~ strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train3, family = binomial)
  preds.1 <- predict(logit.test.1, newdata = test1, type = "response")
  preds.2 <- predict(logit.test.2, newdata = test2, type = "response")
  preds.3 <- predict(logit.test.3, newdata = test3, type = "response")
  names(preds.3) <- test3_inds
  names(preds.2) <- test2_inds
  names(preds.1) <- test1_inds
  preds.round <- c(preds.1, preds.2, preds.3)
  preds.round <- preds.round[order(as.numeric(names(preds.round)))]
  
  
  logit.test.1 <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        household_same +(1|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2 <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        education_difference_abs + average_education +
                        household_same + (1|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3 <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        education_difference_abs + average_education +
                        household_same + (1|village_code_w3),
                              data = train3, family = binomial)
  preds.1 <- predict(logit.test.1, newdata = test1, type = "response")
  preds.2 <- predict(logit.test.2, newdata = test2, type = "response")
  preds.3 <- predict(logit.test.3, newdata = test3, type = "response")
  names(preds.3) <- test3_inds
  names(preds.2) <- test2_inds
  names(preds.1) <- test1_inds
  preds.onlycovars <- c(preds.1, preds.2, preds.3)
  preds.onlycovars <- preds.onlycovars[order(as.numeric(names(preds.onlycovars)))]
  

  logit.test.1.covars <- lme4::glmer(related ~ gender.mm + gender.mf + 
                                     age_difference_abs + average_age  +
                                    strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2.covars <- lme4::glmer(related ~ gender.mm + gender.mf +
                                       age_difference_abs + average_age +
                                       strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3.covars <- lme4::glmer(related ~ gender.mm + gender.mf + 
                                       age_difference_abs + average_age +
                                       strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train3, family = binomial)
  
  preds.1.covars <- predict(logit.test.1.covars, newdata = test1, type = "response")
  preds.2.covars <- predict(logit.test.2.covars, newdata = test2, type = "response")
  preds.3.covars <- predict(logit.test.3.covars, newdata = test3, type = "response")
  names(preds.3.covars) <- test3_inds
  names(preds.2.covars) <- test2_inds
  names(preds.1.covars) <- test1_inds
  preds.round.covars <- c(preds.1.covars, preds.2.covars, preds.3.covars)
  preds.round.covars <- preds.round.covars[order(as.numeric(names(preds.round.covars)))]
  
  
  logit.test.1.all <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2.all <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3.all <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train3, family = binomial)
  
  preds.1.all <- predict(logit.test.1.all, newdata = test1, type = "response")
  preds.2.all <- predict(logit.test.2.all, newdata = test2, type = "response")
  preds.3.all <- predict(logit.test.3.all, newdata = test3, type = "response")
  names(preds.3.all) <- test3_inds
  names(preds.2.all) <- test2_inds
  names(preds.1.all) <- test1_inds
  preds.round.all <- c(preds.1.all, preds.2.all, preds.3.all)
  preds.round.all <- preds.round.all[order(as.numeric(names(preds.round.all)))]
  
  list(preds.round, preds.onlycovars, preds.round.covars, preds.round.all)
} -> preds_all
parallel::stopCluster(cl)

pred.all <- rowSums(sapply(preds_all, `[[`, 1))/5
pred.all.onlycovars <- rowSums(sapply(preds_all, `[[`, 2))/5
pred.all.covars <- rowSums(sapply(preds_all, `[[`, 3))/5
pred.all.all <- rowSums(sapply(preds_all, `[[`, 4))/5

#add on predictions from strain-sharing only model
pred.all.acc <- pred.all.all
pred.all.acc[pred.all.acc >=.5] <- 1
pred.all.acc[pred.all.acc <.5] <- 0

Downsamp_all_reciprocated_ties$predicted <- pred.all.acc

logit.base.roc <- pROC::roc(Downsamp_all_reciprocated_ties$related, pred.all, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")
logit.onlycovars.roc <- pROC::roc(Downsamp_all_reciprocated_ties$related, pred.all.onlycovars, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")
logit.covars.roc <- pROC::roc(Downsamp_all_reciprocated_ties$related, pred.all.covars, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")
logit.all.roc <- pROC::roc(Downsamp_all_reciprocated_ties$related, pred.all.all, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")


roc.list <- list("Strain-Sharing Rate" = logit.base.roc,
                  "Only Socio-Dem Covars" = logit.onlycovars.roc,
                  "SSR + Age + Gender + Household Sharing" = logit.covars.roc,
                  "SSR + All Socio-Dem Covars" = logit.all.roc)

ci.list <- lapply(roc.list, ci.se, specificities = seq(0, 1, l = 25))

dat.ci.list <- lapply(ci.list, function(ciobj) 
  data.frame(x = as.numeric(rownames(ciobj)),
             lower = ciobj[, 1],
             upper = ciobj[, 3]))
```
## All Social or Familial Relationship (on reciprocated ties W3)
```{r plot_roc_all_famililar_social_reciprocated_ties, echo=FALSE, message=FALSE, warning=FALSE}
roc.plot.reciprocated <- ggroc(roc.list) +
  theme_minimal() +
  geom_abline(slope=1, intercept = 1, linetype = "dashed", alpha=1, color = "grey") +
  coord_equal()

for(i in 1:4) {
  roc.plot.reciprocated <- roc.plot.reciprocated + geom_ribbon(
    data = dat.ci.list[[i]],
    aes(x = x, ymin = lower, ymax = upper),
    fill = i + 1,
    alpha = 0.2,
    inherit.aes = F) 
  } 

roc.plot.reciprocated <- roc.plot.reciprocated + 
theme(plot.title = element_text(size=10, hjust = .5),
      legend.position = c(0.6, 0.2),
      legend.text = element_text(size=10),
      text = element_text(size = 10)) + 
      scale_color_manual(name = "AUCs:",
                         values=c("Red", 'Orange', "Green", "Blue"), 
                       labels=c(paste("Strain-Sharing Rate: ", round(logit.base.roc$auc,2),
                                      " ± ",round(sqrt(pROC::var(logit.base.roc)),3), sep=""), 
                                paste("Only Socio-Dem Covars:", round(logit.onlycovars.roc$auc,2 ),
                                      " ± ",round(sqrt(pROC::var(logit.onlycovars.roc)),3), sep=""), 
                                paste("SSR + Age + Gender: ", round(logit.covars.roc$auc, 2),
                                      " ± ",round(sqrt(pROC::var(logit.covars.roc)),3), sep=""),
                                paste("SSR + All Socio-Dem Covars: ", round(logit.all.roc$auc, 2),
                                      " ± ",round(sqrt(pROC::var(logit.all.roc)),3), sep="")
                                )
                       ) +
  ggtitle("All Social or Familial Relationships (reciprocated ties)")

roc.plot.reciprocated <- roc.plot.reciprocated +
  theme_pubr(legend =  c(0.6, 0.2)) +
  theme(plot.title = element_text(size = 10, hjust = .5)) +
  labs_pubr()

roc.plot.reciprocated <- roc.plot.reciprocated + theme(
  plot.title = element_text(size = 10, hjust = .5),
  legend.title = element_text(size = 8),
  legend.text = element_text(size = 8)
)

svglite('Figures/plot_roc_all_famililar_social_reciprocated_ties.svg')
roc.plot.reciprocated
dev.off()

png('Figures/plot_roc_all_famililar_social_reciprocated_ties.png', width = 900, height = 900)
roc.plot.reciprocated
dev.off()
```


```{r ood_5cv_nonkin, message=FALSE, warning=FALSE, include=FALSE}
cl <- parallel::makeCluster(5)
doParallel::registerDoParallel(cl)
foreach(i = c(1:5)) %dopar% {
  inds <- c(1:nrow(Downsamp_all_reciprocated_ties_nkh))
  test1_inds <- sample(inds, size = floor(.33*nrow(Downsamp_all_reciprocated_ties_nkh)), replace = F)
  train1_inds <- inds[!(inds %in% test1_inds)]
  test2_inds <- sample(train1_inds, size = floor(.5*length(train1_inds)), replace = F)
  train2_inds <- inds[!(inds %in% test2_inds)]
  test3_inds <- train2_inds[!(train2_inds %in% test1_inds)]
  train3_inds <- inds[!(inds %in% test3_inds)]
  
  train1 <- Downsamp_all_reciprocated_ties_nkh[train1_inds,]
  test1 <- Downsamp_all_reciprocated_ties_nkh[test1_inds,]
  train2 <- Downsamp_all_reciprocated_ties_nkh[train2_inds,]
  test2 <- Downsamp_all_reciprocated_ties_nkh[test2_inds,]
  train3 <- Downsamp_all_reciprocated_ties_nkh[train3_inds,]
  test3 <- Downsamp_all_reciprocated_ties_nkh[test3_inds,]
  
  logit.test.1 <- lme4::glmer(related ~ strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2 <- lme4::glmer(related ~ strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3 <- lme4::glmer(related ~ strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train3, family = binomial)
  preds.1 <- predict(logit.test.1, newdata = test1, type = "response")
  preds.2 <- predict(logit.test.2, newdata = test2, type = "response")
  preds.3 <- predict(logit.test.3, newdata = test3, type = "response")
  names(preds.3) <- test3_inds
  names(preds.2) <- test2_inds
  names(preds.1) <- test1_inds
  preds.round <- c(preds.1, preds.2, preds.3)
  preds.round <- preds.round[order(as.numeric(names(preds.round)))]
  
    
  logit.test.1 <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        household_same +(1|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2 <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        household_same + (1|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3 <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        household_same + (1|village_code_w3),
                              data = train3, family = binomial)
  preds.1 <- predict(logit.test.1, newdata = test1, type = "response")
  preds.2 <- predict(logit.test.2, newdata = test2, type = "response")
  preds.3 <- predict(logit.test.3, newdata = test3, type = "response")
  names(preds.3) <- test3_inds
  names(preds.2) <- test2_inds
  names(preds.1) <- test1_inds
  preds.onlycovars <- c(preds.1, preds.2, preds.3)
  preds.onlycovars <- preds.onlycovars[order(as.numeric(names(preds.onlycovars)))]
  

  logit.test.1.covars <- lme4::glmer(related ~ gender.mm + gender.mf + 
                                     age_difference_abs + average_age  +
                                    strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2.covars <- lme4::glmer(related ~ gender.mm + gender.mf +
                                       age_difference_abs + average_age +
                                       strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3.covars <- lme4::glmer(related ~ gender.mm + gender.mf + 
                                       age_difference_abs + average_age +
                                       strain_sharing_rate + household_same + (0+strain_sharing_rate|village_code_w3),
                              data = train3, family = binomial)
  
  preds.1.covars <- predict(logit.test.1.covars, newdata = test1, type = "response")
  preds.2.covars <- predict(logit.test.2.covars, newdata = test2, type = "response")
  preds.3.covars <- predict(logit.test.3.covars, newdata = test3, type = "response")
  names(preds.3.covars) <- test3_inds
  names(preds.2.covars) <- test2_inds
  names(preds.1.covars) <- test1_inds
  preds.round.covars <- c(preds.1.covars, preds.2.covars, preds.3.covars)
  preds.round.covars <- preds.round.covars[order(as.numeric(names(preds.round.covars)))]
  
  logit.test.1.all <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train1, family = binomial)
  logit.test.2.all <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train2, family = binomial)
  logit.test.3.all <- lme4::glmer(related ~ gender.mm + gender.mf + indigenous.both +
                        indigenous.one + religion.same + age_difference_abs +
                        average_age + wealth_difference_abs + average_wealth +
                        education_difference_abs + average_education +
                        average_diet + different_diet + average_bristol + different_bristol + watersource_same + 
                        same_usage_painkillers + same_usage_antibiotics + same_usage_antidiarrheal + same_usage_antiparasite + same_usage_vitamins + same_usage_zinc + same_usage_antifungal + same_usage_antihypertensives + same_usage_antidiabetics + same_usage_antiacids + same_usage_laxatives +
                        strain_sharing_rate + (0+strain_sharing_rate|village_code_w3),
                              data = train3, family = binomial)
  
  preds.1.all <- predict(logit.test.1.all, newdata = test1, type = "response")
  preds.2.all <- predict(logit.test.2.all, newdata = test2, type = "response")
  preds.3.all <- predict(logit.test.3.all, newdata = test3, type = "response")
  names(preds.3.all) <- test3_inds
  names(preds.2.all) <- test2_inds
  names(preds.1.all) <- test1_inds
  preds.round.all <- c(preds.1.all, preds.2.all, preds.3.all)
  preds.round.all <- preds.round.all[order(as.numeric(names(preds.round.all)))]
  
  list(preds.round, preds.onlycovars, preds.round.covars, preds.round.all)
} -> preds_nkh
parallel::stopCluster(cl)

pred.all <- rowSums(sapply(preds_nkh, `[[`, 1))/5
pred.all.onlycovars <- rowSums(sapply(preds_nkh, `[[`, 2))/5
pred.all.covars <- rowSums(sapply(preds_nkh, `[[`, 3))/5
pred.all.all <- rowSums(sapply(preds_nkh, `[[`, 4))/5

pred.all.acc <- pred.all.all
pred.all.acc[pred.all.acc >=.5] <- 1
pred.all.acc[pred.all.acc <.5] <- 0

Downsamp_all_reciprocated_ties_nkh$predicted <- pred.all.acc

logit.base.roc <- pROC::roc(Downsamp_all_reciprocated_ties_nkh$related, pred.all, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")
logit.onlycovars.roc <- pROC::roc(Downsamp_all_reciprocated_ties_nkh$related, pred.all.onlycovars, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")
logit.covars.roc <- pROC::roc(Downsamp_all_reciprocated_ties_nkh$related, pred.all.covars, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")
logit.all.roc <- pROC::roc(Downsamp_all_reciprocated_ties_nkh$related, pred.all.all, ci=TRUE, plot=FALSE,
                levels=c("0","1"), direction = "<")


roc.list <- list("Strain-Sharing Rate (reciprocated ties)" = logit.base.roc,
                  "Only Socio-Dem Covars" = logit.onlycovars.roc,
                  "SSR + Age + Gender + Household Sharing" = logit.covars.roc,
                  "SSR + All Socio-Dem Covars (reciprocated ties)" = logit.all.roc)

ci.list <- lapply(roc.list, ci.se, specificities = seq(0, 1, l = 25))

dat.ci.list <- lapply(ci.list, function(ciobj) 
  data.frame(x = as.numeric(rownames(ciobj)),
             lower = ciobj[, 1],
             upper = ciobj[, 3]))
```

## Non-Kin different house Relationships (on reciprocated ties W3) model
```{r plot_roc_non_kin_different_house, echo=FALSE, message=FALSE, warning=FALSE}
roc.plot.nkh.reciprocated <- ggroc(roc.list) +
  theme_minimal() +
  geom_abline(slope=1, intercept = 1, linetype = "dashed", alpha=1, color = "grey") +
  coord_equal()

for(i in 1:4) {
  roc.plot.nkh.reciprocated <- roc.plot.nkh.reciprocated + geom_ribbon(
    data = dat.ci.list[[i]],
    aes(x = x, ymin = lower, ymax = upper),
    fill = i + 1,
    alpha = 0.2,
    inherit.aes = F) 
  } 

roc.plot.nkh.reciprocated <- roc.plot.nkh.reciprocated + 
theme(plot.title = element_text(size=10, hjust = .5),
      legend.position = c(0.6, 0.2),
      legend.text = element_text(size=10),
      text = element_text(size = 10)) + 
      scale_color_manual(name = "AUCs:",
                         values=c("Red", 'Orange', "Green", "Blue"), 
                       labels=c(paste("Strain-Sharing Rate: ", round(logit.base.roc$auc,2),
                                      " ± ",round(sqrt(pROC::var(logit.base.roc)),3), sep=""), 
                                paste("Only Socio-Dem Covars:", round(logit.onlycovars.roc$auc,2 ),
                                      " ± ",round(sqrt(pROC::var(logit.onlycovars.roc)),3), sep=""), 
                                paste("SSR + Age + Gender: ", round(logit.covars.roc$auc, 2),
                                      " ± ",round(sqrt(pROC::var(logit.covars.roc)),3), sep=""),
                                paste("SSR + All Socio-Dem Covars: ", round(logit.all.roc$auc, 2),
                                      " ± ",round(sqrt(pROC::var(logit.all.roc)),3), sep="")
                                )
                       ) +
  ggtitle("Non-Kin Different House Relationships (reciprocated ties)")

roc.plot.nkh.reciprocated <- roc.plot.nkh.reciprocated +
  theme_pubr(legend =  c(0.6, 0.2)) +
  theme(plot.title = element_text(size = 10, hjust = .5)) +
  labs_pubr()

roc.plot.nkh.reciprocated <- roc.plot.nkh.reciprocated + theme(
  plot.title = element_text(size = 10, hjust = .5),
  legend.title = element_text(size = 8),
  legend.text = element_text(size = 8)
)

svglite('Figures/plot_roc_non_kin_different_house_reciprocated.svg')
roc.plot.nkh.reciprocated
dev.off()

rocs_1 <-ggarrange(roc.plot.reciprocated, roc.plot.nkh.reciprocated, nrow = 2, labels = c("A", "B"), widths = 1)
svglite('Figures/plot_allroc_reciprocated_ties.svg', fix_text_size = FALSE)
rocs_1
dev.off()
```


#Predicted network of all social or familial relationships
```{r}
metrics_all_predictions_reciprocated_ties <- Downsamp_all_reciprocated_ties |> 
  inner_join(village_map, by = join_by(village_code_w3==village_code)) |> 
  group_by(village_name_deid) |>  
  summarize( `Percent predicted` = mean(predicted),
             `Percent real` = mean(related),
             TP = sum(predicted == 1 & related == 1 ),
             TN = sum(predicted == 0 & related == 0 ),
             FP = sum(predicted == 1 & related == 0 ),
             FN = sum(predicted == 0 & related == 1 ),
             TPR = TP/(TP+FN),
             FNR = FN/(TP+FN),
             TNR = TN/(TN+FP),
             Precision = TP/(TP+FP),
             Accuracy = (TP+TN)/(TP+TN+FP+FN),
             F1 = 2*TP / (2*TP + FP + FN)
          ) |> 
  rename(`Village name` = village_name_deid)
metrics_all_predictions_reciprocated_ties
```
```{r}
write_tsv(metrics_all_predictions_reciprocated_ties, 'tables/metrics_all_relationships_network_villages_prediction_reciprocated_ties.tsv')
writexl::write_xlsx(metrics_all_predictions_reciprocated_ties, 'tables/metrics_all_relationships_network_villages_prediction_reciprocated_ties.xlsx')
```

#Predicted network of non-kin relationships
```{r}
metrics_non_kin_prediction_reciprocated_ties <- Downsamp_all_reciprocated_ties_nkh |> 
  inner_join(village_map, by = join_by(village_code_w3==village_code)) |> 
  group_by(village_name_deid) |>  
  summarize( `Percent predicted` = mean(predicted),
             `Percent real` = mean(related),
             TP = sum(predicted == 1 & related == 1 ),
             TN = sum(predicted == 0 & related == 0 ),
             FP = sum(predicted == 1 & related == 0 ),
             FN = sum(predicted == 0 & related == 1 ),
             TPR = TP/(TP+FN),
             FNR = FN/(TP+FN),
             TNR = TN/(TN+FP),
             Precision = TP/(TP+FP),
             Accuracy = (TP+TN)/(TP+TN+FP+FN),
             F1 = 2*TP / (2*TP + FP + FN)
          ) |> 
      rename(`Village name` = village_name_deid)

metrics_non_kin_prediction_reciprocated_ties
```

```{r}
write_tsv(metrics_non_kin_prediction_reciprocated_ties, 'tables/metrics_non_kin_network_villages_prediction_reciprocated_ties.tsv')
writexl::write_xlsx(metrics_non_kin_prediction_reciprocated_ties, 'tables/metrics_non_kin_network_villages_prediction_reciprocated_ties.xlsx')
```