To copy the code, click the button in the upper right corner of the code-chunks.
rm(list = ls())
gc()
fpackage.check
: Check if packages are installed (and
install if not) in Rfsave
: Function to save data with time stamp in correct
directoryfload
: Function to load R-objects under new namesftheme
: pretty ggplot2 themefshowdf
: Print objects (tibble
/
data.frame
) nicely on screen in .Rmd
.ffit
: fit a series of (here, generalized linear
mixed-effects) modelsfpackage.check <- function(packages) {
lapply(packages, FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)
}
})
}
fsave <- function(x, file, location = "./data/processed/", ...) {
if (!dir.exists(location))
dir.create(location)
datename <- substr(gsub("[:-]", "", Sys.time()), 1, 8)
totalname <- paste(location, datename, file, sep = "")
print(paste("SAVED: ", totalname, sep = ""))
save(x, file = totalname)
}
fload <- function(fileName) {
load(fileName)
get(ls()[ls() != "fileName"])
}
# extrafont::font_import(paths = c('C:/Users/u244147/Downloads/Jost/', prompt = FALSE))
ftheme <- function() {
# download font at https://fonts.google.com/specimen/Jost/
theme_minimal(base_family = "Jost") + theme(panel.grid.minor = element_blank(), plot.title = element_text(family = "Jost",
face = "bold"), axis.title = element_text(family = "Jost Medium"), axis.title.x = element_text(hjust = 0),
axis.title.y = element_text(hjust = 1), strip.text = element_text(family = "Jost", face = "bold",
size = rel(0.75), hjust = 0), strip.background = element_rect(fill = "grey90", color = NA),
legend.position = "bottom")
}
fshowdf <- function(x, digits = 2, ...) {
knitr::kable(x, digits = digits, "html", ...) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover")) %>%
kableExtra::scroll_box(width = "100%", height = "300px")
}
ffit <- function(formula, data) {
tryCatch({
model <- lme4::glmer(formula, data = data, family = binomial(link = "logit"), control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 1e+05)))
cat("Fitting model:", as.character(formula), "\n")
summary(model)
cat("\n")
return(model)
}, error = function(e) {
cat("Error fitting model:", as.character(formula), "\n")
cat("Error message:", conditionMessage(e), "\n")
return(NULL)
})
}
tidyverse
lme4
: fitting random effects modelslmtest
: diagnostics test (likelihood ratio test)car
: companion applied regression (calculate VIF)texreg
: output to HTML tableggplot2
ggpubr
: format ggplot2 plotsggh4x
: hacks for ggplot2ggtext
: text renderingparallel
: parallel computingpackages = c("tidyverse", "lme4", "lmtest", "car", "texreg", "ggplot2", "ggpubr", "ggh4x", "ggtext",
"parallel")
fpackage.check(packages)
rm(packages)
Load the replicated data-set. To load these file, adjust the filename
in the following code so that it matches the most recent version of the
.RDa
file you have in your ./data/processed/
folder.
You may also obtain them by downloading: Download networkdata.Rda
# list files in processed data folder
list.files("./data/processed/")
# get todays date:
today <- gsub("-", "", Sys.Date())
# use fload
df <- fload(paste0("./data/processed/", today, "networkdata.Rda"))
#subset sports partners at t
df <- df[df$csn == 1,]
#outcome
df$y <- ifelse(df$Ycsn == 1, 1, 0)
#calculate rpelacement candidates: no. of other sports partners next to alter j, at time t, with whom ego does the same activity type
df$nreplace <- NA
for (i in unique(df$ego)) {
for (t in unique(df$period[df$ego == i])) {
for (j in unique(df$alterid[df$ego == i & df$period == t])) {
#get activity type alter j does with ego at t
activity <- df$sporttogether[df$ego == i & df$alterid == j & df$period == t]
#get number of *other* alters with whom ego does the same activity at t (ie, replacement candidates)
df$nreplace[df$ego == i & df$alterid == j & df$period == t] <- length(which(df$sporttogether[df$ego == i & !df$alterid == j & df$period == t] == activity))
}
}
}
df$proximity <- factor(df$proximity, levels = c("far","close","roommate"))
df$close <- ifelse(df$proximity == "close",1,0)
df$roommate <- ifelse(df$proximity == "roommate",1,0)
df$far <- ifelse(df$proximity == "far",1,0)
df$ego_context[is.na(df$ego_context)] <- "missing"
df$ego_context <- factor(df$ego_context, levels = c("club", "informal", "gym", "alone", "missing"))
df$club <- ifelse(df$ego_context == "club",1,0)
df$informal <- ifelse(df$ego_context == "informal",1,0)
df$gym <- ifelse(df$ego_context == "gym",1,0)
df$alone <- ifelse(df$ego_context == "alone",1,0)
df$missing <- ifelse(df$ego_context == "missing",1,0)
df$HH <- ifelse(df$ego_grade > 7 & df$alter_grade > 7, 1, 0)
df$HM <- ifelse( ((df$ego_grade > 7 & df$alter_grade > 5 & df$alter_grade < 8) | (df$alter_grade > 7 & df$ego_grade > 5 & df$ego_grade < 8)), 1, 0)
df$HL <- ifelse( ((df$ego_grade > 7 & df$alter_grade < 6) | (df$alter_grade > 7 & df$ego_grade < 6)), 1, 0)
df$MM <- ifelse( ((df$ego_grade > 5 & df$ego_grade < 8 & df$alter_grade > 5 & df$alter_grade < 8)), 1, 0)
df$ML <- ifelse( ((df$ego_grade > 5 & df$ego_grade < 8 & df$alter_grade < 6) | (df$ego_grade < 6 & df$alter_grade < 8 & df$alter_grade > 5)), 1, 0)
df$LL <- ifelse( ((df$ego_grade < 6 & df$alter_grade < 6)), 1, 0)
df$skills <- factor(
1 * (((df$ego_grade > 5 & df$ego_grade < 8 & df$alter_grade > 5 & df$alter_grade < 8))) + # MM
2 * (df$ego_grade > 7 & df$alter_grade > 7) + # HH
3 * (((df$ego_grade > 7 & df$alter_grade > 5 & df$alter_grade < 8) | (df$alter_grade > 7 & df$ego_grade > 5 & df$ego_grade < 8))) + # HM
4 * (((df$ego_grade > 7 & df$alter_grade < 6) | (df$alter_grade > 7 & df$ego_grade < 6))) + # HL
5 * (((df$ego_grade > 5 & df$ego_grade < 8 & df$alter_grade < 6) | (df$ego_grade < 6 & df$alter_grade < 8 & df$alter_grade > 5))) + # ML
6 * (((df$ego_grade < 6 & df$alter_grade < 6))), # LL
levels = c(1, 2, 3, 4, 5, 6),
labels = c("MM", "HH", "HM", "HL", "ML", "LL")
)
df$dif_skill <- abs(df$ego_grade - df$alter_grade)
df$alter_freq2 <- ifelse(df$alter_freq == "1 keer per week", 1,
ifelse(df$alter_freq == "2 of 3 keer per week", 2.5,
ifelse(df$alter_freq == "4 of 5 keer per week", 4.5,
ifelse(df$alter_freq == "6 keer per week of vaker", 6.5,
ifelse(df$alter_freq == "Minder dan 1 keer per maand", 0.125,
ifelse(df$alter_freq == "1 of 2 keer per maand", 0.375,
ifelse(df$alter_freq == "1 of 2 keer per maand", 0.375, NA)))))))
df$ego_freq <- as.numeric(df$ego_freq)
df$ego_activew1 <- ifelse(df$ego_meanfreq >0, 1, 0)
df <- df %>%
mutate(mean_freq = rowMeans(select(., alter_freq2, ego_freq), na.rm = TRUE))
#df$ego_meanskill[is.na(df$ego_meanskill)] <- mean(df$ego_meanskill[which(!duplicated(df$ego))], na.rm=TRUE )
df$alter_age <- as.numeric(df$alter_age)
df$ff <- ifelse(df$ego_female == 1 & df$alter_female == 1, 1, 0)
df$fm <- ifelse( ((df$ego_female == 1 & df$alter_female == 0) | (df$ego_female == 0 & df$alter_female == 1)), 1, 0)
df$mm <- ifelse(df$ego_female == 0 & df$alter_female == 0, 1, 0)
df$gender <- factor(
1 * (df$ego_female == 0 & df$alter_female == 0) + # MM
2 * ((df$ego_female == 1 & df$alter_female == 0) | (df$ego_female == 0 & df$alter_female == 1)) + # FM
3 * (df$ego_female == 1 & df$alter_female == 1), # FF
levels = c(1, 2, 3),
labels = c("MM", "FM", "FF")
)
df$ego_quit[is.na(df$ego_quit)] <- 0
#exclude kin
df[df$kin == "0",] -> df
#types of sport
df$fitness <- ifelse(df$sporttogether == 1,1,0)
df$endurance <- ifelse(df$sporttogether %in% c(2,5,7), 1, 0)
df$team <- ifelse(df$sporttogether %in% c(3,10,12), 1, 0)
df$misc <- ifelse(!df$sporttogether %in% c(1,2,5,7,3,10,12),1,0)
df$sporttype <- factor(
1 * (df$fitness == 1) +
2 * (df$endurance == 1) +
3 * (df$team == 1) +
4 * (df$misc == 1),
levels = c(1:4),
labels = c("fitness", "endurance", "team", "miscellaneous")
)
# listwise deletion
df %>%
select(proximity, roommate, close, far, frequency.t, closeness.t, duration, csn, bff, study, cdn,
gender, mm, fm, ff, period, y, ego, alterid, ego_activew1, mean_freq, skills, HH, HM, HL, MM,
ML, LL, ego_grade, dif_skill, ego_context, club, informal, gym, alone, missing, sporttogether,
nreplace, sporttype, fitness, endurance, team, misc) %>%
filter(complete.cases(.)) -> df
# describe
df %>%
select(-c(proximity, gender, skills, ego_context, ego, alterid, ego_activew1, sporttogether, sporttype)) %>%
psych::describe() %>%
fshowdf(caption = "descriptive statistics of ego's (non-kin) social relations WITHIN SPORTS at time t")
vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
roommate | 1 | 1426 | 0.12 | 0.33 | 0 | 0.03 | 0.00 | 0.00 | 1.00 | 1.00 | 2.29 | 3.23 | 0.01 |
close | 2 | 1426 | 0.70 | 0.46 | 1 | 0.75 | 0.00 | 0.00 | 1.00 | 1.00 | -0.87 | -1.25 | 0.01 |
far | 3 | 1426 | 0.18 | 0.38 | 0 | 0.10 | 0.00 | 0.00 | 1.00 | 1.00 | 1.69 | 0.85 | 0.01 |
frequency.t | 4 | 1426 | 6.04 | 1.04 | 6 | 6.20 | 1.48 | 1.00 | 7.00 | 6.00 | -1.78 | 4.87 | 0.03 |
closeness.t | 5 | 1426 | 3.00 | 0.93 | 3 | 3.09 | 1.48 | 1.00 | 4.00 | 3.00 | -0.52 | -0.74 | 0.02 |
duration | 6 | 1426 | 3.91 | 3.95 | 2 | 3.25 | 2.97 | 0.00 | 15.00 | 15.00 | 1.31 | 0.95 | 0.10 |
csn | 7 | 1426 | 1.00 | 0.00 | 1 | 1.00 | 0.00 | 1.00 | 1.00 | 0.00 | NaN | NaN | 0.00 |
bff | 8 | 1426 | 0.41 | 0.49 | 0 | 0.39 | 0.00 | 0.00 | 1.00 | 1.00 | 0.35 | -1.88 | 0.01 |
study | 9 | 1426 | 0.16 | 0.37 | 0 | 0.08 | 0.00 | 0.00 | 1.00 | 1.00 | 1.81 | 1.28 | 0.01 |
cdn | 10 | 1426 | 0.33 | 0.47 | 0 | 0.29 | 0.00 | 0.00 | 1.00 | 1.00 | 0.70 | -1.51 | 0.01 |
mm | 11 | 1426 | 0.15 | 0.36 | 0 | 0.07 | 0.00 | 0.00 | 1.00 | 1.00 | 1.93 | 1.74 | 0.01 |
fm | 12 | 1426 | 0.27 | 0.44 | 0 | 0.21 | 0.00 | 0.00 | 1.00 | 1.00 | 1.06 | -0.88 | 0.01 |
ff | 13 | 1426 | 0.58 | 0.49 | 1 | 0.60 | 0.00 | 0.00 | 1.00 | 1.00 | -0.33 | -1.89 | 0.01 |
period | 14 | 1426 | 1.33 | 0.47 | 1 | 1.28 | 0.00 | 1.00 | 2.00 | 1.00 | 0.73 | -1.46 | 0.01 |
y | 15 | 1426 | 0.43 | 0.50 | 0 | 0.41 | 0.00 | 0.00 | 1.00 | 1.00 | 0.28 | -1.92 | 0.01 |
mean_freq | 16 | 1426 | 2.04 | 1.34 | 2 | 1.91 | 1.48 | 0.12 | 6.75 | 6.62 | 0.84 | 0.62 | 0.04 |
HH | 17 | 1426 | 0.21 | 0.41 | 0 | 0.14 | 0.00 | 0.00 | 1.00 | 1.00 | 1.42 | 0.03 | 0.01 |
HM | 18 | 1426 | 0.35 | 0.48 | 0 | 0.31 | 0.00 | 0.00 | 1.00 | 1.00 | 0.62 | -1.61 | 0.01 |
HL | 19 | 1426 | 0.05 | 0.23 | 0 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | 3.91 | 13.32 | 0.01 |
MM | 20 | 1426 | 0.24 | 0.43 | 0 | 0.18 | 0.00 | 0.00 | 1.00 | 1.00 | 1.20 | -0.55 | 0.01 |
ML | 21 | 1426 | 0.11 | 0.31 | 0 | 0.01 | 0.00 | 0.00 | 1.00 | 1.00 | 2.50 | 4.25 | 0.01 |
LL | 22 | 1426 | 0.03 | 0.18 | 0 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | 5.23 | 25.33 | 0.00 |
ego_grade | 23 | 1426 | 6.80 | 1.45 | 7 | 6.95 | 1.48 | 1.00 | 10.00 | 9.00 | -0.97 | 1.21 | 0.04 |
dif_skill | 24 | 1426 | 1.25 | 1.18 | 1 | 1.09 | 1.48 | 0.00 | 9.00 | 9.00 | 1.55 | 4.74 | 0.03 |
club | 25 | 1426 | 0.35 | 0.48 | 0 | 0.31 | 0.00 | 0.00 | 1.00 | 1.00 | 0.64 | -1.59 | 0.01 |
informal | 26 | 1426 | 0.12 | 0.32 | 0 | 0.02 | 0.00 | 0.00 | 1.00 | 1.00 | 2.36 | 3.56 | 0.01 |
gym | 27 | 1426 | 0.25 | 0.43 | 0 | 0.18 | 0.00 | 0.00 | 1.00 | 1.00 | 1.18 | -0.61 | 0.01 |
alone | 28 | 1426 | 0.10 | 0.31 | 0 | 0.01 | 0.00 | 0.00 | 1.00 | 1.00 | 2.58 | 4.68 | 0.01 |
missing | 29 | 1426 | 0.18 | 0.39 | 0 | 0.11 | 0.00 | 0.00 | 1.00 | 1.00 | 1.63 | 0.64 | 0.01 |
nreplace | 30 | 1426 | 1.37 | 1.37 | 1 | 1.21 | 1.48 | 0.00 | 6.00 | 6.00 | 0.70 | -0.54 | 0.04 |
fitness | 31 | 1426 | 0.27 | 0.44 | 0 | 0.21 | 0.00 | 0.00 | 1.00 | 1.00 | 1.03 | -0.95 | 0.01 |
endurance | 32 | 1426 | 0.11 | 0.32 | 0 | 0.02 | 0.00 | 0.00 | 1.00 | 1.00 | 2.43 | 3.92 | 0.01 |
team | 33 | 1426 | 0.16 | 0.36 | 0 | 0.07 | 0.00 | 0.00 | 1.00 | 1.00 | 1.88 | 1.55 | 0.01 |
misc | 34 | 1426 | 0.46 | 0.50 | 0 | 0.45 | 0.00 | 0.00 | 1.00 | 1.00 | 0.17 | -1.97 | 0.01 |
length(unique(paste0(df$ego, "X", df$alterid))) #N_alter =1222
#> [1] 1222
nrow(df) #N_observation = 1426
#> [1] 1426
#tie continuation:
formula <- list(
#model 1: main predictors
y ~ 1 + (1 | ego) + scale(ego_grade) + scale(dif_skill) + scale(closeness.t) + ego_context + as.factor(period),
#model 2: control for sports behavior dyad
y ~ 1 + (1 | ego) + scale(ego_grade) + scale(dif_skill) + scale(closeness.t) + ego_context + as.factor(period) + scale(mean_freq),
#model 3: include "traditional" dyadic controls
y ~ 1 + (1 | ego) + scale(ego_grade) + scale(dif_skill) + scale(closeness.t) + ego_context + as.factor(period) + scale(mean_freq) + proximity + scale(frequency.t) + scale(duration) + gender,
#model 4: add other relational dimensions (multiplexity)
y ~ 1 + (1 | ego) + scale(ego_grade) + scale(dif_skill) + scale(closeness.t) + ego_context + as.factor(period) + scale(mean_freq) + proximity + scale(frequency.t) + scale(duration) + gender + bff + cdn + study,
#model 5: add replacement candidates
y ~ 1 + (1 | ego) + scale(ego_grade) + scale(dif_skill) + scale(closeness.t) + ego_context + as.factor(period) + scale(mean_freq) + proximity + scale(frequency.t) + scale(duration) + gender + bff + cdn + study + scale(nreplace)
)
ans <- lapply(formula, ffit, data = df)
lapply(ans, summary)
do.call(lmtest::lrtest,ans)
M1: main predictors | M2: sports beh. dyad | M3: dyadic covars | M4: multiplexity | M5: replacement candidates | |
---|---|---|---|---|---|
(Intercept) | -0.65 (0.12)*** | -0.68 (0.12)*** | -0.80 (0.24)** | -0.93 (0.26)*** | -0.86 (0.26)*** |
Ego skill | 0.09 (0.07) | 0.01 (0.08) | -0.00 (0.08) | -0.00 (0.08) | 0.00 (0.08) |
Ego-alter skill diference | -0.08 (0.07) | -0.09 (0.07) | -0.09 (0.07) | -0.10 (0.07) | -0.09 (0.07) |
Emotional closeness | 0.55 (0.07)*** | 0.56 (0.07)*** | 0.33 (0.09)*** | 0.23 (0.10)* | 0.24 (0.10)* |
Informal group | -0.22 (0.22) | -0.10 (0.22) | -0.19 (0.23) | -0.25 (0.24) | -0.33 (0.24) |
Commercial gym | 0.47 (0.17)** | 0.47 (0.17)** | 0.30 (0.18) | 0.27 (0.19) | 0.20 (0.19) |
Unorganized | 0.13 (0.22) | 0.23 (0.23) | 0.07 (0.24) | 0.05 (0.24) | -0.05 (0.25) |
Missing | 0.03 (0.24) | 0.06 (0.24) | -0.10 (0.25) | 0.00 (0.25) | 0.01 (0.25) |
Period: waves 2-3 | 0.68 (0.18)*** | 0.68 (0.18)*** | 0.66 (0.18)*** | 0.61 (0.19)*** | 0.59 (0.19)** |
Mean sports frequency dyad | 0.21 (0.07)** | 0.16 (0.07)* | 0.16 (0.08)* | 0.18 (0.08)* | |
Same municipality | 0.36 (0.18)* | 0.39 (0.18)* | 0.39 (0.18)* | ||
Roommate | 0.19 (0.25) | 0.20 (0.25) | 0.17 (0.25) | ||
Communication frequency | 0.54 (0.10)*** | 0.50 (0.10)*** | 0.50 (0.10)*** | ||
Years known | -0.04 (0.07) | -0.07 (0.07) | -0.07 (0.07) | ||
Woman-man | -0.07 (0.22) | -0.12 (0.22) | -0.13 (0.22) | ||
Woman-woman | -0.13 (0.20) | -0.15 (0.20) | -0.15 (0.20) | ||
Friendship | 0.03 (0.18) | 0.01 (0.18) | |||
Confidant | 0.49 (0.18)** | 0.48 (0.18)** | |||
Study partner | -0.21 (0.18) | -0.23 (0.18) | |||
No. of replacement candidates | -0.14 (0.08) | ||||
AIC | 1830.18 | 1823.41 | 1790.83 | 1786.93 | 1785.89 |
BIC | 1882.80 | 1881.30 | 1880.30 | 1892.18 | 1896.40 |
Log Likelihood | -905.09 | -900.71 | -878.42 | -873.47 | -871.94 |
Num. obs. | 1426 | 1426 | 1426 | 1426 | 1426 |
Num. groups: ego | 409 | 409 | 409 | 409 | 409 |
Var: ego (Intercept) | 0.29 | 0.29 | 0.35 | 0.37 | 0.36 |
***p < 0.001; **p < 0.01; *p < 0.05 |
We compare results in sports partnerships in vs. outside fitness.
df %>%
select(sporttogether) %>%
table(.) %>%
prop.table(.)
# of all sports partner observations: 27% are in fitness; this is the most popular sports in dyads
# 11% are in endurance sports (i.e., running, cycling, swimming) 15% are in team (ball) sports
# (i.e., football, hockey, volleybal)
# 1. solution for fitness partners:
ans_fitness <- lapply(formula, ffit, data = df[df$sporttogether == 1, ])
lapply(ans_fitness, summary)
# 2. solution excluding fitness partners:
ans_nfitness <- lapply(formula, ffit, data = df[!df$sporttogether == 1, ])
lapply(ans_nfitness, summary)
# 3. solution for endurance sports partners:
ans_endurance <- lapply(formula, ffit, data = df[df$sporttogether %in% c(2, 5, 7), ])
lapply(ans_endurance, summary)
# 4. solution for team (ball) sports partners:
ans_team <- lapply(formula, ffit, data = df[df$sporttogether %in% c(3, 10, 12), ])
lapply(ans_team, summary)
M1: main predictors | M2: sports beh. dyad | M3: dyadic covars | M4: multiplexity | M5: replacement candidates | |
---|---|---|---|---|---|
(Intercept) | -0.14 (0.56) | -0.04 (0.56) | -0.75 (0.68) | -0.92 (0.70) | -0.96 (0.73) |
Ego skill | 0.19 (0.13) | -0.03 (0.15) | 0.02 (0.16) | 0.03 (0.16) | 0.06 (0.17) |
Ego-alter skill diference | 0.06 (0.13) | -0.01 (0.13) | 0.08 (0.14) | 0.09 (0.14) | 0.10 (0.14) |
Emotional closeness | 0.38 (0.12)** | 0.40 (0.12)** | 0.05 (0.16) | -0.05 (0.18) | -0.04 (0.18) |
Informal group | -0.77 (0.70) | -0.80 (0.69) | -0.81 (0.74) | -1.03 (0.74) | -0.88 (0.77) |
Commercial gym | 0.07 (0.58) | -0.08 (0.58) | -0.31 (0.62) | -0.39 (0.61) | -0.25 (0.63) |
Unorganized | -0.07 (0.65) | 0.06 (0.65) | -0.15 (0.70) | -0.26 (0.69) | -0.21 (0.71) |
Missing | -1.03 (0.71) | -1.10 (0.71) | -1.15 (0.76) | -1.21 (0.75) | -0.98 (0.78) |
Period: waves 2-3 | 0.91 (0.32)** | 0.92 (0.33)** | 0.98 (0.35)** | 0.99 (0.35)** | 0.91 (0.36)* |
Mean sports frequency dyad | 0.40 (0.14)** | 0.37 (0.15)* | 0.34 (0.15)* | 0.37 (0.16)* | |
Same municipality | 0.77 (0.38)* | 0.85 (0.38)* | 0.83 (0.40)* | ||
Roommate | 0.36 (0.46) | 0.41 (0.46) | 0.41 (0.47) | ||
Communication frequency | 0.72 (0.19)*** | 0.68 (0.19)*** | 0.72 (0.19)*** | ||
Years known | 0.12 (0.13) | 0.07 (0.13) | 0.07 (0.13) | ||
Woman-man | 0.18 (0.37) | 0.20 (0.37) | 0.16 (0.38) | ||
Woman-woman | 0.41 (0.34) | 0.45 (0.33) | 0.38 (0.34) | ||
Friendship | 0.22 (0.30) | 0.21 (0.30) | |||
Confidant | 0.35 (0.32) | 0.29 (0.33) | |||
Study partner | -0.52 (0.31) | -0.56 (0.31) | |||
No. of replacement candidates | -0.31 (0.14)* | ||||
AIC | 521.56 | 514.48 | 502.91 | 503.96 | 500.53 |
BIC | 561.14 | 558.02 | 570.20 | 583.13 | 583.66 |
Log Likelihood | -250.78 | -246.24 | -234.45 | -231.98 | -229.26 |
Num. obs. | 387 | 387 | 387 | 387 | 387 |
Num. groups: ego | 196 | 196 | 196 | 196 | 196 |
Var: ego (Intercept) | 0.15 | 0.10 | 0.21 | 0.15 | 0.23 |
***p < 0.001; **p < 0.01; *p < 0.05 |
M1: main predictors | M2: sports beh. dyad | M3: dyadic covars | M4: multiplexity | M5: replacement candidates | |
---|---|---|---|---|---|
(Intercept) | -0.70 (0.13)*** | -0.73 (0.13)*** | -0.73 (0.29)* | -0.82 (0.30)** | -0.81 (0.30)** |
Ego skill | 0.04 (0.09) | -0.00 (0.09) | -0.01 (0.10) | -0.02 (0.10) | -0.02 (0.10) |
Ego-alter skill diference | -0.14 (0.08) | -0.14 (0.08) | -0.16 (0.08) | -0.17 (0.09) | -0.16 (0.09) |
Emotional closeness | 0.61 (0.09)*** | 0.61 (0.09)*** | 0.42 (0.11)*** | 0.31 (0.12)** | 0.32 (0.12)** |
Informal group | -0.22 (0.25) | -0.13 (0.25) | -0.24 (0.27) | -0.30 (0.27) | -0.33 (0.28) |
Commercial gym | 0.32 (0.28) | 0.41 (0.28) | 0.29 (0.30) | 0.22 (0.31) | 0.19 (0.32) |
Unorganized | 0.06 (0.26) | 0.12 (0.27) | 0.00 (0.28) | -0.05 (0.29) | -0.09 (0.30) |
Missing | 0.21 (0.27) | 0.23 (0.27) | 0.10 (0.28) | 0.22 (0.29) | 0.23 (0.29) |
Period: waves 2-3 | 0.57 (0.22)** | 0.58 (0.22)** | 0.50 (0.22)* | 0.44 (0.23) | 0.43 (0.23) |
Mean sports frequency dyad | 0.12 (0.09) | 0.08 (0.09) | 0.08 (0.09) | 0.09 (0.09) | |
Same municipality | 0.30 (0.22) | 0.31 (0.22) | 0.31 (0.22) | ||
Roommate | 0.14 (0.31) | 0.12 (0.32) | 0.10 (0.32) | ||
Communication frequency | 0.50 (0.12)*** | 0.45 (0.12)*** | 0.45 (0.12)*** | ||
Years known | -0.11 (0.09) | -0.11 (0.09) | -0.12 (0.09) | ||
Woman-man | -0.14 (0.28) | -0.21 (0.29) | -0.22 (0.29) | ||
Woman-woman | -0.24 (0.25) | -0.28 (0.26) | -0.28 (0.26) | ||
Friendship | -0.07 (0.23) | -0.08 (0.23) | |||
Confidant | 0.63 (0.23)** | 0.62 (0.24)** | |||
Study partner | -0.08 (0.23) | -0.08 (0.23) | |||
No. of replacement candidates | -0.05 (0.10) | ||||
AIC | 1322.63 | 1322.52 | 1304.67 | 1302.64 | 1304.39 |
BIC | 1372.09 | 1376.92 | 1388.75 | 1401.56 | 1408.25 |
Log Likelihood | -651.32 | -650.26 | -635.33 | -631.32 | -631.19 |
Num. obs. | 1039 | 1039 | 1039 | 1039 | 1039 |
Num. groups: ego | 328 | 328 | 328 | 328 | 328 |
Var: ego (Intercept) | 0.36 | 0.35 | 0.46 | 0.49 | 0.49 |
***p < 0.001; **p < 0.01; *p < 0.05 |
For more information on the (numerical) approach to computing AMEs, see https://www.jochemtolsma.nl/tutorials/me/.
# 2. tie maintenance
dfclose1 <- dfclose0 <- df
dfroommate1 <- dfroommate0 <- df
dfclose1$proximity <- "close"
dfclose0$proximity <- "far"
dfroommate1$proximity <- "roommate"
dfroommate0$proximity <- "far"
s <- 0.001
dfclosenessplus <- dfclosenessmin <- df
dfclosenessplus$closeness.t <- df$closeness.t + s
dfclosenessmin$closeness.t <- df$closeness.t - s
dffriend1 <- dffriend0 <- df
dfstudy1 <- dfstudy0 <- df
dfcdn1 <- dfcdn0 <- df
dffriend1$bff <- 1
dffriend0$bff <- 0
dfstudy1$study <- 1
dfstudy0$study <- 0
dfcdn1$cdn <- 1
dfcdn0$cdn <- 0
dfwomen1 <- dfwomen0 <- dfmixed1 <- dfmixed0 <- df
dfwomen1$gender <- "FF"
dfwomen0$gender <- "MM"
dfmixed1$gender <- "FM"
dfmixed0$gender <- "MM"
dfdurationplus <- dfdurationmin <- df
dfdurationplus$duration <- df$duration + s
dfdurationmin$duration <- df$duration - s
dffrequencyplus <- dffrequencymin <- df
dffrequencyplus$frequency.t <- df$frequency.t + s
dffrequencymin$frequency.t <- df$frequency.t - s
dfperiod21 <- dfperiod20 <- df
dfperiod21$period <- 2
dfperiod20$period <- 1
dfmeanfreqplus <- dfmeanfreqmin <- df
dfmeanfreqplus$mean_freq <- df$mean_freq + s
dfmeanfreqmin$mean_freq <- df$mean_freq - s
# dfquit1 <- dfquit0 <- df dfquit1$ego_quit <- 1 dfquit0$ego_quit <- 0
# dfHH1 <- dfHH0 <- dfHM1 <- dfHM0 <- dfHL1 <- dfHL0 <- dfML1 <- dfML0 <- dfLL1 <- dfLL0 <- df
# dfHH1$skills <- 'HH' dfHH0$skills <- 'MM' dfHM1$skills <- 'HM' dfHM0$skills <- 'MM' dfHL1$skills
# <- 'HL' dfHL0$skills <- 'MM' dfML1$skills <- 'ML' dfML0$skills <- 'MM' dfLL1$skills <- 'LL'
# dfLL0$skills <- 'MM'
dfegoskillplus <- dfegoskillmin <- df
dfegoskillplus$ego_grade <- df$ego_grade + s
dfegoskillmin$ego_grade <- df$ego_grade - s
dfskilldifplus <- dfskilldifmin <- df
dfskilldifplus$dif_skill <- df$dif_skill + s
dfskilldifmin$dif_skill <- df$dif_skill - s
dfinformal1 <- dfinformal0 <- dfgym1 <- dfgym0 <- dfalone1 <- dfalone0 <- dfmissing1 <- dfmissing0 <- df
dfinformal1$ego_context <- "informal"
dfinformal0$ego_context <- "club"
dfgym1$ego_context <- "gym"
dfgym0$ego_context <- "club"
dfalone1$ego_context <- "alone"
dfalone0$ego_context <- "club"
dfmissing1$ego_context <- "missing"
dfmissing0$ego_context <- "club"
dfreplaceplus <- dfreplacemin <- df
dfreplaceplus$nreplace <- df$nreplace + s
dfreplacemin$nreplace <- df$nreplace - s
fpred <- function(x) {
me_close <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfclose1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfclose0)
ame_close <- mean(me_close)
me_roommate <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfroommate1) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfroommate0)
ame_roommate <- mean(me_roommate)
me_closeness <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfclosenessplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfclosenessmin))/(2 * s)
ame_closeness <- mean(me_closeness)
me_friend <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dffriend1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dffriend0)
ame_friend <- mean(me_friend)
me_study <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfstudy1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfstudy0)
ame_study <- mean(me_study)
me_cdn <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfcdn1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfcdn0)
ame_cdn <- mean(me_cdn)
me_women <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfwomen1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfwomen0)
ame_women <- mean(me_women)
me_mixed <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfmixed1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfmixed0)
ame_mixed <- mean(me_mixed)
me_duration <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfdurationplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfdurationmin))/(2 * s)
ame_duration <- mean(me_duration)
me_frequency <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dffrequencyplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dffrequencymin))/(2 * s)
ame_frequency <- mean(me_frequency)
me_period2 <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfperiod21) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfperiod20)
ame_period2 <- mean(me_period2)
# me_quit <- lme4:::predict.merMod(x, type = 'response', re.form = NULL, newdata = dfquit1) -
# lme4:::predict.merMod(x, #type = 'response', re.form = NULL, newdata = dfquit0) ame_quit <-
# mean(me_quit)
# me_HH <- lme4:::predict.merMod(x, type = 'response', re.form = NULL, newdata = dfHH1) -
# lme4:::predict.merMod(x, type = #'response', re.form = NULL, newdata = dfHH0) ame_HH <-
# mean(me_HH) me_HM <- lme4:::predict.merMod(x, type = 'response', re.form = NULL, newdata =
# dfHM1) - lme4:::predict.merMod(x, type = #'response', re.form = NULL, newdata = dfHM0) ame_HM
# <- mean(me_HM) me_HL <- lme4:::predict.merMod(x, type = 'response', re.form = NULL, newdata =
# dfHL1) - lme4:::predict.merMod(x, type = #'response', re.form = NULL, newdata = dfHL0) ame_HL
# <- mean(me_HL) me_ML <- lme4:::predict.merMod(x, type = 'response', re.form = NULL, newdata =
# dfML1) - lme4:::predict.merMod(x, type = #'response', re.form = NULL, newdata = dfML0) ame_ML
# <- mean(me_ML) me_LL <- lme4:::predict.merMod(x, type = 'response', re.form = NULL, newdata =
# dfLL1) - lme4:::predict.merMod(x, type = #'response', re.form = NULL, newdata = dfLL0) ame_LL
# <- mean(me_LL)
me_egoskill <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfegoskillplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfegoskillmin))/(2 * s)
ame_egoskill <- mean(me_egoskill)
me_skilldif <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfskilldifplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfskilldifmin))/(2 * s)
ame_skilldif <- mean(me_skilldif)
me_informal <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfinformal1) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfinformal0)
ame_informal <- mean(me_informal)
me_gym <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfgym1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfgym0)
ame_gym <- mean(me_gym)
me_alone <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfalone1) - lme4:::predict.merMod(x,
type = "response", re.form = NULL, newdata = dfalone0)
ame_alone <- mean(me_alone)
me_missing <- lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfmissing1) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfmissing0)
ame_missing <- mean(me_missing)
me_replace <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfreplaceplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfreplacemin))/(2 * s)
ame_replace <- mean(me_replace)
me_meanfreq <- (lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfmeanfreqplus) -
lme4:::predict.merMod(x, type = "response", re.form = NULL, newdata = dfmeanfreqmin))/(2 * s)
ame_meanfreq <- mean(me_meanfreq)
c(ame_close, ame_roommate, ame_closeness, ame_friend, ame_study, ame_cdn, ame_women, ame_mixed, ame_duration,
ame_frequency, ame_period2, ame_egoskill, ame_skilldif, ame_informal, ame_gym, ame_alone, ame_missing,
ame_replace, ame_meanfreq)
}
# fpred(ans[[5]])
seed <- 242523
nIter <- 500
nCore <- detectCores() - 1
mycl <- makeCluster(rep("localhost", nCore))
clusterEvalQ(mycl, library(lme4))
clusterExport(mycl, varlist = c("ans", "s", "df", "dfclose1", "dfclose0", "dfroommate1", "dfroommate0",
"dfclosenessplus", "dfclosenessmin", "dffriend1", "dffriend0", "dfstudy1", "dfstudy0", "dfcdn1",
"dfcdn0", "dfwomen1", "dfwomen0", "dfmixed1", "dfmixed0", "dfdurationplus", "dfdurationmin", "dffrequencyplus",
"dffrequencymin", "dfperiod21", "dfperiod20", "dfmeanfreqplus", "dfmeanfreqmin", "dfegoskillmin",
"dfegoskillplus", "dfskilldifmin", "dfskilldifplus", "dfinformal1", "dfinformal0", "dfgym1", "dfgym0",
"dfalone1", "dfalone0", "dfmissing1", "dfmissing0", "dfreplaceplus", "dfreplacemin"))
system.time(boo_m <- bootMer(ans[[5]], fpred, nsim = nIter, parallel = "snow", ncpus = nCore, cl = mycl,
seed = seed))
stopCluster(mycl)
save(boo_m, file = "./results/bootm.Rda")
load("./results/bootm.Rda")
#tie formation
#plotdata1 <- data.frame(
# pred = c("Outside municipality", "Same municipality", "Same house", "*Emotional closeness*", "Friendship", "Study partnership", "Confidant", "Male dyad", "Female dyad", "Mixed dyad", #"*Years known*", "*Communication frequency*", "Ego residential change", "Ego study transition", "Period 1", "Period 2"),
# Outcome = "Tie formation",
# ame = c(0, booL[[1]]$t0[1:6], 0, booL[[1]]$t0[7:12], 0, booL[[1]]$t0[13]),
# ame_se = c(0, apply(booL[[1]]$t, 2, sd)[1:6], 0, apply(booL[[1]]$t, 2, sd)[7:12], 0, apply(booL[[1]]$t, 2, sd)[13]),
# ref = c(1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0))
#tie maintenance
plotdata <- data.frame(
pred = c("Outside municipality", "Same municipality", "Same house",
"*Emotional closeness*", "Friendship", "Study partnership", "Confidant", "Male dyad", "Female dyad", "Mixed dyad", "*Years known*", "*Communication frequency*", "Period 1", "Period 2", "*Ego skill*", "*Ego-alter skill difference*", "Sports club", "Informal group", "Commercial gym", "Unorganized", "Missing", "*No. of replacement candidates*", "*Sports frequency dyad*" ),
#Outcome = "Tie maintenance",
ame = c(0, boo_m$t0[1:6], 0, boo_m$t0[7:10], 0, boo_m$t0[11:13], 0, boo_m$t0[14:19] ),
ame_se = c(0, apply(boo_m$t, 2, sd)[1:6], 0, apply(boo_m$t, 2, sd)[7:10], 0, apply(boo_m$t, 2, sd)[11:13], 0, apply(boo_m$t, 2, sd)[14:19]),
ref = c(1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0))
plotdata$pred <- factor(plotdata$pred, levels = rev(c(
"*Ego-alter skill difference*", "*Emotional closeness*", "Sports club", "Informal group", "Commercial gym", "Unorganized", "Missing", "Outside municipality", "Same municipality", "Same house", "*Communication frequency*","*Years known*", "Male dyad", "Female dyad", "Mixed dyad", "Friendship", "Confidant", "Study partnership", "*Ego skill*", "*Sports frequency dyad*", "*No. of replacement candidates*", "Period 1", "Period 2")))
plotdata <- plotdata[order(plotdata$pred),]
row.names(plotdata) <- 1:nrow(plotdata)
plotdata$ref <- as.factor(plotdata$ref)
#also include coefficients as labels, but leave out the labels for the reference level
plotdata$label <- ifelse(plotdata$ref == 1, 0, 1)
#in main text, only main predictors..
plotdata2 <- plotdata[plotdata$pred %in% c( "*Ego-alter skill difference*", "*Emotional closeness*","Sports club", "Informal group", "Commercial gym", "Unorganized", "Missing"),]
p <- ggplot(plotdata2, aes(x = ame, y = pred, #color = Outcome,
shape = ref)) +
geom_vline(xintercept = 0, color = "grey") +
geom_point() +
geom_errorbar(data = subset(plotdata2, label == 1), aes(xmin = ame - 1.96*ame_se, xmax = ame + 1.96*ame_se), width = 0, linewidth = 0.3) +
geom_text(data = subset(plotdata2, label == 1), aes(label = sprintf("%0.2f (%0.2f)", ame, ame_se )), size = 3, color = "black", position = position_nudge(y = 0.4)) +
#facet_grid(Outcome ~., switch = "y", scales = "free_y", space = "free_y") +
labs(x = "AME", y = NULL) +
scale_x_continuous(labels = scales::percent) +
scale_shape_manual("", values = c("1" = 1, "0" = 16)) +
theme(axis.line = element_line(),
legend.position = "none",
strip.background = element_blank(),
strip.text.x = element_text(face = "bold"),
strip.text.y = element_blank(),
axis.text.y = element_markdown()) + guides(shape = "none")
ggsave("./figures/ames.png", height = 2.5)
plotdata %>%
arrange(desc(row_number())) %>%
select(-label) %>%
fshowdf(digits = 3, caption = "Average marginal effects on sports partnership maintenance")
pred | ame | ame_se | ref |
---|---|---|---|
Ego-alter skill difference | -0.016 | 0.012 | 0 |
Emotional closeness | 0.053 | 0.021 | 0 |
Sports club | 0.000 | 0.000 | 1 |
Informal group | -0.066 | 0.049 | 0 |
Commercial gym | 0.042 | 0.041 | 0 |
Unorganized | -0.010 | 0.049 | 0 |
Missing | 0.001 | 0.052 | 0 |
Outside municipality | 0.000 | 0.000 | 1 |
Same municipality | 0.080 | 0.040 | 0 |
Same house | 0.035 | 0.053 | 0 |
Communication frequency | 0.100 | 0.018 | 0 |
Years known | -0.004 | 0.004 | 0 |
Male dyad | 0.000 | 0.000 | 1 |
Female dyad | -0.032 | 0.039 | 0 |
Mixed dyad | -0.027 | 0.044 | 0 |
Friendship | 0.002 | 0.035 | 0 |
Confidant | 0.103 | 0.040 | 0 |
Study partnership | -0.047 | 0.034 | 0 |
Ego skill | 0.000 | 0.011 | 0 |
Sports frequency dyad | 0.028 | 0.012 | 0 |
No. of replacement candidates | -0.020 | 0.012 | 0 |
Period 1 | 0.000 | 0.000 | 1 |
Period 2 | 0.124 | 0.036 | 0 |
#now with all controls:
p2 <- ggplot(plotdata, aes(x = ame, y = pred, #color = Outcome,
shape = ref)) +
geom_vline(xintercept = 0, color = "grey") +
geom_point() +
geom_errorbar(data = subset(plotdata, label == 1), aes(xmin = ame - 1.96*ame_se, xmax = ame + 1.96*ame_se), width = 0, linewidth = 0.3) +
geom_text(data = subset(plotdata, label == 1), aes(label = sprintf("%0.2f (%0.2f)", ame, ame_se )), size = 2, color = "black", position = position_nudge(y = 0.4)) +
#facet_grid(Outcome ~., switch = "y", scales = "free_y", space = "free_y") +
labs(x = "AME", y = NULL) +
scale_x_continuous(labels = scales::percent) +
scale_shape_manual("", values = c("1" = 1, "0" = 16)) +
#ftheme() +
theme(axis.line = element_line(),
legend.position = "none",
strip.background = element_blank(),
strip.text.x = element_text(face = "bold"),
strip.text.y = element_blank(),
axis.text.y = element_markdown()) + guides(shape = "none")
ggsave("./figures/ames_incl_control.png", height = 5)
#knitr::include_graphics("./figures/ames_incl_control.png")
print(p2)
Copyright © 2025 Rob Franken