Dependent variables: Revision

Last compiled on May, 2025

---

This is the code with which we make our dependent variables according to the revision of our manuscript.

---

1 Initatiating R environment

Start out with a custom function to load a set of required packages.

# packages and read data
rm(list = ls())

# (c) Jochem Tolsma
fpackage.check <- function(packages) {
    lapply(packages, FUN = function(x) {
        if (!require(x, character.only = TRUE)) {
            install.packages(x, dependencies = TRUE)
            library(x, character.only = TRUE)
        }
    })
}
packages = c("haven", "coda", "matrixStats", "parallel", "MASS", "doParallel", "dplyr", "cowplot", "tidyverse",
    "naniar", "dotwhisker", "gt", "reshape2", "VGAM", "expss", "networkscaleup")
fpackage.check(packages)
rm(packages)
load("data/dutch_netsize_analyses.rda")

Revision dependent variable for network size

dfans <- df[, c("uni", "hbo", "mbo", "dochterzoon", "tweeling", "corona", 
              "elauto", "scooter", "vegan", 
              "Sophie", "Julia","Sanne","Lisa","Laura","Maria","Linda","Johanna","Monique","Ester",
              "Anna","Elisabeth","Cornelia","Wilhelmina","Amira","Samira","Sara","Daan","Sem",
              "Thomas","Max","Kevin","Johannes","Dennis","Jeroen","Jan","Marcel","Cornelis","Hendrik","Petrus",
              "Willem","Ali","Mohammed","Noor")] # take mean of each of the cats

ref1 <- c(84957, 75214, 145600, 168066, 2500, 1558549)
# ORDER:  uni, hbo, mbo, dochter/zoon, tweeling, corona

ref2 <- c(273259, 460618, 261000)
# ORDER:   elecauto, scooter, vegan,

ref3 <- c(15276,16350,27394,21200,25681,334502,29955,266522,39481, 2692,
          136296,110231,112807,98208,1386,2186,11640,22704,13276,
          40543,17024,23167,307032,36411,49182,186746,35973,134956,118610,86500,
          102296,4213,5003,4517)
# ORDER:   Sophie, Julia,Sanne,Lisa,Laura,Maria,Linda,Johanna,Monique,Ester,
#             Anna,Elisabeth,Cornelia,Wilhelmina,Amira,Samira,Sara,Daan,Sem,
#             Thomas,Max,Kevin,Johannes,Dennis,Jeroen,Jan,Marcel,Cornelis,Hendrik,Petrus,
#             Willem,Ali,Mohammed,Noor
pops <- c(ref1, ref2, ref3)

# Analyze an example ard data set using Zheng et al. (2006) models
ard <- as.matrix(dfans)
subpop_sizes <- pops
known_ind <- c(1:43)
N <- 17407585

# first with primary scaling groups ALL
overdisp.est1 <- overdispersed(ard,
                            known_sizes = subpop_sizes[known_ind],
                            known_ind = known_ind,
                            G1_ind = known_ind,
                            N = N,
                            warmup = 500,
                            iter = 5000)

# primary scaling groups all, except some
overdisp.est2 <- overdispersed(ard,
                            known_sizes = subpop_sizes[known_ind],
                            known_ind = known_ind,
                            G2_ind = c(2, 3, 4, 7, 9, 35), #hbo, mbo,dochterzoon, elauto,vegan, jan
                            B2_ind = c(15, 17, 32), # maria, johanna, johannes
                            N = N,
                            warmup = 500,
                            iter = 5000)
 
# regularize subpopulations to < 10k
z <- which(subpop_sizes < 10000)
overdisp.est3 <- overdispersed(ard,
                            known_sizes = subpop_sizes[known_ind],
                            known_ind = known_ind,
                            G1_ind = z,
                            N = N,
                            warmup = 500,
                            iter = 5000)


# only subpopulations < 50k
z <- which(subpop_sizes < 50000)
ard <- ard[, z]
subpop_sizes <- pops[z]
known_ind <- c(1:length(subpop_sizes))
overdisp.est4 <- overdispersed(ard,
                            known_sizes = subpop_sizes[known_ind],
                            known_ind = known_ind,
                            N = N,
                            warmup = 500,
                            iter = 5000)


#without corona, the very top level one
ard <- as.matrix(dfans[,-c(6)]) 
subpop_sizes <- pops[-6]
known_ind <- c(1:42)
overdisp.est5 <- overdispersed(ard,
                            known_sizes = subpop_sizes[known_ind],
                            known_ind = known_ind,
                            G1_ind = known_ind,
                            N = N,
                            warmup = 500,
                            iter = 5000)

# same iterations in one mean, correlate
x1 <- data.frame(colMeans(overdisp.est1$degrees))
x2 <- data.frame(colMeans(overdisp.est2$degrees))
x3 <- data.frame(colMeans(overdisp.est3$degrees))
x4 <- data.frame(colMeans(overdisp.est4$degrees))
x5 <- data.frame(colMeans(overdisp.est5$degrees))
summary(x1[,1])
summary(x2[,1])
summary(x3[,1])
summary(x4[,1])
summary(x5[,1])
x <- cbind(x1[,1],x2[,1], x3[,1], x4[,1], x5[,1])
cor(x)



#safe to data
df$netsover1 <- colMeans(overdisp.est1$degrees)
df$netsover2 <- colMeans(overdisp.est2$degrees)
df$netsover3 <- colMeans(overdisp.est3$degrees)
df$netsover4 <- colMeans(overdisp.est4$degrees)
df$netsover5 <- colMeans(overdisp.est5$degrees)
cor(df[,c("netsover1", "netsover2", "netsover3", "netsover4", "netsover5")])

Homogeneity metrics for gender.

# GENDER
table(df$woman)
df$numwoman <- rowSums(df[, c("Sophie", "Julia", "Sanne", "Lisa", "Laura", "Maria", "Linda", "Johanna",
    "Monique", "Ester", "Anna", "Elisabeth", "Cornelia", "Wilhelmina", "Amira", "Samira", "Sara", "Noor")])
df$numman <- rowSums(df[, c("Daan", "Sem", "Thomas", "Max", "Kevin", "Johannes", "Dennis", "Jeroen",
    "Jan", "Marcel", "Cornelis", "Hendrik", "Petrus", "Willem", "Ali", "Mohammed")])
df$numnames <- df$numwoman + df$numman
df$samegender <- ifelse(df$woman == 1, (df$numwoman/df$numnames) * 100, (df$numman/df$numnames) * 100)
summary(df$samegender)

summary(df[df$woman == 1, c("samegender")])
summary(df[df$woman == 0, c("samegender")])

Homogeneity metrics for education.

# EDUC independent
df$opl2[df$opleiding == 1] <- 0  # prim/sec
df$opl2[df$opleiding == 2] <- 0  # prim/sec
df$opl2[df$opleiding == 3] <- 0  # prim/sec
df$opl2[df$opleiding == 4] <- 0  # prim/sec

df$opl2[df$opleiding == 5] <- 0  # prim/sec

df$opl2[df$opleiding == 6] <- 1  # prim/sec
df$opl2[df$opleiding == 7] <- 1  # prim/sec
table(df$opleiding)
table(df$opl2)

df$numeduchigh <- df$hbo + df$uni
df$numeduclow <- df$mbo

df$numeduc <- rowSums(df[, c("uni", "hbo", "mbo")])

df$sameeduc <- ifelse(df$opl2 == 0, (df$numeduclow/df$numeduc) * 100, NA)
df$sameeduc <- ifelse(df$opl2 == 1, (df$numeduchigh/df$numeduc) * 100, df$sameeduc)

summary(df$sameeduc)
summary(df[df$opl2 == 1, c("sameeduc")])
summary(df[df$opl2 == 0, c("sameeduc")])

# df$sameeduc[df$oplei == 5] <- (df$hbo / df$numeduc) * 100 df$sameeduc[df$oplei == 6] <- (df$hbo /
# df$numeduc) * 100 df$sameeduc[df$oplei == 7] <- (df$uni / df$numeduc) * 100

# df[138, c('opleiding', 'oplei', 'numeduc', 'mbo', 'hbo', 'uni', 'sameeduc')]

# geen onderwijs / basisonderwijs / cursus inburgering / cursus Nederlandse taal 1 LBO / VBO / VMBO
# (kader- of beroepsgerichte leerweg) / MBO 1 (assistentenopleidi 2 MAVO / HAVO of VWO (eerste drie
# jaar) / ULO / MULO / VMBO (TL of GL) / VSO 3 MBO 2, 3, 4 (basisberoeps-, vak-, middenkader- of
# specialistenopleiding) of MBO 4 HAVO of VWO (overgegaan naar de 4e klas) / HBS / MMS / HBO
# propedeuse of WO Prop 5 HBO (behalve HBO-master) / WO-kandidaats -of WO-bachelor 6 WO-doctoraal
# of WO-master of HBO-master / postdoctoraal onderwijs 7

Homogeneity metrics for ethnic background

# # ETHNIC df$numdutch <- rowSums(df[, c('Sophie',
# 'Julia','Sanne','Lisa','Laura','Maria','Linda','Johanna','Monique','Ester',
# 'Anna','Elisabeth','Cornelia','Wilhelmina', 'Sara', 'Noor',
# 'Daan','Sem','Thomas','Max','Kevin','Johannes','Dennis',
# 'Jeroen','Jan','Marcel','Cornelis','Hendrik','Petrus', 'Willem')]) df$numnodutch <- rowSums(df[,
# c('Ali','Mohammed','Amira','Samira')]) table(df$migr3) df$sameethnic1 <- ifelse(df$migr3 == 1,
# (df$numdutch/df$numnames)*100, (df$numnodutch/df$numnames)*100) df$sameethnic2 <- ifelse(df$migr3
# == 1, (df$numdutch/df$numnames)*100, (df$numnodutch/df$numnames)*100) df$sameethnic2 <-
# ifelse(df$migr3 == 2, NA, df$sameethnic1) df[37, c('migr3', 'sameethnic1', 'numdutch',
# 'numnodutch', 'numnames')]

2 Data storage

Finally, we save the data that we use for the analyses.

# save data
save(df, file = "data/dutch_netsize_analyses_revision.rda")