Ego-net deselection Analysis: Robustness analysis

Goal

Replicate robustness analyses

Set up

Packages

#library
library(tidyverse)
library(data.table)
library(kableExtra)
library(patchwork)
library(lme4)
library(viridis)
library(sjPlot)
library(foreach)
library(doParallel)
library(furrr) #for parallel computing
library(future) #for parallel computing

Import data

Import data and data prep

#load prepared data.
load("datafiles/data-processed/disaggregated_data/2022-06-13_dyad-survival-data-imputed.rda")

MyData <-  nonkin_survival_data_lead_dependent_imputed

#change scientific notation
options(scipen = 999)

#rename dropped_lead into dropped, so we can easily use this code.
MyData <- MyData %>% 
  rename(dropped = dropped_lead)

#last dataprep
MyData <- MyData %>%
  mutate(
    educ_ego = educ_ego - 4,
    age = leeftijd - 15,
    age_sq = age * age,
    age_alter = age_alter - 1,
    educ_alter = educ_alter - 4,
    length_rel_member = length_rel_member - 1,
    length_rel_total = length_rel_total - 1,
    size = size - 1,
    ei_alter_gender_rev = (2 - (MyData$ei_alter_gender + 1)) - 1,
    ei_alter_ethnicity_rev = (2 - (MyData$ei_alter_ethnicity + 1)) - 1,
    origin_rec_nar = ifelse(is.na(origin_rec_nar), 3, origin_rec_nar),
    origin_rec_nar_fac = factor(
      origin_rec_nar,
      levels = 0:3,
      labels = c("None", "Non-western", "Western", "Missing")
    ),
    origin_alter_rec = ifelse(is.na(origin_alter_rec), 3, origin_alter_rec),
    origin_alter_rec_fac = factor(
      origin_alter_rec,
      levels = 0:3,
      labels = c("None", "Non-western", "Western", "Missing")
    ),
    length = ifelse(is.na(length), 4, length),
    length_fac = factor(
      length,
      levels = 1:4,
      labels = c("< 3 years",
                 "3 - 6 years",
                 "> 6 years", 
                 "Length missing")
      )
    )

#scale variables for me models
MyData <- MyData %>% 
  mutate(avsim_alter_educ_cen = scale(avsim_alter_educ),
         avsim_alter_age_cen = scale(avsim_alter_age),
         ei_alter_gender_rev_cen = scale(ei_alter_gender_rev),
         ei_alter_ethnicity_rev_cen = scale(ei_alter_ethnicity_rev),
         length_rel_member_cen = scale(length_rel_member),
         length_rel_total_cen = scale(length_rel_total),
         size_cen = scale(size),
         degree_cen = scale(degree))

Custom functions

#parameter name function
relabel_f <- function(x) {
  df <- x %>% 
    mutate(parameter = case_when(
      effect_name == "dyad_educ_sim_cen" ~ 1,
      effect_name == "dyad_age_sim_cen" ~ 2,
      effect_name == "dyad_gender_sim_cen" ~ 3,
      effect_name == "dyad_ethnicity_sim_cen" ~ 4,
      effect_name == "combined_sim_cen" ~ 5,
      effect_name == "avsim_alter_educ_cen" ~ 6,
      effect_name == "avsim_alter_age_cen" ~ 7,
      effect_name == "ei_alter_gender_rev_cen" ~ 8,
      effect_name == "ei_alter_ethnicity_rev_cen" ~ 9,
      effect_name == "avsim_alter_combined_cen" ~ 10,
      effect_name == "as.factor(dear_alter_rec)1" ~ 11,
      effect_name == "as.factor(dear_alter_rec)2" ~ 12,
      effect_name == "degree_cen" ~ 13,
      effect_name == "avsim_alter_educ_cen:dyad_educ_sim_cen" ~ 14,
      effect_name == "avsim_alter_age_cen:dyad_age_sim_cen" ~ 15,
      effect_name == "ei_alter_gender_rev_cen:dyad_gender_sim_cen" ~ 16,
      effect_name == "ei_alter_ethnicity_rev_cen:dyad_ethnicity_sim_cen" ~ 17,
      effect_name == "combined_sim_cen:avsim_alter_combined_cen" ~ 18,
      effect_name == "time:dyad_educ_sim_cen" ~ 19,
      effect_name == "time:dyad_age_sim_cen" ~ 20,
      effect_name == "time:dyad_gender_sim_cen" ~ 21,
      effect_name == "time:dyad_ethnicity_sim_cen" ~ 22,
      effect_name == "time:as.factor(dear_alter_rec)1" ~ 23,
      effect_name == "time:as.factor(dear_alter_rec)2" ~ 24,
      effect_name == "time:degree_cen" ~ 25,
      effect_name == "gender_facFemale:dyad_gender_sim_cen" ~ 26,
      effect_name == "educ_ego_cen:dyad_educ_sim_cen" ~ 27,
      effect_name == "age_cen:dyad_age_sim_cen" ~ 28,
      effect_name == "origin_rec_nar_facNon-western:dyad_ethnicity_sim_cen" ~ 29,
      effect_name == "origin_rec_nar_facWestern:dyad_ethnicity_sim_cen" ~ 30,
      effect_name == "(Intercept)" ~ 31,
      effect_name == "time" ~ 32,
      effect_name == "educ_ego_cen" ~ 33,
      effect_name == "age_cen" ~ 34,
      effect_name == "gender_facFemale" ~ 35,
      effect_name == "gender_facMissing" ~ 36,
      effect_name == "origin_rec_nar_facNon-western" ~ 37,
      effect_name == "origin_rec_nar_facWestern" ~ 38,
      effect_name == "origin_rec_nar_facMissing" ~ 39,
      effect_name == "divorced_facdivorced" ~ 40,
      effect_name == "divorced_facmissing" ~ 41,
      effect_name == "moving_facnew_residence" ~ 42,
      effect_name == "moving_facnew_municipality" ~ 43,
      effect_name == "moving_facmissing" ~ 44,
      effect_name == "first_child_facFirst child born" ~ 45,
      effect_name == "first_child_facMissing" ~ 46,
      effect_name == "educ_alter_cen" ~ 47,
      effect_name == "age_alter_cen" ~ 48,
      effect_name == "gender_alter_facFemale" ~ 49,
      effect_name == "gender_alter_facMissing" ~ 50,
      effect_name == "origin_alter_rec_facNon-western" ~ 51,
      effect_name == "origin_alter_rec_facWestern" ~ 52,
      effect_name == "origin_alter_rec_facMissing" ~ 53,
      effect_name == "rel_alter_recSame group or club" ~ 54,
      effect_name == "rel_alter_recNeighbour" ~ 55,
      effect_name == "rel_alter_recFriend" ~ 56,
      effect_name == "rel_alter_recAdvisor" ~ 57,
      effect_name == "rel_alter_recOther" ~ 58,
      effect_name == "rel_alter_recMissing" ~ 59,
      effect_name == "times_dropped_earlier_cen" ~ 60,
      effect_name == "length_fac3 - 6 years" ~ 61,
      effect_name == "length_fac> 6 years" ~ 62,
      effect_name == "length_facLength missing" ~ 63,
      effect_name == "size_cen" ~ 64,
      effect_name == "net_density_cen" ~ 65,
      effect_name == "size_net_one" ~ 66
    ),
    parameter_fac = factor(
      parameter,
      levels = 1:66,
      labels = c(
      "Dyadic similarity: education",
      "Dyadic similarity: age",
      "Dyadic similarity: gender",
      "Dyadic similarity: ethnicity",
      "Dyadic similarity: combined",
      "Confidant uniqueness: education",
      "Confidant uniqueness: age",
      "Confidant uniqueness: gender",
      "Confidant uniqueness: ethnicity",
      "Confidant uniqueness: combined",
      "Closeness: dear",
      "Closeness: not asked",
      "Alter Embeddedness",
      "Interaction education",
      "Interaction age",
      "Interaction gender",
      "Interaction ethnicity",
      "Interaction combined",
      "Interaction time education",
      "Interaction time age",
      "Interaction time gender",
      "Interaction time ethnicity",
      "Interaction time dear",
      "Interaction time dear not asked",
      "Interaction time degree",
      "Interaction gender X sim",
      "Interaction educ X sim",
      "Interaction age X sim",
      "Interaction non-western X sim",
      "Interaction western X sim",
      "Constant",
      "Linear time",
      "Education",
      "Age",
      "Gender Female",
      "Gender Missing",
      "Non-western migration background",
      "Western migration background",
      "Migration missing",
      "Divorced",
      "Divorced Missing",
      "Moving: new residence",
      "Moving: new municipality",
      "Moving Missing",
      "First child born",
      "First child born Missing",
      "Education alter",
      "Age alter",
      "Alter female",
      "Alter gender missing",
      "Non-western migration background Alter",
      "Western migration background Alter",
      "Migration alter missing",
      "Same group or club",
      "Neighbour",
      "Friend",
      "Advisor",
      "Other",
      "Rel missing",
      "Times dropped earlier",
      "Length 3-6 years",
      "Length > 6 years",
      "Length missing",
      "Network size",
      "Network density",
      "Size == one"
      )
    ))
}

# x <- naq_0_model_list[[10]]
# model_label <- "main_model"
# model_number <- 10

#extract output from a lme4 model
extract_model_output <- function(x,
                                 model_label,
                                 model_number){


df_model <- summary(x)$coefficients %>% 
  as.tibble() %>% 
  mutate(effect_name = attributes(summary(x)$coefficients)$dimnames[[1]]) %>% 
  select(5,1:4) %>% 
  relabel_f(x = .) %>% 
  mutate(m_label = model_label,
         m_number = model_number) %>% 
  rename(p = `Pr(>|z|)`,
         estimate = Estimate,
         SE = `Std. Error`) %>% 
  mutate(significant = case_when(
             p < 0.05 & p > 0.01 ~ "*",
             p < 0.01 & p > 0.001 ~ "**",
             p < 0.001 ~ "***",
             p > 0.05 ~ ""),
         min = estimate - SE*1.96,
         max = estimate + SE*1.96,
         riskratio = exp(estimate),
         rr_min = exp(min),
         rr_max = exp(max),
         model = factor(
           m_number,
           levels = 1:13,
           labels = c(
             "Model 0",
             "Model 1",
             "Model 2",
             "Model 3",
             "Model 4",
             "Model 5",
             "Model 6",
             "Model 7",
             "Model 8",
             "Model 9",
             "Model 10",
             "Model 11",
             "Model 12"
             )
           )
         ) 
}

Robustness analysis: Demographic differences

Create model directories

#create sub file
if(!dir.exists("results/survival_results/robustness/demo_interactions")){
  dir.create("results/survival_results/robustness/demo_interactions")
}

Create model formulas

#set empty list
formula_list <- list()

#create formula's

## Gender differences
formula_list[[1]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
       time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    gender_fac:dyad_gender_sim_cen +
    size_net_one
)

## Education
formula_list[[2]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
       time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    educ_ego_cen:dyad_educ_sim_cen +
    size_net_one
)

## Age
formula_list[[3]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
       time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    age_cen:dyad_age_sim_cen +
    size_net_one
)


## Age
formula_list[[4]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
       time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    origin_rec_nar_fac:dyad_ethnicity_sim_cen +
    size_net_one
)

#create names
names(formula_list) <- foreach(i = 1:4,
                     .combine = c) %do% {
                       paste0("model_0", i)
                     }

Estimate models

#set up of file name and dir
#set file name to use
dir <- "results/survival_results/robustness/demo_interactions/" 
file_name_model <- paste0(dir, "model_demo-int_")

#initiate do parallel session
registerDoParallel(cores = 4)

nr_list <-  c("01",
              "02",
              "03",
              "04")

#create foreach loop to use parallel computation to run fixed effects models
foreach(i = 1:4,
        .final = function(x) NULL) %do% {
  #i = 8
  #create file_name
  file_name <- paste0(file_name_model, nr_list[i], ".rda")
  
  #check whether filename exists
  if (!file.exists(file_name)) {
    model_results <- glmer(
      formula = formula_list[[i]],
      data = MyData,
      family = binomial(link = "cloglog"),
      nAGQ = 0,
      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun =
                                                                    200000))
    )
    
    save(model_results,
         file = file_name)
  }
        }

NULL

registerDoSEQ()

Import model results

#import demo interaction models
dir <- "results/survival_results/robustness/demo_interactions/"

demo_int_model_list <- list.files(dir,
           full.names = T) %>% 
  map(.x = ., 
      .f = ~get(load(file = .x)))

Extract information from model object

#extract info and store in df
demo_int_model_df <- foreach(a = 1:length(demo_int_model_list),
        .combine = rbind) %do%
  extract_model_output(x = demo_int_model_list[[a]],
                       model_label = "naq_0",
                       model_number = a)

Create table df

demo_int_table <- demo_int_model_df %>% 
  select(parameter_fac, riskratio, `z value`, significant, m_number) %>% 
  rename(statistic = `z value`) %>% 
  mutate(riskratio = round(riskratio,3),
         statistic = round(statistic,3)) %>% 
  pivot_wider(names_from = m_number,
              values_from = c("riskratio", "statistic", "significant")) %>% 
  select(parameter_fac,
         riskratio_1,  significant_1, statistic_1,
         riskratio_2,  significant_2, statistic_2,
         riskratio_3,  significant_3, statistic_3,
         riskratio_4,  significant_4, statistic_4) %>% 
  filter(parameter_fac %in% c(
    "Constant",
    "Linear time",
      "Education",
      "Age",
      "Gender Female",
      "Non-western migration background",
      "Western migration background",
    "Interaction gender X sim",
      "Interaction educ X sim",
      "Interaction age X sim",
      "Interaction non-western X sim",
      "Interaction western X sim"
  ) |
    str_detect(parameter_fac, "Dyadic similarity")) %>% 
  arrange(as.numeric(parameter_fac))

Create table with tab-model

tab_model(demo_int_model_list[[1]],
          demo_int_model_list[[2]],
          demo_int_model_list[[3]],
          demo_int_model_list[[4]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped		dropped		dropped		dropped
Predictors	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic
(Intercept)	3.000 ***	18.662	2.929 ***	29.300	2.930 ***	29.311	2.963 ***	28.703
time	0.755 ***	-35.068	0.755 ***	-35.078	0.755 ***	-35.068	0.755 ***	-35.056
avsim_alter_educ_cen	1.012	1.685	1.012	1.657	1.012	1.682	1.012	1.678
avsim_alter_age_cen	1.019 *	2.427	1.019 *	2.422	1.019 *	2.433	1.019 *	2.437
ei_alter_gender_rev_cen	1.009	1.234	1.009	1.226	1.009	1.232	1.009	1.236
ei_alter_ethnicity_rev_cen	1.018	1.789	1.018	1.799	1.018	1.786	1.018	1.782
educ_ego_cen	0.937 ***	-6.957	0.942 ***	-5.614	0.937 ***	-6.954	0.937 ***	-6.942
age_cen	1.031	1.844	1.031	1.869	1.033	1.697	1.031	1.858
gender_fac: Female	0.974	-0.150	0.894 ***	-5.273	0.893 ***	-5.291	0.893 ***	-5.300
gender_fac: Missing	0.710	-1.955	0.677 **	-2.597	0.676 **	-2.603	0.678 **	-2.584
origin_rec_nar_fac: Non-western	1.243 ***	4.923	1.243 ***	4.916	1.243 ***	4.922	1.300 ***	3.768
origin_rec_nar_fac: Western	1.003	0.082	1.003	0.064	1.003	0.080	1.091	1.007
origin_rec_nar_fac: Missing	1.014	0.418	1.013	0.384	1.014	0.413	1.009	0.251
divorced_fac: divorced	0.938	-1.400	0.938	-1.399	0.938	-1.397	0.938	-1.399
divorced_fac: missing	0.966 *	-2.069	0.966 *	-2.076	0.966 *	-2.066	0.966 *	-2.061
moving_fac: new_residence	1.069	1.310	1.069	1.313	1.069	1.310	1.069	1.313
moving_fac: new_municipality	1.119	1.606	1.119	1.603	1.119	1.605	1.119	1.598
moving_fac: missing	0.986	-0.813	0.986	-0.822	0.986	-0.813	0.986	-0.806
first_child_fac: First child born	1.106	0.831	1.105	0.824	1.106	0.831	1.106	0.832
first_child_fac: Missing	1.011	0.541	1.011	0.523	1.012	0.549	1.012	0.550
educ_alter_cen	0.997	-0.390	0.991	-0.904	0.997	-0.399	0.997	-0.396
age_alter_cen	1.008	0.561	1.008	0.568	1.006	0.329	1.008	0.568
gender_alter_fac: Female	0.848	-0.624	0.969	-1.787	0.969	-1.786	0.969	-1.769
gender_alter_fac: Missing	1.611	1.387	1.745	1.786	1.745	1.784	1.753	1.802
origin_alter_rec_fac: Non-western	1.112 **	2.703	1.112 **	2.714	1.112 **	2.701	1.047	0.640
origin_alter_rec_fac: Western	1.075	1.541	1.076	1.552	1.075	1.539	0.996	-0.048
origin_alter_rec_fac: Missing	1.434 **	2.887	1.436 **	2.898	1.433 **	2.883	1.408 **	2.714
as.factor(dear_alter_rec)1	0.695 ***	-14.723	0.695 ***	-14.717	0.695 ***	-14.719	0.695 ***	-14.722
as.factor(dear_alter_rec)2	0.834 ***	-10.518	0.834 ***	-10.515	0.834 ***	-10.519	0.834 ***	-10.509
rel_alter_rec: Same group or club	0.985	-0.391	0.986	-0.373	0.986	-0.375	0.985	-0.402
rel_alter_rec: Neighbour	0.943	-1.677	0.944	-1.660	0.944	-1.661	0.942	-1.699
rel_alter_rec: Friend	0.819 ***	-8.261	0.820 ***	-8.240	0.820 ***	-8.232	0.819 ***	-8.290
rel_alter_rec: Advisor	1.229 **	2.961	1.230 **	2.971	1.230 **	2.966	1.228 **	2.947
rel_alter_rec: Other	1.130 **	3.032	1.131 **	3.044	1.130 **	3.033	1.129 **	3.001
rel_alter_rec: Missing	0.607 *	-2.214	0.606 *	-2.225	0.608 *	-2.206	0.603 *	-2.245
times_dropped_earlier_cen	0.896 ***	-15.687	0.895 ***	-15.693	0.895 ***	-15.689	0.895 ***	-15.702
length_fac: 3 - 6 years	0.894 ***	-4.714	0.894 ***	-4.707	0.894 ***	-4.697	0.894 ***	-4.702
length_fac: > 6 years	0.868 ***	-6.193	0.868 ***	-6.202	0.868 ***	-6.167	0.868 ***	-6.188
length_fac: Length missing	1.204	1.919	1.205	1.927	1.204	1.921	1.204	1.918
degree_cen	0.919 ***	-6.736	0.919 ***	-6.742	0.919 ***	-6.734	0.920 ***	-6.728
dyad_educ_sim_cen	0.992	-1.143	0.993	-0.974	0.992	-1.139	0.992	-1.137
dyad_age_sim_cen	0.964 ***	-4.726	0.964 ***	-4.743	0.964 ***	-4.630	0.964 ***	-4.724
dyad_gender_sim_cen	0.871	-1.374	0.916 ***	-12.376	0.916 ***	-12.374	0.916 ***	-12.386
dyad_ethnicity_sim_cen	0.987	-1.137	0.987	-1.144	0.987	-1.137	0.961	-1.556
size_cen	1.067 ***	5.883	1.067 ***	5.858	1.067 ***	5.882	1.067 ***	5.874
net_density_cen	1.016	1.420	1.016	1.421	1.015	1.419	1.015	1.412
size_net_one	1.331 ***	5.521	1.330 ***	5.508	1.330 ***	5.516	1.329 ***	5.501
gender_facFemale:dyad_gender_sim_cen	1.107	0.505
educ_ego_cen:dyad_educ_sim_cen			1.007	0.917
age_cen:dyad_age_sim_cen					1.002	0.175
origin_rec_nar_facNon-western:dyad_ethnicity_sim_cen							1.048	0.951
origin_rec_nar_facWestern:dyad_ethnicity_sim_cen							1.066	1.205
Random Effects
σ2	1.64		1.64		1.64		1.64
τ00	0.19 nomem_encr		0.19 nomem_encr		0.19 nomem_encr		0.19 nomem_encr
ICC	0.10		0.10		0.10		0.10
N	6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr
Observations	49449		49449		49449		49449
Marginal R2 / Conditional R2	0.116 / 0.208		0.116 / 0.207		0.116 / 0.207		0.116 / 0.208
Deviance	56995.141		56994.518		56995.352		56993.962
log-Likelihood	-28497.570		-28497.259		-28497.676		-28496.981
p<0.05   ** p<0.01   *** p<0.001

Robustness analysis: Combined dyadic similarity

Data preperation

Create a new combined dyadics similarity measure

#create combined similarity measure
MyData <- MyData %>% 
  mutate(combined_sim = (dyad_age_sim + dyad_educ_sim + dyad_gender_sim + dyad_ethnicity_sim)/4,
         combined_sim_cen = scale(combined_sim))

#create uniqueness for combined sim estimate
combined_unique_df <- MyData %>% 
  select(nomem_encr, 
         dyad_id, 
         survey_wave, 
         age_alter, 
         gender_alter_fac,
         educ_alter,
         origin_alter_rec_fac)

combined_unique_list <- combined_unique_df %>% 
  group_split(nomem_encr, survey_wave)

#set range and variables for test function
#ranges <- c(12,12,1,2)
#variable <- c("age", "educ", "gender", "origin")

f_create_combined_sim <- function(df, ranges, variable) {
  #df = test
  
  #create list to store information
  net_results_list <- list()
  
  #create for every alter avsim alter for age, gender, education, and gender
  for (i in 1:nrow(df)) {
    #i= 1
    #select alter and remove alter from net
    alter <- df[i, ]
    net <- df[-i, ]
    
    #create lists
    net_list <- list()
    alter_list <- list()
    
    #extract age
    net_list[[1]] <- as.vector(t(net$age_alter))
    alter_list[[1]] <- alter$age_alter
    
    #extract educ
    net_list[[2]] <- as.vector(t(net$educ_alter))
    alter_list[[2]] <- alter$educ_alter
    
    #extract gender
    net_list[[3]] <- as.vector(t(net$gender_alter_fac))
    alter_list[[3]] <- alter$gender_alter_fac
    
    #extract origin
    net_list[[4]] <- as.vector(t(net$origin_alter_rec_fac))
    alter_list[[4]] <- alter$origin_alter_rec_fac
    
    #results list
    result_list <- list()
    
    #start inner loop
    for (j in 1:4) {
      #j = 2
      
      if (j < 3) {
        #for con sim
        result_list[[j]] <-
          tibble(alter = alter_list[[j]], net = net_list[[j]]) %>%
          filter(!is.na(net) & !is.na(alter)) %>% #filter out missings
          mutate(
            dyad_sim = 1 - (abs(alter - net) / ranges[[j]])
            #dyadic sim
          ) %>% 
          select(dyad_sim) %>% 
          rename(!!paste0("dyad_sim", "_", variable[j]) := dyad_sim)
        
      } else {
        #for cat sim
        result_list[[j]] <-
          tibble(alter = alter_list[[j]], net = net_list[[j]]) %>%
          filter(!is.na(net) & !is.na(alter)) %>% #filter out missings
          mutate(
            dyad_sim = as.numeric(alter == net)
            #dyadic sim
          ) %>%
          select(dyad_sim) %>% 
          rename(!!paste0("dyad_sim", "_", variable[j]) := dyad_sim) 
      }
    }
    
    #combine estimates and create combined sim
    result_df <- result_list %>% 
      bind_cols() %>% 
      mutate(
        combined_sim = (dyad_sim_age + dyad_sim_educ + dyad_sim_gender + dyad_sim_origin)/4,
        dyad_id =   df[i, ] %>%
          pull(dyad_id),
        nomem_encr =   df[i, ] %>%
          pull(nomem_encr),
        survey_wave = df[i, ] %>% 
           pull(survey_wave),
        avsim_alter_combined = mean(combined_sim)
      ) %>% 
      select(nomem_encr, dyad_id, survey_wave, avsim_alter_combined) %>% 
      distinct()
  
    #store results of alter
    net_results_list[[i]] <- result_df
    
  }
  
  #store in df
  results_df <- net_results_list %>% 
    bind_rows()
  
  return(results_df)
}

ranges_c <- c(12, 12, 1, 2)
variables_c <- c("age", "educ", "gender", "origin")

#use future_map and the f_make_net_variables_df
#start parallel session
plan(multisession, workers = 7)

combined_unique_result_list <- combined_unique_list %>%
  future_map(
    .x = .,
    .f = ~ f_create_combined_sim(
      df = .x,
      variable = variables_c,
      ranges = ranges_c
    ),
    .progress = T
  )

plan(sequential)



combined_unique_sim_df <- combined_unique_result_list %>% bind_rows()

MyData <- MyData %>% 
  left_join(combined_unique_sim_df, by = c("nomem_encr", "dyad_id", "survey_wave"))

#create scale
MyData <- MyData %>% 
  mutate(avsim_alter_combined_cen = scale(avsim_alter_combined))

#mean imputation
MyData <- MyData %>% 
  mutate(avsim_alter_combined = ifelse(is.na(avsim_alter_combined),
                                       mean(avsim_alter_combined, na.rm = T),
                                       avsim_alter_combined),
         avsim_alter_combined_cen = scale(avsim_alter_combined))

Model estimation

Create directories

#create sub file
if(!dir.exists("results/survival_results/robustness/combined_sim")){
  dir.create("results/survival_results/robustness/combined_sim")
}

Create model formulas

combined_sim_formula_list <- list()

combined_sim_formula_list[[1]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time 
)

combined_sim_formula_list[[2]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen
)

combined_sim_formula_list[[3]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_combined_cen +
    size_net_one
)

combined_sim_formula_list[[4]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    size_net_one
)

combined_sim_formula_list[[5]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    size_cen +
    net_density_cen +
    size_net_one
    
)

combined_sim_formula_list[[6]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    size_cen +
    net_density_cen +
    size_net_one
    
)

combined_sim_formula_list[[7]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    size_cen +
    net_density_cen +
    size_net_one
)

combined_sim_formula_list[[8]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    size_cen +
    net_density_cen +
    size_net_one
)

combined_sim_formula_list[[9]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    size_cen +
    net_density_cen +
    size_net_one +
    combined_sim_cen:avsim_alter_combined_cen
)

combined_sim_formula_list[[10]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    combined_sim_cen +
    avsim_alter_combined_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    size_cen +
    net_density_cen +
    size_net_one +
    combined_sim_cen:time
)


#create names
names(combined_sim_formula_list) <- foreach(i = 1:10,
                     .combine = c) %do% {
                       paste0("model_0", i)
                     }

Estimate models

#set up of file name and dir
#set file name to use
dir <- "results/survival_results/robustness/combined_sim/"
file_name_model <- paste0(dir, "model_com-sim_")


#initiate do parallel session
registerDoParallel(cores = 6)

nr_list <-  c("01",
              "02",
              "03",
              "04",
              "05",
              "06",
              "07",
              "08",
              "09",
              "10")

#create foreach loop to use parallel computation to run fixed effects models
foreach(i = 1:10,
        .final = function(x) NULL) %do% {
  #i = 8
  #create file_name
  file_name <- paste0(file_name_model, nr_list[i], ".rda")
  
  #check whether filename exists
  if (!file.exists(file_name)) {
    model_results <- glmer(
      formula = combined_sim_formula_list[[i]],
      data = MyData,
      family = binomial(link = "cloglog"),
      nAGQ = 0,
      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun =
                                                                    200000))
    )
    
    save(model_results,
         file = file_name)
  }
        }

NULL

registerDoSEQ()

Import results

#import demo interaction models
dir <- "results/survival_results/robustness/combined_sim/"

com_sim_model_list <- list.files(dir,
           full.names = T) %>% 
  map(.x = ., 
      .f = ~get(load(file = .x)))

Extract results

#extract info and store in df
com_sim_model_df <- foreach(a = 1:length(com_sim_model_list),
        .combine = rbind) %do%
  extract_model_output(x = com_sim_model_list[[a]],
                       model_label = "com_sim",
                       model_number = a)

Create coefficient plot with combined sim estimate

combined_coefplot <- com_sim_model_df %>%
  filter(str_detect(parameter_fac, "Confidant uniqueness: combined") |
           str_detect(parameter_fac, "Dyadic similarity: combined")) %>%
  filter(model %in% c("Model 1", "Model 2", "Model 3","Model 4", "Model 5", "Model 6")) %>% 
  mutate(
    controls = ifelse(model %in% c("Model 1", "Model 2"), 1, NA),
    controls = ifelse(
      model %in% c("Model 3"),
      2,
      controls
    ),
    controls = ifelse(
      model %in% c("Model 4"),
      3,
      controls
    ),
    controls = ifelse(
      model %in% c("Model 5"),
      4,
      controls
    ),
    controls = ifelse(
      model %in% c("Model 6"),
      5,
      controls
    ),
    controls = factor(
      controls,
      levels = 1:5,
      labels = c(
        "Time trends (M2 and M3)",
        "Dyadic similarity or confidant homogeneity (M4)",
        "Ego, confidant, dyad, and network controls (M5)",
        "Emotional closeness (M6)",
        "Alter embeddedness (M7)"
        )
    ),
    sim_type = case_when(
      str_detect(parameter_fac, "Confidant") ~ 2,
      str_detect(parameter_fac, "Dyadic") ~ 1
    ),
    sim_type = factor(
      sim_type,
      levels = 1:2,
      labels = c("Dyadic similarity", "Confidant uniqueness")
    ),
    variable = case_when(
      str_detect(parameter_fac, "age") ~ 1,
      str_detect(parameter_fac, "education") ~ 2,
      str_detect(parameter_fac, "gender") ~ 3,
      str_detect(parameter_fac, "ethnicity") ~ 4
    ),
    variable = factor(
      variable,
      levels = 1:4,
      labels = c("Age", "Education", "Gender", "Ethnicity")
    )
  ) %>%
  ggplot(aes(x = variable, y = riskratio, shape = controls, group = controls)) +
  geom_hline(yintercept = 1) +
  geom_pointrange(aes(ymin = rr_min, ymax = rr_max),
                  position = position_dodge(width = 0.7)) +
  geom_point(aes(colour = controls,
                 fill = controls),
             size = 1.2,
             position = position_dodge(width = 0.7)) +
  facet_wrap(vars(sim_type)) +
  scale_color_manual(values = c("#1b9e77",
                                "#d95f02",
                                "#7570b3",
                                "#e7298a",
                                "#66a61e")
                     ) +
  scale_fill_manual(values = c("#1b9e77",
                                "#d95f02",
                                "#7570b3",
                                "#e7298a",
                                "#66a61e")
                     ) +
  scale_y_continuous(minor_breaks = c(0.95, 1.05),
                     limits = c(0.85,1.15)) +
  scale_shape_manual(values = c(21,22,23,24,25)) +
  theme(panel.background = element_rect(fill = "#FFFFFF",
                                        colour = "black"),
        plot.background = element_rect(fill = "#FFFFFF"),
        panel.grid = element_line(colour = "grey"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_line(colour = "grey"),
        text = element_text(family = "sans", size = 12), 
        axis.title.x = element_text(hjust = 0.9,face = "bold"),
        axis.text.x = element_blank(),
        axis.line = element_blank(),
        axis.title.y = element_text(hjust = 0.9, face = "bold"),
        axis.ticks = element_blank(),
        strip.background = element_rect(fill = "#A9A9A9",
                                        colour = "black"),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.text = element_text(family = "sans", size = 8), 
        legend.title = element_text(family = "sans", size = 8),
        legend.background = element_rect(fill = "#FFFFFF"),
        legend.key = element_rect(fill = "#FFFFFF")) +
  guides(fill=guide_legend(ncol=2),
         group=guide_legend(ncol=2),
         colour=guide_legend(ncol=2),
         shape=guide_legend(ncol=2)) + 
  labs(shape = "Controlled for:", 
       group = "Controlled for:",
       colour = "Controlled for:",
       fill = "Controlled for:",
       y = "Risk Ratio",
       x = "")

combined_coefplot

ggsave(combined_coefplot,
       file = "plots/results/survival/combined_coefplot_robustness_com_sim.jpg",
       dpi = 320,
       width = 7,
       height = 5)

Table combined sim measure

tab_model(com_sim_model_list[[2]],
          com_sim_model_list[[3]],
          com_sim_model_list[[4]],
          com_sim_model_list[[5]],
          com_sim_model_list[[6]],
          com_sim_model_list[[7]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped		dropped		dropped		dropped		dropped		dropped
Predictors	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic
(Intercept)	1.863 ***	46.863	1.871 ***	47.214	1.857 ***	46.526	2.901 ***	31.002	3.147 ***	32.366	2.749 ***	28.858
time	0.735 ***	-39.520	0.730 ***	-40.421	0.734 ***	-39.621	0.747 ***	-36.460	0.752 ***	-35.639	0.749 ***	-36.106
combined_sim_cen	0.885 ***	-17.682			0.891 ***	-15.024	0.916 ***	-10.282	0.913 ***	-10.613	0.915 ***	-10.337
avsim_alter_combined_cen			0.938 ***	-9.231	0.986	-1.825	1.003	0.368	1.001	0.066	1.004	0.445
size_net_one			1.121 **	2.770	1.118 **	2.711	1.124 *	2.505	1.149 **	2.953	1.330 ***	5.548
educ_ego_cen							0.944 ***	-6.227	0.940 ***	-6.662	0.943 ***	-6.361
age_cen							1.042 **	2.597	1.041 *	2.503	1.040 *	2.460
gender_fac: Female							0.891 ***	-5.484	0.889 ***	-5.547	0.890 ***	-5.551
gender_fac: Missing							0.682 *	-2.562	0.677 **	-2.597	0.683 *	-2.552
origin_rec_nar_fac: Non-western							1.186 ***	4.023	1.191 ***	4.087	1.185 ***	3.983
origin_rec_nar_fac: Western							0.911 **	-2.719	0.908 **	-2.782	0.910 **	-2.751
origin_rec_nar_fac: Missing							1.001	0.039	1.004	0.127	1.002	0.045
divorced_fac: divorced							0.922	-1.773	0.941	-1.329	0.924	-1.723
divorced_fac: missing							0.950 **	-3.069	0.964 *	-2.169	0.953 **	-2.862
moving_fac: new_residence							1.064	1.234	1.069	1.325	1.065	1.238
moving_fac: new_municipality							1.127	1.713	1.120	1.614	1.127	1.710
moving_fac: missing							0.986	-0.852	0.985	-0.905	0.987	-0.778
first_child_fac: First child born							1.102	0.807	1.099	0.779	1.111	0.873
first_child_fac: Missing							1.009	0.450	1.011	0.515	1.011	0.512
educ_alter_cen							1.000	0.028	1.000	-0.022	0.998	-0.207
age_alter_cen							1.009	0.612	1.005	0.367	1.006	0.443
gender_alter_fac: Female							0.968	-1.857	0.976	-1.396	0.966	-1.950
gender_alter_fac: Missing							1.787	1.862	1.808	1.897	1.738	1.771
origin_alter_rec_fac: Non-western							1.036	1.012	1.036	1.005	1.037	1.032
origin_alter_rec_fac: Western							0.976	-0.601	0.978	-0.550	0.975	-0.629
origin_alter_rec_fac: Missing							1.439 **	2.948	1.456 **	3.016	1.418 **	2.825
rel_alter_rec: Same group or club							0.941	-1.632	0.963	-1.007	0.965	-0.952
rel_alter_rec: Neighbour							0.907 **	-2.840	0.922 *	-2.330	0.927 *	-2.186
rel_alter_rec: Friend							0.764 ***	-11.565	0.798 ***	-9.556	0.784 ***	-10.355
rel_alter_rec: Advisor							1.230 **	3.006	1.222 **	2.888	1.237 **	3.087
rel_alter_rec: Other							1.102 *	2.442	1.124 **	2.903	1.123 **	2.897
rel_alter_rec: Missing							0.589 *	-2.368	0.590 *	-2.346	0.610 *	-2.209
times_dropped_earlier_cen							0.895 ***	-15.884	0.895 ***	-15.853	0.895 ***	-15.792
length_fac: 3 - 6 years							0.868 ***	-6.015	0.877 ***	-5.521	0.882 ***	-5.280
length_fac: > 6 years							0.824 ***	-8.639	0.844 ***	-7.532	0.847 ***	-7.351
length_fac: Length missing							1.145	1.407	1.176	1.682	1.174	1.665
size_cen							1.013	1.571	1.012	1.458	1.074 ***	6.518
net_density_cen							0.967 ***	-4.556	0.962 ***	-5.190	1.027 *	2.534
as.factor(dear_alter_rec)1									0.692 ***	-14.969
as.factor(dear_alter_rec)2									0.833 ***	-10.649
degree_cen											0.908 ***	-7.852
Random Effects
σ2	1.64		1.64		1.64		1.64		1.64		1.64
τ00	0.19 nomem_encr		0.18 nomem_encr		0.19 nomem_encr		0.18 nomem_encr		0.19 nomem_encr		0.18 nomem_encr
ICC	0.10		0.10		0.10		0.10		0.10		0.10
N	6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr
Observations	49449		49449		49449		49449		49449		49449
Marginal R2 / Conditional R2	0.084 / 0.178		0.080 / 0.172		0.085 / 0.178		0.111 / 0.199		0.115 / 0.205		0.112 / 0.200
Deviance	58474.835		58688.866		58464.269		57349.078		57102.626		57288.908
log-Likelihood	-29237.417		-29344.433		-29232.135		-28674.539		-28551.313		-28644.454
p<0.05   ** p<0.01   *** p<0.001

tab_model(com_sim_model_list[[9]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped
Predictors	Risk Ratios	Statistic
(Intercept)	2.972 ***	30.010
time	0.754 ***	-35.325
combined_sim_cen	0.923 ***	-8.434
avsim_alter_combined_cen	1.006	0.703
educ_ego_cen	0.939 ***	-6.799
age_cen	1.041 *	2.512
gender_fac: Female	0.885 ***	-5.748
gender_fac: Missing	0.679 *	-2.575
origin_rec_nar_fac: Non-western	1.210 ***	4.397
origin_rec_nar_fac: Western	0.925 *	-2.200
origin_rec_nar_fac: Missing	1.009	0.273
divorced_fac: divorced	0.944	-1.249
divorced_fac: missing	0.967 *	-1.973
moving_fac: new_residence	1.069	1.311
moving_fac: new_municipality	1.119	1.611
moving_fac: missing	0.985	-0.856
first_child_fac: First child born	1.105	0.823
first_child_fac: Missing	1.012	0.583
educ_alter_cen	0.998	-0.202
age_alter_cen	1.003	0.226
gender_alter_fac: Female	0.975	-1.447
gender_alter_fac: Missing	1.764	1.818
origin_alter_rec_fac: Non-western	1.023	0.631
origin_alter_rec_fac: Western	0.963	-0.902
origin_alter_rec_fac: Missing	1.437 **	2.899
as.factor(dear_alter_rec)1	0.699 ***	-14.510
as.factor(dear_alter_rec)2	0.839 ***	-10.182
rel_alter_rec: Same group or club	0.983	-0.465
rel_alter_rec: Neighbour	0.939	-1.804
rel_alter_rec: Friend	0.814 ***	-8.633
rel_alter_rec: Advisor	1.228 **	2.958
rel_alter_rec: Other	1.140 **	3.258
rel_alter_rec: Missing	0.605 *	-2.231
times_dropped_earlier_cen	0.895 ***	-15.802
length_fac: 3 - 6 years	0.890 ***	-4.887
length_fac: > 6 years	0.863 ***	-6.476
length_fac: Length missing	1.207	1.947
degree_cen	0.920 ***	-6.666
size_cen	1.063 ***	5.540
net_density_cen	1.013	1.219
size_net_one	1.331 ***	5.535
combined_sim_cen:avsim_alter_combined_cen	1.017 **	2.836
Random Effects
σ2	1.64
τ00 nomem_encr	0.19
ICC	0.10
N nomem_encr	6996
Observations	49449
Marginal R2 / Conditional R2	0.115 / 0.206
Deviance	57051.740
log-Likelihood	-28525.870
p<0.05   ** p<0.01   *** p<0.001

tab_model(com_sim_model_list[[10]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped
Predictors	Risk Ratios	Statistic
(Intercept)	3.002 ***	30.363
time	0.751 ***	-35.197
combined_sim_cen	0.890 ***	-8.424
avsim_alter_combined_cen	1.001	0.173
educ_ego_cen	0.939 ***	-6.741
age_cen	1.039 *	2.389
gender_fac: Female	0.889 ***	-5.578
gender_fac: Missing	0.678 **	-2.584
origin_rec_nar_fac: Non-western	1.188 ***	4.024
origin_rec_nar_fac: Western	0.909 **	-2.774
origin_rec_nar_fac: Missing	1.004	0.130
divorced_fac: divorced	0.941	-1.315
divorced_fac: missing	0.967 *	-2.004
moving_fac: new_residence	1.069	1.321
moving_fac: new_municipality	1.118	1.595
moving_fac: missing	0.985	-0.898
first_child_fac: First child born	1.106	0.834
first_child_fac: Missing	1.012	0.576
educ_alter_cen	0.998	-0.249
age_alter_cen	1.003	0.235
gender_alter_fac: Female	0.973	-1.522
gender_alter_fac: Missing	1.758	1.806
origin_alter_rec_fac: Non-western	1.037	1.032
origin_alter_rec_fac: Western	0.977	-0.571
origin_alter_rec_fac: Missing	1.439 **	2.913
as.factor(dear_alter_rec)1	0.700 ***	-14.478
as.factor(dear_alter_rec)2	0.840 ***	-10.158
rel_alter_rec: Same group or club	0.983	-0.460
rel_alter_rec: Neighbour	0.939	-1.804
rel_alter_rec: Friend	0.815 ***	-8.621
rel_alter_rec: Advisor	1.227 **	2.942
rel_alter_rec: Other	1.139 **	3.231
rel_alter_rec: Missing	0.611 *	-2.188
times_dropped_earlier_cen	0.895 ***	-15.768
length_fac: 3 - 6 years	0.891 ***	-4.859
length_fac: > 6 years	0.864 ***	-6.425
length_fac: Length missing	1.198	1.873
degree_cen	0.921 ***	-6.646
size_cen	1.064 ***	5.602
net_density_cen	1.013	1.220
size_net_one	1.323 ***	5.433
time:combined_sim_cen	1.020 *	2.284
Random Effects
σ2	1.64
τ00 nomem_encr	0.19
ICC	0.10
N nomem_encr	6996
Observations	49449
Marginal R2 / Conditional R2	0.115 / 0.205
Deviance	57053.631
log-Likelihood	-28526.815
p<0.05   ** p<0.01   *** p<0.001

Robustness analysis: Recoded age similarity

Data preperation

fage_rec <- function (x) {
  y <- ifelse(x < 16, 1, x)
  y <- ifelse(x > 15 & x < 21, 2, y)
  y <- ifelse(x > 20 & x < 26, 3, y)
  y <- ifelse(x > 25 & x < 31, 4, y)
  y <- ifelse(x > 30 & x < 36, 5, y)
  y <- ifelse(x > 35 & x < 41, 6, y)
  y <- ifelse(x > 40 & x < 46, 7, y)
  y <- ifelse(x > 45 & x < 51, 8, y)
  y <- ifelse(x > 50 & x < 56, 9, y)
  y <- ifelse(x > 55 & x < 61, 10, y)
  y <- ifelse(x > 60 & x < 66, 11, y)
  y <- ifelse(x > 65 & x < 71, 12, y)
  y <- ifelse(x > 70, 13, y)
  return(y)
}


MyData <- MyData %>% 
  mutate(age_rec = fage_rec((age+15)),
         dyad_age_sim_rec = abs(age_rec - (age_alter+1))/age_rec,
         dyad_age_sim_rec_cen = scale(dyad_age_sim_rec),
         dyad_age_sim_rec_rev = 5.5 - dyad_age_sim_rec,
         dyad_age_sim_rec_rev_cen = scale(dyad_age_sim_rec_rev),
         dyad_age_sim_rec_2 = 1 - ((abs(age_rec - (age_alter+1))/12)/age_rec),
         dyad_age_sim_rec_2_cen = scale(dyad_age_sim_rec_2))

age_dyad_rec_list <- MyData %>% 
  select(nomem_encr, dyad_id, survey_wave, age_alter) %>% 
  mutate(age_alter = age_alter + 1) %>% 
  group_split(nomem_encr, survey_wave) 


f_create_age_rec_sim <- function(df) { #test = age_dyad_rec_list[[4]]
  #df = test
  
  #create list to store information
  net_results_list <- list()
  
  #create for every alter avsim alter for age
  for (i in 1:nrow(df)) {
    #i= 1
    #select alter and remove alter from net
    alter <- df[i,]
    net <- df[-i,]
    
    #extract age
    age_net <- as.vector(t(net$age_alter))
    age_alter <- alter$age_alter
    
    #for con sim
    result_df <-
      tibble(alter = age_net, net = age_alter) %>%
      filter(!is.na(net) &
               !is.na(alter)) %>% #filter out missings
      mutate(
        dyad_age_sim = abs(alter - net) / alter,
        avsim_alter_age_rec = mean(dyad_age_sim)
      )
    
    #combine estimates and create combined sim
    result_df <- result_df %>%
      mutate(
        dyad_id =   df[i,] %>%
          pull(dyad_id),
        nomem_encr =   df[i,] %>%
          pull(nomem_encr),
        survey_wave = df[i,] %>%
          pull(survey_wave)
      ) %>%
      select(nomem_encr, dyad_id, survey_wave, avsim_alter_age_rec) %>%
      distinct()
  
  #store results of alter
  net_results_list[[i]] <- result_df
  }
  
  #store in df
  results_df <- net_results_list %>%
    bind_rows()
  
  return(results_df)
}

  #use future_map and the f_make_net_variables_df
#start parallel session
plan(multisession, workers = 7)

age_dyad_rec_result_list <- age_dyad_rec_list %>%
  future_map(
    .x = .,
    .f = ~ f_create_age_rec_sim(
      df = .x
    ),
    .progress = T
  )

plan(sequential)

#combine results into one df
age_dyad_rec_result_list <- age_dyad_rec_result_list %>% 
  bind_rows()

#add to MyData
MyData <- MyData %>% 
  left_join(age_dyad_rec_result_list, by = c("nomem_encr", "dyad_id", "survey_wave"))


#mean imputation
MyData <- MyData %>% 
  mutate(avsim_alter_age_rec_rev = 1- avsim_alter_age_rec,
    avsim_alter_age_rec_rev = ifelse(is.na(avsim_alter_age_rec_rev),
                                       mean(avsim_alter_age_rec_rev, na.rm = T),
                                       avsim_alter_age_rec_rev),
         avsim_alter_age_rec_rev_cen = scale(avsim_alter_age_rec_rev))

Model estimation

Create formulas

#create empty list
age_rec_formula_list <- list()

#create formula's

#model 1
age_rec_formula_list[[1]] <- as.formula(
  dropped ~ (1 | nomem_encr) 
)

#model 2
age_rec_formula_list[[2]] <- as.formula(dropped ~ (1 | nomem_encr) +
                                  time)



#model 3
age_rec_formula_list[[3]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen
)

#model 4
age_rec_formula_list[[4]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen
)

#model 5
age_rec_formula_list[[5]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen
)

#model 6
age_rec_formula_list[[6]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    size_net_one
)

#model 7
age_rec_formula_list[[7]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    size_net_one +
    as.factor(dear_alter_rec)
)

#model 8
age_rec_formula_list[[8]]  <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen + 
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    size_net_one
)


#model 9
age_rec_formula_list[[9]]  <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen + 
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    as.factor(dear_alter_rec)  +
    size_cen +
    net_density_cen +
    size_net_one
)


#model 10
age_rec_formula_list[[10]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    avsim_alter_age_rec_rev_cen:dyad_age_sim_rec_rev_cen +
    size_net_one
)


#model 11
age_rec_formula_list[[11]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_rec_rev_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    as.factor(dear_alter_rec) +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_rec_rev_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    time:dyad_age_sim_rec_rev_cen +
    size_net_one
)

Create directories

#create sub file
if(!dir.exists("results/survival_results/robustness/age_rec")){
  dir.create("results/survival_results/robustness/age_rec")
}

Estimate models

#set up of file name and dir
#set file name to use
dir <- "results/survival_results/robustness/age_rec/"
file_name_model <- paste0(dir, "model_age-rec_")


#initiate do parallel session
registerDoParallel(cores = 6)

nr_list <-  c("01",
              "02",
              "03",
              "04",
              "05",
              "06",
              "07",
              "08",
              "09",
              "10",
              "11")

#create foreach loop to use parallel computation to run fixed effects models
foreach(i = 1:11,
        .final = function(x) NULL) %do% {
  #i = 8
  #create file_name
  file_name <- paste0(file_name_model, nr_list[i], ".rda")
  
  #check whether filename exists
  if (!file.exists(file_name)) {
    model_results <- glmer(
      formula = age_rec_formula_list[[i]],
      data = MyData,
      family = binomial(link = "cloglog"),
      nAGQ = 0,
      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun =
                                                                    200000))
    )
    
    save(model_results,
         file = file_name)
  }
        }

NULL

registerDoSEQ()

Import results

#import demo interaction models
dir <- "results/survival_results/robustness/age_rec/"

age_rec_model_list <- list.files(dir,
           full.names = T) %>% 
  map(.x = ., 
      .f = ~get(load(file = .x)))

Tables

tab_model(age_rec_model_list[[3]],
          age_rec_model_list[[4]],
          age_rec_model_list[[5]],
          age_rec_model_list[[6]],
          age_rec_model_list[[7]],
          age_rec_model_list[[8]],
          age_rec_model_list[[9]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped		dropped		dropped		dropped		dropped		dropped		dropped
Predictors	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic
(Intercept)	1.855 ***	46.567	1.876 ***	47.501	1.856 ***	46.559	2.823 ***	29.832	3.071 ***	31.282	2.677 ***	27.789	2.926 ***	29.328
time	0.737 ***	-39.212	0.730 ***	-40.380	0.737 ***	-39.182	0.747 ***	-36.339	0.753 ***	-35.494	0.750 ***	-35.978	0.754 ***	-35.198
dyad_educ_sim_cen	0.989	-1.588			0.987	-1.904	0.993	-0.952	0.992	-1.118	0.992	-1.066	0.991	-1.207
dyad_age_sim_rec_rev_cen	0.925 ***	-10.971			0.922 ***	-10.137	0.960 ***	-4.973	0.960 ***	-4.948	0.960 ***	-4.889	0.960 ***	-4.873
dyad_gender_sim_cen	0.904 ***	-15.488			0.903 ***	-15.037	0.922 ***	-11.657	0.917 ***	-12.302	0.921 ***	-11.770	0.916 ***	-12.375
dyad_ethnicity_sim_cen	0.972 ***	-3.829			0.972 **	-3.260	0.988	-1.093	0.987	-1.112	0.987	-1.127	0.987	-1.140
avsim_alter_educ_cen			1.005	0.758	1.009	1.225	1.011	1.508	1.013	1.808	1.012	1.719	1.014 *	1.973
avsim_alter_age_rec_rev_cen			0.963 ***	-5.432	1.007	0.886	1.015	1.860	1.013	1.660	1.016 *	1.983	1.014	1.772
ei_alter_gender_rev_cen			0.977 ***	-3.472	1.005	0.685	1.003	0.397	1.009	1.147	1.005	0.676	1.010	1.355
ei_alter_ethnicity_rev_cen			0.982 **	-2.647	0.999	-0.075	1.018	1.842	1.017	1.715	1.019	1.869	1.018	1.739
educ_ego_cen							0.944 ***	-6.241	0.940 ***	-6.660	0.943 ***	-6.369	0.939 ***	-6.748
age_cen							1.062 ***	3.670	1.062 ***	3.633	1.059 ***	3.495	1.059 ***	3.485
gender_fac: Female							0.896 ***	-5.203	0.895 ***	-5.235	0.895 ***	-5.260	0.894 ***	-5.286
gender_fac: Missing							0.650 **	-2.870	0.646 **	-2.904	0.652 **	-2.850	0.647 **	-2.886
origin_rec_nar_fac: Non-western							1.238 ***	4.881	1.245 ***	4.964	1.236 ***	4.840	1.243 ***	4.924
origin_rec_nar_fac: Western							1.003	0.061	1.003	0.084	1.001	0.032	1.002	0.057
origin_rec_nar_fac: Missing							1.010	0.289	1.013	0.381	1.010	0.289	1.013	0.378
divorced_fac: divorced							0.917	-1.886	0.936	-1.441	0.919	-1.840	0.937	-1.413
divorced_fac: missing							0.948 **	-3.179	0.962 *	-2.277	0.951 **	-2.975	0.965 *	-2.138
moving_fac: new_residence							1.065	1.240	1.070	1.333	1.065	1.248	1.070	1.331
moving_fac: new_municipality							1.125	1.686	1.118	1.584	1.125	1.685	1.118	1.587
moving_fac: missing							0.986	-0.848	0.985	-0.886	0.987	-0.767	0.986	-0.815
first_child_fac: First child born							1.104	0.820	1.099	0.781	1.112	0.884	1.107	0.840
first_child_fac: Missing							1.009	0.417	1.010	0.472	1.010	0.479	1.011	0.522
educ_alter_cen							0.998	-0.269	0.998	-0.313	0.996	-0.499	0.996	-0.507
age_alter_cen							0.994	-0.440	0.990	-0.708	0.992	-0.556	0.988	-0.800
gender_alter_fac: Female							0.963 *	-2.110	0.969	-1.777	0.961 *	-2.238	0.967	-1.897
gender_alter_fac: Missing							1.708	1.700	1.722	1.725	1.662	1.614	1.682	1.651
origin_alter_rec_fac: Non-western							1.110 **	2.666	1.111 **	2.692	1.111 **	2.684	1.112 **	2.704
origin_alter_rec_fac: Western							1.070	1.453	1.074	1.524	1.069	1.426	1.073	1.500
origin_alter_rec_fac: Missing							1.395 **	2.693	1.417 **	2.787	1.376 *	2.575	1.399 **	2.684
rel_alter_rec: Same group or club							0.945	-1.499	0.966	-0.909	0.969	-0.845	0.986	-0.373
rel_alter_rec: Neighbour							0.911 **	-2.692	0.927 *	-2.191	0.931 *	-2.056	0.944	-1.664
rel_alter_rec: Friend							0.767 ***	-11.269	0.802 ***	-9.276	0.786 ***	-10.108	0.818 ***	-8.370
rel_alter_rec: Advisor							1.231 **	3.006	1.222 **	2.878	1.238 **	3.090	1.229 **	2.958
rel_alter_rec: Other							1.087 *	2.089	1.108 *	2.545	1.108 *	2.554	1.125 **	2.923
rel_alter_rec: Missing							0.585 *	-2.386	0.584 *	-2.379	0.605 *	-2.235	0.601 *	-2.250
times_dropped_earlier_cen							0.895 ***	-15.850	0.895 ***	-15.836	0.895 ***	-15.765	0.895 ***	-15.759
length_fac: 3 - 6 years							0.875 ***	-5.654	0.884 ***	-5.154	0.889 ***	-4.931	0.897 ***	-4.554
length_fac: > 6 years							0.833 ***	-8.130	0.853 ***	-7.031	0.856 ***	-6.868	0.872 ***	-5.998
length_fac: Length missing							1.129	1.259	1.159	1.523	1.158	1.517	1.183	1.731
size_cen							1.013	1.627	1.012	1.485	1.074 ***	6.584	1.064 ***	5.658
net_density_cen							0.966 ***	-4.670	0.962 ***	-5.212	1.027 *	2.466	1.013	1.211
size_net_one							1.120 *	2.425	1.145 **	2.872	1.327 ***	5.489	1.321 ***	5.390
as.factor(dear_alter_rec)1									0.691 ***	-15.034			0.698 ***	-14.554
as.factor(dear_alter_rec)2									0.830 ***	-10.840			0.836 ***	-10.378
degree_cen											0.907 ***	-7.882	0.920 ***	-6.682
Random Effects
σ2	1.64		1.64		1.64		1.64		1.64		1.64		1.64
τ00	0.19 nomem_encr		0.18 nomem_encr		0.19 nomem_encr		0.18 nomem_encr		0.19 nomem_encr		0.18 nomem_encr		0.19 nomem_encr
ICC	0.10		0.10		0.10		0.10		0.10		0.10		0.10
N	6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr		6996 nomem_encr
Observations	49449		49449		49449		49449		49449		49449		49449
Marginal R2 / Conditional R2	0.086 / 0.179		0.079 / 0.172		0.086 / 0.179		0.112 / 0.200		0.116 / 0.206		0.113 / 0.201		0.116 / 0.207
Deviance	58353.479		58727.140		58350.465		57287.163		57037.398		57226.577		56993.909
log-Likelihood	-29176.740		-29363.570		-29175.233		-28643.582		-28518.699		-28613.289		-28496.954
p<0.05   ** p<0.01   *** p<0.001

tab_model(age_rec_model_list[[10]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped
Predictors	Risk Ratios	Statistic
(Intercept)	2.918 ***	29.234
time	0.755 ***	-35.157
avsim_alter_educ_cen	1.015 *	2.015
avsim_alter_age_rec_rev_cen	1.005	0.606
ei_alter_gender_rev_cen	1.011	1.413
ei_alter_ethnicity_rev_cen	1.018	1.769
educ_ego_cen	0.939 ***	-6.779
age_cen	1.056 ***	3.292
gender_fac: Female	0.894 ***	-5.269
gender_fac: Missing	0.641 **	-2.950
origin_rec_nar_fac: Non-western	1.241 ***	4.894
origin_rec_nar_fac: Western	1.001	0.031
origin_rec_nar_fac: Missing	1.013	0.394
divorced_fac: divorced	0.938	-1.391
divorced_fac: missing	0.964 *	-2.148
moving_fac: new_residence	1.069	1.316
moving_fac: new_municipality	1.118	1.583
moving_fac: missing	0.987	-0.787
first_child_fac: First child born	1.107	0.842
first_child_fac: Missing	1.011	0.544
educ_alter_cen	0.997	-0.442
age_alter_cen	0.991	-0.612
gender_alter_fac: Female	0.966	-1.929
gender_alter_fac: Missing	1.674	1.627
origin_alter_rec_fac: Non-western	1.111 **	2.685
origin_alter_rec_fac: Western	1.073	1.509
origin_alter_rec_fac: Missing	1.420 **	2.786
as.factor(dear_alter_rec)1	0.698 ***	-14.563
as.factor(dear_alter_rec)2	0.836 ***	-10.367
rel_alter_rec: Same group or club	0.988	-0.312
rel_alter_rec: Neighbour	0.945	-1.619
rel_alter_rec: Friend	0.821 ***	-8.202
rel_alter_rec: Advisor	1.231 **	2.979
rel_alter_rec: Other	1.130 **	3.015
rel_alter_rec: Missing	0.590 *	-2.326
times_dropped_earlier_cen	0.895 ***	-15.705
length_fac: 3 - 6 years	0.897 ***	-4.554
length_fac: > 6 years	0.873 ***	-5.946
length_fac: Length missing	1.177	1.681
degree_cen	0.920 ***	-6.668
dyad_educ_sim_cen	0.991	-1.206
dyad_age_sim_rec_rev_cen	0.951 ***	-5.452
dyad_gender_sim_cen	0.916 ***	-12.374
dyad_ethnicity_sim_cen	0.987	-1.151
size_cen	1.064 ***	5.662
net_density_cen	1.014	1.250
size_net_one	1.320 ***	5.367
avsim_alter_age_rec_rev_cen:dyad_age_sim_rec_rev_cen	0.996 *	-2.355
Random Effects
σ2	1.64
τ00 nomem_encr	0.19
ICC	0.10
N nomem_encr	6996
Observations	49449
Marginal R2 / Conditional R2	0.117 / 0.208
Deviance	56987.980
log-Likelihood	-28493.990
p<0.05   ** p<0.01   *** p<0.001

Robustness analysis: different life events specifications

Model estimation

Create formulas

#create empty list
le_formula_list <- list()

#model 1
le_formula_list[[1]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_fac +
    moving_fac +
    first_child_fac +
    employment_fac +
    unemployment_fac + 
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    size_net_one
)

#model 2
le_formula_list[[2]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_event +
    new_residence_event +
    new_municipality_event +
    first_child +
    employment_event +
    unemployed_event + 
    move_event_missing +
    unemployed_event_missing + 
    divorced_event_missing + 
    first_child_event_missing +
    employment_event_missing +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    size_net_one
)


#model 3
le_formula_list[[3]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    avsim_alter_educ_cen +
    avsim_alter_age_cen +
    ei_alter_gender_rev_cen +
    ei_alter_ethnicity_rev_cen +
    educ_ego_cen +
    age_cen +
    gender_fac +
    origin_rec_nar_fac +
    divorced_transition +
    new_residence_transition +
    new_municipality_transition +
    first_child_transition +
    employment_transition +
    unemployed_transition + 
    move_event_missing +
    unemployed_event_missing + 
    divorced_event_missing + 
    first_child_event_missing +
    employment_event_missing +
    educ_alter_cen +
    age_alter_cen +
    gender_alter_fac +
    origin_alter_rec_fac +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    size_cen +
    net_density_cen +
    size_net_one
)

Create directories

#create sub file
if(!dir.exists("results/survival_results/robustness/life_events")){
  dir.create("results/survival_results/robustness/life_events")
}

Estimate models

#set up of file name and dir
#set file name to use
dir <- "results/survival_results/robustness/life_events/" 
file_name_model <- paste0(dir, "model_le_")

#initiate do parallel session
registerDoParallel(cores = 3)

nr_list <-  c("01",
              "02",
              "03")

#create foreach loop to use parallel computation to run fixed effects models
foreach(i = 1:3,
        .final = function(x) NULL) %do% {
  #i = 8
  #create file_name
  file_name <- paste0(file_name_model, nr_list[i], ".rda")
  
  #check whether filename exists
  if (!file.exists(file_name)) {
    model_results <- glmer(
      formula = le_formula_list[[i]],
      data = MyData,
      family = binomial(link = "cloglog"),
      nAGQ = 0,
      control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun =
                                                                    200000))
    )
    
    save(model_results,
         file = file_name)
  }
        }

NULL

registerDoSEQ()

Import results

#import naq_0 effects models
dir <- "results/survival_results/robustness/life_events"

life_events_model_list <- list.files(dir,
           full.names = T) %>% 
  map(.x = ., 
      .f = ~get(load(file = .x)))

Tables

tab_model(life_events_model_list[[1]],
          life_events_model_list[[2]],
          life_events_model_list[[3]],
          show.p = T,
          show.ci = F,
          p.style = "stars",
          show.se =  F,
          show.stat = T,
          show.dev = T,
          show.loglik = T,
          show.ngroups = T,
          digits = 3)

	dropped		dropped		dropped
Predictors	Risk Ratios	Statistic	Risk Ratios	Statistic	Risk Ratios	Statistic
(Intercept)	2.826 ***	29.828	2.821 ***	31.187	2.769 ***	30.617
time	0.744 ***	-36.633	0.744 ***	-36.787	0.743 ***	-36.857
avsim_alter_educ_cen	1.010	1.386	1.010	1.389	1.010	1.378
avsim_alter_age_cen	1.010	1.381	1.011	1.394	1.011	1.398
ei_alter_gender_rev_cen	1.002	0.223	1.002	0.242	1.002	0.295
ei_alter_ethnicity_rev_cen	1.019	1.879	1.019	1.883	1.019	1.900
educ_ego_cen	0.949 ***	-5.566	0.949 ***	-5.622	0.948 ***	-5.672
age_cen	1.030	1.816	1.031	1.859	1.031	1.900
gender_fac: Female	0.893 ***	-5.391	0.893 ***	-5.406	0.892 ***	-5.436
gender_fac: Missing	0.795	-1.328	0.792	-1.351	0.798	-1.305
origin_rec_nar_fac: Non-western	1.233 ***	4.782	1.231 ***	4.757	1.233 ***	4.785
origin_rec_nar_fac: Western	1.002	0.045	1.002	0.043	1.002	0.058
origin_rec_nar_fac: Missing	1.013	0.384	1.013	0.399	1.011	0.335
divorced_fac: divorced	0.928	-1.622
divorced_fac: missing	0.988	-0.677
moving_fac: new_residence	1.063	1.214
moving_fac: new_municipality	1.128	1.715
moving_fac: missing	1.008	0.459
first_child: first_child_fac: First child born	1.118	0.931
first_child: first_child_fac: Missing	0.993	-0.345
employment_fac: Employed	1.116	1.323
employment_fac: Missing	0.890 ***	-5.093
unemployment_fac: Employed	1.101	1.516
unemployment_fac: Missing	0.916	-1.013
educ_alter_cen	1.001	0.100	1.001	0.076	1.000	0.033
age_alter_cen	1.014	0.954	1.014	0.944	1.014	0.990
gender_alter_fac: Female	0.966 *	-1.974	0.965 *	-1.984	0.966	-1.958
gender_alter_fac: Missing	1.794	1.876	1.790	1.870	1.775	1.844
origin_alter_rec_fac: Non-western	1.108 **	2.624	1.108 **	2.617	1.106 **	2.582
origin_alter_rec_fac: Western	1.071	1.464	1.070	1.458	1.071	1.467
origin_alter_rec_fac: Missing	1.430 **	2.906	1.430 **	2.902	1.424 **	2.868
rel_alter_rec: Same group or club	0.934	-1.818	0.933	-1.827	0.932	-1.860
rel_alter_rec: Neighbour	0.902 **	-2.970	0.902 **	-2.989	0.901 **	-3.019
rel_alter_rec: Friend	0.761 ***	-11.538	0.760 ***	-11.572	0.759 ***	-11.636
rel_alter_rec: Advisor	1.222 **	2.901	1.222 **	2.900	1.223 **	2.913
rel_alter_rec: Other	1.081	1.945	1.081	1.936	1.080	1.914
rel_alter_rec: Missing	0.581 *	-2.430	0.582 *	-2.418	0.589 *	-2.372
times_dropped_earlier_cen	0.893 ***	-16.096	0.893 ***	-16.209	0.890 ***	-16.506
length_fac: 3 - 6 years	0.874 ***	-5.691	0.875 ***	-5.668	0.876 ***	-5.610
length_fac: > 6 years	0.836 ***	-8.000	0.836 ***	-7.974	0.838 ***	-7.879
length_fac: Length missing	1.153	1.485	1.154	1.487	1.154	1.484
dyad_educ_sim_cen	0.993	-0.918	0.993	-0.919	0.994	-0.883
dyad_age_sim_cen	0.967 ***	-4.380	0.967 ***	-4.362	0.968 ***	-4.335
dyad_gender_sim_cen	0.922 ***	-11.533	0.922 ***	-11.529	0.922 ***	-11.536
dyad_ethnicity_sim_cen	0.987	-1.123	0.987	-1.132	0.987	-1.147
size_cen	1.015	1.822	1.015	1.839	1.015	1.834
net_density_cen	0.968 ***	-4.430	0.968 ***	-4.397	0.968 ***	-4.359
size_net_one	1.117 *	2.347	1.116 *	2.346	1.116 *	2.337
divorced_event			0.941	-1.323
new_residence_event			1.085	1.614
new_municipality_event			1.014	0.190
first_child			1.111	0.878
employment_event			1.106	1.270
unemployed_event			1.137 *	2.109
unemployed_event_missing			0.918	-0.994	0.910	-1.084
divorced_event_missing			0.992	-0.439	0.997	-0.171
first_child_event_missing			0.803	-0.915	0.992	-0.399
employment_event_missing			0.891 ***	-5.210	0.890 ***	-5.208
divorced_transition					0.985	-0.439
new_residence_transition					1.131 ***	3.303
new_municipality_transition					1.093	1.766
first_child_transition					1.090	0.925
employment_transition					1.003	0.040
unemployed_transition					1.074	1.467
Random Effects
σ2	1.64		1.64		1.64
τ00	0.18 nomem_encr		0.18 nomem_encr		0.18 nomem_encr
ICC	0.10		0.10		0.10
N	6996 nomem_encr		6996 nomem_encr		6996 nomem_encr
Observations	49449		49449		49449
Marginal R2 / Conditional R2	0.114 / 0.201		0.114 / 0.201		0.114 / 0.201
Deviance	57263.625		57261.262		57242.704
log-Likelihood	-28631.813		-28630.631		-28621.352
p<0.05   ** p<0.01   *** p<0.001

Robustness analysis: Hybrid Model specification

Data preperation and model estimation

#Create ego fixed effects. 
MyData <- MyData %>% 
  group_by(nomem_encr) %>% 
  mutate(mean_dyad_educ = mean(dyad_educ_sim, na.rm = T),
         mean_dyad_age = mean(dyad_age_sim, na.rm = T),
         mean_dyad_gender = mean(dyad_gender_sim, na.rm = T),
         mean_dyad_ethnicity = mean(dyad_ethnicity_sim, na.rm = T),
         mean_net_educ = mean(avsim_alter_educ, na.rm = T),
         mean_net_age = mean(avsim_alter_age, na.rm = T),
         mean_net_gender = mean(ei_alter_gender_rev, na.rm = T),
         mean_net_ethnicity = mean(ei_alter_ethnicity_rev, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(dev_dyad_educ = dyad_educ_sim - mean_dyad_educ,
         dev_dyad_age = dyad_age_sim - mean_dyad_age,
         dev_dyad_gender =  dyad_gender_sim - mean_dyad_gender,
         dev_dyad_ethnicity = dyad_ethnicity_sim - mean_dyad_ethnicity,
         dev_net_educ = avsim_alter_educ - mean_net_educ,
         dev_net_age = avsim_alter_age - mean_net_age,
         dev_net_gender = ei_alter_gender_rev - mean_net_gender,
         dev_net_ethnicity = ei_alter_ethnicity_rev - mean_net_ethnicity)


#Crete eogo-surv wave fixed effects.
#model 1 with only main effects
M1 <- glmer(
  formula = dropped ~ (1|nomem_encr) + time + 
    mean_dyad_educ +
    mean_dyad_age +
    mean_dyad_gender +
    mean_dyad_ethnicity +
    dev_dyad_educ +
    dev_dyad_age + 
    dev_dyad_gender +
    dev_dyad_ethnicity +
    mean_net_educ +
    mean_net_age +
    mean_net_gender +
    mean_net_ethnicity +
    dev_net_educ +
    dev_net_age +
    dev_net_gender +
    dev_net_ethnicity,
  data = MyData,
  family = binomial(link = "cloglog"),
  nAGQ = 0,
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun =
                                                                200000))
)

tab_model(
  M1,
  show.p = T,
  show.ci = F,
  p.style = "stars",
  show.se =  F,
  show.stat = T,
  show.dev = T,
  show.loglik = T,
  show.ngroups = T,
  digits = 3
)

	dropped
Predictors	Risk Ratios	Statistic
(Intercept)	7.414 ***	14.939
time	0.737 ***	-39.077
mean_dyad_educ	0.919	-1.032
mean_dyad_age	0.362 ***	-8.421
mean_dyad_gender	0.790 ***	-6.026
mean_dyad_ethnicity	0.911 **	-2.815
dev_dyad_educ	0.923	-1.679
dev_dyad_age	0.525 ***	-8.530
dev_dyad_gender	0.764 ***	-13.265
dev_dyad_ethnicity	0.949	-1.309
mean_net_educ	1.099	0.928
mean_net_age	0.814 *	-2.018
mean_net_gender	0.929 ***	-3.468
mean_net_ethnicity	1.005	0.170
dev_net_educ	1.069	1.031
dev_net_age	1.423 ***	4.302
dev_net_gender	1.054 ***	3.616
dev_net_ethnicity	0.987	-0.575
Random Effects
σ2	1.64
τ00 nomem_encr	0.19
ICC	0.10
N nomem_encr	6996
Observations	49449
Marginal R2 / Conditional R2	0.089 / 0.182
Deviance	58276.470
log-Likelihood	-29138.235
p<0.05   ** p<0.01   *** p<0.001

---

Copyright © 2023 Jeroense Thijmen