I used imputation to fill CMOI and EMOI. The missing data is not random. And the CMOI and EMOI is not normally distributed, so I used predictive mean matching (pmm) to impute the missing data. The data is imputated from 121 to 180 subjects. And EWR can be predicted by alpha power (p = .02). For EDICT the prediction is marginally significant (p = .05). For the mediation model, the indirect effect of alpha power on EWR through CWR as mediator is significant (p = .008). However, the direct effect and total effect are insignificant. Similar indirect effect can also be found on EDICT through CDICT as mediator (p = .04), where the direct effect and total effect are insignificant.
library(bruceR)## Warning: package 'bruceR' was built under R version 4.3.3##
## bruceR (v2024.6)
## Broadly Useful Convenient and Efficient R functions
##
## Packages also loaded:
## ✔ data.table ✔ emmeans
## ✔ dplyr ✔ lmerTest
## ✔ tidyr ✔ effectsize
## ✔ stringr ✔ performance
## ✔ ggplot2 ✔ interactions
##
## Main functions of `bruceR`:
## cc() Describe() TTEST()
## add() Freq() MANOVA()
## .mean() Corr() EMMEANS()
## set.wd() Alpha() PROCESS()
## import() EFA() model_summary()
## print_table() CFA() lavaan_summary()
##
## For full functionality, please install all dependencies:
## install.packages("bruceR", dep=TRUE)
##
## Online documentation:
## https://psychbruce.github.io/bruceR
##
## To use this package in publications, please cite:
## Bao, H.-W.-S. (2024). bruceR: Broadly useful convenient and efficient R functions (Version 2024.6) [Computer software]. https://CRAN.R-project.org/package=bruceRset.wd()## ✔ Set working directory to "C:/Users/psyuser/Desktop/ARWA analysis"set.seed(6)
rm(list = ls())
library(tidyverse)## Warning: package 'readr' was built under R version 4.3.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ data.table::between() masks dplyr::between()
## ✖ Matrix::expand() masks tidyr::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ data.table::first() masks dplyr::first()
## ✖ lubridate::hour() masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag() masks stats::lag()
## ✖ data.table::last() masks dplyr::last()
## ✖ lubridate::mday() masks data.table::mday()
## ✖ lubridate::minute() masks data.table::minute()
## ✖ lubridate::month() masks data.table::month()
## ✖ Matrix::pack() masks tidyr::pack()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second() masks data.table::second()
## ✖ purrr::transpose() masks data.table::transpose()
## ✖ Matrix::unpack() masks tidyr::unpack()
## ✖ lubridate::wday() masks data.table::wday()
## ✖ lubridate::week() masks data.table::week()
## ✖ lubridate::yday() masks data.table::yday()
## ✖ lubridate::year() masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(mice)## Warning: package 'mice' was built under R version 4.3.3## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
##
## Attaching package: 'mice'
##
## The following object is masked from 'package:bruceR':
##
## cc
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## cbind, rbindlibrary(naniar)## Warning: package 'naniar' was built under R version 4.3.3data <- read.csv("./data/EEG_RS_lanExp.csv")
# Merge the PCA scores back into the original data
data_income <- data %>%
mutate(Age_EEG_yr = Age_EEG / 12) %>% # convert Age_EEG from months to years
# keep unique EEGID values
distinct(EEGID, .keep_all = TRUE) %>%
dplyr::select(CWR, CDICT, EWR, EDICT, PW, CMOI, EMOI, Age_EEG_yr, FID, EEGID, Age_EEG_yr)
# Test if the missing data is random using Little's MCAR test
mcar_test <- naniar::mcar_test(data_income[, c("CMOI", "EMOI")])
print(mcar_test)## # A tibble: 1 × 4
## statistic df p.value missing.patterns
## <dbl> <dbl> <dbl> <int>
## 1 8.30e-31 0 0 2# If the missing data is random, perform imputation using the mice package
if (mcar_test$p.value > 0.05) {
print("The missing data is at random.")
imputed_data <- mice(data_income[, c("CMOI", "EMOI")], method = "norm.predict", m = 5, maxit = 50, seed = 500)
completed_data <- complete(imputed_data, 1)
# Replace the original columns with the imputed data
data_income$CMOI <- completed_data$CMOI
data_income$EMOI <- completed_data$EMOI
} else {
print("The missing data is not completely at random.")
}## [1] "The missing data is not completely at random."# Check if CMOI and EMOI are normally distributed
shapiro_test_cmo <- shapiro.test(data_income$CMOI)
shapiro_test_emo <- shapiro.test(data_income$EMOI)
print(shapiro_test_cmo)##
## Shapiro-Wilk normality test
##
## data: data_income$CMOI
## W = 0.73803, p-value = 9.87e-16print(shapiro_test_emo)##
## Shapiro-Wilk normality test
##
## data: data_income$EMOI
## W = 0.70746, p-value < 2.2e-16# Analyze the pattern of missing data
missing_pattern <- naniar::miss_var_summary(data_income)
print(missing_pattern)## # A tibble: 10 × 3
## variable n_miss pct_miss
## <chr> <int> <num>
## 1 CMOI 39 19.3
## 2 EMOI 39 19.3
## 3 PW 16 7.92
## 4 CWR 5 2.48
## 5 EDICT 4 1.98
## 6 CDICT 3 1.49
## 7 EWR 3 1.49
## 8 Age_EEG_yr 3 1.49
## 9 FID 3 1.49
## 10 EEGID 0 0# Print the results of the normality tests
if (shapiro_test_cmo$p.value < 0.05) {
print("CMOI is not normally distributed.")
} else {
print("CMOI is normally distributed.")
}## [1] "CMOI is not normally distributed."if (shapiro_test_emo$p.value < 0.05) {
print("EMOI is not normally distributed.")
} else {
print("EMOI is normally distributed.")
}## [1] "EMOI is not normally distributed."# Choose imputation method based on normality test results
imputation_method <- ifelse(shapiro_test_cmo$p.value < 0.05 | shapiro_test_emo$p.value < 0.05, "pmm", "norm.predict")
# Perform imputation using the mice package
imputed_data <- mice(data_income[, c("CMOI", "EMOI")], method = imputation_method, m = 5, maxit = 50, seed = 500)##
## iter imp variable
## 1 1 CMOI EMOI
## 1 2 CMOI EMOI
## 1 3 CMOI EMOI
## 1 4 CMOI EMOI
## 1 5 CMOI EMOI
## 2 1 CMOI EMOI
## 2 2 CMOI EMOI
## 2 3 CMOI EMOI
## 2 4 CMOI EMOI
## 2 5 CMOI EMOI
## 3 1 CMOI EMOI
## 3 2 CMOI EMOI
## 3 3 CMOI EMOI
## 3 4 CMOI EMOI
## 3 5 CMOI EMOI
## 4 1 CMOI EMOI
## 4 2 CMOI EMOI
## 4 3 CMOI EMOI
## 4 4 CMOI EMOI
## 4 5 CMOI EMOI
## 5 1 CMOI EMOI
## 5 2 CMOI EMOI
## 5 3 CMOI EMOI
## 5 4 CMOI EMOI
## 5 5 CMOI EMOI
## 6 1 CMOI EMOI
## 6 2 CMOI EMOI
## 6 3 CMOI EMOI
## 6 4 CMOI EMOI
## 6 5 CMOI EMOI
## 7 1 CMOI EMOI
## 7 2 CMOI EMOI
## 7 3 CMOI EMOI
## 7 4 CMOI EMOI
## 7 5 CMOI EMOI
## 8 1 CMOI EMOI
## 8 2 CMOI EMOI
## 8 3 CMOI EMOI
## 8 4 CMOI EMOI
## 8 5 CMOI EMOI
## 9 1 CMOI EMOI
## 9 2 CMOI EMOI
## 9 3 CMOI EMOI
## 9 4 CMOI EMOI
## 9 5 CMOI EMOI
## 10 1 CMOI EMOI
## 10 2 CMOI EMOI
## 10 3 CMOI EMOI
## 10 4 CMOI EMOI
## 10 5 CMOI EMOI
## 11 1 CMOI EMOI
## 11 2 CMOI EMOI
## 11 3 CMOI EMOI
## 11 4 CMOI EMOI
## 11 5 CMOI EMOI
## 12 1 CMOI EMOI
## 12 2 CMOI EMOI
## 12 3 CMOI EMOI
## 12 4 CMOI EMOI
## 12 5 CMOI EMOI
## 13 1 CMOI EMOI
## 13 2 CMOI EMOI
## 13 3 CMOI EMOI
## 13 4 CMOI EMOI
## 13 5 CMOI EMOI
## 14 1 CMOI EMOI
## 14 2 CMOI EMOI
## 14 3 CMOI EMOI
## 14 4 CMOI EMOI
## 14 5 CMOI EMOI
## 15 1 CMOI EMOI
## 15 2 CMOI EMOI
## 15 3 CMOI EMOI
## 15 4 CMOI EMOI
## 15 5 CMOI EMOI
## 16 1 CMOI EMOI
## 16 2 CMOI EMOI
## 16 3 CMOI EMOI
## 16 4 CMOI EMOI
## 16 5 CMOI EMOI
## 17 1 CMOI EMOI
## 17 2 CMOI EMOI
## 17 3 CMOI EMOI
## 17 4 CMOI EMOI
## 17 5 CMOI EMOI
## 18 1 CMOI EMOI
## 18 2 CMOI EMOI
## 18 3 CMOI EMOI
## 18 4 CMOI EMOI
## 18 5 CMOI EMOI
## 19 1 CMOI EMOI
## 19 2 CMOI EMOI
## 19 3 CMOI EMOI
## 19 4 CMOI EMOI
## 19 5 CMOI EMOI
## 20 1 CMOI EMOI
## 20 2 CMOI EMOI
## 20 3 CMOI EMOI
## 20 4 CMOI EMOI
## 20 5 CMOI EMOI
## 21 1 CMOI EMOI
## 21 2 CMOI EMOI
## 21 3 CMOI EMOI
## 21 4 CMOI EMOI
## 21 5 CMOI EMOI
## 22 1 CMOI EMOI
## 22 2 CMOI EMOI
## 22 3 CMOI EMOI
## 22 4 CMOI EMOI
## 22 5 CMOI EMOI
## 23 1 CMOI EMOI
## 23 2 CMOI EMOI
## 23 3 CMOI EMOI
## 23 4 CMOI EMOI
## 23 5 CMOI EMOI
## 24 1 CMOI EMOI
## 24 2 CMOI EMOI
## 24 3 CMOI EMOI
## 24 4 CMOI EMOI
## 24 5 CMOI EMOI
## 25 1 CMOI EMOI
## 25 2 CMOI EMOI
## 25 3 CMOI EMOI
## 25 4 CMOI EMOI
## 25 5 CMOI EMOI
## 26 1 CMOI EMOI
## 26 2 CMOI EMOI
## 26 3 CMOI EMOI
## 26 4 CMOI EMOI
## 26 5 CMOI EMOI
## 27 1 CMOI EMOI
## 27 2 CMOI EMOI
## 27 3 CMOI EMOI
## 27 4 CMOI EMOI
## 27 5 CMOI EMOI
## 28 1 CMOI EMOI
## 28 2 CMOI EMOI
## 28 3 CMOI EMOI
## 28 4 CMOI EMOI
## 28 5 CMOI EMOI
## 29 1 CMOI EMOI
## 29 2 CMOI EMOI
## 29 3 CMOI EMOI
## 29 4 CMOI EMOI
## 29 5 CMOI EMOI
## 30 1 CMOI EMOI
## 30 2 CMOI EMOI
## 30 3 CMOI EMOI
## 30 4 CMOI EMOI
## 30 5 CMOI EMOI
## 31 1 CMOI EMOI
## 31 2 CMOI EMOI
## 31 3 CMOI EMOI
## 31 4 CMOI EMOI
## 31 5 CMOI EMOI
## 32 1 CMOI EMOI
## 32 2 CMOI EMOI
## 32 3 CMOI EMOI
## 32 4 CMOI EMOI
## 32 5 CMOI EMOI
## 33 1 CMOI EMOI
## 33 2 CMOI EMOI
## 33 3 CMOI EMOI
## 33 4 CMOI EMOI
## 33 5 CMOI EMOI
## 34 1 CMOI EMOI
## 34 2 CMOI EMOI
## 34 3 CMOI EMOI
## 34 4 CMOI EMOI
## 34 5 CMOI EMOI
## 35 1 CMOI EMOI
## 35 2 CMOI EMOI
## 35 3 CMOI EMOI
## 35 4 CMOI EMOI
## 35 5 CMOI EMOI
## 36 1 CMOI EMOI
## 36 2 CMOI EMOI
## 36 3 CMOI EMOI
## 36 4 CMOI EMOI
## 36 5 CMOI EMOI
## 37 1 CMOI EMOI
## 37 2 CMOI EMOI
## 37 3 CMOI EMOI
## 37 4 CMOI EMOI
## 37 5 CMOI EMOI
## 38 1 CMOI EMOI
## 38 2 CMOI EMOI
## 38 3 CMOI EMOI
## 38 4 CMOI EMOI
## 38 5 CMOI EMOI
## 39 1 CMOI EMOI
## 39 2 CMOI EMOI
## 39 3 CMOI EMOI
## 39 4 CMOI EMOI
## 39 5 CMOI EMOI
## 40 1 CMOI EMOI
## 40 2 CMOI EMOI
## 40 3 CMOI EMOI
## 40 4 CMOI EMOI
## 40 5 CMOI EMOI
## 41 1 CMOI EMOI
## 41 2 CMOI EMOI
## 41 3 CMOI EMOI
## 41 4 CMOI EMOI
## 41 5 CMOI EMOI
## 42 1 CMOI EMOI
## 42 2 CMOI EMOI
## 42 3 CMOI EMOI
## 42 4 CMOI EMOI
## 42 5 CMOI EMOI
## 43 1 CMOI EMOI
## 43 2 CMOI EMOI
## 43 3 CMOI EMOI
## 43 4 CMOI EMOI
## 43 5 CMOI EMOI
## 44 1 CMOI EMOI
## 44 2 CMOI EMOI
## 44 3 CMOI EMOI
## 44 4 CMOI EMOI
## 44 5 CMOI EMOI
## 45 1 CMOI EMOI
## 45 2 CMOI EMOI
## 45 3 CMOI EMOI
## 45 4 CMOI EMOI
## 45 5 CMOI EMOI
## 46 1 CMOI EMOI
## 46 2 CMOI EMOI
## 46 3 CMOI EMOI
## 46 4 CMOI EMOI
## 46 5 CMOI EMOI
## 47 1 CMOI EMOI
## 47 2 CMOI EMOI
## 47 3 CMOI EMOI
## 47 4 CMOI EMOI
## 47 5 CMOI EMOI
## 48 1 CMOI EMOI
## 48 2 CMOI EMOI
## 48 3 CMOI EMOI
## 48 4 CMOI EMOI
## 48 5 CMOI EMOI
## 49 1 CMOI EMOI
## 49 2 CMOI EMOI
## 49 3 CMOI EMOI
## 49 4 CMOI EMOI
## 49 5 CMOI EMOI
## 50 1 CMOI EMOI
## 50 2 CMOI EMOI
## 50 3 CMOI EMOI
## 50 4 CMOI EMOI
## 50 5 CMOI EMOIcompleted_data <- complete(imputed_data, 1)
# Replace the original columns with the imputed data
data_income$CMOI <- completed_data$CMOI
data_income$EMOI <- completed_data$EMOI
# output the data_income to csv
write.csv(data_income, file = "./data/EEG_RS_lanExp_income_imputed.csv", row.names = FALSE)# library(lme4)
library(tidyverse)
library(lmerTest)
library(conflicted)
data_income <- read.csv("./data/EEG_RS_lanExp_income_imputed.csv")
# Center the variables
data_income <- data_income %>%
mutate(
PW_c = scale(PW, center = TRUE, scale = FALSE),
CMOI_c = scale(CMOI, center = TRUE, scale = FALSE),
EMOI_c = scale(EMOI, center = TRUE, scale = FALSE),
Age_EEG_yr_c = scale(Age_EEG_yr, center = TRUE, scale = FALSE)
# EducationIncomeFactor_c = scale(EducationIncomeFactor, center = TRUE, scale = FALSE)
)
# keep only the columns data for the analysis
data_income <- data_income %>%
dplyr::select(CWR, CDICT, EWR, EDICT, PW_c, CMOI_c, EMOI_c, Age_EEG_yr_c, FID, EEGID, Age_EEG_yr) %>%
drop_na()
demo_info <- read.csv("./data/assessment-resting-state-v8-200-drop-duplicated.csv")
# left join data_income with demo_info by selected column age gender
demo_data_summary <- data_income %>%
left_join(demo_info %>% dplyr::select(EEGID, Gender), by = "EEGID", suffix = c("", "")) %>%
summarise(count = n(),
mean_age = mean(Age_EEG_yr, na.rm = TRUE),
gender_count = list(table(Gender))
)
knitr::kable(demo_data_summary)| count | mean_age | gender_count |
|---|---|---|
| 180 | 8.135648 | 73, 107 |
# Fit the models with centered variables
model_CWR_MOI <- lmerTest::lmer(CWR ~ PW_c + CMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
model_CDICT_MOI <- lmerTest::lmer(CDICT ~ PW_c + CMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
model_EWR_MOI <- lmerTest::lmer(EWR ~ PW_c + EMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
model_EDICT_MOI <- lmerTest::lmer(EDICT ~ PW_c + EMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
# Print the model summaries
summary(model_CWR_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CWR ~ PW_c + CMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1616.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.64311 -0.45677 -0.00036 0.39047 2.41967
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 563.5 23.74
## Residual 198.7 14.10
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 75.340 2.602 95.084 28.955 < 2e-16 ***
## PW_c 196.584 66.171 152.027 2.971 0.00345 **
## CMOI_c 1.528 3.014 167.998 0.507 0.61280
## Age_EEG_yr_c 21.505 2.625 102.325 8.194 7.61e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c CMOI_c
## PW_c 0.014
## CMOI_c -0.017 -0.084
## Ag_EEG_yr_c -0.011 0.017 0.138summary(model_CDICT_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDICT ~ PW_c + CMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1177.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.27085 -0.48138 0.07815 0.51146 2.03674
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 37.39 6.115
## Residual 19.41 4.406
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.3868 0.6970 94.7228 36.423 < 2e-16 ***
## PW_c 45.0838 19.4624 164.3635 2.316 0.0218 *
## CMOI_c -0.4782 0.8471 157.9740 -0.564 0.5733
## Age_EEG_yr_c 3.6772 0.7066 100.8590 5.204 1.03e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c CMOI_c
## PW_c 0.013
## CMOI_c -0.020 -0.060
## Ag_EEG_yr_c -0.016 0.021 0.145summary(model_EWR_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EWR ~ PW_c + EMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1332.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3215 -0.3717 -0.0737 0.3993 3.4203
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 134.62 11.603
## Residual 33.46 5.784
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 18.836 1.239 98.544 15.206 < 2e-16 ***
## PW_c 67.213 28.547 141.271 2.354 0.0199 *
## EMOI_c 5.502 1.303 174.494 4.221 3.90e-05 ***
## Age_EEG_yr_c 6.064 1.249 108.292 4.856 4.06e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c EMOI_c
## PW_c 0.015
## EMOI_c 0.013 0.109
## Ag_EEG_yr_c -0.007 0.009 -0.168summary(model_EDICT_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EDICT ~ PW_c + EMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1072
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3277 -0.4047 -0.1276 0.4202 2.4060
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 26.25 5.123
## Residual 8.78 2.963
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.4923 0.5591 98.9687 13.401 < 2e-16 ***
## PW_c 27.4859 14.0215 152.1534 1.960 0.0518 .
## EMOI_c 2.5207 0.6164 168.9476 4.090 6.67e-05 ***
## Age_EEG_yr_c 2.9810 0.5668 107.2404 5.259 7.44e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c EMOI_c
## PW_c 0.014
## EMOI_c 0.015 0.088
## Ag_EEG_yr_c -0.010 0.013 -0.174# Fit the models with centered variables
moderate_CWR_MOI <- lmerTest::lmer(CWR ~ PW_c * CMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
moderate_CDICT_MOI <- lmerTest::lmer(CDICT ~ PW_c * CMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
moderate_EWR_MOI <- lmerTest::lmer(EWR ~ PW_c * EMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
moderate_EDICT_MOI <- lmerTest::lmer(EDICT ~ PW_c * EMOI_c + Age_EEG_yr_c + (1 | FID), data = data_income)
# Print the model summaries
summary(moderate_CWR_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CWR ~ PW_c * CMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1602.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.68492 -0.47054 0.00031 0.39122 2.38725
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 546.5 23.38
## Residual 199.0 14.11
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 75.443 2.570 95.256 29.351 < 2e-16 ***
## PW_c 188.799 66.048 152.940 2.859 0.00485 **
## CMOI_c 1.953 2.996 164.500 0.652 0.51533
## Age_EEG_yr_c 21.746 2.597 101.921 8.375 3.14e-13 ***
## PW_c:CMOI_c 154.709 84.692 129.529 1.827 0.07004 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c CMOI_c A_EEG_
## PW_c 0.013
## CMOI_c -0.016 -0.085
## Ag_EEG_yr_c -0.010 0.014 0.141
## PW_c:CMOI_c 0.020 -0.066 0.061 0.051summary(moderate_CDICT_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDICT ~ PW_c * CMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1168.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.26309 -0.48142 0.07873 0.50976 2.02896
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 37.48 6.122
## Residual 19.57 4.424
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.3865 0.6985 93.0030 36.345 < 2e-16 ***
## PW_c 45.1748 19.5571 163.4170 2.310 0.0221 *
## CMOI_c -0.4763 0.8496 155.3345 -0.561 0.5759
## Age_EEG_yr_c 3.6766 0.7087 98.8035 5.187 1.14e-06 ***
## PW_c:CMOI_c -1.0060 25.4751 138.7221 -0.039 0.9686
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c CMOI_c A_EEG_
## PW_c 0.011
## CMOI_c -0.020 -0.061
## Ag_EEG_yr_c -0.015 0.019 0.146
## PW_c:CMOI_c 0.022 -0.059 0.030 0.046summary(moderate_EWR_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EWR ~ PW_c * EMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1323.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3193 -0.3671 -0.0676 0.4013 3.4171
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 134.90 11.615
## Residual 33.71 5.806
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 18.845 1.241 98.513 15.189 < 2e-16 ***
## PW_c 66.554 28.681 140.661 2.321 0.0218 *
## EMOI_c 5.457 1.311 172.455 4.163 4.95e-05 ***
## Age_EEG_yr_c 6.100 1.254 107.715 4.866 3.91e-06 ***
## PW_c:EMOI_c -13.509 34.048 112.805 -0.397 0.6923
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c EMOI_c A_EEG_
## PW_c 0.014
## EMOI_c 0.012 0.113
## Ag_EEG_yr_c -0.006 0.005 -0.172
## PW_c:EMOI_c -0.019 0.058 0.082 -0.067summary(moderate_EDICT_MOI)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EDICT ~ PW_c * EMOI_c + Age_EEG_yr_c + (1 | FID)
## Data: data_income
##
## REML criterion at convergence: 1064.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3198 -0.4039 -0.1283 0.4053 2.3973
##
## Random effects:
## Groups Name Variance Std.Dev.
## FID (Intercept) 26.33 5.131
## Residual 8.84 2.973
## Number of obs: 180, groups: FID, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.4884 0.5602 98.8000 13.368 < 2e-16 ***
## PW_c 27.7741 14.0823 151.2725 1.972 0.0504 .
## EMOI_c 2.5326 0.6189 166.6479 4.092 6.64e-05 ***
## Age_EEG_yr_c 2.9686 0.5690 106.6256 5.217 9.00e-07 ***
## PW_c:EMOI_c 5.9725 16.9979 120.5521 0.351 0.7259
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PW_c EMOI_c A_EEG_
## PW_c 0.013
## EMOI_c 0.014 0.091
## Ag_EEG_yr_c -0.009 0.010 -0.177
## PW_c:EMOI_c -0.020 0.055 0.058 -0.063# Perform the likelihood ratio test
comparison_CWR <- anova(model_CWR_MOI, moderate_CWR_MOI)## refitting model(s) with ML (instead of REML)# Print the comparison results
print(comparison_CWR)## Data: data_income
## Models:
## model_CWR_MOI: CWR ~ PW_c + CMOI_c + Age_EEG_yr_c + (1 | FID)
## moderate_CWR_MOI: CWR ~ PW_c * CMOI_c + Age_EEG_yr_c + (1 | FID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_CWR_MOI 6 1650.3 1669.5 -819.15 1638.3
## moderate_CWR_MOI 7 1648.9 1671.2 -817.45 1634.9 3.4013 1 0.06514 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform the likelihood ratio test
comparison_CDICT <- anova(model_CDICT_MOI, moderate_CDICT_MOI)## refitting model(s) with ML (instead of REML)# Print the comparison results
print(comparison_CDICT)## Data: data_income
## Models:
## model_CDICT_MOI: CDICT ~ PW_c + CMOI_c + Age_EEG_yr_c + (1 | FID)
## moderate_CDICT_MOI: CDICT ~ PW_c * CMOI_c + Age_EEG_yr_c + (1 | FID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_CDICT_MOI 6 1200.6 1219.8 -594.3 1188.6
## moderate_CDICT_MOI 7 1202.6 1224.9 -594.3 1188.6 6e-04 1 0.9799# Perform the likelihood ratio test
comparison_EWR <- anova(model_EWR_MOI, moderate_EWR_MOI)## refitting model(s) with ML (instead of REML)# Print the comparison results
print(comparison_EWR)## Data: data_income
## Models:
## model_EWR_MOI: EWR ~ PW_c + EMOI_c + Age_EEG_yr_c + (1 | FID)
## moderate_EWR_MOI: EWR ~ PW_c * EMOI_c + Age_EEG_yr_c + (1 | FID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_EWR_MOI 6 1359.6 1378.8 -673.80 1347.6
## moderate_EWR_MOI 7 1361.4 1383.8 -673.72 1347.4 0.1605 1 0.6887# Perform the likelihood ratio test
comparison_EDICT <- anova(model_EDICT_MOI, moderate_EDICT_MOI)## refitting model(s) with ML (instead of REML)# Print the comparison results
print(comparison_EDICT)## Data: data_income
## Models:
## model_EDICT_MOI: EDICT ~ PW_c + EMOI_c + Age_EEG_yr_c + (1 | FID)
## moderate_EDICT_MOI: EDICT ~ PW_c * EMOI_c + Age_EEG_yr_c + (1 | FID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model_EDICT_MOI 6 1093.2 1112.4 -540.63 1081.2
## moderate_EDICT_MOI 7 1095.1 1117.5 -540.56 1081.1 0.1248 1 0.7239# library(lme4)
# library(mediation)
# # Fit the mediator model (M ~ X)
# mediator_CWR <- lme4::lmer(CWR ~ PW_c + Age_EEG_yr_c + (1 | FID), data = data_income)
#
# # Fit the outcome model (Y ~ X + M)
# outcome_EWR <- lme4::lmer(EWR ~ PW_c + CWR + Age_EEG_yr_c + (1 | FID), data = data_income)
#
# # Perform mediation analysis
# mediation_EWR <- mediate(mediator_CWR, outcome_EWR, treat = "PW_c", mediator = "CWR", boot = FALSE, sims = 1000)
#
# # Summary of mediation analysis
# summary(mediation_EWR)
#
# # also for EDICT
# # Fit the mediator model (M ~ X)
# mediator_CDICT <- lme4::lmer(CDICT ~ PW_c + Age_EEG_yr_c + (1 | FID), data = data_income)
#
# # Fit the outcome model (Y ~ X + M)
# outcome_EDICT <- lme4::lmer(EDICT ~ PW_c + Age_EEG_yr_c + CDICT + (1 | FID), data = data_income)
#
# # Perform mediation analysis
# mediation_EDICT <- mediate(mediator_CDICT, outcome_EDICT, treat = "PW_c", mediator = "CDICT", boot = FALSE, sims = 1000)
#
# # Summary of mediation analysis
# summary(mediation_EDICT)
#
# ###################
# # take age as moderator
# ###################
# print('age')
# # Fit the mediator model (M ~ X)
# mediator_CWR <- lme4::lmer(CWR ~ PW_c + (1 | FID), data = data_income)
#
# # Fit the outcome model (Y ~ X + M)
# outcome_EWR <- lme4::lmer(EWR ~ PW_c + CWR * Age_EEG_yr_c + (1 | FID), data = data_income)
#
# # Perform mediation analysis
# mediation_EWR <- mediate(mediator_CWR, outcome_EWR, treat = "PW_c", mediator = "CWR", boot = FALSE, sims = 1000)
#
# # Summary of mediation analysis
# summary(mediation_EWR)
#
# # also for EDICT
# # Fit the mediator model (M ~ X)
# mediator_CDICT <- lme4::lmer(CDICT ~ PW_c + (1 | FID), data = data_income)
#
# # Fit the outcome model (Y ~ X + M)
# outcome_EDICT <- lme4::lmer(EDICT ~ PW_c * Age_EEG_yr_c + CDICT + (1 | FID), data = data_income)
#
# # Perform mediation analysis
# mediation_EDICT <- mediate(mediator_CDICT, outcome_EDICT, treat = "PW_c", mediator = "CDICT", boot = FALSE, sims = 1000)
#
# # Summary of mediation analysis
# summary(mediation_EDICT)# Function to perform regression on each quartile
perform_regression <- function(data, response_var) {
# Create quartile groups based on the response variable
data <- data %>%
mutate(quartile = cut(data[[response_var]], breaks = quantile(data[[response_var]], probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), include.lowest = TRUE, labels = FALSE))
# Initialize a list to store regression results
regression_results <- list()
# Loop through each quartile and perform regression
for (i in 1:4) {
# Subset data for the current quartile
quartile_data <- data %>%
dplyr::filter(quartile == i)
# Perform regression
model <- lm(as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")), data = quartile_data)
# Store the summary of the regression model
regression_results[[paste0("Q", i)]] <- summary(model)
}
return(regression_results)
}
# Perform regression for EWR
ewr_regression_results <- perform_regression(data_income, "EWR")
# Perform regression for EDICT
edict_regression_results <- perform_regression(data_income, "EDICT")
# Print the regression results
print("EWR Regression Results:")## [1] "EWR Regression Results:"print(ewr_regression_results)## $Q1
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8731 -0.8415 -0.1513 0.7495 2.2370
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.4174 1.6596 -1.457 0.1520
## PW_c 4.6054 5.8386 0.789 0.4343
## Age_EEG_yr 0.4989 0.2212 2.256 0.0289 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.001 on 46 degrees of freedom
## Multiple R-squared: 0.1143, Adjusted R-squared: 0.07576
## F-statistic: 2.967 on 2 and 46 DF, p-value: 0.06136
##
##
## $Q2
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7013 -3.1856 0.6701 3.5988 5.6699
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9783 6.1975 0.965 0.341
## PW_c -14.3804 20.9254 -0.687 0.496
## Age_EEG_yr 0.6426 0.7631 0.842 0.405
##
## Residual standard error: 4.052 on 38 degrees of freedom
## Multiple R-squared: 0.03666, Adjusted R-squared: -0.01404
## F-statistic: 0.723 on 2 and 38 DF, p-value: 0.4918
##
##
## $Q3
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7499 -4.1240 -0.9662 4.6029 9.2597
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.9915 7.5162 5.188 5.18e-06 ***
## PW_c -69.9848 32.5055 -2.153 0.0368 *
## Age_EEG_yr -1.5852 0.8978 -1.766 0.0844 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.845 on 44 degrees of freedom
## Multiple R-squared: 0.1369, Adjusted R-squared: 0.09767
## F-statistic: 3.489 on 2 and 44 DF, p-value: 0.03921
##
##
## $Q4
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2938 -4.1984 0.2211 3.6262 9.1810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.8162 6.0470 4.765 2.49e-05 ***
## PW_c -13.8419 31.0029 -0.446 0.6577
## Age_EEG_yr 1.3343 0.6818 1.957 0.0573 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.949 on 40 degrees of freedom
## Multiple R-squared: 0.1016, Adjusted R-squared: 0.05664
## F-statistic: 2.261 on 2 and 40 DF, p-value: 0.1174print("EDICT Regression Results:")## [1] "EDICT Regression Results:"print(edict_regression_results)## $Q1
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5495 -0.4091 -0.1989 0.5300 0.9389
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.82131 0.73382 2.482 0.0167 *
## PW_c 0.76685 2.50853 0.306 0.7612
## Age_EEG_yr -0.19202 0.09748 -1.970 0.0548 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4806 on 47 degrees of freedom
## Multiple R-squared: 0.0785, Adjusted R-squared: 0.03929
## F-statistic: 2.002 on 2 and 47 DF, p-value: 0.1464
##
##
## $Q2
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2975 -1.1405 0.5065 0.9521 1.6792
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2583 1.8822 0.137 0.8916
## PW_c -0.1608 7.9414 -0.020 0.9839
## Age_EEG_yr 0.4329 0.2363 1.832 0.0748 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.225 on 38 degrees of freedom
## Multiple R-squared: 0.08746, Adjusted R-squared: 0.03943
## F-statistic: 1.821 on 2 and 38 DF, p-value: 0.1757
##
##
## $Q3
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2910 -1.4494 -0.4522 1.6680 3.8155
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.55487 3.67053 2.331 0.0245 *
## PW_c -2.65947 13.76562 -0.193 0.8477
## Age_EEG_yr 0.09144 0.43135 0.212 0.8331
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.219 on 43 degrees of freedom
## Multiple R-squared: 0.001947, Adjusted R-squared: -0.04447
## F-statistic: 0.04195 on 2 and 43 DF, p-value: 0.959
##
##
## $Q4
##
## Call:
## lm(formula = as.formula(paste(response_var, "~ PW_c + Age_EEG_yr")),
## data = quartile_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6904 -2.4715 -0.1817 2.2104 5.1399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.6860 3.5325 3.308 0.00199 **
## PW_c -19.1471 17.7460 -1.079 0.28707
## Age_EEG_yr 0.8077 0.4009 2.015 0.05068 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.914 on 40 degrees of freedom
## Multiple R-squared: 0.1309, Adjusted R-squared: 0.08743
## F-statistic: 3.012 on 2 and 40 DF, p-value: 0.06047