See Matlab script get_FNC_corrs.m under ~/restingstate/data/derivatives/gift/analyses-dissertation/batch-jul-19.
library(tidyverse)
Registered S3 method overwritten by 'rvest':
method from
read_xml.response xml2
[30m── [1mAttaching packages[22m ─────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──[39m
[30m[32m✔[30m [34mggplot2[30m 3.2.0 [32m✔[30m [34mpurrr [30m 0.3.2
[32m✔[30m [34mtibble [30m 2.1.3 [32m✔[30m [34mdplyr [30m 0.8.1
[32m✔[30m [34mtidyr [30m 0.8.3 [32m✔[30m [34mstringr[30m 1.4.0
[32m✔[30m [34mreadr [30m 1.3.1 [32m✔[30m [34mforcats[30m 0.4.0[39m
[30m── [1mConflicts[22m ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31m✖[30m [34m.GlobalEnv[30m::[32mfilter()[30m masks [34mdplyr[30m::filter(), [34mstats[30m::filter()
[31m✖[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(psych)
Attaching package: ‘psych’
The following objects are masked from ‘package:ggplot2’:
%+%, alpha
library(sjPlot)
Learn more about sjPlot with 'browseVignettes("sjPlot")'.
library(sjmisc)
Attaching package: ‘sjmisc’
The following object is masked from ‘package:purrr’:
is_empty
The following object is masked from ‘package:tidyr’:
replace_na
The following object is masked from ‘package:tibble’:
add_case
library(nlme)
Attaching package: ‘nlme’
The following object is masked from ‘package:dplyr’:
collapse
library(jtools)
Attaching package: ‘jtools’
The following objects are masked from ‘package:sjmisc’:
%nin%, center
library(mediation)
Loading required package: MASS
Attaching package: ‘MASS’
The following object is masked _by_ ‘.GlobalEnv’:
select
The following object is masked from ‘package:dplyr’:
select
Loading required package: Matrix
Attaching package: ‘Matrix’
The following object is masked from ‘package:tidyr’:
expand
Loading required package: mvtnorm
Loading required package: sandwich
mediation: Causal Mediation Analysis
Version: 4.4.7
Attaching package: ‘mediation’
The following object is masked from ‘package:psych’:
mediate
library(effects)
Loading required package: carData
lattice theme set by effectsTheme()
See ?effectsTheme for details.
library(reghelper)
Attaching package: ‘reghelper’
The following object is masked from ‘package:psych’:
ICC
The following object is masked from ‘package:base’:
beta
select <- dplyr::select
filter <- dplyr::filter
# read in static FNC stats
## FNC = correlation between a pair of ICs
fnc_pl <- read.csv("/Users/sarenseeley/Dropbox/Dissertation/fMRI data analysis/data/fnc_corrs_networks_pl.csv")
fnc_ot <- read.csv("/Users/sarenseeley/Dropbox/Dissertation/fMRI data analysis/data/fnc_corrs_networks_ot.csv")
# read in subject IDs
subj <- read.csv("/Users/sarenseeley/Dropbox/Dissertation/fMRI data analysis/data/subj.csv", header=F)
subj <- subj %>% transmute(ID = as.character(subj$V1))
subj$ID <- paste0('D', subj$ID) # add the "D" to make it consistent with the master dataset and other data
# drop the "pl/ot" suffix and add a new column "tx" for the long dataframe
fnc_pl <- fnc_pl %>%
rename_all(.funs = funs(sub("\\_pl", "", names(fnc_pl))))
fnc_pl$tx <- "placebo"
names(fnc_pl)
fnc_pl <- bind_cols(subj, fnc_pl) # add the "ID" variable
# drop the "pl/ot" suffix and add a new column "tx" for the long dataframe
fnc_ot <- fnc_ot %>%
rename_all(.funs = funs(sub("\\_ot", "", names(fnc_ot))))
fnc_ot$tx <- "oxytocin"
names(fnc_ot)
fnc_ot <- bind_cols(subj, fnc_ot) # add the "ID" variable
# create the long FNC dataframe
long_fnc <- bind_rows(fnc_pl, fnc_ot)
long_fnc$X <- NULL # remove the "X" variable
# should have 76 (38 subjects*2 sessions) observations of 30 variables (28 IC pairs + ID + tx)
These include mean framewise displacement, time elapsed from spray administration to when the resting state sequence began, PANAS, STAI-S, “subjective rating of experience” post-scan questions.
# read in MRQC group stats to get mean FD
meanfd <- read_tsv("~/Desktop/restingstate/data/derivatives/mriqc/group_bold.tsv", col_names = TRUE, col_types=cols_only(
bids_name = col_character(),
fd_mean = col_number()))
# split up meanfd by session
meanfd_pl <- dplyr::filter(meanfd, grepl("ses-txA",bids_name))
meanfd_ot <- dplyr::filter(meanfd, grepl("ses-txB",bids_name))
# then rename for merging with master dataset variables
meanfd_pl <- rename(meanfd_pl, ID = bids_name)
meanfd_pl$tx <- "placebo"
meanfd_ot <- rename(meanfd_ot, ID = bids_name)
meanfd_ot$tx <- "oxytocin"
meanfd_pl <- mutate_if(tibble::as_tibble(meanfd_pl),
is.character,
stringr::str_replace_all, pattern = "sub-", replacement = "D")
meanfd_pl <- mutate_if(tibble::as_tibble(meanfd_pl),
is.character,
stringr::str_replace_all, pattern = "_ses-txA_task-rest_bold", replacement = "")
meanfd_ot <- mutate_if(tibble::as_tibble(meanfd_ot),
is.character,
stringr::str_replace_all, pattern = "sub-", replacement = "D")
meanfd_ot <- mutate_if(tibble::as_tibble(meanfd_ot),
is.character,
stringr::str_replace_all, pattern = "_ses-txB_task-rest_bold", replacement = "")
long_meanfd <- bind_rows(meanfd_ot, meanfd_pl)
# read in master dataset
data <- readRDS("~/Dropbox/GLASS Lab/OT Study/data/master-dataset/ot-fmri_master-dataset_020719.rds")
# note that these data are in wide format
# in particular, need to pull out the state variables and combine into one long variable as with meanfd above
# placebo session:
pl_vars <- subset(data, select=c(ID,
tot_pre_panas_pa_A,
tot_pre_panas_na_A,
tot_pre_stai_A,
tot_post_panas_pa_A,
tot_post_panas_na_A,
tot_post_stai_A,
timetorest_A))
pl_vars$tx <- "placebo"
pl_vars <- pl_vars %>%
rename_all(.funs = funs(sub("\\_A", "", names(pl_vars)))) # remove the "_A" from column names
# oxytocin session:
ot_vars <- subset(data, select=c(ID,
tot_pre_panas_pa_B,
tot_pre_panas_na_B,
tot_pre_stai_B,
tot_post_panas_pa_B,
tot_post_panas_na_B,
tot_post_stai_B,
timetorest_B))
ot_vars$tx <- "oxytocin"
ot_vars <- ot_vars %>%
rename_all(.funs = funs(sub("\\_B", "", names(ot_vars)))) # remove the "_B" from column names
# bind the oxytocin and placebo session data together into long format
long_state <- bind_rows(pl_vars, ot_vars)
# merge the state variables with mean FD by ID and tx
long_state <- left_join(long_state, long_meanfd, by=c("ID", "tx"))
long_state_fnc <- left_join(long_fnc, long_state, by=c("ID", "tx"))
long_state_fnc <- long_state_fnc %>% mutate(tx = recode_factor(tx, "placebo" = "placebo", "oxytocin" = "oxytocin", .ordered = TRUE))
### Self-reported post-scan thoughts
# using the % of content that was spouse/death-related from their post scan reports as a very imperfect measure of ruminative content for Aim 1 H3
# (this will probably change)
rumPL <- read.csv("/Users/sarenseeley/Dropbox/Grants\ and\ Applications/NRSA\ F31\ 2018-2019/Resubmission\ Reviews/power-analysis/rumination_txA.csv") # placebo
rumOT <- read.csv("/Users/sarenseeley/Dropbox/Grants\ and\ Applications/NRSA\ F31\ 2018-2019/Resubmission\ Reviews/power-analysis/rumination_txB.csv") # oxytocin
# rename post_sroe variables for binding rows and make factors into characters
# (otherwise they get all screwed up, w/warning "Unequal factor levels: coercing to characterbinding character and factor vector")
rumPL <- rumPL %>% rename(post_sroe_1 = post_sroe_1_A,
post_sroe_2 = post_sroe_2_A) %>% mutate(
ID = as.character(ID),
tx = as.character(tx),
post_sroe_1 = as.character(post_sroe_1),
post_sroe_2 = as.character(post_sroe_2)
) %>% mutate(ID = str_replace(ID, "sub-", "D")) # replace "sub-" prefix in ID to "D" so that the ID variables match across dataframes to merge later
rumPL$X <- NULL
rumOT <- rumOT %>% rename(post_sroe_1 = post_sroe_1_B,
post_sroe_2 = post_sroe_2_B) %>% mutate(
ID = as.character(ID),
tx = as.character(tx),
post_sroe_1 = as.character(post_sroe_1),
post_sroe_2 = as.character(post_sroe_2)
) %>%
mutate(ID = str_replace(ID, "sub-", "D")) # replace "sub-" prefix in ID to "D" so that the ID variables match across dataframes to merge later
rum <- bind_rows(rumPL, rumOT)
# recode "tx" variable
rum <- rum %>% mutate(tx = recode_factor(tx, "txA" = "placebo", "txB" = "oxytocin", .ordered = TRUE))
## this will only be 75 observations, not 76, because of the one participant (D115) who received the pre-scan Qualtrics instead of the post-scan Qualtrics at one visit, so we have post PANAS and STAI but not SROE responses from him at the oxytocin session
# create composite variables for (1) focus on spouse in post-scan responses, and (2) of all of the things they reported thinking about, what % of what they mentioned was related to their spouse/death?
rum$focus_spouse <- (rum$spouse_1+rum$spouse_2)
rum$focus_task <- (rum$task_1+rum$task_2)
rum$focus_whoto <- (rum$whoto_1+rum$whoto_2)
rum$focus_otherpics <- (rum$otherpics_1+rum$otherpics_2)
rum$focus_all <- (rum$spouse_1+rum$spouse_2+rum$task_1+rum$task_2+rum$whoto_1+rum$whoto_2+rum$otherpics_1+rum$otherpics_2+rum$somatic_1+rum$somatic_2+rum$metacog_1+rum$metacog_2)
rum$focus_spouse_perc <- (rum$focus_spouse/rum$focus_all) # of the total number of things they were thinking about, what proportion were related to their spouse in some way (positive or negative)?
rum$focus_nonspouse_perc <- (1-(rum$focus_spouse_perc)) # of all the total number of things they were thinking about, what proportion of thought content was NOT related to their spouse?
# drop extraneous variables
rum <- rum %>% dplyr::select(-c(post_sroe_1, spouse_1, task_1, whoto_1, otherpics_1, somatic_1, metacog_1, post_sroe_2, spouse_2, task_2, whoto_2, otherpics_2, somatic_2, metacog_2))
# replace NaNs (generated by dividing zero by zero, for D129 placebo session) with 0
fix_nan <- function(x){
x[is.nan(x)] <- 0
x
}
rum$focus_spouse_perc <- fix_nan(rum$focus_spouse_perc)
rum$focus_nonspouse_perc <- fix_nan(rum$focus_nonspouse_perc)
# merge with data from master dataset
long_state_fnc_sroe <- left_join(long_state_fnc, rum, by=c("tx", "ID"))
These include participant demographics and characteristics, and trait measures.
# read in master dataset
data <- readRDS("~/Dropbox/GLASS Lab/OT Study/data/master-dataset/ot-fmri_master-dataset_020719.rds")
# add the variables from the master dataset to the long dataset containing the state and FNC variables,
# dropping the wide-format state variables as well as startdate (when they did the Qualtrics), DOB, and DOD
data <- subset(data, select=-c(startdate,
dob,
dateofdeath,
pre_panas_1_v1,
pre_panas_2_v1,
pre_panas_3_v1,
pre_panas_4_v1,
pre_panas_5_v1,
pre_panas_6_v1,
pre_panas_7_v1,
pre_panas_8_v1,
pre_panas_9_v1,
pre_panas_10_v1,
pre_panas_11_v1,
pre_panas_12_v1,
pre_panas_13_v1,
pre_panas_14_v1,
pre_panas_15_v1,
pre_panas_16_v1,
pre_panas_17_v1,
pre_panas_18_v1,
pre_panas_19_v1,
pre_panas_20_v1,
pre_stai_1_v1,
pre_stai_2_v1,
pre_stai_3_v1,
pre_stai_4_v1,
pre_stai_5_v1,
pre_stai_6_v1,
pre_stai_7_v1,
pre_stai_8_v1,
pre_stai_9_v1,
pre_stai_10_v1,
pre_stai_11_v1,
pre_stai_12_v1,
pre_stai_13_v1,
pre_stai_14_v1,
pre_stai_15_v1,
pre_stai_16_v1,
pre_stai_17_v1,
pre_stai_18_v1,
pre_stai_19_v1,
pre_stai_20_v1,
tot_pre_panas_pa_v1,
tot_pre_panas_na_v1,
pre_stai_1r_v1,
pre_stai_2r_v1,
pre_stai_5r_v1,
pre_stai_8r_v1,
pre_stai_10r_v1,
pre_stai_11r_v1,
pre_stai_15r_v1,
pre_stai_16r_v1,
pre_stai_19r_v1,
pre_stai_20r_v1,
tot_pre_stai_v1,
pre_panas_1_v2,
pre_panas_2_v2,
pre_panas_3_v2,
pre_panas_4_v2,
pre_panas_5_v2,
pre_panas_6_v2,
pre_panas_7_v2,
pre_panas_8_v2,
pre_panas_9_v2,
pre_panas_10_v2,
pre_panas_11_v2,
pre_panas_12_v2,
pre_panas_13_v2,
pre_panas_14_v2,
pre_panas_15_v2,
pre_panas_16_v2,
pre_panas_17_v2,
pre_panas_18_v2,
pre_panas_19_v2,
pre_panas_20_v2,
pre_stai_1_v2,
pre_stai_2_v2,
pre_stai_3_v2,
pre_stai_4_v2,
pre_stai_5_v2,
pre_stai_6_v2,
pre_stai_7_v2,
pre_stai_8_v2,
pre_stai_9_v2,
pre_stai_10_v2,
pre_stai_11_v2,
pre_stai_12_v2,
pre_stai_13_v2,
pre_stai_14_v2,
pre_stai_15_v2,
pre_stai_16_v2,
pre_stai_17_v2,
pre_stai_18_v2,
pre_stai_19_v2,
pre_stai_20_v2,
tot_pre_panas_pa_v2,
tot_pre_panas_na_v2,
pre_stai_1r_v2,
pre_stai_2r_v2,
pre_stai_5r_v2,
pre_stai_8r_v2,
pre_stai_10r_v2,
pre_stai_11r_v2,
pre_stai_15r_v2,
pre_stai_16r_v2,
pre_stai_19r_v2,
pre_stai_20r_v2,
tot_pre_stai_v2,
post_panas_1_v1,
post_panas_2_v1,
post_panas_3_v1,
post_panas_4_v1,
post_panas_5_v1,
post_panas_6_v1,
post_panas_7_v1,
post_panas_8_v1,
post_panas_9_v1,
post_panas_10_v1,
post_panas_11_v1,
post_panas_12_v1,
post_panas_13_v1,
post_panas_14_v1,
post_panas_15_v1,
post_panas_16_v1,
post_panas_17_v1,
post_panas_18_v1,
post_panas_19_v1,
post_panas_20_v1,
post_stai_1_v1,
post_stai_2_v1,
post_stai_3_v1,
post_stai_4_v1,
post_stai_5_v1,
post_stai_6_v1,
post_stai_7_v1,
post_stai_8_v1,
post_stai_9_v1,
post_stai_10_v1,
post_stai_11_v1,
post_stai_12_v1,
post_stai_13_v1,
post_stai_14_v1,
post_stai_15_v1,
post_stai_16_v1,
post_stai_17_v1,
post_stai_18_v1,
post_stai_19_v1,
post_stai_20_v1,
post_sroe_1_v1,
post_sroe_2_v1,
post_sroe_3_v1,
tot_post_panas_pa_v1,
tot_post_panas_na_v1,
post_stai_1r_v1,
post_stai_2r_v1,
post_stai_5r_v1,
post_stai_8r_v1,
post_stai_10r_v1,
post_stai_11r_v1,
post_stai_15r_v1,
post_stai_16r_v1,
post_stai_19r_v1,
post_stai_20r_v1,
tot_post_stai_v1,
post_panas_1_v2,
post_panas_2_v2,
post_panas_3_v2,
post_panas_4_v2,
post_panas_5_v2,
post_panas_6_v2,
post_panas_7_v2,
post_panas_8_v2,
post_panas_9_v2,
post_panas_10_v2,
post_panas_11_v2,
post_panas_12_v2,
post_panas_13_v2,
post_panas_14_v2,
post_panas_15_v2,
post_panas_16_v2,
post_panas_17_v2,
post_panas_18_v2,
post_panas_19_v2,
post_panas_20_v2,
post_stai_1_v2,
post_stai_2_v2,
post_stai_3_v2,
post_stai_4_v2,
post_stai_5_v2,
post_stai_6_v2,
post_stai_7_v2,
post_stai_8_v2,
post_stai_9_v2,
post_stai_10_v2,
post_stai_11_v2,
post_stai_12_v2,
post_stai_13_v2,
post_stai_14_v2,
post_stai_15_v2,
post_stai_16_v2,
post_stai_17_v2,
post_stai_18_v2,
post_stai_19_v2,
post_stai_20_v2,
tot_post_panas_pa_v2,
tot_post_panas_na_v2,
post_stai_1r_v2,
post_stai_2r_v2,
post_stai_5r_v2,
post_stai_8r_v2,
post_stai_10r_v2,
post_stai_11r_v2,
post_stai_15r_v2,
post_stai_16r_v2,
post_stai_19r_v2,
post_stai_20r_v2,
tot_post_stai_v2,
post_sroe_1_v2,
post_sroe_2_v2,
post_sroe_3_v2,
timetorest_A,
timetorest_B,
gAAT_RT_spouse_pull_A,
gAAT_RT_spouse_pull_B,
gAAT_RT_spouse_push_A,
gAAT_RT_spouse_push_B,
gAAT_RT_stranger_pull_A,
gAAT_RT_stranger_pull_B,
gAAT_RT_stranger_push_A,
gAAT_RT_stranger_push_B,
gAAT_RT_whoto_pull_A,
gAAT_RT_whoto_pull_B,
gAAT_RT_whoto_push_A,
gAAT_RT_whoto_push_B,
gAAT_RT_neutral_pull_A,
gAAT_RT_neutral_pull_B,
gAAT_RT_neutral_push_A,
gAAT_RT_neutral_push_B,
gAAT_RT_death_pull_A,
gAAT_RT_death_pull_B,
gAAT_RT_death_push_A,
gAAT_RT_death_push_B,
gAAT_push_v1,
gAAT_push_v2,
tot_pre_panas_pa_A,
tot_pre_panas_na_A,
tot_pre_stai_A,
tot_post_panas_pa_A,
tot_post_panas_na_A,
tot_post_stai_A,
post_sroe_1_A,
post_sroe_2_A,
post_sroe_3_A,
cortisol_BL_A,
cortisol_PR_A,
cortisol_PO_A,
tot_pre_panas_pa_B,
tot_pre_panas_na_B,
tot_pre_stai_B,
tot_post_panas_pa_B,
tot_post_panas_na_B,
tot_post_stai_B,
post_sroe_1_B,
post_sroe_2_B,
post_sroe_3_B,
cortisol_BL_B,
cortisol_PR_B,
cortisol_PO_B))
# drop D142 & D147, who were not included in the resting state fMRI data analysis due to poor data quality
data <- filter(data, !grepl("D142|D147",ID))
# join the data from the master dataset with the long dataset w/FNC and state variables
data <- left_join(data, long_state_fnc_sroe, by=c("ID"))
min(data$tx) # check that levels were set correctly (console should display Levels: placebo < oxytocin)
# make variables for retired/non-retired; post-grad ed. vs. no post-grad ed.
data <- data %>% mutate(employment_retired = as.factor(ifelse(employment == "retired", 1, 0)),
education_postgrad = recode_factor(education, "less than high school" = "less than postgrad", "high school grad" = "less than postgrad", "some college" = "less than postgrad", "college grad" = "less than postgrad", "1 year grad school" = "postgrad", "2 years grad school" = "postgrad", "3 years grad school" = "postgrad", "4+ years grad school" = "postgrad", .ordered = TRUE),
education_collegegrad = recode_factor(education, "less than high school" = "< college grad", "high school grad" = "< college grad", "some college" = "< college grad", "college grad" = "college grad", "1 year grad school" = "college grad", "2 years grad school" = "college grad", "3 years grad school" = "college grad", "4+ years grad school" = "college grad", .ordered = TRUE))
levels(data$education_collegegrad)
# make a "time since death in months" variable
data <- data %>% mutate(timesincedeath_mon = timesincedeath/30.417)
# mean-center the FNCs (any variable starting with "fnc")
# and add the suffix "_cm"
meancent <- function(x) { x - mean(x, na.rm=TRUE) } #simple worker function to mean center a variable
data <- data %>% mutate_at(vars(starts_with("fnc")), funs(cm=meancent))
# also mean-center any continuous variables that will be used in interactions
data <- data %>% mutate_at(vars("tot_bdi", "tot_icg", "age_yrs", "fd_mean", "yrs_together", "timesincedeath_mon"), funs(cm=meancent))
# save as an .rds file and a csv file
saveRDS(data, "/Users/sarenseeley/Dropbox/Dissertation/Non-fMRI data analysis/master_with_static_FNC.rds")
write_csv(data, "/Users/sarenseeley/Dropbox/Dissertation/Non-fMRI data analysis/master_with_static_FNC.csv")
# rm(list=ls()) clears the global environment
# split up by session
data_pl <- data %>% filter(tx == "placebo")
data_ot <- data %>% filter(tx == "oxytocin")
Figure 2. Hypothesized major network interactions in complicated grief, showing the flow and content of spontaneous thought as potentially more automatically constrained and less variable due to (a) lesser influence from executive control regions over default mode influence26, (b) greater influence from default mode core subsystem regions48, and (c) greater influence of the salience network – so thoughts are intrusive, inflexible, and fixated on the loss.
Based on my model, the connections that I should care most about are DMNcore-DMNmtl, DMNcore-SN, SN-dorsal attention, SN-FPN, DMNcore-FPN, and FPN-dorsal attention.
There was no dorsal attention network component in the gICA results, leaving (A) DMNcore-DMNmtl, (B) DMNcore-SN, (C) SN-FPN, and (D) DMNcore-FPN as the connections of interest that I can look at in this dataset.
gICA results yielded 4 DMN components, with peaks in the PCC (27), precuneus (13), mOFC/ventral striatum (21), and retrosplenial cortex (10). 27 and 13 are most representative of DMNcore, and 10 is most representative of DMNmtl. 21 is interesting…peak could actually be in nucleus accumbens, and associated Neurosynth topic keywords are mostly stuff in the vein of self-referential, social, and value/reward.
For the final DMN components, I selected 27 as most clearly DMNcore and 10 as most clearly DMNmtl.
gICA results yielded one salience network component, primarily centered on the dACC (26). There is also another component that’s basically all insula (posterior & anterior) with the peak in the anterior insula. I classified this as SN as well. However, the fALFF is pretty low (2.65), whereas all the other solidly “good” components have fALFF >4.2 (the “fALFF > 6.0”" rule of thumb for artifact cutoff is different given that ICA-AROMA already removed a lot of high frequency noise).
gICA results yielded two FPN components that both look strongly related to the visuospatial network. I selected 6 (Right) vs. 27 (Left), because 27 has some eyeball movement in it and maps more strongly to the Stanford RSN template.
(look @ this article later: https://academic.oup.com/cercor/article/28/2/726/4637602#116373137)
Relevant to model:
Intra-network connectivity:
fnc_dmn27_dmn10)Inter-network connectivity:
fnc_dmn10_sn26)fnc_dmn27_sn26)fnc_fpn6_sn26)fnc_dmn27_fpn6)Aim 1. To test whether dynamic functional connectivity under placebo differs among widowed older adults who are adjusting well, versus those who are adjusting poorly (i.e., complicated grief).
Question 1. Is grief severity associated with functional connectivity between (a) DMNcore-FPN, (b) DMNmtl-DMNcore, or (c) DMNcore-SN?
Static FNC Aim 1. To test whether static functional connectivity between model-relevant IC pairs is predicted by/predicts (???) complicated grief symptom severity.
- H1: Complicated grief symptom severity will predict lower variability in spontaneous thought flow over time, as indexed by fewer functional connectivity state transitions over the resting state scan.
Aim 1, H1: Not relevant to static FNC.
- H2: Complicated grief symptom severity will predict greater automatic constraints on thought content, as indexed by more dwell time in states of default mode-salience network interconnectivity.
Static Aim 1, H2: Complicated grief symptom severity will predict/be predicted by stronger correlation between DMN and SN IC pair(s).
- H3: Ruminative content in self-reported pre-resting state thought content will mediate the relationship between grief symptom severity and automatic constraints as indexed by dwell time in states of default mode-salience network interconnectivity.
Question 2. Does spouse-related thought content mediate the relationship between grief severity complicated grief symptom severity and DMN-SN FNC?
Static Aim 1, H2: Higher proportion of spouse-related thought content will mediate the relationship between complicated grief symptom severity and DMN-SN FNC.
# First, test within-default network FNC
a1h2a <- lm(tot_icg ~ fnc_dmn27_dmn10, data=data_pl)
summary(a1h2a)
Call:
lm(formula = tot_icg ~ fnc_dmn27_dmn10, data = data_pl)
Residuals:
Min 1Q Median 3Q Max
-18.785 -9.235 -1.199 8.367 24.430
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.611 2.199 9.828 9.85e-12 ***
fnc_dmn27_dmn10 8.985 6.749 1.331 0.191
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12.41 on 36 degrees of freedom
Multiple R-squared: 0.04693, Adjusted R-squared: 0.02045
F-statistic: 1.773 on 1 and 36 DF, p-value: 0.1914
summ(a1h2a, confint=TRUE) # use jtools `summ` function to get confidence intervals
[4mMODEL INFO:[24m
[3mObservations:[23m 38
[3mDependent Variable:[23m tot_icg
[3mType:[23m OLS linear regression
[4mMODEL FIT:[24m
[3mF[23m(1,36) = 1.77, [3mp[23m = 0.19
[3mR² = [23m0.05
[3mAdj. R² = [23m0.02
[3mStandard errors: OLS[23m
-------------------------------------------------------------
Est. 2.5% 97.5% t val. p
--------------------- ------- ------- ------- -------- ------
(Intercept) 21.61 17.15 26.07 9.83 0.00
fnc_dmn27_dmn10 8.99 -4.70 22.67 1.33 0.19
-------------------------------------------------------------
fnc_dmn27_dmn10 did not predict tot_icg (95% CI’s -4.70, 33.67). The overall model fit was not significant, F(1,36) = 1.77, p = .19.
# then, test between-network FNC
a1h2b <- lm(tot_icg ~ fnc_dmn27_sn26 + fnc_dmn10_sn26 + fnc_dmn27_fpn6 + fnc_fpn6_sn26, data=data_pl)
summary(a1h2b)
Call:
lm(formula = tot_icg ~ fnc_dmn27_sn26 + fnc_dmn10_sn26 + fnc_dmn27_fpn6 +
fnc_fpn6_sn26, data = data_pl)
Residuals:
Min 1Q Median 3Q Max
-20.818 -7.106 -1.099 9.211 21.129
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 29.897 4.381 6.824 8.73e-08 ***
fnc_dmn27_sn26 17.018 6.767 2.515 0.017 *
fnc_dmn10_sn26 7.093 8.076 0.878 0.386
fnc_dmn27_fpn6 -7.724 6.773 -1.140 0.262
fnc_fpn6_sn26 -5.483 7.900 -0.694 0.493
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11.88 on 33 degrees of freedom
Multiple R-squared: 0.1997, Adjusted R-squared: 0.1027
F-statistic: 2.058 on 4 and 33 DF, p-value: 0.1089
summ(a1h2b, confint=TRUE) # use jtools `summ` function to get confidence intervals
[4mMODEL INFO:[24m
[3mObservations:[23m 38
[3mDependent Variable:[23m tot_icg
[3mType:[23m OLS linear regression
[4mMODEL FIT:[24m
[3mF[23m(4,33) = 2.06, [3mp[23m = 0.11
[3mR² = [23m0.20
[3mAdj. R² = [23m0.10
[3mStandard errors: OLS[23m
-------------------------------------------------------------
Est. 2.5% 97.5% t val. p
-------------------- ------- -------- ------- -------- ------
(Intercept) 29.90 20.98 38.81 6.82 0.00
fnc_dmn27_sn26 17.02 3.25 30.79 2.51 0.02
fnc_dmn10_sn26 7.09 -9.34 23.52 0.88 0.39
fnc_dmn27_fpn6 -7.72 -21.50 6.06 -1.14 0.26
fnc_fpn6_sn26 -5.48 -21.56 10.59 -0.69 0.49
-------------------------------------------------------------
In the model testing all between-network connections, fnc_dmn27_sn26 was the only component pair that was a significant predictor of grief severity, b = 17.02, t = 2.51, p = .017, 95% CI’s 3.25, 30.79). The overall model fit was not significant, F(4,33) = 2.06, p = .11.
A basic scatterplot:
plot_model(a1h2b, type = "pred", terms = c("fnc_dmn27_sn26"), show.data=TRUE) + xlab("FNC between salience and default network core") + ylab("Grief severity (ICG)") + labs(title="Predicted ICG scores") + theme_apa()
plot(data_pl$fnc_dmn27_sn26, data_pl$tot_icg)
Next, I tested fnc_dmn27_sn26 in a model that included age, sex, and time since death.
a1h2c <- lm(tot_icg ~ fnc_dmn27_sn26*(age_yrs_cm + sex_m + timesincedeath_mon_cm), data=data_pl) # look at two-way interactions only
summary(a1h2c)
Call:
lm(formula = tot_icg ~ fnc_dmn27_sn26 * (age_yrs_cm + sex_m +
timesincedeath_mon_cm), data = data_pl)
Residuals:
Min 1Q Median 3Q Max
-27.174 -5.283 1.129 5.979 15.334
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.1556 2.6133 11.922 6.57e-13 ***
fnc_dmn27_sn26 21.7524 6.2272 3.493 0.001504 **
age_yrs_cm -0.0799 0.4002 -0.200 0.843092
sex_m.L 13.1018 3.5813 3.658 0.000967 ***
timesincedeath_mon_cm 0.1561 0.3100 0.503 0.618349
fnc_dmn27_sn26:age_yrs_cm -0.5569 0.8257 -0.674 0.505194
fnc_dmn27_sn26:sex_m.L 22.5152 8.2273 2.737 0.010329 *
fnc_dmn27_sn26:timesincedeath_mon_cm 0.7677 0.7562 1.015 0.318147
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.3 on 30 degrees of freedom
Multiple R-squared: 0.4528, Adjusted R-squared: 0.3251
F-statistic: 3.546 on 7 and 30 DF, p-value: 0.006798
summ(a1h2c, confint=TRUE)
[4mMODEL INFO:[24m
[3mObservations:[23m 38
[3mDependent Variable:[23m tot_icg
[3mType:[23m OLS linear regression
[4mMODEL FIT:[24m
[3mF[23m(7,30) = 3.55, [3mp[23m = 0.01
[3mR² = [23m0.45
[3mAdj. R² = [23m0.33
[3mStandard errors: OLS[23m
----------------------------------------------------------------------------------
Est. 2.5% 97.5% t val. p
------------------------------------------ ------- ------- ------- -------- ------
(Intercept) 31.16 25.82 36.49 11.92 0.00
fnc_dmn27_sn26 21.75 9.03 34.47 3.49 0.00
age_yrs_cm -0.08 -0.90 0.74 -0.20 0.84
sex_m.L 13.10 5.79 20.42 3.66 0.00
timesincedeath_mon_cm 0.16 -0.48 0.79 0.50 0.62
fnc_dmn27_sn26:age_yrs_cm -0.56 -2.24 1.13 -0.67 0.51
fnc_dmn27_sn26:sex_m.L 22.52 5.71 39.32 2.74 0.01
fnc_dmn27_sn26:timesincedeath_mon_cm 0.77 -0.78 2.31 1.02 0.32
----------------------------------------------------------------------------------
DMNcore-SN FNC remained a significant predictor of grief severity in interaction with sex, when age, sex, and time since death were included in the model. The overall model fit was significant, F(7,30) = 3.55, p = .009, adjusted R2 = .33.
The following plot shows the interaction:
plot_model(a1h2c, type = "pred", terms = c("fnc_dmn27_sn26", "sex_m"), show.data=TRUE) + xlab("FNC between salience and default network core") + ylab("Grief severity (ICG)") + labs(title="Predicted ICG scores")
There was a significant interaction between FNC and sex, such that DMN-SN FNC was positively associated with grief severity in men but not women. (Note that grief severity and sex are pretty confounded in this dataset, with higher mean ICG among the men than the women…)
Aim 2. To identify how intranasal oxytocin alters dynamic functional connectivity in older adults, and whether effects of oxytocin are moderated by grief severity.
H1: Under oxytocin, the sample as a whole will show increased time in salience network-dominant states, relative to placebo.
Paired t-tests will be used to compare oxytocin and placebo conditions to test whether oxytocin (A) increases dwell time in salience-network dominant states and/or (B) reduces n state transitions in the sample.
H2: Complicated grief symptom severity will moderate oxytocin effects on DFNC.
Julia said:
Response: I actually would recommend doing a paired t-test (equally, when no data are missing, a mixed model with oxytocin state as a fixed predictor and a random intercept per participant) to compare the differences in dwell times between oxytocin states. The reason is that when you also have an interaction with grief symptom severity in the model, comparing the two oxytocin states must be done at a particular level of grief symptom severity; the difference between these states will depend on the particular level of grief symptom severity chosen. Even if the interaction is not statistically significant, the model likely did not estimate the interaction as exactly zero. Thus, you’re having to utilize a non-zero estimate of a parameter that you actually think is zero in your calculation of the oxytocin differences. I don’t love that.
Of course, if the interaction is statistically significant, then you definitely should be taking it into account when you’re estimating the difference between oxytocin states. In this case, the grief symptom severity really does moderate that difference between states. Given all this, I would propose looking at the interaction model first. If the interaction is not statistically significant, look at the model with just oxytocin states (again, assuming no missing data, this would be the paired t-test).
Following Julia’s advice, test Q4/H2 (interaction model) first, then follow up with LME (instead of t-test, so things like sex can be added to the model) if there is no effect of ICG or interaction:
# fnc_dmn27_sn26 / a
# fnc_dmn10_sn26 / b
# fnc_fpn6_sn26 / c
# fnc_dmn10_dmn27 /d
# from MFO meeting 6/20:
# run it with just icg*tx (+/- other stuff), then see if BDI changes results (+ tot_bdi_cm)
# run 3 models: icg*tx alone, icg*tx + other stuff, icg*tx + other stuff + BDI
# use ML over REML so that models with different fixed effects can be compared
## Models with tot_icg*tx ONLY
# SN component pairs
a2h2a <- lme(fnc_dmn27_sn26 ~ tot_icg_cm*tx, random = ~ 1|ID, data, method="ML") # DMN core-SN dACC
a2h2b <- lme(fnc_dmn10_sn26 ~ tot_icg_cm*tx, random = ~ 1|ID, data, method="ML") # DMN mtl-SN dACC
a2h2c <- lme(fnc_fpn6_sn26 ~ tot_icg_cm*tx, random = ~ 1|ID, data, method="ML") # R FPN-SN dACC
# within-DMN pair
a2h2d <- lme(fnc_dmn27_dmn10 ~ tot_icg_cm*tx, random = ~ 1|ID, data, method="ML") # DMN core-SN dACC
## Models with other covariates
# SN component pairs
a2h2a1 <- lme(fnc_dmn27_sn26 ~ tot_icg_cm*tx + sex_m + age_yrs_cm + fd_mean + meds_hrt + timetorest + meds_psychoactive, random = ~ 1|ID, data, method="ML") # DMN core-SN dACC
a2h2b1 <- lme(fnc_dmn10_sn26 ~ tot_icg_cm*tx + sex_m + age_yrs_cm + fd_mean + meds_hrt + timetorest + meds_psychoactive, random = ~ 1|ID, data, method="ML") # DMN mtl-SN dACC
a2h2c1 <- lme(fnc_fpn6_sn26 ~ tot_icg_cm*tx + sex_m + age_yrs_cm + fd_mean + meds_hrt + timetorest + meds_psychoactive, random = ~ 1|ID, data, method="ML") # R FPN-SN dACC
# within-DMN pair
a2h2d1 <- lme(fnc_dmn27_dmn10 ~ tot_icg_cm*tx + sex_m + age_yrs_cm + fd_mean + meds_hrt + timetorest + meds_psychoactive, random = ~ 1|ID, data, method="ML") # DMN core-SN dACC
# Models w/BDI
# SN component pairs
a2h2a2 <- lme(fnc_dmn27_sn26 ~ tot_icg_cm*tx*tot_bdi_cm, random = ~ 1|ID, data, method="ML") # DMN core-SN dACC
a2h2b2 <- lme(fnc_dmn10_sn26 ~ tot_icg_cm*tx*tot_bdi_cm, random = ~ 1|ID, data, method="ML") # DMN mtl-SN dACC
a2h2c2 <- lme(fnc_fpn6_sn26 ~ tot_icg_cm*tx*tot_bdi_cm, random = ~ 1|ID, data, method="ML") # R FPN-SN dACC
# within-DMN pair
a2h2d2 <- lme(fnc_dmn27_dmn10 ~ tot_icg_cm*tx*tot_bdi_cm, random = ~ 1|ID, data, method="ML") # DMN core-SN dACC
Note that sex_m and tot_icg are confounded - the mean ICG score for men is ~31 while for women it’s ~20, p = .06.
anova(a2h2a) # DMN core-SN dACC
numDF denDF F-value p-value
(Intercept) 1 36 30.595264 <.0001
tot_icg_cm 1 36 4.531186 0.0402
tx 1 36 2.472479 0.1246
tot_icg_cm:tx 1 36 0.677543 0.4159
anova(a2h2a1) # w/covariates
numDF denDF F-value p-value
(Intercept) 1 34 29.165981 <.0001
tot_icg_cm 1 32 4.319508 0.0458
tx 1 34 2.381523 0.1320
sex_m 1 32 0.008121 0.9288
age_yrs_cm 1 32 0.031568 0.8601
fd_mean 1 34 0.354832 0.5553
meds_hrt 1 32 0.383940 0.5399
timetorest 1 34 2.714524 0.1087
meds_psychoactive 1 32 0.010532 0.9189
tot_icg_cm:tx 1 34 0.142984 0.7077
anova(a2h2a2) # w/BDI
numDF denDF F-value p-value
(Intercept) 1 34 31.545943 <.0001
tot_icg_cm 1 34 4.671982 0.0378
tx 1 34 2.517478 0.1218
tot_bdi_cm 1 34 0.117518 0.7339
tot_icg_cm:tx 1 34 0.689875 0.4120
tot_icg_cm:tot_bdi_cm 1 34 3.001101 0.0923
tx:tot_bdi_cm 1 34 2.269676 0.1412
tot_icg_cm:tx:tot_bdi_cm 1 34 0.385522 0.5388
anova(a2h2a, a2h2a1, a2h2a2) # model comparison
Model df AIC BIC logLik Test L.Ratio p-value
a2h2a 1 6 61.44989 75.43429 -24.72495
a2h2a1 2 12 70.07932 98.04812 -23.03966 1 vs 2 3.370575 0.7611
a2h2a2 3 10 63.25768 86.56502 -21.62884 2 vs 3 2.821635 0.2439
Grief severity, but not tx or their interaction, has a significant effect on DMN core-SN ACC FNC - similar to the result with just the placebo data. This effect holds when sex, age, mean FD and hormone replacement treatment are included in the model.
Model summary and confidence intervals for the estimate:
summary(a2h2a)
Linear mixed-effects model fit by maximum likelihood
Data: data
AIC BIC logLik
61.44989 75.43429 -24.72495
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 0.2260328 0.268737
Fixed effects: fnc_dmn27_sn26 ~ tot_icg_cm * tx
Value Std.Error DF t-value p-value
(Intercept) -0.27222976 0.04921625 36 -5.531299 0.0000
tot_icg_cm 0.00846803 0.00397811 36 2.128658 0.0402
tx.L 0.07042754 0.04478950 36 1.572412 0.1246
tot_icg_cm:tx.L -0.00297997 0.00362030 36 -0.823130 0.4159
Correlation:
(Intr) tt_cg_ tx.L
tot_icg_cm 0
tx.L 0 0
tot_icg_cm:tx.L 0 0 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.72884794 -0.64406077 0.02956914 0.64004024 1.98658252
Number of Observations: 76
Number of Groups: 38
round(intervals(a2h2a)[["fixed"]],4) # 95% CIs
lower est. upper
(Intercept) -0.3694 -0.2722 -0.1751
tot_icg_cm 0.0006 0.0085 0.0163
tx.L -0.0180 0.0704 0.1588
tot_icg_cm:tx.L -0.0101 -0.0030 0.0042
attr(,"label")
[1] "Fixed effects:"
# variance explained by the entire model
fitted_a2h2a <- fitted(a2h2a)
data$fitted_a2h2a <- as.vector(fitted_a2h2a)
forR2 <- lm(fnc_dmn27_sn26 ~ fitted_a2h2a, data=data)
summary(forR2)
Call:
lm(formula = fnc_dmn27_sn26 ~ fitted_a2h2a, data = data)
Residuals:
Min 1Q Median 3Q Max
-0.42141 -0.15045 0.02023 0.17715 0.47552
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1307 0.0385 3.394 0.00111 **
fitted_a2h2a 1.4801 0.1120 13.215 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.205 on 74 degrees of freedom
Multiple R-squared: 0.7024, Adjusted R-squared: 0.6984
F-statistic: 174.6 on 1 and 74 DF, p-value: < 2.2e-16
# ICC (variance explained by between-person variance)
# function from the `reghelper` package
round(ICC(a2h2a),2)
[1] 0.41
# panel by tx
ggplot(data, aes(fnc_dmn27_sn26, tot_icg, fill=tx)) + geom_point(aes(color = tx)) +
xlab("DMN core-SN connectivity") + ylab("Grief severity") + theme_bw(base_size=12) + scale_color_discrete() + facet_grid(. ~ tx) + geom_smooth(method=lm, aes(color = tx)) + theme_apa() + scale_fill_grey(start=.8, end=.8)
The effect is being driven by the association between grief severity and FNC in the placebo condition.
anova(a2h2b) # DMN mtl-SN dACC
numDF denDF F-value p-value
(Intercept) 1 36 12.086395 0.0013
tot_icg_cm 1 36 1.102735 0.3007
tx 1 36 7.015436 0.0119
tot_icg_cm:tx 1 36 0.086318 0.7706
anova(a2h2b1) # w/covariates
numDF denDF F-value p-value
(Intercept) 1 34 12.608381 0.0011
tot_icg_cm 1 32 1.150360 0.2915
tx 1 34 6.996639 0.0123
sex_m 1 32 0.051012 0.8227
age_yrs_cm 1 32 1.566661 0.2198
fd_mean 1 34 2.157875 0.1510
meds_hrt 1 32 1.621930 0.2120
timetorest 1 34 1.006385 0.3229
meds_psychoactive 1 32 1.123969 0.2970
tot_icg_cm:tx 1 34 0.016566 0.8983
anova(a2h2b2) # w/BDI
numDF denDF F-value p-value
(Intercept) 1 34 13.267921 0.0009
tot_icg_cm 1 34 1.210535 0.2790
tx 1 34 7.468969 0.0099
tot_bdi_cm 1 34 3.459973 0.0715
tot_icg_cm:tx 1 34 0.091898 0.7636
tot_icg_cm:tot_bdi_cm 1 34 2.059271 0.1604
tx:tot_bdi_cm 1 34 3.456061 0.0717
tot_icg_cm:tx:tot_bdi_cm 1 34 0.871266 0.3572
anova(a2h2b, a2h2b1, a2h2b2) # model comparison
Model df AIC BIC logLik Test L.Ratio p-value
a2h2b 1 6 -5.338165 8.646235 8.669082
a2h2b1 2 12 -1.455771 26.513029 12.727885 1 vs 2 8.117606 0.2296
a2h2b2 3 10 -7.606898 15.700435 13.803449 2 vs 3 2.151127 0.3411
# variance explained by the entire model
fitted_a2h2b <- fitted(a2h2b)
data$fitted_a2h2b <- as.vector(fitted_a2h2b)
forR2 <- lm(fnc_dmn10_sn26 ~ fitted_a2h2b, data=data)
summary(forR2)
Call:
lm(formula = fnc_dmn10_sn26 ~ fitted_a2h2b, data = data)
Residuals:
Min 1Q Median 3Q Max
-0.302106 -0.082537 0.002776 0.077691 0.277115
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.04444 0.01708 -2.602 0.0112 *
fitted_a2h2b 1.37446 0.08764 15.684 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1181 on 74 degrees of freedom
Multiple R-squared: 0.7687, Adjusted R-squared: 0.7656
F-statistic: 246 on 1 and 74 DF, p-value: < 2.2e-16
# ICC (variance explained by between-person variance)
# function from the `reghelper` package
round(ICC(a2h2b),2)
[1] 0.53
There was a significant main effect of tx (but not tot_icg or their interaction) on the DMN mtl-SN dACC pair, in which these components were more strongly correlated in the oxytocin condition. This effect holds when sex, age, mean FD and hormone replacement treatment are included in the model.
t.test(data$fnc_dmn10_sn26 ~ data$tx, paired=TRUE)
Paired t-test
data: data$fnc_dmn10_sn26 by data$tx
t = -2.682, df = 37, p-value = 0.01087
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.17632769 -0.02456044
sample estimates:
mean of the differences
-0.1004441
# create a bar plot with error bars (95% CI)
p <- ggplot(data, aes(tx, fnc_dmn10_sn26, fill=tx)) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_se, position=position_dodge(.9), width = 0.25, color="black") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) +
scale_fill_manual(values=c("#FFC107", "#1E88E5"))
p
m <- lme(fnc_dmn10_sn26 ~ tx, random = ~ 1|ID, data, method="ML") # use lme() instead of paired t-test to get effect estimate to plot (this gives same result as paired t-test)
ef <- effect("tx", m)
ef <- as.data.frame(ef)
ggplot(ef, aes(tx, fit, color=tx)) + geom_line() + geom_point() + geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.2) +
scale_color_discrete() + theme_bw(base_size=12) + xlab("Session") + ylab("Fitted values for DMN mtl (IC10) - SN (IC26)") + labs(color='Session')
anova(a2h2c) # R FPN-SN dACC
numDF denDF F-value p-value
(Intercept) 1 36 27.681369 <.0001
tot_icg_cm 1 36 4.967929 0.0322
tx 1 36 2.487758 0.1235
tot_icg_cm:tx 1 36 2.178898 0.1486
anova(a2h2c1) # w/covariates
numDF denDF F-value p-value
(Intercept) 1 34 25.952482 <.0001
tot_icg_cm 1 32 4.657648 0.0385
tx 1 34 2.607650 0.1156
sex_m 1 32 0.204197 0.6544
age_yrs_cm 1 32 0.002682 0.9590
fd_mean 1 34 2.998659 0.0924
meds_hrt 1 32 0.048902 0.8264
timetorest 1 34 1.146024 0.2919
meds_psychoactive 1 32 0.454598 0.5050
tot_icg_cm:tx 1 34 2.915342 0.0969
anova(a2h2c2) # w/BDI
numDF denDF F-value p-value
(Intercept) 1 34 26.811864 <.0001
tot_icg_cm 1 34 4.811880 0.0352
tx 1 34 2.465044 0.1257
tot_bdi_cm 1 34 0.812129 0.3738
tot_icg_cm:tx 1 34 2.159004 0.1509
tot_icg_cm:tot_bdi_cm 1 34 0.057067 0.8126
tx:tot_bdi_cm 1 34 0.918632 0.3446
tot_icg_cm:tx:tot_bdi_cm 1 34 0.752677 0.3917
anova(a2h2c, a2h2c1, a2h2c2) # model comparison
Model df AIC BIC logLik Test L.Ratio p-value
a2h2c 1 6 30.96778 44.95218 -9.483889
a2h2c1 2 12 37.01705 64.98585 -6.508527 1 vs 2 5.950723 0.4287
a2h2c2 3 10 36.18506 59.49239 -8.092529 2 vs 3 3.168004 0.2052
Grief severity, but not tx or their interaction, has a significant effect on R FPN-SN ACC FNC. This effect holds when sex, age, mean FD and hormone replacement treatment are included in the model.
Model summary and confidence intervals for the estimate:
summary(a2h2c)
Linear mixed-effects model fit by maximum likelihood
Data: data
AIC BIC logLik
30.96778 44.95218 -9.483889
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 0.1362308 0.2425831
Fixed effects: fnc_fpn6_sn26 ~ tot_icg_cm * tx
Value Std.Error DF t-value p-value
(Intercept) 0.19208002 0.03650803 36 5.261309 0.0000
tot_icg_cm 0.00657724 0.00295091 36 2.228885 0.0322
tx.L -0.06376955 0.04043051 36 -1.577263 0.1235
tot_icg_cm:tx.L 0.00482387 0.00326796 36 1.476109 0.1486
Correlation:
(Intr) tt_cg_ tx.L
tot_icg_cm 0
tx.L 0 0
tot_icg_cm:tx.L 0 0 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.3547201 -0.4245518 0.1132649 0.6106873 1.8400883
Number of Observations: 76
Number of Groups: 38
round(intervals(a2h2c)[["fixed"]],4) # 95% CIs
lower est. upper
(Intercept) 0.1200 0.1921 0.2641
tot_icg_cm 0.0008 0.0066 0.0124
tx.L -0.1436 -0.0638 0.0160
tot_icg_cm:tx.L -0.0016 0.0048 0.0113
attr(,"label")
[1] "Fixed effects:"
# panel by tx
ggplot(data, aes(fnc_fpn6_sn26, tot_icg, fill=tx)) + geom_point(aes(color = tx)) +
xlab("FPN-SN connectivity") + ylab("Grief severity") + theme_bw(base_size=12) + scale_color_discrete() + facet_grid(. ~ tx) + geom_smooth(method=lm, aes(color = tx)) + theme_apa() + scale_fill_grey(start=.8, end=.8)
# variance explained by the entire model
fitted_a2h2c <- fitted(a2h2c)
data$fitted_a2h2c <- as.vector(fitted_a2h2c)
forR2 <- lm(fnc_fpn6_sn26 ~ fitted_a2h2c, data=data)
summary(forR2)
Call:
lm(formula = fnc_fpn6_sn26 ~ fitted_a2h2c, data = data)
Residuals:
Min 1Q Median 3Q Max
-0.50451 -0.10227 0.01775 0.12225 0.43372
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.12411 0.04096 -3.030 0.00337 **
fitted_a2h2c 1.64611 0.17543 9.383 3.15e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.203 on 74 degrees of freedom
Multiple R-squared: 0.5433, Adjusted R-squared: 0.5372
F-statistic: 88.05 on 1 and 74 DF, p-value: 3.15e-14
# ICC (variance explained by between-person variance)
# function from the `reghelper` package
round(ICC(a2h2c),2)
[1] 0.24
The effect is being driven by the association between grief severity and FNC in the oxytocin condition.
# to visualize effects (sig. and ns) of oxytocin, create a long dataset with one column for FNC (all kinds), one column for IC pair name, one column for `tx`, one column for `ID`
# first, subset columns to select only FNC pairs of interest
fnc_bar_data <- data %>% select(contains("sn26"), contains("dmn10"), contains("dmn27"), -contains("_cm"), -contains("dmn13"), -contains("dmn10_fpn6"), -contains("dmn21"), -contains("sn12"), -contains("fpn17"), "ID", "tx", "tot_icg", "group", "sex_m")
fnc_bar_data <- gather(fnc_bar_data, fnc,value,-tot_icg, -tx, -ID, -group, -sex_m) %>% mutate(tx = factor(ifelse(tx == "placebo", "placebo", "oxytocin")))
# show as a bar chart w/error bars (mean SE):
ggplot(fnc_bar_data, aes(tx, value, fill=tx)) + stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_se, position=position_dodge(.9), width = 0.25, color="black") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.25) +
ylab("Correlation") + facet_grid(. ~ fnc) + scale_fill_discrete() + theme_apa(facet.title.size=7.75)
# panel by group
ggplot(fnc_bar_data, aes(group, value, fill=tx)) + stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_se, position=position_dodge(.9), width = 0.25, color="black") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) +
ylab("Correlation") + facet_grid(. ~ fnc) + theme_apa(facet.title.size=7.75)
anova(a2h2d)
numDF denDF F-value p-value
(Intercept) 1 36 14.128882 0.0006
tot_icg_cm 1 36 1.838644 0.1836
tx 1 36 0.155134 0.6960
tot_icg_cm:tx 1 36 0.189230 0.6662
anova(a2h2d1)
numDF denDF F-value p-value
(Intercept) 1 34 14.724719 0.0005
tot_icg_cm 1 32 1.916183 0.1759
tx 1 34 0.140201 0.7104
sex_m 1 32 1.073780 0.3079
age_yrs_cm 1 32 2.236282 0.1446
fd_mean 1 34 0.194583 0.6619
meds_hrt 1 32 0.178177 0.6758
timetorest 1 34 0.076769 0.7834
meds_psychoactive 1 32 0.340334 0.5637
tot_icg_cm:tx 1 34 0.123834 0.7271
anova(a2h2d2)
numDF denDF F-value p-value
(Intercept) 1 34 13.563076 0.0008
tot_icg_cm 1 34 1.765014 0.1928
tx 1 34 0.159221 0.6924
tot_bdi_cm 1 34 0.015113 0.9029
tot_icg_cm:tx 1 34 0.194215 0.6622
tot_icg_cm:tot_bdi_cm 1 34 0.543228 0.4662
tx:tot_bdi_cm 1 34 2.935086 0.0958
tot_icg_cm:tx:tot_bdi_cm 1 34 0.013287 0.9089
anova(a2h2d, a2h2d1, a2h2d2)
Model df AIC BIC logLik Test L.Ratio p-value
a2h2d 1 6 37.73946 51.72386 -12.86973
a2h2d1 2 12 45.40314 73.37194 -10.70157 1 vs 2 4.336315 0.6313
a2h2d2 3 10 41.96037 65.26771 -10.98019 2 vs 3 0.557232 0.7568
There was no effect of tot_icg or oxytocin on within-DMN connectivity.
dependent.vars <- cbind(data$fnc_dmn10_sn26, data$fnc_dmn27_sn26, data$fnc_fpn6_sn26)
summary(aov(dependent.vars ~ data$tot_icg * data$tx))
Response 1 :
Df Sum Sq Mean Sq F value Pr(>F)
data$tot_icg 1 0.0977 0.097675 1.6855 0.19834
data$tx 1 0.1917 0.191691 3.3079 0.07311 .
data$tot_icg:data$tx 1 0.0024 0.002359 0.0407 0.84069
Residuals 72 4.1724 0.057950
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response 2 :
Df Sum Sq Mean Sq F value Pr(>F)
data$tot_icg 1 0.8341 0.83415 6.4086 0.01354 *
data$tx 1 0.1885 0.18848 1.4481 0.23278
data$tot_icg:data$tx 1 0.0517 0.05165 0.3968 0.53073
Residuals 72 9.3716 0.13016
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response 3 :
Df Sum Sq Mean Sq F value Pr(>F)
data$tot_icg 1 0.5032 0.50323 6.1590 0.01541 *
data$tx 1 0.1545 0.15453 1.8913 0.17332
data$tot_icg:data$tx 1 0.1353 0.13534 1.6565 0.20220
Residuals 72 5.8828 0.08171
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(aov(dependent.vars ~ data$tot_icg * data$tx * data$tot_bdi))
Response 1 :
Df Sum Sq Mean Sq F value Pr(>F)
data$tot_icg 1 0.0977 0.097675 1.8368 0.17981
data$tx 1 0.1917 0.191691 3.6048 0.06186 .
data$tot_bdi 1 0.2792 0.279176 5.2500 0.02505 *
data$tot_icg:data$tx 1 0.0024 0.002359 0.0444 0.83383
data$tot_icg:data$tot_bdi 1 0.1662 0.166157 3.1247 0.08160 .
data$tx:data$tot_bdi 1 0.0887 0.088700 1.6680 0.20089
data$tot_icg:data$tx:data$tot_bdi 1 0.0224 0.022361 0.4205 0.51887
Residuals 68 3.6160 0.053176
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response 2 :
Df Sum Sq Mean Sq F value Pr(>F)
data$tot_icg 1 0.8341 0.83415 6.5833 0.01250 *
data$tx 1 0.1885 0.18848 1.4876 0.22681
data$tot_bdi 1 0.0210 0.02098 0.1656 0.68533
data$tot_icg:data$tx 1 0.0517 0.05165 0.4076 0.52531
data$tot_icg:data$tot_bdi 1 0.5358 0.53582 4.2289 0.04358 *
data$tx:data$tot_bdi 1 0.1699 0.16993 1.3411 0.25089
data$tot_icg:data$tx:data$tot_bdi 1 0.0289 0.02886 0.2278 0.63469
Residuals 68 8.6160 0.12671
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response 3 :
Df Sum Sq Mean Sq F value Pr(>F)
data$tot_icg 1 0.5032 0.50323 6.0170 0.01674 *
data$tx 1 0.1545 0.15453 1.8477 0.17855
data$tot_bdi 1 0.0849 0.08493 1.0155 0.31715
data$tot_icg:data$tx 1 0.1353 0.13534 1.6183 0.20766
data$tot_icg:data$tot_bdi 1 0.0060 0.00597 0.0714 0.79018
data$tx:data$tot_bdi 1 0.0576 0.05759 0.6886 0.40956
data$tot_icg:data$tx:data$tot_bdi 1 0.0472 0.04718 0.5642 0.45518
Residuals 68 5.6871 0.08363
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
YES, functional connectivity between DMNcore and SNACC (fnc_dmn27_sn26) was predicted of grief severity (tot_icg). There was a significant interaction with sex, such that DMN-SN FNC was positively associated with grief severity in men but not women (sex also predictive of grief severity, though…)
plot_model(a1h2c, type = "pred", terms = c("fnc_dmn27_sn26", "sex_m"), show.data=TRUE) + xlab("FNC between salience and default network core") + ylab("Grief severity (ICG)") + labs(title="Predicted ICG scores")
# this is the model that includes main and interaction effects of age, sex, and mean FD
YES, oxytocin increased connectivity between DMNMTL-SNdACC across the sample (no effect of grief severity). This effect held when age, sex, mean FD, and HRT were included as covariates in the model.
NO, not in any of the component pairs I looked at.
In this aim, I wanted to test whether static FNC would predict change in grief severity from study to 20-24 months post-study.
t.test(data_pl$tot_icg, data_pl$tot_icg_20mFU,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
Paired t-test
data: data_pl$tot_icg and data_pl$tot_icg_20mFU
t = 4.4218, df = 30, p-value = 0.0001185
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
2.586522 7.026381
sample estimates:
mean of the differences
4.806452
describe(data_pl$tot_icg)
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 38 22.79 12.54 21.5 22.06 10.38 4 51 47 0.48 -0.58 2.03
describe(data_pl$tot_icg_20mFU)
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 31 15.45 10.55 13 14.44 13.34 2 45 43 0.75 0.05 1.89
# yes, ICG on average is significantly decreased at the follow-up
# compute ICG change and time elapsed between death and follow-up
data_pl$tot_icg_change <- data_pl$tot_icg_20mFU - data_pl$tot_icg
describe(data_pl$tot_icg_change)
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 31 -4.81 6.05 -4 -4.68 4.45 -19 7 26 -0.31 -0.44 1.09
hist(data_pl$tot_icg_change)
data_pl$time_death_20mFU_mos <- data_pl$timesincedeath_mon + data_pl$timeto20mFU_mos
describe(data_pl$time_death_20mFU_mos)
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 31 36.43 7.24 35.67 35.39 6.39 26.63 55.56 28.93 1.03 0.41 1.3
hist(data_pl$time_death_20mFU_mos)
# remove the cases with no follow-up data (7 out of the 38)
data_pl_fu <- data_pl %>% filter(!is.na(tot_icg_change))
# Predicting ICG change from DMN mtl-SN dACC FNC, controlling for time to follow-up
a3h3a <- lm(tot_icg_change ~ fnc_dmn10_sn26 + timeto20mFU, data=data_pl_fu)
summary(a3h3a)
Call:
lm(formula = tot_icg_change ~ fnc_dmn10_sn26 + timeto20mFU, data = data_pl_fu)
Residuals:
Min 1Q Median 3Q Max
-10.3949 -3.6184 0.1642 2.4244 11.2999
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 86.34634 29.73199 2.904 0.00711 **
fnc_dmn10_sn26 -8.61529 3.63124 -2.373 0.02477 *
timeto20mFU -0.13692 0.04487 -3.051 0.00495 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.166 on 28 degrees of freedom
Multiple R-squared: 0.3199, Adjusted R-squared: 0.2714
F-statistic: 6.586 on 2 and 28 DF, p-value: 0.004526
summ(a3h3a, confint=TRUE) # use jtools `summ` function to get confidence intervals
[4mMODEL INFO:[24m
[3mObservations:[23m 31
[3mDependent Variable:[23m tot_icg_change
[3mType:[23m OLS linear regression
[4mMODEL FIT:[24m
[3mF[23m(2,28) = 6.59, [3mp[23m = 0.00
[3mR² = [23m0.32
[3mAdj. R² = [23m0.27
[3mStandard errors: OLS[23m
--------------------------------------------------------------
Est. 2.5% 97.5% t val. p
-------------------- ------- -------- -------- -------- ------
(Intercept) 86.35 25.44 147.25 2.90 0.01
fnc_dmn10_sn26 -8.62 -16.05 -1.18 -2.37 0.02
timeto20mFU -0.14 -0.23 -0.05 -3.05 0.00
--------------------------------------------------------------
# model diagnostics/assumption-checking:
plot_model(a3h3a, type = "diag")
[[1]]
[[2]]
[[3]]
[[4]]
plot(data_pl_fu$tot_icg_change ~ data_pl_fu$fnc_dmn10_sn26)
plot_model(a3h3a, type = "pred", terms = c("fnc_dmn10_sn26"), show.data=TRUE) + xlab("DMN mtl-SN dACC connectivity") + ylab("ICG change") + labs(title="Predicted ICG change scores")
DMN mtl-SN dACC FNC does continue to predict ICG scores at follow-up when time to follow up is included in the model.
icg <- subset(data_pl_fu, select=c(ID, tot_icg, timeto20mFU, tot_icg_change, group)) %>% mutate(time = rep("T1: study"))
icg_fu <- subset(data_pl_fu, select=c(ID, tot_icg_20mFU,tot_icg_change, group)) %>% rename(tot_icg = tot_icg_20mFU) %>% mutate(time = rep("T2: follow-up"))
change <- bind_rows(icg, icg_fu) %>% mutate(ID=as.factor(ID))
# Plot trajectories
ggplot(data = change, aes(x = time, y = tot_icg, group = ID, color=group)) +
geom_hline(yintercept=25, linetype="dashed", color="light blue 4", size=.5) +
geom_hline(yintercept=22.79, linetype="dotted", color="dark grey", size=.5) +
geom_hline(yintercept=15.45, linetype="dotted", color="dark grey", size=.5) +
geom_point() + geom_line() +
annotate("text", label = "Clinical cutoff", x = .55, y = 27, size = 3, colour = "light blue 4", fontface = "italic") +
annotate("text", label = "Study M = 22.78", x = .74, y = 21, size = 3, colour = "dark grey") +
annotate("text", label = "Follow-up M = 15.45", x = .74, y = 14, size = 3, colour = "dark grey")
mean(data_pl$tot_icg)
[1] 22.78947
mean(data_pl_fu$tot_icg)
[1] 20.25806
mean(data_pl_fu$tot_icg_20mFU)
[1] 15.45161
# read in data
data <- readRDS("/Users/sarenseeley/Dropbox/Dissertation/Non-fMRI data analysis/master_with_static_FNC.rds")
# make a placebo session-only dataset to get rid of duplicate observations per person,
# which change the descriptive stats
data_pl <- data %>% filter(tx == "placebo")
# also make an oxytocin session dataset while we're at it
data_ot <- data %>% filter(tx == "oxytocin")
count(data_pl,group) # n in each group
count(data_pl_fu,group_20mFU) # n in each group at follow-up
# t-tests
## age
t.test(data_pl$age_yrs ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$age_yrs, group=data_pl$group)
## time since death
t.test(data_pl$timesincedeath ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$timesincedeath, group=data_pl$group)
## relationship length
t.test(data_pl$yrs_together ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$yrs_together, group=data_pl$group)
## total n prescription medications
t.test(data_pl$meds_total ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$meds_total, group=data_pl$group)
## BDI
t.test(data_pl$tot_icg ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_icg, group=data_pl$group)
## ICG
t.test(data_pl$tot_bdi ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_bdi, group=data_pl$group)
## ICG at follow-up
t.test(data_pl$tot_icg_20mFU ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_icg, group=data_pl$group)
describeBy(data_pl$tot_icg_20mFU, group=data_pl$group)
bi.bars(data_pl, "tot_icg_20mFU", "group")
No significant difference in age, length of relationship, time since death, or total number of prescription meds. Unsurprisingly, the CG group has higher mean BDI and ICG scores.
# chi-square tests
## sex
chisq.test(data_pl$group, y = data_pl$sex_m, correct = FALSE, simulate.p.value = TRUE)
# use the simulated p value because if you don't, with small cell sizes, many of the expected values will be very small
# and therefore the approximations of p may not be correct (R will warn you about this)
data_pl %>% group_by(group) %>%
count(sex_m) %>%
mutate(`(\\%)` = prop.table(n)*100)
## race
chisq.test(data_pl$group, y = data_pl$race, correct = FALSE, simulate.p.value = TRUE)
data_pl %>% group_by(group) %>%
count(race) %>%
mutate(`(\\%)` = prop.table(n)*100)
## ethnicity
chisq.test(data_pl$group, y = data_pl$ethnicity_hisp, correct = FALSE, simulate.p.value = TRUE)
data_pl %>% group_by(group) %>%
count(ethnicity_hisp) %>%
mutate(`(\\%)` = prop.table(n)*100)
## education (college grad vs. less than college degree)
chisq.test(data_pl$group, y = data_pl$education_collegegrad, correct = FALSE)
data_pl %>% group_by(group) %>%
count(education_collegegrad) %>%
mutate(`(\\%)` = prop.table(n)*100)
## employment (retired vs. still working full/part-time or looking for work)
chisq.test(data_pl$group, y = data_pl$employment_retired, correct = FALSE, simulate.p.value = TRUE)
data_pl %>% group_by(group) %>%
count(employment_retired) %>%
mutate(`(\\%)` = prop.table(n)*100)
## psychoactive meds (yes/no)
chisq.test(data_pl$group, y = data_pl$meds_psychoactive, correct = FALSE, simulate.p.value = TRUE)
data_pl %>% group_by(group) %>%
count(meds_psychoactive) %>%
mutate(`(\\%)` = prop.table(n)*100)
## hormone or opiate meds (yes/no)
chisq.test(data_pl$group, y = data_pl$meds_hormone_opiate, correct = FALSE)
data_pl %>% group_by(group) %>%
count(meds_hormone_opiate) %>%
mutate(`(\\%)` = prop.table(n)*100)
## order of sessions (B first vs. A first)
chisq.test(data_pl$group, y = data_pl$tx_v1, correct = FALSE)
data_pl %>% group_by(group) %>%
count(tx_v1) %>%
mutate(`(\\%)` = prop.table(n)*100)
No significant group differences in the frequencies, though the CG group does have a higher proportion of men and the NCG has a higher proportion of people taking hormones (not HRT - this category includes things like levothyroxine, insulin, etc.) or opiate medications (n=1). I suspect this is confounded with the imbalance in sex between the two groups (i.e., more women being treated for low thyroid).
bi.bars(data_pl, "meds_hormone_opiate", "group", main="Non-HRT hormone & opiate meds by group (NCG, CG)", zero=FALSE)
# on the y-axis, 1 = No, 2 = Yes, taking
table(data_pl$meds_hormone_opiate, data_pl$sex_m, data_pl$group)
# t-tests
### PANAS
## pre PANAS NA (placebo)
t.test(data_pl$tot_pre_panas_na ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_pre_panas_na, group=data_pl$group)
## pre PANAS NA (oxytocin)
t.test(data_ot$tot_pre_panas_na ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_ot$tot_pre_panas_na, group=data_ot$group)
## pre PANAS PA (placebo)
t.test(data_pl$tot_pre_panas_pa ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_pre_panas_pa, group=data_pl$group)
## pre PANAS PA (oxytocin)
t.test(data_ot$tot_pre_panas_pa ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_ot$tot_pre_panas_pa, group=data_ot$group)
### STAI-S
## pre STAI (placebo)
t.test(data_pl$tot_pre_stai ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_pre_stai, group=data_pl$group)
## pre STAI (oxytocin)
t.test(data_ot$tot_pre_stai ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_ot$tot_pre_stai, group=data_ot$group)
## post STAI (placebo)
t.test(data_pl$tot_post_stai ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$tot_post_stai, group=data_pl$group)
## post STAI (oxytocin)
t.test(data_ot$tot_post_stai ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_ot$tot_post_stai, group=data_ot$group)
### Framewise displacement
## mean framewise displacement (placebo session)
t.test(data_pl$fd_mean ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$fd_mean, group=data_pl$group)
## mean framewise displacement (oxytocin session)
t.test(data_ot$fd_mean ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_ot$fd_mean, group=data_ot$group)
### Time to rest
## mean time to resting state start (placebo session)
t.test(data_pl$timetorest ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_pl$timetorest, group=data_pl$group)
## mean time to resting state start (oxytocin session)
t.test(data_ot$timetorest ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
describeBy(data_ot$timetorest, group=data_ot$group)
### Thought content
# spouse related (%), non-spouse related (%)
t.test(data_pl$focus_spouse_perc ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
t.test(data_pl$focus_nonspouse_perc ~ data_pl$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
t.test(data_ot$focus_spouse_perc ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
t.test(data_ot$focus_nonspouse_perc ~ data_ot$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$focus_spouse_perc ~ data$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$focus_nonspouse_perc ~ data$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
At the placebo session, the groups differ on tot_pre_stai & tot_pre_panas_na (CG > NCG).
At the oxytocin session, the groups differ on tot_pre_stai & tot_pre_panas_na (CG > NCG).
Mean framewise displacement and time elapsed from spray admin. to resting state do not significantly differ at either session. However, the NCG group overall has more movement than the CG group during the oxytocin session (p = .077), and it may make sense to control for this.
Percent of spouse-related or non-spouse-related thought content doesn’t differ by group at either session, but when collapsing across session, the CG group tends to report more spouse-related content, p = .08.
# by treatment
describe(data$timetorest)
t.test(data$timetorest ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$timetorest ~ data$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$fd_mean ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$tot_pre_panas_na ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$tot_pre_panas_pa ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$tot_pre_stai ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$tot_post_panas_na ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$tot_post_panas_pa ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data$tot_post_stai ~ data$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
# D115 is missing post-scan SROE at one session, so test the SROE variables after removing their case
data_no115 <- data %>% filter(ID != "D115")
t.test(data_no115$focus_spouse_perc ~ data_no115$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
t.test(data_no115$focus_nonspouse_perc ~ data_no115$tx,
alternative = c("two.sided"),
mu = 0, paired = TRUE, var.equal = FALSE,
conf.level = 0.95)
Only pre-scan STAI differs by tx, as shown below. People are generally more anxious before the OT session (although from the descriptives below, it doesn’t look like this should be a significant difference. Hmm.)
describeBy(data$tot_pre_stai, group=data$tx)
histBy(data, "tot_pre_stai", "tx")
bi.bars(data, "tot_pre_stai", "tx", main="pre-scan STAI by session (PL, OT)", zero=FALSE)
Models should definitely include tot_bdi as a predictor, and possibly pre-scan negative affect and anxiety (tot_pre_panas_na, tot_pre_stai), as the two groups differ in these. May also want to consider including fd_mean, as the NCG group on average moved a bit more during the OT session compared to the CG group (p < .077).
Maybe age as well, since that and framewise displacement/motion overall are known to impact FNC?
For the analyses involving data from the oxytocin session, include age and sex based on the Ebner paper?
But, my analyses aren’t using group, so…?