Ego-net deselection Analysis: Main analysis

Goal

Recreate the results.

Set up

Packages

#library
library(tidyverse)
library(data.table)
library(kableExtra)
library(patchwork)
library(ggpubr)
library(lme4)
library(viridis)
library(survival)
library(sjPlot)
library(foreach)
library(doParallel)
library(ggeffects)

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 == "avsim_alter_educ_cen" ~ 5,
      effect_name == "avsim_alter_age_cen" ~ 6,
      effect_name == "ei_alter_gender_rev_cen" ~ 7,
      effect_name == "ei_alter_ethnicity_rev_cen" ~ 8,
      effect_name == "as.factor(dear_alter_rec)1" ~ 9,
      effect_name == "as.factor(dear_alter_rec)2" ~ 10,
      effect_name == "degree_cen" ~ 11,
      effect_name == "avsim_alter_educ_cen:dyad_educ_sim_cen" ~ 12,
      effect_name == "avsim_alter_age_cen:dyad_age_sim_cen" ~ 13,
      effect_name == "ei_alter_gender_rev_cen:dyad_gender_sim_cen" ~ 14,
      effect_name == "ei_alter_ethnicity_rev_cen:dyad_ethnicity_sim_cen" ~ 15,
      effect_name == "time:dyad_educ_sim_cen" ~ 16,
      effect_name == "time:dyad_age_sim_cen" ~ 17,
      effect_name == "time:dyad_gender_sim_cen" ~ 18,
      effect_name == "time:dyad_ethnicity_sim_cen" ~ 19,
      effect_name == "time:as.factor(dear_alter_rec)1" ~ 20,
      effect_name == "time:as.factor(dear_alter_rec)2" ~ 21,
      effect_name == "time:degree_cen" ~ 22,
      effect_name == "(Intercept)" ~ 23,
      effect_name == "time" ~ 24,
      effect_name == "educ_ego_cen" ~ 25,
      effect_name == "age_cen" ~ 26,
      effect_name == "gender_facFemale" ~ 27,
      effect_name == "gender_facMissing" ~ 28,
      effect_name == "origin_rec_nar_facNon-western" ~ 29,
      effect_name == "origin_rec_nar_facWestern" ~ 30,
      effect_name == "origin_rec_nar_facMissing" ~ 31,
      effect_name == "divorced_facdivorced" ~ 32,
      effect_name == "divorced_facmissing" ~ 33,
      effect_name == "moving_facnew_residence" ~ 34,
      effect_name == "moving_facnew_municipality" ~ 35,
      effect_name == "moving_facmissing" ~ 36,
      effect_name == "first_child_facFirst child born" ~ 37,
      effect_name == "first_child_facMissing" ~ 38,
      effect_name == "educ_alter_cen" ~ 39,
      effect_name == "age_alter_cen" ~ 40,
      effect_name == "gender_alter_facFemale" ~ 41,
      effect_name == "gender_alter_facMissing" ~ 42,
      effect_name == "origin_alter_rec_facNon-western" ~ 43,
      effect_name == "origin_alter_rec_facWestern" ~ 44,
      effect_name == "origin_alter_rec_facMissing" ~ 45,
      effect_name == "rel_alter_recSame group or club" ~ 46,
      effect_name == "rel_alter_recNeighbour" ~ 47,
      effect_name == "rel_alter_recFriend" ~ 48,
      effect_name == "rel_alter_recAdvisor" ~ 49,
      effect_name == "rel_alter_recOther" ~ 50,
      effect_name == "rel_alter_recMissing" ~ 51,
      effect_name == "times_dropped_earlier_cen" ~ 52,
      effect_name == "length_fac3 - 6 years" ~ 53,
      effect_name == "length_fac> 6 years" ~ 54,
      effect_name == "length_facLength missing" ~ 55,
      effect_name == "size_cen" ~ 56,
      effect_name == "net_density_cen" ~ 57,
      effect_name == "size_net_one" ~ 58
    ),
    parameter_fac = factor(
      parameter,
      levels = 1:58,
      labels = c(
      "Dyadic similarity: education",
      "Dyadic similarity: age",
      "Dyadic similarity: gender",
      "Dyadic similarity: ethnicity",
      "Confidant uniqueness: education",
      "Confidant uniqueness: age",
      "Confidant uniqueness: gender",
      "Confidant uniqueness: ethnicity",
      "Closeness: dear",
      "Closeness: not asked",
      "Alter Embeddedness",
      "Interaction education",
      "Interaction age",
      "Interaction gender",
      "Interaction ethnicity",
      "Interaction time education",
      "Interaction time age",
      "Interaction time gender",
      "Interaction time ethnicity",
      "Interaction time dear",
      "Interaction time dear not asked",
      "Interaction time degree",
      "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"
             )
           )
         ) 
}

Import data

#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))

Analyses

Time specification

#linear time
ml_m0.1 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) + time,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))

#^2 time
ml_m0.2 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))


#^3 time
ml_m0.3 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))



#^4 time
ml_m0.4 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3 +
        time_4,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))



#^5 time
ml_m0.5 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3 +
        time_4 +
        time_5,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))


#^6 time
ml_m0.6 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3 +
        time_4 +
        time_5 +
        time_6,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))


#^7 time
ml_m0.7 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3 +
        time_4 + 
        time_5 +
        time_6 +
        time_7,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))



#^8 time
ml_m0.8 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3 +
        time_4 +
        time_5 +
        time_6 +
        time_7 +
        time_8,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))


#^9 time
ml_m0.9 <- MyData %>%
  glmer(formula = dropped ~ (1|nomem_encr) +
        time + 
        time_2 +
        time_3 +
        time_4 + 
        time_5 +
        time_6 +
        time_7 +
        time_8 +
        time_9,
      family = binomial(link = 'cloglog'),
      nAGQ = 0,
      control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=200000)))


#store in list
time_models <- list(ml_m0.1,ml_m0.2, ml_m0.3, ml_m0.4, ml_m0.5, ml_m0.6,
                    ml_m0.7, ml_m0.8, ml_m0.9)

#extract fit info
#create fit list
fit_time_list <- list()

#loop
for(i in 1:length(time_models)) {#i = 1
  model <- time_models[[i]] #extract correct model
  fit <- llikAIC(model)
  fit_table <- as_tibble(fit$AICtab) %>% 
    transpose() %>% 
    rename(AIC = V1,
           BIC = V2,
           LogLikelihood = V3,
           Deviance = V4,
           DF = V5) %>% 
    mutate(model = paste0("Model ", i)) #create model collumn
  
  fit_time_list[[i]] <- fit_table %>%
    select(model,AIC,BIC,LogLikelihood, Deviance, DF) #export fit model
}


#run lrt tests
models_anova_time <- anova(ml_m0.1,ml_m0.2, ml_m0.3, ml_m0.4, ml_m0.5, ml_m0.6,
                    ml_m0.7, ml_m0.8, ml_m0.9, 
              test = "LRT")

#combine fit with lrt tests
fit_time_list <- fit_time_list %>% 
  rbindlist() %>%
  cbind(.,as.data.frame(models_anova_time[,c(1,5,6,7,8)])) %>% 
  as_tibble(.) %>% 
  select(model, AIC, BIC, LogLikelihood, Deviance, DF, Df, Chisq, `Pr(>Chisq)`)

#rename collumns
names(fit_time_list) <- c("Model", "AIC", "BIC","LogLikelihood", "Deviance", "DF","DF difference", "Deviance difference","Pr(>Chi)")

#fit table
fit_time_list %>%
  kbl(caption = "Different time specifications", digits = 3) %>%
  kable_classic(full_width = F, bootstrap_options = c("hover", "condensed"), fixed_thead = T)

Different time specifications

Model	AIC	BIC	LogLikelihood	Deviance	DF	DF difference	Deviance difference	Pr(>Chi)
Model 1	58785.90	58812.33	-29389.95	58779.90	49446	NA	NA	NA
Model 2	58507.80	58543.04	-29249.90	58499.80	49445	1	280.099	0.000
Model 3	58466.57	58510.61	-29228.28	58456.57	49444	1	43.238	0.000
Model 4	58463.73	58516.58	-29225.86	58451.73	49443	1	4.839	0.028
Model 5	58461.24	58522.90	-29223.62	58447.24	49442	1	4.484	0.034
Model 6	58463.02	58533.49	-29223.51	58447.02	49441	1	0.222	0.638
Model 7	58464.67	58543.95	-29223.33	58446.67	49440	1	0.352	0.553
Model 8	58466.29	58554.37	-29223.14	58446.29	49439	1	0.382	0.537
Model 9	58468.24	58565.13	-29223.12	58446.24	49438	1	0.050	0.822

ML models

First create directories for storing results of ME models

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

if(!dir.exists("results/survival_results/main/mixed_effects/")){
  dir.create("results/survival_results/main/mixed_effects/")
}

if(!dir.exists("results/survival_results/main/time_interaction/")){
  dir.create("results/survival_results/main/time_interaction/")
}

Create a formula list. We can use this to estimate models with parallel computation.

#set empty list
formula_list <- list()

#create formula's

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

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

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

#model 4
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
)

#model 5
formula_list[[5]] <- 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 +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen
)

#model 6
formula_list[[6]] <- 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 +
    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 7
formula_list[[7]] <- 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 +
    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 +
    as.factor(dear_alter_rec)
)

#model 8
formula_list[[8]]  <- 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 +
    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 +
    size_net_one
)


#model 9
formula_list[[9]]  <- 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 +
    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 +
    as.factor(dear_alter_rec)  +
    size_cen +
    net_density_cen +
    size_net_one
)

#model 10
formula_list[[10]] <- 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 +
    avsim_alter_educ_cen:dyad_educ_sim_cen +
    size_net_one
)


#model 11
formula_list[[11]] <- 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 +
    avsim_alter_age_cen:dyad_age_sim_cen +
    size_net_one
)


#model 12
formula_list[[12]] <- 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 +
    ei_alter_gender_rev_cen:dyad_gender_sim_cen +
    size_net_one
)

#model 13
formula_list[[13]] <- 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 +
    ei_alter_ethnicity_rev_cen:dyad_ethnicity_sim_cen +
    size_net_one
)

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

Estimate models in parallel and store results in list.

#set up of file name and dir
#set file name to use
dir <- "results/survival_results/main/mixed_effects/"
file_name_model <- paste0(dir, "model_me_")

#initiate do parallel session
registerDoParallel(cores = 7)

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

#create foreach loop to use parallel computation to run fixed effects models
foreach(i = 1:13,
        .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()

Interaction between time and different forms of similarity and hurdles to tie maintenance.

#set up of file name and dir
#set file name to use
dir <- "results/survival_results/main/time_interaction/" 
file_name_model <- paste0(dir, "model_time_int_")

Create formulas and store these in a list

formula_list_time <- list() 

formula_list_time[[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 +
    size_cen +
    net_density_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    time:dyad_ethnicity_sim_cen +
    size_net_one
)

formula_list_time[[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 +
    size_cen +
    net_density_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    time:dyad_gender_sim_cen +
    size_net_one
)

formula_list_time[[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 +
    size_cen +
    net_density_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    time:dyad_age_sim_cen +
    size_net_one
)

formula_list_time[[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 +
    size_cen +
    net_density_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    time:dyad_educ_sim_cen +
    size_net_one
)

formula_list_time[[5]] <- 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 +
    size_cen +
    net_density_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    time:degree_cen +
    size_net_one
)


formula_list_time[[6]] <- 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 +
    rel_alter_rec +
    times_dropped_earlier_cen +
    length_fac +
    degree_cen +
    size_cen +
    net_density_cen +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    as.factor(dear_alter_rec) +
    time:as.factor(dear_alter_rec) +
    size_net_one
)

formula_list_time[[7]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    as.factor(dear_alter_rec) +
    length_fac +
    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 +
    rel_alter_rec +
    times_dropped_earlier_cen +
    degree_cen +
    size_cen +
    net_density_cen +
    length_fac:as.factor(dear_alter_rec) +
    size_net_one
)

formula_list_time[[8]] <- as.formula(
  dropped ~ (1 | nomem_encr) +
    time +
    dyad_educ_sim_cen +
    dyad_age_sim_cen +
    dyad_gender_sim_cen +
    dyad_ethnicity_sim_cen +
    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 +
    size_cen +
    net_density_cen +
    length_fac:degree_cen +
    size_net_one
)

#create names
names(formula_list_time) <- foreach(i = 1:8,
                     .combine = c) %do% {
                       paste0("time_het_model_", i)
                     }

Estimate models and store results

#run models
registerDoParallel(cores = 7)

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

#create foreach loop to use parallel computation to run fixed effects models
foreach(i = 1:8,
        .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_time[[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 and add to model list

#import naq_0 effects models
dir <- "results/survival_results/main/mixed_effects/"

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

#import naq_0 effects models
dir <- "results/survival_results/main/time_interaction/"

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

Extract model info and store in separate tibbles

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


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

Output: tables and Plots

Coefficient plot main results

combined_coefplot <- me_model_df %>%
  filter(str_detect(parameter_fac, "Confidant uniqueness") |
           str_detect(parameter_fac, "Dyadic similarity")) %>%
  filter(model %in% c("Model 2", "Model 3", "Model 4","Model 5", "Model 6", "Model 7")) %>% 
  mutate(
    controls = ifelse(model %in% c("Model 2", "Model 3"), 1, NA),
    controls = ifelse(
      model %in% c("Model 4"),
      2,
      controls
    ),
    controls = ifelse(
      model %in% c("Model 5"),
      3,
      controls
    ),
    controls = ifelse(
      model %in% c("Model 6"),
      4,
      controls
    ),
    controls = ifelse(
      model %in% c("Model 7"),
      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", "Migration background")
    )
  ) %>%
  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_x_discrete(guide = guide_axis(n.dodge = 2)) +
  scale_x_discrete(label = scales::label_wrap(9)) +
  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_text(),
        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.jpg",
       dpi = 320,
       width = 7,
       height = 5)

Coefficient plot covariates

#control plot
control_coefplot_df <- me_model_df %>% 
  filter(m_number == 6) %>%
  filter(!(
    str_detect(parameter_fac, "Confidant uniqueness") |
      str_detect(parameter_fac, "Dyadic similarity") |
      str_detect(parameter_fac, "issing") |
      str_detect(parameter_fac, "Intercept") |
      str_detect(parameter_fac, "Size ==") |
      str_detect(parameter_fac, "Constant") |
      str_detect(parameter_fac, "time")
  )) %>%
  mutate(row = row_number()) %>% 
  mutate(
    group = case_when(
      row_number() < 10 ~ 1,
      row_number() > 9 & row_number() < 15 ~ 2,
      row_number() > 14 & row_number() < 23  ~ 3,
      row_number() > 22 ~ 4
    ),
    group = factor(
      group,
      levels = 1:4,
      labels = c( "Ego",
                 "Confidant",
                 "Dyad",
                 "CDN")
    ),
    sign = ifelse(p > 0.05, 0, 1),
    sign = factor(
           sign, 
           levels = 0:1, 
           labels = c("Not signficant at P < 0.05", "Significant at P < 0.05")
           ),
    parameter_fac = fct_rev(parameter_fac)
  )

control_coefplot <- control_coefplot_df %>% 
  ggplot(aes(x = parameter_fac, 
             y = riskratio,
             shape = sign)) +
  geom_hline(yintercept = 1) +
  geom_pointrange(aes(ymin = rr_min, ymax = rr_max),
                  position = position_dodge(width = 0.4)) +
  geom_point(aes(colour = sign,
                 fill = sign),
             size = 1,
             position = position_dodge(width = 0.4)) +
  facet_grid(vars(group),  scales = "free", switch = "y", space = "free_y") +
  coord_flip() +
  scale_x_discrete(position = "top") +
  scale_fill_manual(values = c("#1b9e77",
                               "#d95f02")) +
  scale_colour_manual(values = c("#1b9e77",
                               "#d95f02")) +
  scale_y_continuous(limits = c(0.25,1.75)) +
  scale_shape_manual(values = c(21,22)) +
  theme(panel.background = element_rect(fill = "#FFFFFF"),
        plot.background = element_rect(fill = "#FFFFFF"),
        panel.grid = 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_text(),
        axis.line = element_blank(),
        axis.title.y = element_text(hjust = 0.9, face = "bold"),
        axis.ticks = element_blank(),
        strip.background = element_rect(fill = "#A9A9A9"),
        strip.text = element_text(family = "sans", size = 8),
        panel.grid.minor = element_blank(),
        legend.position = "top",
        legend.direction = "horizontal",
        legend.text = element_text(family = "sans", size = 8), 
        legend.title = element_blank(),
        legend.background = element_rect(fill = "#FFFFFF"),
        legend.key = element_rect(fill = "#FFFFFF")) +
  labs(y = "Risk Ratio", x = "")

control_coefplot

ggsave(control_coefplot,
       file = "plots/results/survival/control_coefplot.jpg",
       dpi = 320,
       width = 5.5,
       height = 6)

Marginal effect plots

Uniqueness and dyadic similarity interaction

#create average marginal effects
me_effects_gender <- ggpredict(me_model_list[[12]],
                               terms = c("dyad_gender_sim_cen[-2.18351765957899,0.459267619659605]",
                                         "ei_alter_gender_rev_cen[-2.0842322544561,1.25777983594881]"),
                               ci.lvl = 0.95)



me_plot_gender <- me_effects_gender %>%
  as_tibble() %>%
  mutate(
    group = factor(
      as.numeric(group),
      levels = 1:2,
      labels = c("Unique",
                 "Common")
    ),
    x = ifelse(x < 0, 1, 2),
    sim = factor(
      as.numeric(x),
      levels = 1:2,
      labels = c("Ego-dyad different gender",
                 "Ego-dyad identical gender")
    )
  )

#create plot
me_plot_gender <- me_effects_gender %>%
  as_tibble() %>%
  mutate(
    group = factor(
      as.numeric(group),
      levels = 1:2,
      labels = c("Unique",
                 "Common")
    ),
    x = ifelse(x < 0, 1, 2),
    sim = factor(
      as.numeric(x),
      levels = 1:2,
      labels = c("Dissimilar gender",
                 "Similar gender")
    )
  ) %>% 
  ggplot(aes(x = sim,
             y = predicted, 
             group = group)) +
  geom_col(aes(colour = group,
             fill = group),
           position = position_dodge(0.7),
           width = 0.5,
           alpha = 0.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(0.7),
                width = 0.2) +
  geom_label(aes(label = round(predicted, 
                               3)),
             position = position_dodge(0.7),
             vjust = -0.33) +
  scale_y_continuous(limits = c(0,1)) +
  scale_alpha(guide = "none") +
  scale_fill_manual(values = c("#1b9e77",
                               "#d95f02")) +
  scale_colour_manual(values = c("#1b9e77",
                               "#d95f02")) +
  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(),
        text = element_text(family = "sans", size = 12), 
        axis.title.x = element_text(hjust = 0.9,face = "bold"),
        axis.text.x = element_text(),
        axis.line = element_blank(),
        axis.title.y = element_text( face = "bold"),
        axis.ticks = element_blank(),
        strip.background = element_rect(fill = "#A9A9A9"),
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        legend.background = element_rect(fill = "#FFFFFF"),
        legend.key = element_rect(fill = "#FFFFFF")) +
  labs(y = "Predicted probability",
       x = "",
       title = "",
       fill = "Confidant uniqueness: Gender",
       group = "Confidant uniqueness: Gender",
       colour = "Confidant uniqueness: Gender")

me_plot_gender

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

Time interactions

#create average marginal effect plots for m1, m2, and m3
if (!file.exists("results/survival_results/2023-06-14_me_effects.rda")) {
  m1 <- time_int_model_list[[2]]
  m2 <- time_int_model_list[[5]]
  m3 <- time_int_model_list[[6]]
  
  
  me_effects_age <- ggeffect(m1,
                             terms = c("dyad_age_sim_cen[-1,1]", "time[1:10]"),
                             ci.lvl = 0.95)
  
  me_effects_degree <- ggeffect(m2,
                                terms = c("degree_cen[-1,1]", "time[1:10]"),
                                ci.lvl = 0.95)
  
  
  me_effects_dear <- ggeffect(m3,
                              terms = c("as.factor(dear_alter)", "time[1:10]"),
                              ci.lvl = 0.95)
  
  me_effects_list <- list(me_effects_age,
                          me_effects_degree,
                          me_effects_dear)
  save(me_effects_list,
       file = "results/survival_results/2023-06-14_me_effects.rda")
} else {
  load("results/survival_results/2023-06-14_me_effects.rda")
}

#create plot extract function
me_effects_plot_function <- function(x,
                                     model){ #x = me_effects_age1
  df <- x %>% 
  as_tibble() %>% 
  mutate(time = as.numeric(group),
         model = model)
}

#use extract function to get information from me_effects object
me_effects_df <- foreach(a = 1:3,
        .combine = rbind) %do% {
          me_effects_plot_function(x = me_effects_list[[a]],
                                   model = a)
        }


#1b9e77
#d95f02

#create plot for degree and for dyadic similarity
age_me_effects_plot <- me_effects_df %>%
  filter(model == 1) %>%
  mutate(
    model = case_when(
      model == 1 ~ "Dyadic similarity: age",
      model == 2 ~ "Alter embeddedness"
    ),
    x = factor(x,
               levels = c(-1,1),
               labels = c("-1 SD",
                          "+1 SD"))
  ) %>% 
  ggplot(aes(
    x = time,
    y = predicted,
    group = x,
    colour = x,
    fill = x
  )) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.4,
              linetype = "blank") +
  facet_wrap(vars(model),
            ncol = 1) +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(limits = c(0,1)) +
  scale_fill_manual(values = c("#1b9e77",
                               "#d95f02")) +
  scale_colour_manual(values = c("#1b9e77",
                               "#d95f02")) +
  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 = "right",
        legend.direction = "vertical",
        legend.text = element_text(family = "sans", size = 8), 
        legend.title = element_blank(),
        legend.background = element_rect(fill = "#FFFFFF"),
        legend.key = element_rect(fill = "#FFFFFF")) +
  labs(y = "Predicted probability",
       x = "")

degree_me_effects_plot <- me_effects_df %>%
  filter(model == 2) %>%
  mutate(
    model = case_when(
      model == 1 ~ "Dyadic similarity: age",
      model == 2 ~ "Alter embeddedness"
    ),
    x = factor(x,
               levels = c(-1,1),
               labels = c("-1 SD",
                          "+1 SD"))
  ) %>% 
  ggplot(aes(
    x = time,
    y = predicted,
    group = x,
    colour = x,
    fill = x
  )) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.4,
              linetype = "blank") +
  facet_wrap(vars(model),
            ncol = 1) +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(limits = c(0,1)) +
  scale_fill_manual(values = c("#1b9e77",
                               "#d95f02")) +
  scale_colour_manual(values = c("#1b9e77",
                               "#d95f02")) +
  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 = "right",
        legend.direction = "vertical",
        legend.text = element_text(family = "sans", size = 8), 
        legend.title = element_blank(),
        legend.background = element_rect(fill = "#FFFFFF"),
        legend.key = element_rect(fill = "#FFFFFF")) +
  labs(y = "Predicted probability",
       x = "")

#create plot for emotional closeness
discr_me_effects_plot <- me_effects_df %>%
  filter(model  == 3) %>%
  mutate(
    model = case_when(
      model == 3 ~ "Emotional closeness"
    ),
    x = factor(x,
               levels = c(0,1,2),
               labels = c("Not dear",
                          "Dear",
                          "Not asked"))
  ) %>% 
  ggplot(aes(
    x = time,
    y = predicted,
    group = x,
    colour = x,
    fill = x
  )) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.4,
              linetype = "blank") +
  facet_wrap(vars(model)) +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(limits = c(0,1)) +
   scale_fill_manual(values = c("#1b9e77",
                               "#d95f02",
                               "#7570b3")) +
  scale_colour_manual(values = c("#1b9e77",
                               "#d95f02",
                               "#7570b3")) +
  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_text(),
        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 = "right",
        legend.direction = "vertical",
        legend.text = element_text(family = "sans", size = 8), 
        legend.title = element_blank(),
        legend.background = element_rect(fill = "#FFFFFF"),
        legend.key = element_rect(fill = "#FFFFFF")) +
  labs(y = "Predicted probability",
       x = "Period")


me_effects_plot <- age_me_effects_plot +
  degree_me_effects_plot +
  discr_me_effects_plot +
  plot_layout(ncol = 1)

me_effects_plot

ggsave(me_effects_plot,
       file = "plots/results/survival/me_effects_plot.jpg",
       dpi = 320,
       width = 5,
       height = 8)

#1b9e77
#d95f02
#7570b3
#e7298a
#66a61e

Tables

3, 5, 6, 7, 8

Table 2

tab_model(me_model_list[[3]],
          me_model_list[[4]],
          me_model_list[[5]],
          me_model_list[[6]],
          me_model_list[[7]],
          me_model_list[[8]],
          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.859 ***	46.716	1.877 ***	47.520	1.860 ***	46.711	2.825 ***	29.808	3.078 ***	31.283	2.679 ***	27.754
time	0.736 ***	-39.278	0.729 ***	-40.485	0.737 ***	-39.212	0.748 ***	-36.248	0.753 ***	-35.368	0.750 ***	-35.881
dyad_educ_sim_cen	0.990	-1.571			0.987	-1.880	0.993	-0.934	0.992	-1.059	0.993	-1.039
dyad_age_sim_cen	0.923 ***	-11.889			0.920 ***	-11.829	0.968 ***	-4.281	0.964 ***	-4.792	0.968 ***	-4.229
dyad_gender_sim_cen	0.905 ***	-15.339			0.904 ***	-14.880	0.921 ***	-11.681	0.917 ***	-12.307	0.921 ***	-11.789
dyad_ethnicity_sim_cen	0.973 ***	-3.728			0.973 **	-3.134	0.988	-1.101	0.987	-1.106	0.987	-1.136
avsim_alter_educ_cen			1.005	0.798	1.008	1.174	1.010	1.391	1.011	1.538	1.011	1.572
avsim_alter_age_cen			0.985 *	-2.105	1.011	1.499	1.010	1.272	1.018 *	2.279	1.011	1.485
ei_alter_gender_rev_cen			0.976 ***	-3.491	1.005	0.649	1.003	0.367	1.008	1.032	1.005	0.631
ei_alter_ethnicity_rev_cen			0.981 **	-2.810	0.999	-0.157	1.019	1.902	1.018	1.763	1.019	1.928
educ_ego_cen							0.942 ***	-6.416	0.938 ***	-6.867	0.941 ***	-6.542
age_cen							1.034 *	2.069	1.032	1.935	1.033	1.957
gender_fac: Female							0.895 ***	-5.249	0.894 ***	-5.251	0.894 ***	-5.297
gender_fac: Missing							0.679 **	-2.584	0.675 **	-2.611	0.681 *	-2.571
origin_rec_nar_fac: Non-western							1.238 ***	4.884	1.245 ***	4.967	1.236 ***	4.841
origin_rec_nar_fac: Western							1.003	0.076	1.005	0.109	1.002	0.047
origin_rec_nar_fac: Missing							1.011	0.313	1.014	0.415	1.011	0.314
divorced_fac: divorced							0.917	-1.885	0.937	-1.424	0.919	-1.840
divorced_fac: missing							0.949 **	-3.156	0.963 *	-2.209	0.952 **	-2.947
moving_fac: new_residence							1.064	1.229	1.069	1.314	1.064	1.234
moving_fac: new_municipality							1.126	1.693	1.119	1.601	1.126	1.694
moving_fac: missing							0.986	-0.820	0.985	-0.880	0.987	-0.744
first_child_fac: First child born							1.104	0.819	1.098	0.774	1.112	0.882
first_child_fac: Missing							1.010	0.470	1.010	0.501	1.011	0.522
educ_alter_cen							0.999	-0.136	0.998	-0.197	0.997	-0.374
age_alter_cen							1.013	0.868	1.010	0.708	1.010	0.692
gender_alter_fac: Female							0.965 *	-2.022	0.971	-1.662	0.963 *	-2.149
gender_alter_fac: Missing							1.778	1.848	1.791	1.868	1.728	1.755
origin_alter_rec_fac: Non-western							1.110 **	2.682	1.111 **	2.695	1.111 **	2.697
origin_alter_rec_fac: Western							1.072	1.496	1.076	1.566	1.071	1.468
origin_alter_rec_fac: Missing							1.432 **	2.914	1.453 **	2.997	1.411 **	2.787
rel_alter_rec: Same group or club							0.945	-1.513	0.966	-0.924	0.968	-0.858
rel_alter_rec: Neighbour							0.910 **	-2.717	0.926 *	-2.212	0.931 *	-2.074
rel_alter_rec: Friend							0.768 ***	-11.185	0.803 ***	-9.175	0.787 ***	-10.019
rel_alter_rec: Advisor							1.232 **	3.023	1.222 **	2.885	1.239 **	3.103
rel_alter_rec: Other							1.092 *	2.209	1.113 **	2.654	1.113 **	2.668
rel_alter_rec: Missing							0.590 *	-2.357	0.589 *	-2.348	0.611 *	-2.199
times_dropped_earlier_cen							0.895 ***	-15.778	0.895 ***	-15.765	0.896 ***	-15.694
length_fac: 3 - 6 years							0.871 ***	-5.829	0.881 ***	-5.325	0.886 ***	-5.091
length_fac: > 6 years							0.829 ***	-8.356	0.849 ***	-7.242	0.852 ***	-7.074
length_fac: Length missing							1.149	1.449	1.180	1.712	1.178	1.704
size_cen							1.014	1.770	1.014	1.749	1.076 ***	6.711
net_density_cen							0.968 ***	-4.472	0.964 ***	-4.976	1.029 **	2.642
size_net_one							1.124 *	2.484	1.151 **	2.985	1.332 ***	5.565
as.factor(dear_alter_rec)1									0.688 ***	-15.193
as.factor(dear_alter_rec)2									0.828 ***	-10.979
degree_cen											0.907 ***	-7.925
Random Effects
σ2	1.64		1.64		1.64		1.64		1.64		1.64
τ00	0.19 nomem_encr		0.19 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.086 / 0.179		0.078 / 0.172		0.086 / 0.180		0.112 / 0.200		0.115 / 0.207		0.113 / 0.201
Deviance	58349.556		58754.084		58344.430		57294.872		57039.621		57233.590
log-Likelihood	-29174.778		-29377.042		-29172.215		-28647.436		-28519.811		-28616.795
p<0.05   ** p<0.01   *** p<0.001

Table 3

tab_model(
          me_model_list[[10]],
          me_model_list[[11]],
          me_model_list[[12]],
          me_model_list[[13]],
          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.511 ***	27.676	3.511 ***	27.677	3.502 ***	27.611	3.518 ***	27.436
time	0.755 ***	-35.067	0.755 ***	-35.068	0.755 ***	-35.047	0.755 ***	-35.068
avsim_alter_educ_cen	1.012	1.622	1.012	1.670	1.012	1.631	1.012	1.683
avsim_alter_age_cen	1.019 *	2.432	1.019 *	2.455	1.016 *	2.116	1.019 *	2.419
ei_alter_gender_rev_cen	1.009	1.231	1.009	1.217	1.012	1.615	1.009	1.230
ei_alter_ethnicity_rev_cen	1.018	1.788	1.018	1.793	1.018	1.798	1.017	1.463
educ_ego	0.978 ***	-6.909	0.978 ***	-6.958	0.978 ***	-6.868	0.978 ***	-6.954
age_cen	1.031	1.846	1.031	1.871	1.031	1.872	1.031	1.845
gender_fac: Female	0.893 ***	-5.294	0.893 ***	-5.301	0.882 ***	-5.830	0.894 ***	-5.289
gender_fac: Missing	0.676 **	-2.600	0.677 **	-2.595	0.677 **	-2.598	0.677 **	-2.597
origin_rec_nar_fac: Non-western	1.243 ***	4.925	1.244 ***	4.931	1.240 ***	4.871	1.237 ***	4.421
origin_rec_nar_fac: Western	1.003	0.081	1.004	0.089	1.005	0.113	0.995	-0.098
origin_rec_nar_fac: Missing	1.014	0.411	1.014	0.414	1.013	0.366	1.012	0.346
divorced_fac: divorced	0.938	-1.397	0.938	-1.393	0.938	-1.379	0.938	-1.398
divorced_fac: missing	0.966 *	-2.068	0.966 *	-2.061	0.965 *	-2.128	0.966 *	-2.067
moving_fac: new_residence	1.069	1.310	1.069	1.313	1.069	1.319	1.069	1.310
moving_fac: new_municipality	1.119	1.607	1.119	1.606	1.123	1.646	1.119	1.603
moving_fac: missing	0.986	-0.812	0.986	-0.809	0.988	-0.732	0.986	-0.811
first_child_fac: First child born	1.106	0.831	1.106	0.830	1.103	0.808	1.106	0.832
first_child_fac: Missing	1.011	0.544	1.011	0.534	1.011	0.535	1.011	0.540
educ_alter_cen	0.997	-0.407	0.997	-0.393	0.997	-0.337	0.997	-0.397
age_alter_cen	1.008	0.561	1.008	0.554	1.008	0.532	1.008	0.558
gender_alter_fac: Female	0.969	-1.783	0.969	-1.773	0.977	-1.313	0.969	-1.779
gender_alter_fac: Missing	1.747	1.788	1.746	1.786	1.781	1.848	1.746	1.786
origin_alter_rec_fac: Non-western	1.112 **	2.704	1.112 **	2.700	1.110 **	2.654	1.115 **	2.700
origin_alter_rec_fac: Western	1.075	1.541	1.075	1.547	1.075	1.536	1.078	1.570
origin_alter_rec_fac: Missing	1.434 **	2.887	1.433 **	2.883	1.429 **	2.857	1.433 **	2.885
as.factor(dear_alter_rec)1	0.695 ***	-14.719	0.695 ***	-14.728	0.693 ***	-14.824	0.695 ***	-14.722
as.factor(dear_alter_rec)2	0.834 ***	-10.516	0.834 ***	-10.527	0.831 ***	-10.686	0.834 ***	-10.517
rel_alter_rec: Same group or club	0.985	-0.386	0.985	-0.385	0.982	-0.483	0.985	-0.389
rel_alter_rec: Neighbour	0.943	-1.678	0.943	-1.680	0.941	-1.741	0.943	-1.680
rel_alter_rec: Friend	0.819 ***	-8.262	0.819 ***	-8.266	0.819 ***	-8.296	0.819 ***	-8.266
rel_alter_rec: Advisor	1.229 **	2.962	1.229 **	2.958	1.232 **	3.002	1.229 **	2.960
rel_alter_rec: Other	1.130 **	3.030	1.130 **	3.018	1.129 **	3.016	1.130 **	3.028
rel_alter_rec: Missing	0.607 *	-2.213	0.608 *	-2.207	0.601 *	-2.261	0.608 *	-2.210
times_dropped_earlier_cen	0.895 ***	-15.689	0.895 ***	-15.690	0.895 ***	-15.748	0.896 ***	-15.688
length_fac: 3 - 6 years	0.894 ***	-4.710	0.894 ***	-4.705	0.892 ***	-4.762	0.894 ***	-4.709
length_fac: > 6 years	0.868 ***	-6.189	0.868 ***	-6.190	0.867 ***	-6.238	0.868 ***	-6.187
length_fac: Length missing	1.204	1.917	1.204	1.920	1.215 *	2.016	1.204	1.915
degree_cen	0.919 ***	-6.738	0.919 ***	-6.742	0.918 ***	-6.832	0.919 ***	-6.742
dyad_educ_sim_cen	0.992	-1.121	0.992	-1.134	0.992	-1.095	0.992	-1.142
dyad_age_sim_cen	0.964 ***	-4.728	0.966 ***	-3.917	0.965 ***	-4.665	0.965 ***	-4.721
dyad_gender_sim_cen	0.916 ***	-12.379	0.916 ***	-12.378	0.932 ***	-8.314	0.916 ***	-12.377
dyad_ethnicity_sim_cen	0.987	-1.136	0.987	-1.134	0.987	-1.168	0.984	-1.014
size_cen	1.067 ***	5.883	1.067 ***	5.895	1.067 ***	5.834	1.067 ***	5.886
net_density_cen	1.015	1.419	1.016	1.427	1.018	1.633	1.016	1.423
size_net_one	1.330 ***	5.509	1.332 ***	5.532	1.340 ***	5.650	1.329 ***	5.494
avsim_alter_educ_cen:dyad_educ_sim_cen	0.999	-0.121
avsim_alter_age_cen:dyad_age_sim_cen			1.002	0.410
ei_alter_gender_rev_cen:dyad_gender_sim_cen					1.026 ***	3.499
ei_alter_ethnicity_rev_cen:dyad_ethnicity_sim_cen							0.998	-0.281
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.208		0.116 / 0.208
Deviance	56995.366		56995.230		56983.839		56995.301
log-Likelihood	-28497.683		-28497.615		-28491.919		-28497.650
p<0.05   ** p<0.01   *** p<0.001

Table 4

time_int_table <- time_int_model_df %>% 
  filter(m_number %in% c(1:6)) %>% 
  mutate(m_number = case_when(
    m_number == 1 ~ 4,
    m_number == 2 ~ 3,
    m_number == 3 ~ 2,
    m_number == 4 ~ 1,
    m_number == 5 ~ 5,
    m_number == 6 ~ 6,
  )
  ) %>% 
  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,
         riskratio_6,  significant_6,statistic_6,
         riskratio_5,  significant_5,statistic_5) %>% 
  filter(parameter_fac %in% c(
    "Constant",
    "Linear time",
    "Closeness: dear",
      "Closeness: not asked",
      "Alter Embeddedness"
  ) |
    str_detect(parameter_fac, "Interaction")|
    str_detect(parameter_fac, "Dyad")) %>% 
  arrange(as.numeric(parameter_fac))
  
write_csv(time_int_table,
          file = "results/survival_results/tables/2023-06-14_table4.csv")

tab_model(time_int_model_list[[4]],
          time_int_model_list[[3]],
          time_int_model_list[[2]],
          time_int_model_list[[1]],
          time_int_model_list[[5]],
          time_int_model_list[[6]],
          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)	3.513 ***	27.695	3.516 ***	27.705	3.516 ***	27.709	3.516 ***	27.707	3.503 ***	27.618	3.342 ***	24.939
time	0.755 ***	-35.075	0.753 ***	-35.087	0.753 ***	-34.903	0.754 ***	-35.037	0.754 ***	-35.109	0.784 ***	-16.114
avsim_alter_educ_cen	1.012	1.681	1.012	1.674	1.012	1.684	1.012	1.684	1.013	1.706	1.012	1.694
avsim_alter_age_cen	1.019 *	2.430	1.019 *	2.463	1.019 *	2.425	1.019 *	2.441	1.019 *	2.457	1.019 *	2.424
ei_alter_gender_rev_cen	1.009	1.230	1.009	1.249	1.009	1.240	1.009	1.220	1.009	1.245	1.009	1.245
ei_alter_ethnicity_rev_cen	1.018	1.787	1.018	1.797	1.018	1.793	1.018	1.808	1.018	1.803	1.018	1.791
educ_ego	0.978 ***	-6.972	0.978 ***	-6.963	0.978 ***	-6.937	0.978 ***	-6.965	0.978 ***	-6.959	0.978 ***	-6.933
age_cen	1.031	1.844	1.031	1.831	1.031	1.834	1.031	1.855	1.031	1.863	1.031	1.861
gender_fac: Female	0.893 ***	-5.294	0.894 ***	-5.289	0.894 ***	-5.280	0.894 ***	-5.283	0.893 ***	-5.294	0.894 ***	-5.288
gender_fac: Missing	0.677 **	-2.595	0.677 **	-2.598	0.677 **	-2.594	0.677 **	-2.596	0.676 **	-2.604	0.677 **	-2.596
origin_rec_nar_fac: Non-western	1.243 ***	4.921	1.243 ***	4.923	1.243 ***	4.925	1.242 ***	4.895	1.242 ***	4.908	1.242 ***	4.908
origin_rec_nar_fac: Western	1.003	0.079	1.004	0.089	1.003	0.078	1.004	0.104	1.002	0.058	1.004	0.088
origin_rec_nar_fac: Missing	1.014	0.409	1.014	0.410	1.014	0.404	1.014	0.417	1.014	0.413	1.014	0.410
divorced_fac: divorced	0.938	-1.398	0.937	-1.408	0.937	-1.403	0.937	-1.407	0.937	-1.406	0.937	-1.412
divorced_fac: missing	0.966 *	-2.079	0.966 *	-2.031	0.966 *	-2.042	0.966 *	-2.078	0.966 *	-2.036	0.966 *	-2.072
moving_fac: new_residence	1.069	1.315	1.068	1.287	1.069	1.311	1.069	1.319	1.068	1.287	1.069	1.310
moving_fac: new_municipality	1.119	1.602	1.118	1.587	1.118	1.586	1.119	1.599	1.118	1.587	1.121	1.626
moving_fac: missing	0.986	-0.810	0.985	-0.854	0.986	-0.849	0.986	-0.833	0.986	-0.817	0.986	-0.839
first_child_fac: First child born	1.106	0.837	1.106	0.834	1.106	0.833	1.106	0.830	1.105	0.828	1.107	0.843
first_child_fac: Missing	1.011	0.543	1.012	0.566	1.012	0.551	1.011	0.544	1.012	0.564	1.011	0.538
educ_alter_cen	0.997	-0.381	0.997	-0.442	0.997	-0.411	0.997	-0.393	0.997	-0.434	0.997	-0.393
age_alter_cen	1.008	0.554	1.008	0.567	1.008	0.572	1.008	0.547	1.008	0.522	1.008	0.545
gender_alter_fac: Female	0.969	-1.785	0.968	-1.794	0.968	-1.805	0.969	-1.789	0.969	-1.787	0.969	-1.789
gender_alter_fac: Missing	1.747	1.789	1.746	1.786	1.741	1.777	1.746	1.786	1.742	1.778	1.741	1.782
origin_alter_rec_fac: Non-western	1.112 **	2.707	1.113 **	2.731	1.112 **	2.699	1.113 **	2.720	1.113 **	2.725	1.112 **	2.699
origin_alter_rec_fac: Western	1.075	1.543	1.076	1.552	1.075	1.535	1.075	1.549	1.075	1.548	1.075	1.541
origin_alter_rec_fac: Missing	1.434 **	2.889	1.431 **	2.875	1.435 **	2.893	1.434 **	2.888	1.431 **	2.876	1.436 **	2.897
as.factor(dear_alter_rec)1	0.695 ***	-14.718	0.695 ***	-14.700	0.695 ***	-14.710	0.695 ***	-14.726	0.696 ***	-14.683	0.757 ***	-6.177
as.factor(dear_alter_rec)2	0.834 ***	-10.519	0.834 ***	-10.514	0.834 ***	-10.518	0.834 ***	-10.510	0.834 ***	-10.474	0.890 ***	-3.985
rel_alter_rec: Same group or club	0.986	-0.380	0.985	-0.397	0.985	-0.387	0.985	-0.387	0.986	-0.364	0.985	-0.386
rel_alter_rec: Neighbour	0.943	-1.669	0.944	-1.657	0.943	-1.686	0.943	-1.684	0.944	-1.650	0.943	-1.682
rel_alter_rec: Friend	0.819 ***	-8.255	0.820 ***	-8.244	0.819 ***	-8.262	0.819 ***	-8.268	0.821 ***	-8.192	0.819 ***	-8.303
rel_alter_rec: Advisor	1.230 **	2.970	1.226 **	2.920	1.228 **	2.951	1.229 **	2.960	1.226 **	2.931	1.231 **	2.988
rel_alter_rec: Other	1.130 **	3.040	1.127 **	2.964	1.130 **	3.023	1.129 **	3.015	1.131 **	3.043	1.130 **	3.030
rel_alter_rec: Missing	0.607 *	-2.214	0.610 *	-2.195	0.610 *	-2.194	0.608 *	-2.209	0.612 *	-2.178	0.606 *	-2.226
times_dropped_earlier_cen	0.895 ***	-15.692	0.896 ***	-15.686	0.896 ***	-15.677	0.895 ***	-15.697	0.896 ***	-15.684	0.895 ***	-15.697
length_fac: 3 - 6 years	0.894 ***	-4.712	0.894 ***	-4.671	0.894 ***	-4.668	0.894 ***	-4.699	0.896 ***	-4.603	0.893 ***	-4.737
length_fac: > 6 years	0.868 ***	-6.184	0.869 ***	-6.146	0.869 ***	-6.150	0.869 ***	-6.171	0.871 ***	-6.056	0.867 ***	-6.229
length_fac: Length missing	1.204	1.922	1.203	1.914	1.203	1.911	1.203	1.914	1.206	1.937	1.201	1.890
degree_cen	0.919 ***	-6.734	0.919 ***	-6.747	0.919 ***	-6.737	0.919 ***	-6.744	0.893 ***	-6.984	0.919 ***	-6.802
size_cen	1.067 ***	5.876	1.067 ***	5.892	1.067 ***	5.865	1.067 ***	5.888	1.065 ***	5.714	1.068 ***	5.933
net_density_cen	1.015	1.413	1.016	1.425	1.016	1.424	1.016	1.424	1.015	1.398	1.016	1.426
dyad_educ_sim_cen	1.000	0.006	0.992	-1.129	0.992	-1.135	0.992	-1.155	0.992	-1.137	0.992	-1.131
dyad_age_sim_cen	0.964 ***	-4.732	0.938 ***	-4.785	0.965 ***	-4.716	0.964 ***	-4.732	0.964 ***	-4.744	0.965 ***	-4.718
dyad_gender_sim_cen	0.916 ***	-12.384	0.916 ***	-12.377	0.898 ***	-8.156	0.916 ***	-12.373	0.916 ***	-12.387	0.916 ***	-12.378
dyad_ethnicity_sim_cen	0.987	-1.135	0.987	-1.135	0.987	-1.139	0.973	-1.754	0.987	-1.159	0.987	-1.130
size_net_one	1.331 ***	5.519	1.330 ***	5.506	1.330 ***	5.510	1.331 ***	5.520	1.325 ***	5.431	1.333 ***	5.552
time:dyad_educ_sim_cen	0.994	-0.822
time:dyad_age_sim_cen			1.022 *	2.520
time:dyad_gender_sim_cen					1.016	1.796
time:dyad_ethnicity_sim_cen							1.011	1.347
time:degree_cen									1.023 **	2.809
time:as.factor(dear_alter_rec)1											0.940 *	-2.292
time:as.factor(dear_alter_rec)2											0.953 **	-2.706
Random Effects
σ2	1.64		1.64		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		0.19 nomem_encr		0.19 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.116 / 0.208		0.116 / 0.207		0.115 / 0.207		0.116 / 0.207		0.116 / 0.207		0.117 / 0.209
Deviance	56994.757		56988.163		56991.423		56993.411		56986.896		56988.199
log-Likelihood	-28497.379		-28494.081		-28495.712		-28496.706		-28493.448		-28494.100
p<0.05   ** p<0.01   *** p<0.001

Table B.1

tab_model(me_model_list[[6]],
          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.825 ***	29.808
time	0.748 ***	-36.248
avsim_alter_educ_cen	1.010	1.391
avsim_alter_age_cen	1.010	1.272
ei_alter_gender_rev_cen	1.003	0.367
ei_alter_ethnicity_rev_cen	1.019	1.902
educ_ego_cen	0.942 ***	-6.416
age_cen	1.034 *	2.069
gender_fac: Female	0.895 ***	-5.249
gender_fac: Missing	0.679 **	-2.584
origin_rec_nar_fac: Non-western	1.238 ***	4.884
origin_rec_nar_fac: Western	1.003	0.076
origin_rec_nar_fac: Missing	1.011	0.313
divorced_fac: divorced	0.917	-1.885
divorced_fac: missing	0.949 **	-3.156
moving_fac: new_residence	1.064	1.229
moving_fac: new_municipality	1.126	1.693
moving_fac: missing	0.986	-0.820
first_child_fac: First child born	1.104	0.819
first_child_fac: Missing	1.010	0.470
educ_alter_cen	0.999	-0.136
age_alter_cen	1.013	0.868
gender_alter_fac: Female	0.965 *	-2.022
gender_alter_fac: Missing	1.778	1.848
origin_alter_rec_fac: Non-western	1.110 **	2.682
origin_alter_rec_fac: Western	1.072	1.496
origin_alter_rec_fac: Missing	1.432 **	2.914
rel_alter_rec: Same group or club	0.945	-1.513
rel_alter_rec: Neighbour	0.910 **	-2.717
rel_alter_rec: Friend	0.768 ***	-11.185
rel_alter_rec: Advisor	1.232 **	3.023
rel_alter_rec: Other	1.092 *	2.209
rel_alter_rec: Missing	0.590 *	-2.357
times_dropped_earlier_cen	0.895 ***	-15.778
length_fac: 3 - 6 years	0.871 ***	-5.829
length_fac: > 6 years	0.829 ***	-8.356
length_fac: Length missing	1.149	1.449
dyad_educ_sim_cen	0.993	-0.934
dyad_age_sim_cen	0.968 ***	-4.281
dyad_gender_sim_cen	0.921 ***	-11.681
dyad_ethnicity_sim_cen	0.988	-1.101
size_cen	1.014	1.770
net_density_cen	0.968 ***	-4.472
size_net_one	1.124 *	2.484
Random Effects
σ2	1.64
τ00 nomem_encr	0.18
ICC	0.10
N nomem_encr	6996
Observations	49449
Marginal R2 / Conditional R2	0.112 / 0.200
Deviance	57294.872
log-Likelihood	-28647.436
p<0.05   ** p<0.01   *** p<0.001

GOF fit table B.2

#create list with models.
models_list <- list(
  me_model_list[[1]],
  me_model_list[[2]],
  me_model_list[[3]],
  me_model_list[[4]],
  me_model_list[[5]],
  me_model_list[[6]],
  me_model_list[[7]],
  me_model_list[[8]],
  me_model_list[[9]],
  me_model_list[[10]],
  me_model_list[[11]],
  time_int_model_list[[4]],
  time_int_model_list[[3]],
  time_int_model_list[[2]],
  time_int_model_list[[1]],
  time_int_model_list[[5]],
  time_int_model_list[[6]]
)


#create fit list
fit_list <- list()

#loop
for(i in 1:length(models_list)) {#i = 1
  model <- models_list[[i]] #extract correct model
  fit <- llikAIC(model)
  fit_table <- as_tibble(fit$AICtab) %>% 
    transpose() %>% 
    rename(AIC = V1,
           BIC = V2,
           LogLikelihood = V3,
           Deviance = V4,
           DF = V5) %>% 
    mutate(model = paste0("Model ", i)) #create model collumn
  
  fit_list[[i]] <- fit_table %>%
    select(model,AIC,BIC,LogLikelihood, Deviance, DF) #export fit model
}

gof_table <- fit_list %>%
  rbindlist()

#fit table
fit_list %>%
  rbindlist() %>%
  kbl(caption = "Table 3. Fit statistics of discrete time hazard models 1-11", digits = 3) %>%
  kable_classic(full_width = F, bootstrap_options = c("hover", "condensed"), fixed_thead = T)

Table 3. Fit statistics of discrete time hazard models 1-11

model	AIC	BIC	LogLikelihood	Deviance	DF
Model 1	60518.26	60535.88	-30257.13	60514.26	49447
Model 2	58785.90	58812.33	-29389.95	58779.90	49446
Model 3	58363.56	58425.22	-29174.78	58349.56	49442
Model 4	58768.08	58829.75	-29377.04	58754.08	49442
Model 5	58366.43	58463.33	-29172.22	58344.43	49438
Model 6	57384.87	57781.26	-28647.44	57294.87	49404
Model 7	57133.62	57547.63	-28519.81	57039.62	49402
Model 8	57325.59	57730.79	-28616.79	57233.59	49403
Model 9	57091.39	57514.20	-28497.69	56995.39	49401
Model 10	57093.37	57524.99	-28497.68	56995.37	49400
Model 11	57093.23	57524.86	-28497.62	56995.23	49400
Model 12	57092.76	57524.38	-28497.38	56994.76	49400
Model 13	57086.16	57517.79	-28494.08	56988.16	49400
Model 14	57089.42	57521.05	-28495.71	56991.42	49400
Model 15	57091.41	57523.04	-28496.71	56993.41	49400
Model 16	57084.90	57516.52	-28493.45	56986.90	49400
Model 17	57088.20	57528.63	-28494.10	56988.20	49399

#export GOF fit list. 
write_csv(gof_table, file = "results/survival_results/tables/2023-07-07_goftable.csv")

Correlation plot

Dyadic similarity and Confidant uniqueness

All variables

#create df
dataframe <- MyData %>% 
  select(
    dyad_educ_sim_cen,
    dyad_age_sim_cen,
    dyad_gender_sim_cen,
    dyad_ethnicity_sim_cen,
    avsim_alter_educ_cen,
    avsim_alter_age_cen,
    ei_alter_gender_rev_cen,
    ei_alter_ethnicity_rev_cen,
    degree_cen,
    dear_alter_rec,
    educ_ego_cen,
    age_cen,
    gender_fac,
    origin_rec_nar,
    divorced_fac,
    moving_fac,
    first_child,
    educ_alter_cen,
    age_alter_cen,
    gender_alter_fac,
    origin_alter_rec,
    rel_alter_rec,
    length,
    times_dropped_earlier_cen,
    size_cen,
    net_density_cen) %>% 
 mutate(
    no_mig = ifelse(origin_rec_nar == 0, 1, 0),
    nonwestern_mig = ifelse(origin_rec_nar == 1, 1, 0),
    western_mig =  ifelse(origin_rec_nar == 2, 1, 0),
    no_mig_a =  ifelse(origin_alter_rec == 0, 1, 0),
    nonwestern_mig_a =  ifelse(origin_alter_rec == 1, 1, 0),
    western_mig_a =  ifelse(origin_alter_rec == 2, 1, 0),
    rel_a_b3 = ifelse(length == 1, 1, 0),
    rel_a_36 = ifelse(length == 2, 1, 0),
    rel_a_o6 = ifelse(length == 3, 1, 0),
    rel_missing = ifelse(is.na(length), 1, 0),
    coll = ifelse(as.numeric(rel_alter_rec) == 4, 1, 0),
    samegroup = ifelse(as.numeric(rel_alter_rec) == 5, 1, 0),
    neighbour = ifelse(as.numeric(rel_alter_rec) == 6, 1, 0),
    friend = ifelse(as.numeric(rel_alter_rec) == 7, 1, 0),
    advisor = ifelse(as.numeric(rel_alter_rec) == 8, 1, 0),
    otherrel = ifelse(as.numeric(rel_alter_rec) == 9, 1, 0),
    not_dear = ifelse(as.numeric(dear_alter_rec) == 0, 1, 0),
    dear = ifelse(as.numeric(dear_alter_rec) == 1, 1, 0),
    not_asked = ifelse(as.numeric(dear_alter_rec) == 2, 1, 0),
    not_divorced = ifelse(as.numeric(divorced_fac) == 1, 1, 0),
    divorced = ifelse(as.numeric(divorced_fac) == 2, 1, 0),
    divorced_missing = ifelse(as.numeric(divorced_fac) == 3, 1, 0),
    no_move = ifelse(as.numeric(moving_fac) == 1, 1, 0),
    new_res = ifelse(as.numeric(moving_fac) == 2, 1, 0),
    new_mun = ifelse(as.numeric(moving_fac) == 3, 1, 0),
    move_mis = ifelse(as.numeric(moving_fac) == 4, 1, 0),
    female = ifelse(as.numeric(gender_fac) == 2, 1, 0),
    first_child_rec = ifelse(first_child == 1, 1, 0),
    first_child_missing = ifelse(first_child == 2, 1, 0),
    female_alter = ifelse(gender_alter_fac == "Female", 1, 0),
    male_alter = ifelse(gender_alter_fac == "Male", 1, 0),
    missing_gender_alter = ifelse(gender_alter_fac == "Missing", 1, 0)
  ) %>% 
  select(
    dyad_educ_sim_cen,
    dyad_gender_sim_cen,
    dyad_age_sim_cen,
    dyad_ethnicity_sim_cen,
    avsim_alter_educ_cen,
    avsim_alter_age_cen,
    ei_alter_gender_rev_cen,
    ei_alter_ethnicity_rev_cen,
    dear,
    not_asked,
    degree_cen,
    educ_ego_cen,
    age_cen,
    female,
    nonwestern_mig,
    western_mig,
    divorced,
    new_res,
    new_mun,
    first_child_rec,
    educ_alter_cen,
    age_alter_cen,
    female_alter,
    nonwestern_mig_a,
    western_mig_a,
    samegroup,
    neighbour,
    friend,
    advisor,
    otherrel,
    times_dropped_earlier_cen,
    rel_a_36,
    rel_a_o6,
    net_density_cen,
    size_cen
  )

#crete lower tri identifier
lower_tri_identifier <- dataframe %>%
  cor(., use = "complete.obs") %>% 
  lower.tri() %>%
  as_tibble() %>%
  mutate(variable_x = 1:35) %>%
  pivot_longer(1:35,
               names_to = "variable_y",
               values_to = "identifier") %>%
  select(identifier)


#create plot df
plot_df <- dataframe %>%
  cor(., use = "complete.obs") %>%
  as_tibble() %>%
  mutate(variable_x = factor(
    1:35,
    levels = 1:35,
    labels = c(
      "Dyadic similarity: education",
      "Dyadic similarity: gender",
      "Dyadic similarity: age",
      "Dyadic similarity: ethnicity",
      "Confidant uniqueness education",
      "Confidant uniquenessy age",
      "Confidant uniqueness gender",
      "Confidant uniqueness ethnicity",
      "Confidant is dear",
      "Closeness not asked",
      "Embeddedness",
      "Education ego",
      "Age",
      "Female",
      "Non-western migration background (ref. Dutch)",
      "Western migration background (ref. Dutch)",
      "Divorced (ref. not)",
      "New residence (ref. no move)",
      "New municipality (ref. no move)",
      "First child born",
      "Education alter",
      "Age alter",
      "Alter is female",
      "Non-western migration background confidant (ref. Dutch)",
      "Western migration background confidant (ref. Dutch)",
      "Same group or club",
      "Neighbour",
      "Friend",
      "Advisor",
      "Other relation",
      "Times dropped earlier",
      "Knows confidant for 3-6 years",
      "Knows confidant  > 6 years",
      "Net density",
      "Net size"
    )
  )) %>%
  tidyr::pivot_longer(1:35,
                      names_to = "variable_y",
                      values_to = "value") %>%
  mutate(
    variable_y = case_when(
      variable_y == "dyad_educ_sim_cen" ~ 1,
      variable_y == "dyad_gender_sim_cen" ~ 2,
      variable_y == "dyad_age_sim_cen" ~ 3,
      variable_y == "dyad_ethnicity_sim_cen" ~ 4,
      variable_y == "avsim_alter_educ_cen" ~ 5,
      variable_y == "avsim_alter_age_cen" ~ 6,
      variable_y == "ei_alter_gender_rev_cen" ~ 7,
      variable_y == "ei_alter_ethnicity_rev_cen" ~ 8,
      variable_y == "dear" ~ 9,
      variable_y == "not_asked" ~ 10,
      variable_y == "degree_cen" ~ 11,
      variable_y == "educ_ego_cen" ~ 12,
      variable_y == "age_cen" ~ 13,
      variable_y == "female" ~ 14,
      variable_y == "nonwestern_mig" ~ 15,
      variable_y == "western_mig" ~ 16,
      variable_y == "divorced" ~ 17,
      variable_y == "new_res" ~ 18,
      variable_y == "new_mun" ~ 19,
      variable_y == "first_child_rec" ~ 20,
      variable_y == "educ_alter_cen" ~ 21,
      variable_y == "age_alter_cen" ~ 22,
      variable_y == "female_alter" ~ 23,
      variable_y == "nonwestern_mig_a" ~ 24,
      variable_y == "western_mig_a" ~ 25,
      variable_y == "samegroup" ~ 26,
      variable_y == "neighbour" ~ 27,
      variable_y == "friend" ~ 28,
      variable_y == "advisor" ~ 29,
      variable_y == "otherrel" ~ 30,
      variable_y == "times_dropped_earlier_cen" ~ 31,
      variable_y == "rel_a_36" ~ 32,
      variable_y == "rel_a_o6" ~ 33,
      variable_y == "net_density_cen" ~ 34,
      variable_y == "size_cen" ~ 35
    ),
    variable_y = factor(
      variable_y,
      levels = 1:35,
      labels = c(
        "Dyadic similarity: education",
        "Dyadic similarity: gender",
        "Dyadic similarity: age",
        "Dyadic similarity: ethnicity",
        "Confidant uniqueness education",
        "Confidant uniquenessy age",
        "Confidant uniqueness gender",
        "Confidant uniqueness ethnicity",
        "Confidant is dear",
        "Closeness not asked",
        "Embeddedness",
        "Education ego",
        "Age",
        "Female",
        "Non-western migration background (ref. Dutch)",
        "Western migration background (ref. Dutch)",
        "Divorced (ref. not)",
        "New residence (ref. no move)",
        "New municipality (ref. no move)",
        "First child born",
        "Education alter",
        "Age alter",
        "Alter is female",
        "Non-western migration background confidant (ref. Dutch)",
        "Western migration background confidant (ref. Dutch)",
        "Same group or club",
        "Neighbour",
        "Friend",
        "Advisor",
        "Other relation",
        "Times dropped earlier",
        "Knows confidant for 3-6 years",
        "Knows confidant  > 6 years",
        "Net density",
        "Net size"
      )
    )
  ) %>% 
    bind_cols(lower_tri_identifier) %>%
    mutate(value = ifelse(identifier == FALSE, NA, value))

#correlation plot for all variables
correlation_plot_all <- plot_df %>% 
  mutate(value_abs = abs(value)) %>% 
ggplot(aes(x = variable_x, y = variable_y, 
           colour = value,
           size = value_abs)) +
  geom_point(alpha = 0.8,
            na.rm = T,) +
  theme_minimal() +
  scale_colour_viridis(
    option = "D",
    na.value = "white",
    begin = 0,
    end = 1,
    direction = -1
  ) +
  theme(
    axis.text.x = element_text(angle = 90,
                               vjust = 0.4),
    panel.background = element_rect(colour = "black"),
    axis.title.y = element_text(hjust = 1),
    axis.title.x = element_text(hjust = 1),
    text = element_text(colour = "black"),
    strip.text = element_text(colour = "black"),
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_line(linetype = "dotted"),
    panel.grid.major.x = element_line(linetype = "dotted"),
    legend.position = "right"
  ) +
  scale_size(guide = 'none') +
  labs(x = "",
       y = "",
       size = "Correlation",
       colour = "Correlation")

correlation_plot_all

ggsave(correlation_plot_all,
       file = "plots/results/survival/correlation_plot_all.jpg",
       dpi = 320,
       width = 10,
       height = 10)

Dyadic similarity and confidant uniqueness

#create df
dataframe_indep <- dataframe %>% 
  select(
     dyad_educ_sim_cen,
    dyad_gender_sim_cen,
    dyad_age_sim_cen,
    dyad_ethnicity_sim_cen,
    avsim_alter_educ_cen,
    avsim_alter_age_cen,
    ei_alter_gender_rev_cen,
    ei_alter_ethnicity_rev_cen,
    dear,
    not_asked,
    degree_cen)


#crete lower tri identifier
lower_tri_identifier <- dataframe_indep %>%
  cor(., use = "complete.obs") %>% 
  lower.tri() %>%
  as_tibble() %>%
  mutate(variable_x = 1:11) %>%
  pivot_longer(1:11,
               names_to = "variable_y",
               values_to = "identifier") %>%
  select(identifier)


#create plot df
indep_plot_df <- dataframe_indep %>%
  cor(., use = "complete.obs") %>%
  as_tibble() %>%
  mutate(variable_x = factor(
    1:11,
    levels = 1:11,
    labels = c(
      "Dyadic similarity: education",
      "Dyadic similarity: gender",
      "Dyadic similarity: age",
      "Dyadic similarity: ethnicity",
      "Confidant uniqueness education",
      "Confidant uniquenessy age",
      "Confidant uniqueness gender",
      "Confidant uniqueness ethnicity",
      "Confidant is dear",
      "Closeness not asked",
      "Embeddedness"
    )
  )) %>%
  tidyr::pivot_longer(1:11,
                      names_to = "variable_y",
                      values_to = "value") %>%
  mutate(
    variable_y = case_when(
      variable_y == "dyad_educ_sim_cen" ~ 1,
      variable_y == "dyad_gender_sim_cen" ~ 2,
      variable_y == "dyad_age_sim_cen" ~ 3,
      variable_y == "dyad_ethnicity_sim_cen" ~ 4,
      variable_y == "avsim_alter_educ_cen" ~ 5,
      variable_y == "avsim_alter_age_cen" ~ 6,
      variable_y == "ei_alter_gender_rev_cen" ~ 7,
      variable_y == "ei_alter_ethnicity_rev_cen" ~ 8,
      variable_y == "dear" ~ 9,
      variable_y == "not_asked" ~ 10,
      variable_y == "degree_cen" ~ 11
    ),
    variable_y = factor(
      variable_y,
      levels = 1:11,
      labels = c(
        "Dyadic similarity: education",
        "Dyadic similarity: gender",
        "Dyadic similarity: age",
        "Dyadic similarity: ethnicity",
        "Confidant uniqueness education",
        "Confidant uniquenessy age",
        "Confidant uniqueness gender",
        "Confidant uniqueness ethnicity",
        "Confidant is dear",
        "Closeness not asked",
        "Embeddedness"
      )
    )
  ) %>% 
    bind_cols(lower_tri_identifier) %>%
    mutate(value = ifelse(identifier == FALSE, NA, value))

#correlation plot for all variables
correlation_plot_indep <- indep_plot_df %>% 
  mutate(value_abs = abs(value)) %>% 
ggplot(aes(x = variable_x, y = variable_y, 
           colour = value,
           size = value_abs)) +
  geom_point(alpha = 0.8,
            na.rm = T,) +
  geom_text(aes(label = round(value, 2),
                na.rm = T),
            size = 2,
            colour = "black") +
  theme_minimal() +
  scale_colour_viridis(
    option = "D",
    na.value = "white",
    begin = 0,
    end = 1,
    direction = -1
  ) +
  theme(
    axis.line = element_blank(),
    axis.text.x = element_text(angle = 90,
                               vjust = 0.4),
    panel.background = element_rect(colour = "black"),
    axis.title.y = element_text(hjust = 1),
    axis.title.x = element_text(hjust = 1),
    text = element_text(colour = "black"),
    strip.text = element_text(colour = "black"),
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_line(linetype = "dotted"),
    panel.grid.major.x = element_line(linetype = "dotted"),
    legend.position = "right"
  ) +
  scale_size(guide = 'none') +
  labs(x = "",
       y = "",
       size = "Correlation",
       colour = "Correlation")

correlation_plot_indep

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

---

Copyright © 2023 Jeroense Thijmen