SR AMY analysis¶
import data¶
In [1]:
library(psych)
library(tidyverse)
In [2]:
## import sorted, wrangled data
df <- read_csv ("../data/all_wrangled.csv", na = "NA")
wrangle¶
In [3]:
## shorten amgydala variable names
df <- df %>% rename(L_amy = AngerFear_Tyszka_L_ALL_Clust,
R_amy = AngerFear_Tyszka_R_ALL_Clust)
## create composite amygdala variable
df <- df %>% rowwise () %>% mutate(Bilat_amy = mean(c(L_amy,R_amy), na.rm=T))
df[df=="NaN"]<-NA
## create SR displacement variable
#df <- df %>% rowwise () %>% mutate(SRdisp = abs(50-SRamy))
## create composite MASQ variable
df <- df %>% rowwise () %>% mutate(MASQ = sum(c(MASQDep,MASQAnx,MASQAa,MASQAd), na.rm=T))
## calculate pre - post anxiety
df <- df %>% rowwise () %>% mutate(STAI_dif = (STAI_S_post- STAI_S))
demographics¶
In [4]:
## Filter participants who have completed data (scan data & 5 questionnaires)
df <- df %>% filter(!is.na(L_amy)& !is.na(R_amy) & !is.na(NEON)& !is.na(PSSTOT)& !is.na(STAITOT)& !is.na(CESDTOT) & !is.na(MASQ))
nrow(df)
df %>% group_by(gender) %>% tally
## Ethnicity key
#1 white (non hispanic caucasian)
#2 black (african americans)
#3 asian
#4 latino/a
#5 Pacific Islander
#6 multi/other
df %>% group_by(ethnicity) %>% tally
describe(df$age)
### SRamy Subsample
df.mdot <- df %>% filter(!is.na(SRamy))
nrow(df.mdot)
df.mdot %>% group_by(gender) %>% tally
describe(df.mdot$age)
df.mdot %>% group_by(ethnicity) %>% tally
QA¶
In [5]:
### LOOK FOR OUTLIERS
#### OUTLIER DETECTION: SRAMY (1 found: MDOT 103)
## Get Descriptive statistics
stats<-describe(df$SRamy)
##### identify outliers (+/- 3 SD from mean)
no_outlier = df %>% filter(SRamy < (stats$mean + 3*stats$sd) & SRamy > (stats$mean - 3*stats$sd))
## Get Descriptive statistics without outlier
stats_no_outlier<-describe(no_outlier$SRamy)
## Create table to compare descriptive stats with/without outlier
describe <- matrix(c(stats$n, stats$mean, stats$sd, stats$min, stats$max, stats$range, stats_no_outlier$n, stats_no_outlier$mean, stats_no_outlier$sd, stats_no_outlier$min, stats_no_outlier$max, stats_no_outlier$range), ncol = 6, byrow=TRUE)
rownames(describe) <- c("Original","Removed")
colnames(describe) <- c("n", "mean","sd","min","max","range")
describe
In [6]:
### OPTIONAL: Remove SRAmy outlier here (MDOT_ID = 103, SRAmy = 4) by replacing value with NA
#### Note that participant will be included in imaging analyses but NOT SRAmy analyses
df$SRamy[df$SRamy == 4] <- NA
attach(df)
brain¶
In [7]:
## Fit simple linear regression models for all associations with SRamy, FDR correct
NEON.fitL <- lm(NEON ~L_amy)
NEON.fitR <- lm(NEON ~R_amy)
PSSTOT.fitL <- lm(PSSTOT ~ L_amy)
PSSTOT.fitR <- lm(PSSTOT ~ R_amy)
STAITOT.fitL <- lm(STAITOT ~ L_amy)
STAITOT.fitR <- lm(STAITOT ~ R_amy)
CESDTOT.fitL <- lm(CESDTOT ~ L_amy)
CESDTOT.fitR <- lm(CESDTOT ~ R_amy)
MASQ.fitL <- lm(MASQ ~ L_amy)
MASQ.fitR <- lm(MASQ ~ R_amy)
## FDR correction on FIT L & R models
a1 <- coef(summary(NEON.fitL))[,4]
a2 <- coef(summary(NEON.fitR))[,4]
b1 <- coef(summary(PSSTOT.fitL))[,4]
b2 <- coef(summary(PSSTOT.fitR))[,4]
c1 <- coef(summary(STAITOT.fitL))[,4]
c2 <- coef(summary(STAITOT.fitR))[,4]
d1 <- coef(summary(CESDTOT.fitL))[,4]
d2 <- coef(summary(CESDTOT.fitR))[,4]
e1 <- coef(summary(MASQ.fitL))[,4]
e2 <- coef(summary(MASQ.fitR))[,4]
## index p value for L & R amy
p <- c(a1[2], a2[2],
b1[2], b2[2],
c1[2], c2[2],
d1[2], d2[2],
e1[2], e2[2])
p
#### FDR correction for multiple comparison
#https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
p.adjust(p, method = "fdr", n = length(p))
SRamy¶
In [8]:
## Fit simple linear regression models for all associations with SRamy, FDR correct
summary(L_amy.fit1 <- lm(L_amy ~ SRamy))
summary(R_amy.fit1 <- lm(R_amy ~ SRamy))
## FDR correction on FIT 1 models
a <- coef(summary(L_amy.fit1))[,4]
b <- coef(summary(R_amy.fit1))[,4]
## index p value
p <- c(a[2],
b[2])
p
#### FDR correction for multiple comparison
#https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
p.adjust(p, method = "fdr", n = length(p))
test for model improvement¶
In [9]:
df.SR <- df %>% filter(!is.na(SRamy))
nrow(df.SR)
In [10]:
L_amy.fit1 <- lm(L_amy ~ SRamy, df.SR)
L_amy.fit2 <- lm(L_amy ~ poly(SRamy,2),df.SR)
anova(L_amy.fit1,L_amy.fit2) # doesn't
R_amy.fit1 <- lm(R_amy ~ SRamy, df.SR)
R_amy.fit2 <- lm(R_amy ~ poly(SRamy,2),df.SR)
anova(R_amy.fit1,R_amy.fit2) # doesn't
cor.test(L_amy, R_amy)
SRamy + brain¶
In [11]:
## Fit simple linear regression models for all associations with SRamy, FDR correct
NEON.fit1 <- lm(NEON ~ SRamy +Bilat_amy)
PSSTOT.fit1 <- lm(PSSTOT ~ SRamy +Bilat_amy)
STAITOT.fit1 <- lm(STAITOT ~ SRamy +Bilat_amy)
CESDTOT.fit1 <- lm(CESDTOT ~ SRamy +Bilat_amy)
MASQ.fit1 <- lm(MASQ ~ SRamy +Bilat_amy)
## FDR correction on Betas
a1 <- coef(summary(NEON.fit1))[,4]
b1 <- coef(summary(PSSTOT.fit1))[,4]
c1 <- coef(summary(STAITOT.fit1))[,4]
d1 <- coef(summary(CESDTOT.fit1))[,4]
e1 <- coef(summary(MASQ.fit1))[,4]
## index p value for SRamy, SRamy^2
p <- c(a1[2],
b1[2],
c1[2],
d1[2],
e1[2])
p
#### FDR correction for multiple comparison
#https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
p.adjust(p, method = "fdr", n = length(p))
Full model stats¶
In [12]:
## get p value for full model
lmp <- function (modelobject) {
if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
f <- summary(modelobject)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=F)
attributes(p) <- NULL
return(p)
}
## Get other full model stats
## get adjusted R squared
summary(MASQ.fit1)$adj.r.squared
## get F statistic
summary(MASQ.fit1)$fstatistic
# get p value for
lmp(MASQ.fit1)
In [13]:
## FDR correction on full model stats (FIT 1)
a1 <- lmp(NEON.fit1)
b1 <- lmp(PSSTOT.fit1)
c1 <- lmp(STAITOT.fit1)
e1 <- lmp(CESDTOT.fit1)
f1 <- lmp(MASQ.fit1)
## index p value for SRamy, SRamy^2
p <- c(a1[1],
b1[1],
c1[1],
e1[1],
f1[1])
p
#### FDR correction for multiple comparison
#https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
p.adjust(p, method = "fdr", n = length(p))
test for model improvement¶
In [14]:
## Fit simple linear regression models for all associations with SRamy, FDR correct
NEON.fit1 <- lm(NEON ~ SRamy + Bilat_amy, df.SR)
NEON.fit2 <- lm(NEON ~ poly(SRamy,2) + Bilat_amy, df.SR)
PSSTOT.fit1 <- lm(PSSTOT ~ SRamy + Bilat_amy, df.SR)
PSSTOT.fit2 <- lm(PSSTOT ~ poly(SRamy,2) + Bilat_amy, df.SR)
STAITOT.fit1 <- lm(STAITOT ~ SRamy + Bilat_amy, df.SR)
STAITOT.fit2 <- lm(STAITOT ~ poly(SRamy,2) + Bilat_amy, df.SR)
CESDTOT.fit1 <- lm(CESDTOT ~ SRamy + Bilat_amy, df.SR)
CESDTOT.fit2 <- lm(CESDTOT ~ poly(SRamy,2) + Bilat_amy, df.SR)
MASQ.fit1 <- lm(MASQ ~ SRamy + Bilat_amy, df.SR)
MASQ.fit2 <- lm(MASQ ~ poly(SRamy,2) + Bilat_amy, df.SR)
## Test for improvement of model fit
anova(NEON.fit1,NEON.fit2) #improves
anova(PSSTOT.fit1,PSSTOT.fit2)
anova(STAITOT.fit1,STAITOT.fit2) #improves
anova(CESDTOT.fit1,CESDTOT.fit2)
anova(MASQ.fit1,MASQ.fit2) #trend
In [15]:
## FDR correction on full model stats (FIT 2)
a1 <- lmp(NEON.fit2)
b1 <- lmp(PSSTOT.fit2)
c1 <- lmp(STAITOT.fit2)
e1 <- lmp(CESDTOT.fit2)
f1 <- lmp(MASQ.fit2)
## index p value for SRamy, SRamy^2
p <- c(a1[1],
b1[1],
c1[1],
e1[1],
f1[1])
p
#### FDR correction for multiple comparison
#https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html
p.adjust(p, method = "fdr", n = length(p))
In [16]:
options(scipen = 999)
In [17]:
summary(NEON.fit2)
plot¶
Scale variables¶
In [18]:
# SCALE VARS
df.SR$NEON <- scale(df.SR$NEON, center = TRUE, scale = TRUE)
df.SR$PSSTOT <- scale(df.SR$PSSTOT, center = TRUE, scale = TRUE)
df.SR$STAITOT <- scale(df.SR$STAITOT, center = TRUE, scale = TRUE)
df.SR$CESDTOT <- scale(df.SR$CESDTOT, center = TRUE, scale = TRUE)
df.SR$MASQ <- scale(df.SR$MASQ, center = TRUE, scale = TRUE)
df.SR$Bilat_amy <- scale(df.SR$Bilat_amy, center = TRUE, scale = TRUE)
In [20]:
# rerun regressions with scaled vars
NEON.fit1 <- lm(NEON ~ SRamy + Bilat_amy, df.SR)
NEON.fit2 <- lm(NEON ~ poly(SRamy,2) + Bilat_amy, df.SR)
PSSTOT.fit1 <- lm(PSSTOT ~ SRamy + Bilat_amy, df.SR)
PSSTOT.fit2 <- lm(PSSTOT ~ poly(SRamy,2) + Bilat_amy, df.SR)
STAITOT.fit1 <- lm(STAITOT ~ SRamy + Bilat_amy, df.SR)
STAITOT.fit2 <- lm(STAITOT ~ poly(SRamy,2) + Bilat_amy, df.SR)
CESDTOT.fit1 <- lm(CESDTOT ~ SRamy + Bilat_amy, df.SR)
CESDTOT.fit2 <- lm(CESDTOT ~ poly(SRamy,2) + Bilat_amy, df.SR)
MASQ.fit1 <- lm(MASQ ~ SRamy + Bilat_amy, df.SR)
MASQ.fit2 <- lm(MASQ ~ poly(SRamy,2) + Bilat_amy, df.SR)
Dot-whisker plot¶
In [21]:
library(dotwhisker)
In [22]:
MASQ_1 <- broom::tidy(MASQ.fit1) %>% filter(term != "Bilat_amy") %>% mutate(model = "MASQ 1")
PSS_1 <- broom::tidy(PSSTOT.fit1) %>% filter(term != "Bilat_amy") %>% mutate(model = "PSS 1")
CESD_1 <- broom::tidy(CESDTOT.fit1) %>% filter(term != "Bilat_amy") %>% mutate(model = "CESD 1")
NEON_1 <- broom::tidy(NEON.fit1) %>% filter(term != "Bilat_amy") %>% mutate(model = "NEON 1")
STAI_T_1 <- broom::tidy(STAITOT.fit1) %>% filter(term != "Bilat_amy") %>% mutate(model = "STAI_T 1")
NEON <- broom::tidy (NEON.fit2) %>% filter(term != "Bilat_amy") %>% mutate(model = "NEON")
STAI_T <- broom::tidy(STAITOT.fit2) %>% filter(term != "Bilat_amy") %>% mutate(model = "STAI-T")
MASQ <- broom::tidy(MASQ.fit2) %>% filter(term != "Bilat_amy") %>% mutate(model = "MASQ")
PSS <- broom::tidy(PSSTOT.fit2) %>% filter(term != "Bilat_amy") %>% mutate(model = "PSS")
CESD <- broom::tidy(CESDTOT.fit2) %>% filter(term != "Bilat_amy") %>% mutate(model = "CESD")
In [23]:
models_1 <- rbind(NEON_1,PSS_1, STAI_T_1, CESD_1,MASQ_1)
models_2 <- rbind(NEON,PSS, STAI_T, CESD,MASQ)
In [24]:
#dwplot(models_1,
# vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))
dwplot(models_2,
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))
Scatterplots¶
In [25]:
# NEON plot
fit <- lm(NEON ~ poly(SRamy,2), data = df.SR)
prd <- data.frame(SRamy = seq(from = range(df.SR$SRamy)[1], to = range(df.SR$SRamy)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
ggplot(prd, aes(x = SRamy, y = fit)) +
geom_smooth(aes(ymin = lci, ymax = uci, fill = ""),colour="darkgoldenrod1", stat = "identity") +
geom_point(data = df.SR, aes(x = SRamy, y = NEON, alpha = .7, size = 7)) +
ggtitle(paste("")) +
theme_minimal(base_size = 16) +
labs (x ="SRA", y = "NEON") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "gray"),
legend.position = "none") +
scale_fill_manual(values=c("darkgoldenrod1"))+
geom_vline(xintercept = 50, colour = "grey60", linetype = 2) +
scale_x_continuous (limits = c(10, 85), breaks = c(10,20,30,40,50,60,70,80)) +
scale_y_continuous(limits=c(-2, 3.7),breaks=c(-1.5,0,1.5,3))
In [26]:
# PSSTOT plot
fit <- lm(PSSTOT ~ SRamy , data = df.SR)
prd <- data.frame(SRamy = seq(from = range(df.SR$SRamy)[1], to = range(df.SR$SRamy)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
ggplot(prd, aes(x = SRamy, y = fit)) +
geom_smooth(aes(ymin = lci, ymax = uci, fill = ""),colour="darkseagreen", stat = "identity") +
geom_point(data = df.SR, aes(x = SRamy, y = PSSTOT, alpha = .7, size = 7)) +
ggtitle(paste("")) +
theme_minimal(base_size = 16) +
labs (x ="SRA", y = "PSS") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "gray"),
legend.position = "none") +
scale_fill_manual(values=c("darkseagreen"))+
geom_vline(xintercept = 50, colour = "grey60", linetype = 2)+
scale_x_continuous (limits = c(10, 85), breaks = c(10,20,30,40,50,60,70,80)) +
scale_y_continuous(limits=c(-2, 3.7),breaks=c(-1.5,0,1.5,3))
In [27]:
# CESDTOT plot
fit <- lm(CESDTOT ~ SRamy , data = df.SR)
prd <- data.frame(SRamy = seq(from = range(df.SR$SRamy)[1], to = range(df.SR$SRamy)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
ggplot(prd, aes(x = SRamy, y = fit)) +
geom_smooth(aes(ymin = lci, ymax = uci, fill = ""),colour="darkseagreen", stat = "identity") +
geom_point(data = df.SR, aes(x = SRamy, y = CESDTOT, alpha = .7, size = 7)) +
ggtitle(paste(" ")) +
theme_minimal(base_size = 16) +
labs (x ="SRA", y = "CESD") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "gray"),
legend.position = "none") +
scale_fill_manual(values=c("darkseagreen"))+
geom_vline(xintercept = 50, colour = "grey60", linetype = 2)+
scale_x_continuous (limits = c(10, 85), breaks = c(10,20,30,40,50,60,70,80)) +
scale_y_continuous(limits=c(-2, 3.7),breaks=c(-1.5,0,1.5,3))
In [28]:
# MASQ plot
fit <- lm(MASQ ~ poly(SRamy,2), data = df.SR)
prd <- data.frame(SRamy = seq(from = range(df.SR$SRamy)[1], to = range(df.SR$SRamy)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
ggplot(prd, aes(x = SRamy, y = fit)) +
geom_smooth(aes(ymin = lci, ymax = uci, fill = ""),color="darkgoldenrod1", stat = "identity",linetype=5) +
geom_point(data = df.SR, aes(x = SRamy, y = MASQ, alpha = .7, size = 7)) +
ggtitle(paste(" ")) +
theme_minimal(base_size = 16) +
labs (x ="SRA", y = "MASQ") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "gray"),
legend.position = "none") +
scale_fill_manual(values=c("darkgoldenrod1"))+
geom_vline(xintercept = 50, colour = "grey60", linetype = 2)+
scale_x_continuous (limits = c(10, 85), breaks = c(10,20,30,40,50,60,70,80)) +
scale_y_continuous(limits=c(-2, 3.7),breaks=c(-1.5,0,1.5,3))
In [29]:
# STAITOT plot
fit <- lm(STAITOT ~ poly(SRamy,2) , data = df.SR)
prd <- data.frame(SRamy = seq(from = range(df.SR$SRamy)[1], to = range(df.SR$SRamy)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
ggplot(prd, aes(x = SRamy, y = fit)) +
geom_smooth(aes(ymin = lci, ymax = uci, fill = ""),colour="darkgoldenrod1", stat = "identity") +
geom_point(data = df.SR, aes(x = SRamy, y = STAITOT, alpha = .7, size = 7)) +
ggtitle(paste("")) +
theme_minimal(base_size = 16) +
labs (x ="SRA", y = "STAI-T") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "gray"),
legend.position = "none") +
scale_fill_manual(values=c("darkgoldenrod1"))+
geom_vline(xintercept = 50, colour = "grey60", linetype = 2)+
scale_x_continuous (limits = c(10, 85), breaks = c(10,20,30,40,50,60,70,80)) +
scale_y_continuous(limits=c(-2, 3.7),breaks=c(-1.5,0,1.5,3))
In [30]:
# CESDTOT plot
fit <- lm(Bilat_amy ~ SRamy , data = df.SR)
prd <- data.frame(SRamy = seq(from = range(df.SR$SRamy)[1], to = range(df.SR$SRamy)[2], length.out = 100))
err <- predict(fit, newdata = prd, se.fit = TRUE)
prd$lci <- err$fit - 1.96 * err$se.fit
prd$fit <- err$fit
prd$uci <- err$fit + 1.96 * err$se.fit
ggplot(prd, aes(x = SRamy, y = fit)) +
geom_smooth(aes(ymin = lci, ymax = uci, fill = ""),colour="darkseagreen", stat = "identity", linetype = 5) +
geom_point(data = df.SR, aes(x = SRamy, y = Bilat_amy, alpha = .7, size = 7)) +
ggtitle(paste(" ")) +
theme_minimal(base_size = 16) +
labs (x ="SRA", y = "Bilateral TRA") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "gray"),
legend.position = "none") +
scale_fill_manual(values=c("darkseagreen"))+
geom_vline(xintercept = 50, colour = "grey60", linetype = 2)+
scale_x_continuous (limits = c(10, 85), breaks = c(10,20,30,40,50,60,70,80)) +
scale_y_continuous(limits=c(-2.5, 3.7),breaks=c(-1.5,0,1.5,3))
