PSILAUT- Analysis
Yididya
2025-06-16
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(readxl)## Warning: package 'readxl' was built under R version 4.3.3library(ggpubr)
library(dplyr)
library(afex)## Warning: package 'afex' was built under R version 4.3.3## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
##
## The following object is masked from 'package:lme4':
##
## lmerlibrary(emmeans)## Warning: package 'emmeans' was built under R version 4.3.3## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(effsize)## Warning: package 'effsize' was built under R version 4.3.3library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(rstatix)##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filterlibrary(ggplot2)
library(broom)## Warning: package 'broom' was built under R version 4.3.3library(Matrix)
library(lme4)
library(lmerTest)##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step# load data
data <- read_csv("/Users/edie/Downloads/PSILAUT_DMN_data_YK_20250427_TW.csv")## Rows: 125 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): ID, DMNDMN, Mean_FD
## dbl (4): Dose, Group, Age, P_Naive
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.str(data)## spc_tbl_ [125 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ID : chr [1:125] "PSF001" "PSF001" "PSF002" "PSF002" ...
## $ Dose : num [1:125] 0 5 0 5 0 5 5 0 0 5 ...
## $ Group : num [1:125] 0 0 0 0 0 0 0 0 0 0 ...
## $ Age : num [1:125] 21 21 23 23 19 19 29 29 26 26 ...
## $ P_Naive: num [1:125] 1 NA 1 NA 1 NA 0 NA 1 NA ...
## $ DMNDMN : chr [1:125] "0.452838799" "0.157142269" "0.156038518" "0.254695207" ...
## $ Mean_FD: chr [1:125] "0.096" "0.27" "0.134" "0.092" ...
## - attr(*, "spec")=
## .. cols(
## .. ID = col_character(),
## .. Dose = col_double(),
## .. Group = col_double(),
## .. Age = col_double(),
## .. P_Naive = col_double(),
## .. DMNDMN = col_character(),
## .. Mean_FD = col_character()
## .. )
## - attr(*, "problems")=<externalptr>summary(data)## ID Dose Group Age
## Length:125 Min. :0.00 Min. :0.00 Min. :18.00
## Class :character 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:23.00
## Mode :character Median :5.00 Median :1.00 Median :28.00
## Mean :2.52 Mean :0.56 Mean :29.32
## 3rd Qu.:5.00 3rd Qu.:1.00 3rd Qu.:35.00
## Max. :5.00 Max. :1.00 Max. :59.00
##
## P_Naive DMNDMN Mean_FD
## Min. :0.0000 Length:125 Length:125
## 1st Qu.:1.0000 Class :character Class :character
## Median :1.0000 Mode :character Mode :character
## Mean :0.7903
## 3rd Qu.:1.0000
## Max. :1.0000
## NA's :63# Convert character columns to numeric
data$DMNDMN <- as.numeric(data$DMNDMN)## Warning: NAs introduced by coerciondata$Mean_FD <- as.numeric(data$Mean_FD)## Warning: NAs introduced by coercion# Standardize missing values before conversion
data$DMNDMN[data$DMNDMN %in% c("", "NA", "NaN")] <- NA
data$Mean_FD[data$Mean_FD %in% c("", "NA", "NaN")] <- NA
#Convert from character to numeric
data$DMNDMN <- as.numeric(data$DMNDMN)
data$Mean_FD <- as.numeric(data$Mean_FD)
# Check how many missing values
cat("DMNDMN missing:", sum(is.na(data$DMNDMN)), "\n")## DMNDMN missing: 14cat("Mean_FD missing:", sum(is.na(data$Mean_FD)), "\n")## Mean_FD missing: 14both DMNDMD and Mean_FD have 14 missing values
# Convert Dose to factor: 0 = Placebo, 5 = Treatment
data$Dose <- factor(data$Dose,
levels = c(0, 5),
labels = c("Placebo", "Treatment"))
# Convert Group to factor: 0 = Control, 1 = ASC
data$Group <- factor(data$Group,
levels = c(0, 1),
labels = c("Control", "ASC"))
# Convert P_Naive to factor: 1 = Yes, 0 = No
data$P_Naive <- factor(data$P_Naive,
levels = c(0, 1),
labels = c("No", "Yes"))## Pairing the data
# Extract demo variables (one row per ID)
demo_data <- data %>%
select(ID, Group, Age, P_Naive) %>%
distinct()
# Pivot DMNDMN wide by Dose
wide_dmndmn <- data %>%
select(ID, Dose, DMNDMN) %>%
pivot_wider(names_from = Dose, values_from = DMNDMN)
# Pivot Mean_FD wide by Dose
wide_meanfd <- data %>%
select(ID, Dose, Mean_FD) %>%
pivot_wider(names_from = Dose, values_from = Mean_FD)
# Join all together: demo + DMNDMN + Mean_FD
wide_data <- demo_data %>%
left_join(wide_dmndmn, by = "ID") %>%
left_join(wide_meanfd, by = "ID", suffix = c("_DMNDMN", "_Mean_FD"))
# Inspect the wide data
head(wide_data)## # A tibble: 6 × 8
## ID Group Age P_Naive Placebo_DMNDMN Treatment_DMNDMN Placebo_Mean_FD
## <chr> <fct> <dbl> <fct> <dbl> <dbl> <dbl>
## 1 PSF001 Control 21 Yes 0.453 0.157 0.096
## 2 PSF001 Control 21 <NA> 0.453 0.157 0.096
## 3 PSF002 Control 23 Yes 0.156 0.255 0.134
## 4 PSF002 Control 23 <NA> 0.156 0.255 0.134
## 5 PSF003 Control 19 Yes 0.136 0.116 0.084
## 6 PSF003 Control 19 <NA> 0.136 0.116 0.084
## # ℹ 1 more variable: Treatment_Mean_FD <dbl># Summary to see missing values, etc.
summary(wide_data)## ID Group Age P_Naive Placebo_DMNDMN
## Length:124 Control:55 Min. :18.0 No :13 Min. :0.0901
## Class :character ASC :69 1st Qu.:23.0 Yes :49 1st Qu.:0.1645
## Mode :character Median :28.0 NA's:62 Median :0.2099
## Mean :29.4 Mean :0.2237
## 3rd Qu.:35.0 3rd Qu.:0.2585
## Max. :59.0 Max. :0.4860
## NA's :9
## Treatment_DMNDMN Placebo_Mean_FD Treatment_Mean_FD
## Min. :0.08252 Min. :0.0620 Min. :0.0580
## 1st Qu.:0.16169 1st Qu.:0.0920 1st Qu.:0.0920
## Median :0.22041 Median :0.1090 Median :0.1250
## Mean :0.21578 Mean :0.1223 Mean :0.1422
## 3rd Qu.:0.25997 3rd Qu.:0.1365 3rd Qu.:0.1720
## Max. :0.38655 Max. :0.3690 Max. :0.3720
## NA's :20 NA's :9 NA's :20# **Hypothesis I**
# Testing normality of DMN under placepo for each group
ggqqplot(wide_data, x = "Placebo_DMNDMN", facet.by = "Group",
title = "QQ Plots of Placebo DMN Connectivity by Group")## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_qq()`).## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_qq_line()`).
## Removed 9 rows containing non-finite outside the scale range
## (`stat_qq_line()`).
ggplot(wide_data, aes(x = Placebo_DMNDMN, fill = Group)) +
geom_histogram(aes(y = ..density..), bins = 30, alpha = 0.4, position = 'identity') +
stat_density(aes(color = Group), geom = "line", position = "identity", size = 1) +
facet_wrap(~Group) +
labs(title = "Placebo DMN Connectivity Histograms with Density Curves by Group",
x = "Placebo DMN Connectivity (DMNDMN)",
y = "Density") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_bin()`).## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_density()`).
# Subset data by Group and run Shapiro-Wilk test for normality
control_dmn <- wide_data$Placebo_DMNDMN[wide_data$Group == "Control"]
asc_dmn <- wide_data$Placebo_DMNDMN[wide_data$Group == "ASC"]
# Shapiro-Wilk tests
shapiro_control <- shapiro.test(control_dmn)
shapiro_asc <- shapiro.test(asc_dmn)
# Print results
cat("Shapiro-Wilk test for Control group:\n")## Shapiro-Wilk test for Control group:print(shapiro_control)##
## Shapiro-Wilk normality test
##
## data: control_dmn
## W = 0.84929, p-value = 2.946e-05cat("\nShapiro-Wilk test for ASC group:\n")##
## Shapiro-Wilk test for ASC group:print(shapiro_asc)##
## Shapiro-Wilk normality test
##
## data: asc_dmn
## W = 0.87037, p-value = 3.608e-06Shapiro-Wilk tests were conducted to assess the normality of placebo DMN connectivity values within each group. The test indicated a significant deviation from normality in both the Control group (W = 0.85, p < 0.001) and the ASC group (W = 0.87, p < 0.001). Visual inspection of the distribution plots further confirmed that the data do not follow a normal bell-shaped curve in either group. Therefore, the assumption of normality was violated for both groups.
# Wilcoxon rank-sum test
wilcox.test(Placebo_DMNDMN ~ Group, data = wide_data)##
## Wilcoxon rank sum test with continuity correction
##
## data: Placebo_DMNDMN by Group
## W = 1886, p-value = 0.08832
## alternative hypothesis: true location shift is not equal to 0Hypothesis I
A Wilcoxon rank-sum test was conducted to compare placebo DMN connectivity between ASC and Control groups. The results showed no statistically significant difference between the groups (W = 1886, p = 0.088), indicating that baseline DMN connectivity does not differ significantly between ASC and non-ASC individuals under the placebo condition.
## DMN connectivity differs between ASC and Control groups while controlling for age
# Filter data for placebo only
placebo_data <- data %>% filter(Dose == "Placebo")
# Fit linear model controlling for Age
model_placebo <- lm(DMNDMN ~ Group + Age, data = placebo_data)
# Check summary
summary(model_placebo)##
## Call:
## lm(formula = DMNDMN ~ Group + Age, data = placebo_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.138572 -0.054398 -0.003636 0.037384 0.266765
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.275400 0.040955 6.724 1.07e-08 ***
## GroupASC -0.022389 0.022538 -0.993 0.325
## Age -0.001352 0.001347 -1.003 0.320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08335 on 55 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.0396, Adjusted R-squared: 0.004676
## F-statistic: 1.134 on 2 and 55 DF, p-value: 0.3292A linear regression was conducted to examine whether baseline (placebo) DMN connectivity differed between ASC and Control groups, controlling for age. The model was not statistically significant, \(F(2, 55) = 1.13, p = 0.33\), explaining only about 4% of the variance in DMN connectivity.
Neither group status (\(\beta = -0.022, SE = 0.023, p = 0.33\)) nor age (\(\beta = -0.0014, SE = 0.0013, p = 0.32\)) significantly predicted DMN connectivity under placebo.
These results suggest no significant difference in baseline DMN connectivity between ASC and non-ASC individuals, even after controlling for age.
## Fit linear model controlling for Group, Age, and Mean_FD under Placebo
model_placebo_fd <- lm(DMNDMN ~ Group + Age + Mean_FD, data = placebo_data)
# View summary
summary(model_placebo_fd)##
## Call:
## lm(formula = DMNDMN ~ Group + Age + Mean_FD, data = placebo_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.149322 -0.051625 -0.007762 0.038627 0.257798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.297470 0.046775 6.360 4.5e-08 ***
## GroupASC -0.014980 0.023786 -0.630 0.531
## Age -0.001404 0.001349 -1.041 0.303
## Mean_FD -0.202255 0.206813 -0.978 0.332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08338 on 54 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.05631, Adjusted R-squared: 0.003886
## F-statistic: 1.074 on 3 and 54 DF, p-value: 0.3678A linear regression was conducted to examine whether baseline (placebo) DMN connectivity differed between ASC and Control groups, controlling for age and Mean_FD. The overall model was not significant, \(F(3, 54) = 1.07, p = 0.37\), explaining approximately 6% of the variance.
None of the predictors significantly influenced DMN connectivity at baseline:
- Group status (\(\beta = -0.015, SE = 0.024, p = 0.53\))
- Age (\(\beta = -0.0014, SE = 0.0013, p = 0.30\))
- Mean_FD (\(\beta = -0.20, SE = 0.21, p = 0.33\)).
These results suggest no significant baseline difference in DMN connectivity between ASC and non-ASC individuals when controlling for age and head motion (Mean_FD).
## Violin plot comparing DMN connectivity under placebo.
# Filter placebo data
placebo_data <- data %>% filter(Dose == "Placebo")
ggplot(placebo_data, aes(x = Group, y = DMNDMN, fill = Group)) +
geom_violin(trim = FALSE, alpha = 0.5) +
geom_jitter(width = 0.1, size = 1, alpha = 0.7) +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "white") +
labs(title = "Baseline DMN Connectivity by Group under Placebo",
y = "DMN Connectivity (DMNDMN)") +
theme_minimal() +
theme(legend.position = "none")## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_ydensity()`).## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_summary()`).## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).
# **Hypothesis II**
# Fit the model
model <- lmer(DMNDMN ~ Dose * Group + (1 | ID), data = data)
# Summary with p-values
summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DMNDMN ~ Dose * Group + (1 | ID)
## Data: data
##
## REML criterion at convergence: -242.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4762 -0.5653 -0.1193 0.4839 3.0589
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.001826 0.04273
## Residual 0.003818 0.06179
## Number of obs: 111, groups: ID, 63
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.23813 0.01552 102.88206 15.341 <2e-16 ***
## DoseTreatment -0.04045 0.01862 54.56298 -2.173 0.0341 *
## GroupASC -0.02507 0.02006 100.86246 -1.250 0.2142
## DoseTreatment:GroupASC 0.05550 0.02438 52.01098 2.276 0.0270 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DsTrtm GrpASC
## DoseTretmnt -0.616
## GroupASC -0.774 0.476
## DsTrtm:GASC 0.470 -0.763 -0.587# tests differences within each group/doseA linear mixed-effects model was conducted to examine the interaction between drug condition (Dose: Placebo vs Treatment) and group (ASC vs Control) on DMN connectivity, accounting for repeated measures by including a random intercept for participants.
The fixed effect of DoseTreatment was significant (estimate = -0.040, p = 0.034), indicating that treatment with psilocybin led to a decrease in DMN connectivity across groups compared to placebo.
The fixed effect of GroupASC was not significant (p = 0.214), suggesting no baseline difference in DMN connectivity between ASC and Control groups.
Importantly, the DoseTreatment × GroupASC interaction was significant (estimate = 0.056, p = 0.027), indicating that the effect of psilocybin on DMN connectivity differed between the ASC and Control groups.
This significant interaction supports the hypothesis of a drug-by-group effect, meaning the treatment influenced DMN connectivity differently depending on group status.
# fixed effects
emm <- emmeans(model, ~ Dose * Group)
pairs(emm, simple = "each") ## $`simple contrasts for Dose`
## Group = Control:
## contrast estimate SE df t.ratio p.value
## Placebo - Treatment 0.0405 0.0187 57.8 2.163 0.0347
##
## Group = ASC:
## contrast estimate SE df t.ratio p.value
## Placebo - Treatment -0.0150 0.0158 51.8 -0.953 0.3450
##
## Degrees-of-freedom method: kenward-roger
##
## $`simple contrasts for Group`
## Dose = Placebo:
## contrast estimate SE df t.ratio p.value
## Control - ASC 0.0251 0.0201 102 1.248 0.2150
##
## Dose = Treatment:
## contrast estimate SE df t.ratio p.value
## Control - ASC -0.0304 0.0206 103 -1.476 0.1431
##
## Degrees-of-freedom method: kenward-rogerThe fixed effects table shows:
DoseTreatment: Significant main effect (estimate = -0.0405, p = 0.034). On average, DMN connectivity decreases under treatment compared to placebo.
GroupASC: No significant main effect (estimate = -0.0251, p = 0.214), so no baseline difference between ASC and Control.
DoseTreatment × GroupASC interaction: Significant (estimate = 0.0555, p = 0.027), meaning the treatment effect differs between groups.
The simple contrasts clarify the interaction:
In Control group, treatment significantly reduces DMN connectivity compared to placebo (estimate = 0.0405, p = 0.035).
In ASC group, the treatment effect is not significant (estimate = -0.0150, p = 0.345).
Comparing groups at each dose shows no significant difference in either placebo or treatment condition.
## Refit model including Age and Mean_FD
model_cov <- lmer(DMNDMN ~ Dose * Group + Age + Mean_FD + (1 | ID), data = data)
summary(model_cov)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DMNDMN ~ Dose * Group + Age + Mean_FD + (1 | ID)
## Data: data
##
## REML criterion at convergence: -230.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.50712 -0.57042 -0.08804 0.49365 2.99113
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.001886 0.04343
## Residual 0.003733 0.06109
## Number of obs: 111, groups: ID, 63
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.2720495 0.0338640 68.9082124 8.034 1.75e-11 ***
## DoseTreatment -0.0348656 0.0187062 54.8737296 -1.864 0.0677 .
## GroupASC -0.0169857 0.0205571 97.9175352 -0.826 0.4107
## Age -0.0004844 0.0009635 56.6667166 -0.503 0.6171
## Mean_FD -0.2003545 0.1215827 97.1987954 -1.648 0.1026
## DoseTreatment:GroupASC 0.0536345 0.0241556 50.4952034 2.220 0.0309 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DsTrtm GrpASC Age Men_FD
## DoseTretmnt -0.171
## GroupASC -0.199 0.493
## Age -0.811 -0.056 -0.081
## Mean_FD -0.397 -0.164 -0.219 0.041
## DsTrtm:GASC 0.178 -0.759 -0.575 0.029 0.034# Get estimated marginal means for Dose * Group, controlling for Age and Mean_FD at their mean
emm_cov <- emmeans(model_cov, ~ Dose * Group)
# Plot adjusted means
plot(emm_cov, comparison = TRUE)
The interaction between Dose and Group is significant (p = 0.031), meaning the effect of Treatment on DMN connectivity depends on whether someone is ASC or Control, even after controlling for age and motion.
The main effect of Dose alone shows a trend but is not significant.
Group alone (ASC vs Control at placebo) and covariates (Age, Mean_FD) are not significant predictors on their own.
Random effects show variability between subjects (ID), as expected for paired data.