This is the custom code that was used for the main analyses for the manuscript “Lonely individuals process the world in idiosyncratic ways” under revision at Psychological Science. We include a small, simulated dataset is provided to efficiently demo the code.
This code uses simulated data with 10 participants, 5 brain regions, and 5 videos. There are two input files:
1. Loneliness values of 10 participants
2. Dyad-level inter-subject correlation for each region of the brain
Read in all the files (please change the .csv paths to the correct paths on your local machine):
lone <- read.csv("~/Desktop/code/loneliness_simulated.csv")
isc <- read.csv("~/Desktop/code/dyad_isc_simulated.csv")First, we create binary variables for loneliness. If a participant has a loneliness score that is greater-than-median, then they are labeled as having “high” loneliness. Otherwise, they are labeled as having “low” loneliness.
lone$loneliness_binary <- ifelse(lone$loneliness > median(lone$loneliness, na.rm = TRUE), "high", "low")Next, we transform the dataframe with the loneliness values into a dyad-level dataframe.
dyad_list <- unique(isc[,c("sub1", "sub2")])
lone_dyad <- merge(dyad_list, lone, by.x = c("sub1"), by.y = c("pID"))
colnames(lone_dyad)[c(3:4)] <- paste("sub1", colnames(lone_dyad)[c(3:4)], sep = "_")
lone_dyad <- merge(lone_dyad, lone, by.x = c("sub2"), by.y = c("pID"))
colnames(lone_dyad)[5:6] <- paste("sub2", colnames(lone_dyad)[5:6], sep = "_")
lone_dyad$dyad_loneliness_binary <- ifelse(lone_dyad$sub1_loneliness_binary==c("low") &
lone_dyad$sub2_loneliness_binary==c("low"), "low_low", ifelse(
lone_dyad$sub1_loneliness_binary==c("high") &
lone_dyad$sub2_loneliness_binary==c("high"), "high_high", "low_high"))
lone_dyad$dyad_loneliness_max <- ifelse(lone_dyad$sub1_loneliness < lone_dyad$sub2_loneliness, lone_dyad$sub2_loneliness, lone_dyad$sub1_loneliness)Merge dyad-level dataframes together:
lone_isc_dyad <- merge(lone_dyad, isc, by = c("sub1", "sub2"))Create doubled dataframes (adding redundancy) to allow fully-crossed random effects and Fisher’s z transform ISCs:
library(DescTools)
lone_isc_dyad_2 <- lone_isc_dyad
colnames(lone_isc_dyad_2)[1:2] <- c("sub2", "sub1")
lone_isc_dyad_double <- rbind(lone_isc_dyad, lone_isc_dyad_2)
lone_isc_dyad_double$fisher_z <- FisherZ(lone_isc_dyad_double$cor)Relate ISC with binarized loneliness
n_unique_dyad <- 45 #input the number of unique dyads in dataset; here, it is 45
library(emmeans)## Welcome to emmeans.
## NOTE -- Important change from versions <= 1.41:
## Indicator predictors are now treated as 2-level factors by default.
## To revert to old behavior, use emm_options(cov.keep = character(0))library(lmerTest)## Loading required package: lme4## Warning: package 'lme4' was built under R version 3.6.2## Loading required package: Matrix##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(plyr)
# Set contrasts:
Contrasts <- list(HHvsLL = c(1, 0, -1),
HHvsLH = c(1, -1, 0),
LHvsLL = c(0, 1, -1))
dyad_ISC_lone_bin <- function(data)
{
md <- (lmer(scale(fisher_z) ~ dyad_loneliness_binary + (1|sub1) + (1|sub2), data))
emm_md <- emmeans(md, ~ dyad_loneliness_binary, contr = Contrasts, adjust = "none")
HHvsLL_B <- summary(emm_md$contrasts)[1,2]
HHvsLH_B <- summary(emm_md$contrasts)[2,2]
LHvsLL_B <- summary(emm_md$contrasts)[3,2]
HHvsLL_p <- summary(emm_md$contrasts)[1,6]
HHvsLH_p <- summary(emm_md$contrasts)[2,6]
LHvsLL_p <- summary(emm_md$contrasts)[3,6]
HHvsLL_t <- summary(emm_md$contrasts)[1,5]
HHvsLH_t <- summary(emm_md$contrasts)[2,5]
LHvsLL_t <- summary(emm_md$contrasts)[3,5]
HHvsLL_se <- summary(emm_md$contrasts)[1,3]
HHvsLH_se <- summary(emm_md$contrasts)[2,3]
LHvsLL_se <- summary(emm_md$contrasts)[2,3]
HHvsLL_se <- HHvsLL_se * (sqrt(2*(n_unique_dyad-1)-1)/sqrt((n_unique_dyad-1)-1))
HHvsLH_se <- HHvsLH_se * (sqrt(2*(n_unique_dyad-1)-1)/sqrt((n_unique_dyad-1)-1))
LHvsLL_se <- LHvsLL_se * (sqrt(2*(n_unique_dyad-1)-1)/sqrt((n_unique_dyad-1)-1))
HHvsLL_p_penalized <- 2 * pt(-abs(HHvsLL_t), (n_unique_dyad-1))
HHvsLH_p_penalized <- 2 * pt(-abs(HHvsLH_t), (n_unique_dyad-1))
LHvsLL_p_penalized <- 2 * pt(-abs(LHvsLL_t), (n_unique_dyad-1))
return(data.frame(HHvsLL_B, HHvsLL_t, HHvsLL_se, HHvsLL_p, HHvsLL_p_penalized,
HHvsLH_B, HHvsLH_t, HHvsLH_se, HHvsLH_p, HHvsLH_p_penalized,
LHvsLL_B, LHvsLL_t, LHvsLL_se, LHvsLL_p, LHvsLL_p_penalized))
}
dyad_ISC_lone_results_bin <- ddply(lone_isc_dyad_double, .(parcel), dyad_ISC_lone_bin)## boundary (singular) fit: see ?isSingular## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## boundary (singular) fit: see ?isSingular## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## boundary (singular) fit: see ?isSingular## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimatesFunction to rearrange data:
library(reshape2)
func_rearrange <- function(data){
melt_B <- melt(data[,c("parcel","HHvsLL_B", "HHvsLH_B", "LHvsLL_B")], id.vars = c("parcel"), measure.vars = c("HHvsLL_B", "HHvsLH_B", "LHvsLL_B"))
melt_B$variable <- gsub("_B","",melt_B$variable)
colnames(melt_B) <- c("parcel", "contrast", "B")
melt_p <- melt(data[,c("parcel","HHvsLL_p_penalized", "HHvsLH_p_penalized", "LHvsLL_p_penalized")], id.vars = c("parcel"), measure.vars =
c("HHvsLL_p_penalized", "HHvsLH_p_penalized", "LHvsLL_p_penalized"))
melt_p$variable <- gsub("_p_penalized","",melt_p$variable)
colnames(melt_p) <- c("parcel", "contrast", "p_penalized")
melt_se <- melt(data[,c("parcel","HHvsLL_se", "HHvsLH_se", "LHvsLL_se")], id.vars = c("parcel"), measure.vars =
c("HHvsLL_se", "HHvsLH_se", "LHvsLL_se"))
melt_se$variable <- gsub("_se","",melt_se$variable)
colnames(melt_se) <- c("parcel", "contrast", "se")
new_data <- merge(melt_B, melt_p, by = c("parcel", "contrast"))
new_data <- merge(new_data, melt_se, by = c("parcel", "contrast"))
new_data$corrected_p <- p.adjust(new_data$p_penalized, method = "fdr")
return(new_data)
}
dyad_ISC_lone_results_bin <- func_rearrange(dyad_ISC_lone_results_bin)
dyad_ISC_lone_results_bin## parcel contrast B p_penalized se corrected_p
## 1 brain_region1 HHvsLH -0.42145245 0.12391156 0.3821808 0.2828294
## 2 brain_region1 HHvsLL -0.81399503 0.02244255 0.4893324 0.2828294
## 3 brain_region1 LHvsLL -0.39254258 0.15112667 0.3821808 0.2828294
## 4 brain_region2 HHvsLH -0.40131511 0.16969765 0.4088872 0.2828294
## 5 brain_region2 HHvsLL -0.59980039 0.15770443 0.5935843 0.2828294
## 6 brain_region2 LHvsLL -0.19848528 0.49351972 0.4088872 0.7188224
## 7 brain_region3 HHvsLH 0.38663484 0.16510904 0.3896041 0.2828294
## 8 brain_region3 HHvsLL 0.56123428 0.11668089 0.4988370 0.2828294
## 9 brain_region3 LHvsLL 0.17459944 0.52713646 0.3896041 0.7188224
## 10 brain_region4 HHvsLH 0.55526673 0.05355366 0.3981597 0.2828294
## 11 brain_region4 HHvsLL 0.57689981 0.14527611 0.5534345 0.2828294
## 12 brain_region4 LHvsLL 0.02163308 0.93874843 0.3981597 0.9387484
## 13 brain_region5 HHvsLH -0.10889220 0.69798828 0.3965573 0.8053711
## 14 brain_region5 HHvsLL -0.17518722 0.62602097 0.5077396 0.7825262
## 15 brain_region5 LHvsLL -0.06629501 0.81314497 0.3965573 0.8712268Dyad-level analysis: Relate ISC with non-binarized loneliness
dyad_ISC_lone_nobin <- function(data){
md <- summary(lmer(scale(fisher_z) ~ scale(dyad_loneliness_max) + (1|sub1) + (1|sub2), data))
B <- md$coefficients[[2,1]]
df <- md$coefficients[[2,3]]
p <- md$coefficients[[2,5]]
t <- md$coefficients[[2,4]]
p_penalized <- 2 * pt(-abs(t), df/2)
return(data.frame(B, df, p, p_penalized))
}
dyad_ISC_lone_results_nobin <- ddply(lone_isc_dyad_double, .(parcel), dyad_ISC_lone_nobin)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingulardyad_ISC_lone_results_nobin$corrected_p <- p.adjust(dyad_ISC_lone_results_nobin$p_penalized, method = "fdr")
dyad_ISC_lone_results_nobin## parcel B df p p_penalized corrected_p
## 1 brain_region1 -0.14720519 22.96618 0.17842037 0.1914425 0.4473701
## 2 brain_region2 -0.14351814 30.07026 0.25952715 0.2684220 0.4473701
## 3 brain_region3 0.09377832 88.00000 0.37931303 0.3817081 0.4771351
## 4 brain_region4 0.18752021 22.06556 0.09490692 0.1087459 0.4473701
## 5 brain_region5 -0.07356318 88.00000 0.49078504 0.4925988 0.4925988