PassiveActive
Parids Miri
3/4/2021
Include libraries
#adding message=F and warning=F supresses the comments
library(forcats) # for modifying the factors of a column
library(Hmisc)
library("knitr") # for knitting RMarkdown
library("lme4") # for mixed effects models
library("tidyverse") # for wrangling, plotting, etc.
library("broom") #for using augment
library("effects")
library("naniar") # for using the replace_with_na function.
library("emmeans") # for planned contrast in R
# library("lmerTest") #for seeing the p value of lmer
library("glmmTMB")
library("tidybayes")
library(plyr) # to use functiosn like rename
#library(dplyr)
library(tidyr)
library("kableExtra") # for making nice tables
library("janitor") # for cleaning column names
library("modelr") # for doing modeling stuff
library("brms") # Bayesian regression models with Stan
library("rstanarm") # for Bayesian models
library("cowplot") # for making figure panels
library("ggrepel") # for labels in ggplots
library("gganimate") # for animations
library("GGally") # for pairs plot
library("bayesplot") # for visualization of Bayesian model fits
library("stats")
library("gridExtra")
#library(simr) # for poewr analysis
library("boot") # for bootstrapping
#library("simpleboot") # for lm.booth
library("visreg") # for visualization of the mixed model output
#library("bootpredictlme4") # for using the function predict
library("car") # for Anova(type = "3")
library("afex") # for aov_ez and then fitting the result into Anova
library(sjPlot) # for plot_model function
library(e1071) # for skewness calculation
library(psych) # for data cleaning
library(ggplot2) # for plotting data
library(lsr) # cohen's d
library("ggpubr") # t-test
library("pracma") #strcmp
library("stringr")
library("plotly")
library(lsr) # cohen's d
library("ggpubr") # t-test filter out the participants
24 exclusion cases for the physiology: either the box went out of battery, they were wearig loution or experiment didn’t go as planned. the viberations didn’t start or they picked up the phone in the middle of the study. unexpected events.
for stressor 1 and 2, there were 8 cases that need further investigation. For some reason, the number of CDA.nSCR is detected as zero which is inacurate.
# filter list for skin conductance as well as problematic trials
# ids of 1 to 8 are from the pilots
problematic_skin_conductance_ids = c("s001","s002", "s003", "s004", "s005", "s006", "s007","s008", "s011", "s012", "s019", "s029" , "s030","s037","s039","s043", "s050","s052", "s057","s071","s089", "s095", "s108", "s120")
#----------------------------------------------------------------------------------------------------------------------
# filter list
problematic_self_report_ids = c( "s001","s002", "s003", "s004", "s005", "s006", "s007","s008", "s011", "s012", "s019","s029","s030", "s037", "s039", "s043", "s050", "s057", "s071", "s089", "s095", "s108", "s120")#passive & actives
passive = c('s034', 's045', 's048', 's054', 's059', 's060', 's063', 's075',
's077', 's083', 's098', 's102', 's107', 's110', 's113', 's119')
active = c('s016', 's023', 's025', 's027', 's032', 's040', 's044', 's067',
's068', 's078', 's080', 's085', 's090', 's096', 's109', 's111',
's118')Read the data
# read the STAI, PA, and NA dataset (all self-report) with (5 conditions)
df.qualtrics.STAI.scores.corrected.ids.excluded = read.csv("../data/qualtrics/HPC-STAI-Y1_with_individual_difference_scores.csv") %>%
dplyr::select(-X, -X1)
# read the Ledalab and Model free ways of Skin Conductance physio (10 conditions )
ledalab_results = read.csv("../data/physio_csv/ledalab_results_physio_with_missing_data.csv")
ledalab_driver_integral = read.csv("../data/physio_csv/SkinCondutance_driver_240_integral.csv")
ledalab_group_id = read.csv("../data/physio_csv/id_map_to_group.csv")
# passive active
df.breathing_act_pas = read.csv("../data/Breathing/all_passive_active_ratio.csv")
df.breathing_act_pas = df.breathing_act_pas %>%
rename(condition = `block`)
df.antegral.sc.driver = left_join(ledalab_driver_integral,ledalab_group_id, by = "id")## Warning: Column `id` joining factors with different levels, coercing to
## character vectordf.antegral.sc.driver = left_join(df.antegral.sc.driver,df.breathing_act_pas, by = c("id","condition"))## Warning: Column `id` joining character vector and factor, coercing into
## character vector#————————- March1————————— ## SC.Driver (active, passive) [show horieh]
#remove bad and outlier cases.
#stick with the T group
temp = df.antegral.sc.driver %>%
rename( block = condition) %>%
dplyr::select(id, block, antegral, group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>% # outliers
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>% #bad antegral
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2")) %>%
filter(antegral < 40000)
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining character vector and factor, coercing into
## character vector#---------------- hist passive & active--------------
hist.active.data = temp2 %>%
filter(group == "Active")
hist(hist.active.data$antegral, breaks = 100) 
hist.passive.data = temp2 %>%
filter(group == "Passive")
hist(hist.passive.data$antegral, breaks = 100) 
#---------------is active breathing predictor of SC.Driver---------
#bocketize the active and passive groups
df.antegral= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active")))## Warning: Column `id` joining character vector and factor, coercing into
## character vectorfit.lmer.1 = lmer(data = df.antegral, antegral ~ active_perc + (1|id))
summary(fit.lmer.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: antegral ~ active_perc + (1 | id)
## Data: df.antegral
##
## REML criterion at convergence: 1540.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9745 -0.5337 -0.0260 0.5947 2.1370
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 45685666 6759
## Residual 18556162 4308
## Number of obs: 77, groups: id, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 15823.69 2148.44 36.67 7.365 9.81e-09 ***
## active_perc -3865.22 3049.84 36.43 -1.267 0.213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## active_perc -0.833p1 = fit.lmer.1 %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = active_perc, y = antegral
, col = id)) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("Sum SC Driver") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
# OR (which i think is the correct one)
df.antegral= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
spread(block,antegral) %>%
mutate( diff_antegral = stressor2 - stressor1)## Warning: Column `id` joining character vector and factor, coercing into
## character vectorfit.lm = lm(data = df.antegral, diff_antegral ~ active_perc)
summary(fit.lm)##
## Call:
## lm(formula = diff_antegral ~ active_perc, data = df.antegral)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14594.1 -3299.0 -432.6 4064.2 10748.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3541.9 1575.2 -2.248 0.0308 *
## active_perc 468.8 2211.4 0.212 0.8333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5237 on 36 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.001247, Adjusted R-squared: -0.0265
## F-statistic: 0.04494 on 1 and 36 DF, p-value: 0.8333#----------------------if i bocketize passive and active
#---------------------will i see a signifciant difference?--
fit.lmer = lmer(formula = antegral ~ block * group + active_perc+ (1 | id), data = temp2)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: antegral ~ block * group + active_perc + (1 | id)
## Data: temp2
##
## REML criterion at convergence: 1145.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.80710 -0.56314 -0.03184 0.40812 1.92052
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 35357585 5946
## Residual 16228295 4028
## Number of obs: 61, groups: id, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 24719.41 12042.17 26.85 2.053 0.0500 *
## blockstressor2 -2694.26 1424.27 26.43 -1.892 0.0695 .
## groupPassive -7130.12 10501.35 26.89 -0.679 0.5030
## active_perc -11251.05 12196.92 26.66 -0.922 0.3646
## blockstressor2:groupPassive -260.81 2079.42 26.69 -0.125 0.9011
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blcks2 grpPss actv_p
## blckstrssr2 -0.059
## groupPassiv -0.983 0.068
## active_perc -0.989 0.000 0.969
## blckstrs2:P 0.022 -0.685 -0.084 0.018# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 5 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 24719. -660. 11832. 1510. 46167.
## 2 blockstressor2 -2694. 61.0 1441. -5576. 66.9
## 3 groupPassive -7130. 502. 10375. -26273. 13786.
## 4 active_perc -11251. 626. 11972. -33221. 11822.
## 5 blockstressor2:groupPassive -261. -1.16 2038. -4323. 3874.p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = block, y = antegral
, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("Sum SC Driver") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
Global mean (active passive) [show horiah]
temp = ledalab_results %>%
dplyr::select(id, block, Global.Mean..muS., group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>%
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>%
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2"))
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining factors with different levels, coercing to
## character vector#---------------- hist passive & active--------------
hist.active.data = temp2 %>%
filter(group == "Active")
hist(hist.active.data$Global.Mean..muS., breaks = 100) 
hist.passive.data = temp2 %>%
filter(group == "Passive")
hist(hist.passive.data$Global.Mean..muS., breaks = 100) 
#---------------is active breathing predictor of globalmean ---------
df.Global.Mean..muS.= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active")))## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lmer.1 = lmer(data = df.Global.Mean..muS., Global.Mean..muS. ~ active_perc + (1|id))
summary(fit.lmer.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Global.Mean..muS. ~ active_perc + (1 | id)
## Data: df.Global.Mean..muS.
##
## REML criterion at convergence: 377.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.71099 -0.44930 -0.00098 0.46592 1.58080
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 28.893 5.375
## Residual 1.429 1.195
## Number of obs: 76, groups: id, 38
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.621 1.637 36.000 8.933 1.16e-10 ***
## active_perc -4.418 2.298 36.000 -1.923 0.0624 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## active_perc -0.842#----------------
# OR which i think is more accurate
df.Global.Mean..muS.= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,block, Global.Mean..muS.,active_perc) %>%
spread(block,Global.Mean..muS.) %>%
mutate( diff_Global.Mean..muS. = stressor2 - stressor1)## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lm = lm(data = df.Global.Mean..muS., diff_Global.Mean..muS. ~ active_perc)
summary(fit.lm)##
## Call:
## lm(formula = diff_Global.Mean..muS. ~ active_perc, data = df.Global.Mean..muS.)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9152 -1.0255 0.0739 1.1762 3.6421
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5512 0.4454 -1.238 0.224
## active_perc -0.5169 0.6253 -0.827 0.414
##
## Residual standard error: 1.481 on 36 degrees of freedom
## Multiple R-squared: 0.01863, Adjusted R-squared: -0.008632
## F-statistic: 0.6833 on 1 and 36 DF, p-value: 0.4139#----------------------if i bocketize passive and active
#---------------------will i see a signifciant difference?--
fit.lmer = lmer(formula = Global.Mean..muS. ~ block * group + active_perc +(1 | id), data = temp2)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Global.Mean..muS. ~ block * group + active_perc + (1 | id)
## Data: temp2
##
## REML criterion at convergence: 270.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.70741 -0.52949 -0.08499 0.45355 1.67335
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 19.379 4.402
## Residual 1.175 1.084
## Number of obs: 60, groups: id, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.9948 8.3995 27.0281 1.785 0.0855 .
## blockstressor2 -0.8506 0.3833 28.0000 -2.219 0.0347 *
## groupPassive -0.6592 7.2446 27.0810 -0.091 0.9282
## active_perc -4.2331 8.5249 26.9999 -0.497 0.6235
## blockstressor2:groupPassive 0.2502 0.5611 28.0000 0.446 0.6590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blcks2 grpPss actv_p
## blckstrssr2 -0.023
## groupPassiv -0.986 0.026
## active_perc -0.991 0.000 0.973
## blckstrs2:P 0.016 -0.683 -0.039 0.000# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 5 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 15.0 -0.292 8.40 -2.11 31.8
## 2 blockstressor2 -0.851 0.00693 0.377 -1.59 -0.103
## 3 groupPassive -0.659 0.207 7.30 -15.2 14.4
## 4 active_perc -4.23 0.364 8.53 -20.8 13.2
## 5 blockstressor2:groupPassive 0.250 -0.0135 0.572 -0.877 1.33p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = block, y = global_mean_mu_s, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("Global Mean") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
##nSCR *SCR (active passive) [show horiah]
temp = ledalab_results %>%
dplyr::select(id, block, CDA.SCR..muS.,CDA.nSCR, group2) %>%
mutate(product_nSCR_and_SCR = CDA.SCR..muS.*CDA.nSCR) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>%
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>%
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2"))
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining factors with different levels, coercing to
## character vector#---------------- hist passive & active--------------
hist.active.data = temp2 %>%
filter(group == "Active")
hist(hist.active.data$product_nSCR_and_SCR, breaks = 100) 
hist.passive.data = temp2 %>%
filter(group == "Passive")
hist(hist.passive.data$product_nSCR_and_SCR, breaks = 100) 
#---------------is active breathing predictor of nSCR*SCR ---------
df.product_nSCR_and_SCR= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active")))## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lmer.1 = lmer(data = df.product_nSCR_and_SCR, product_nSCR_and_SCR ~ active_perc + (1|id))
summary(fit.lmer.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: product_nSCR_and_SCR ~ active_perc + (1 | id)
## Data: df.product_nSCR_and_SCR
##
## REML criterion at convergence: 699.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3707 -0.3779 -0.1150 0.2912 4.7980
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 543.7 23.32
## Residual 340.1 18.44
## Number of obs: 76, groups: id, 38
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 36.559 8.036 36.000 4.549 5.9e-05 ***
## active_perc -16.888 11.282 36.000 -1.497 0.143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## active_perc -0.842#----------------
# OR which i think is more accurate
df.product_nSCR_and_SCR= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,block, product_nSCR_and_SCR,active_perc) %>%
spread(block,product_nSCR_and_SCR) %>%
mutate( diff_product_nSCR_and_SCR = stressor2 - stressor1)## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lm = lm(data = df.product_nSCR_and_SCR, diff_product_nSCR_and_SCR ~ active_perc)
summary(fit.lm)##
## Call:
## lm(formula = diff_product_nSCR_and_SCR ~ active_perc, data = df.product_nSCR_and_SCR)
##
## Residuals:
## Min 1Q Median 3Q Max
## -115.112 -7.171 2.669 11.365 34.087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.096 7.176 -2.382 0.0226 *
## active_perc 9.695 10.074 0.962 0.3422
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.86 on 36 degrees of freedom
## Multiple R-squared: 0.02509, Adjusted R-squared: -0.001996
## F-statistic: 0.9263 on 1 and 36 DF, p-value: 0.3422#----------------------if i bocketize passive and active
#---------------------will i see a signifciant difference?--
fit.lmer = lmer(formula = product_nSCR_and_SCR ~ block * group + (1 | id), data = temp2)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: product_nSCR_and_SCR ~ block * group + (1 | id)
## Data: temp2
##
## REML criterion at convergence: 532.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7227 -0.3653 -0.0875 0.1928 4.6386
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 462.4 21.50
## Residual 338.3 18.39
## Number of obs: 60, groups: id, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.943 7.074 41.996 3.667 0.000684 ***
## blockstressor2 -8.441 6.503 28.000 -1.298 0.204891
## groupPassive 12.047 10.355 41.996 1.163 0.251247
## blockstressor2:groupPassive -6.763 9.519 28.000 -0.710 0.483296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blcks2 grpPss
## blckstrssr2 -0.460
## groupPassiv -0.683 0.314
## blckstrs2:P 0.314 -0.683 -0.460# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 4 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 25.9 0.367 7.00 12.1 40.4
## 2 blockstressor2 -8.44 0.0152 6.59 -21.6 4.90
## 3 groupPassive 12.0 -0.209 10.6 -9.19 34.0
## 4 blockstressor2:groupPassive -6.76 0.143 9.58 -26.1 11.6p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = block, y = product_n_scr_and_scr, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("product.") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
## CDA.nSCR (active,passive) [show horieh] filter(! id %in% c(“s086”, “s040”, “s048”)) outliers filter( ! id %in% c(“s020”, “s036”, “s094”)) visual driver issues 13 (passive) vs 15 (active)
temp = ledalab_results %>%
dplyr::select(id, block, CDA.nSCR, group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>%
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>%
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2"))
#filter( product_nSCR_and_SCR < 120)
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining factors with different levels, coercing to
## character vector#---------------- hist passive & active--------------
hist.active.data = temp2 %>%
filter(group == "Active")
hist(hist.active.data$CDA.nSCR, breaks = 100) 
hist.passive.data = temp2 %>%
filter(group == "Passive")
hist(hist.passive.data$CDA.nSCR, breaks = 100) 
#---------------is active breathing predictor of nSCR---------
df.nscR= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active")))## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lmer.1 = lmer(data = df.nscR, CDA.nSCR ~ active_perc + (1|id))
summary(fit.lmer.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDA.nSCR ~ active_perc + (1 | id)
## Data: df.nscR
##
## REML criterion at convergence: 673.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.68182 -0.57005 0.02598 0.55387 1.85518
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 425.7 20.63
## Residual 222.2 14.91
## Number of obs: 76, groups: id, 38
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 62.631 6.970 36.000 8.986 9.96e-11 ***
## active_perc -17.284 9.784 36.000 -1.767 0.0858 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## active_perc -0.842 #OR
df.CDA.nSCR= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,block, CDA.nSCR,active_perc) %>%
spread(block,CDA.nSCR) %>%
mutate( diff_CDA.nSCR = stressor2 - stressor1)## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lm = lm(data = df.CDA.nSCR, diff_CDA.nSCR ~ active_perc)
summary(fit.lm)##
## Call:
## lm(formula = diff_CDA.nSCR ~ active_perc, data = df.CDA.nSCR)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.514 -7.464 -0.163 6.853 29.294
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.486 4.721 -2.009 0.0521 .
## active_perc -7.876 6.628 -1.188 0.2425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.7 on 36 degrees of freedom
## Multiple R-squared: 0.03775, Adjusted R-squared: 0.01102
## F-statistic: 1.412 on 1 and 36 DF, p-value: 0.2425#----------------------if i bocketize passive and active
#---------------------will i see a signifciant difference?--
fit.lmer = lmer(formula = CDA.nSCR ~ block * group + (1 | id), data = temp2)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDA.nSCR ~ block * group + (1 | id)
## Data: temp2
##
## REML criterion at convergence: 502.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52161 -0.43140 0.04456 0.35999 1.73891
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 476.4 21.83
## Residual 131.8 11.48
## Number of obs: 60, groups: id, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 56.562 6.165 34.706 9.174 8.33e-11 ***
## blockstressor2 -15.937 4.059 28.000 -3.927 0.000511 ***
## groupPassive 7.152 9.025 34.706 0.792 0.433491
## blockstressor2:groupPassive 4.580 5.941 28.000 0.771 0.447212
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blcks2 grpPss
## blckstrssr2 -0.329
## groupPassiv -0.683 0.225
## blckstrs2:P 0.225 -0.683 -0.329# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 4 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 56.6 0.310 6.15 44.2 68.8
## 2 blockstressor2 -15.9 -0.0643 4.05 -23.9 -8.21
## 3 groupPassive 7.15 -0.160 9.00 -10.4 24.4
## 4 blockstressor2:groupPassive 4.58 -0.0511 6.17 -7.25 16.9p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = block, y = cda_n_scr, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("nSCR") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
## CDA.SCR (active,passive) [show horieh]
temp = ledalab_results %>%
dplyr::select(id, block, CDA.SCR..muS., group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>%
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>%
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2"))
# filter(CDA.SCR..muS. <1)
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining factors with different levels, coercing to
## character vector#---------------- hist passive & active--------------
hist.active.data = temp2 %>%
filter(group == "Active")
hist(hist.active.data$CDA.SCR..muS., breaks = 100) 
hist.passive.data = temp2 %>%
filter(group == "Passive")
hist(hist.passive.data$CDA.SCR..muS., breaks = 100) 
#---------------is active breathing predictor of nSCR---------
df.nscR= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active")))## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lmer.1 = lmer(data = df.nscR, CDA.SCR..muS. ~ active_perc + (1|id))
summary(fit.lmer.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDA.SCR..muS. ~ active_perc + (1 | id)
## Data: df.nscR
##
## REML criterion at convergence: 6.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8134 -0.4185 -0.0730 0.3915 3.5406
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05845 0.2418
## Residual 0.02521 0.1588
## Number of obs: 76, groups: id, 38
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.50026 0.08018 36.00000 6.239 3.33e-07 ***
## active_perc -0.15226 0.11256 36.00000 -1.353 0.185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## active_perc -0.842 #OR
df.CDA.SCR..muS.= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,block, CDA.SCR..muS.,active_perc) %>%
spread(block,CDA.SCR..muS.) %>%
mutate( diff_CDA.nSCR = stressor2 - stressor1)## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lm = lm(data = df.CDA.SCR..muS., diff_CDA.nSCR ~ active_perc)
summary(fit.lm)##
## Call:
## lm(formula = diff_CDA.nSCR ~ active_perc, data = df.CDA.SCR..muS.)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.67021 -0.08760 0.01544 0.13585 0.30545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.14141 0.05968 -2.369 0.0233 *
## active_perc 0.04723 0.08379 0.564 0.5765
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1984 on 36 degrees of freedom
## Multiple R-squared: 0.008748, Adjusted R-squared: -0.01879
## F-statistic: 0.3177 on 1 and 36 DF, p-value: 0.5765#----------------------if i bocketize passive and active
#---------------------will i see a signifciant difference?--
fit.lmer = lmer(formula = CDA.SCR..muS. ~ block * group + (1 | id), data = temp2)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CDA.SCR..muS. ~ block * group + (1 | id)
## Data: temp2
##
## REML criterion at convergence: 7.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5680 -0.4528 -0.0854 0.2849 3.3542
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.05163 0.2272
## Residual 0.02391 0.1546
## Number of obs: 60, groups: id, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.40966 0.06871 38.16826 5.962 6.32e-07 ***
## blockstressor2 -0.08998 0.05467 28.00000 -1.646 0.111
## groupPassive 0.09884 0.10058 38.16826 0.983 0.332
## blockstressor2:groupPassive -0.03636 0.08002 28.00000 -0.454 0.653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blcks2 grpPss
## blckstrssr2 -0.398
## groupPassiv -0.683 0.272
## blckstrs2:P 0.272 -0.683 -0.398# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 4 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.410 0.00214 0.0695 0.277 0.548
## 2 blockstressor2 -0.0900 0.00241 0.0547 -0.198 0.0213
## 3 groupPassive 0.0988 -0.00653 0.101 -0.0966 0.294
## 4 blockstressor2:groupPassive -0.0364 0.0000401 0.0826 -0.196 0.132p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = block, y = cda_scr_mu_s, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("CDA.SCR") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
## STAI (active, passive) [show horieh]
df.STAI = df.qualtrics.STAI.scores.corrected.ids.excluded %>%
dplyr::select(baseline_STAI_score,practice_STAI_score,prestressor1_STAI_score,poststressor1_STAI_score,poststressor2_STAI_score, id, group2) %>%
gather(STAI, STAI_value, baseline_STAI_score,practice_STAI_score,prestressor1_STAI_score,poststressor1_STAI_score,poststressor2_STAI_score) %>%
mutate(STAI = factor(STAI, levels = c("baseline_STAI_score", "practice_STAI_score", "prestressor1_STAI_score", "poststressor1_STAI_score", "poststressor2_STAI_score"), labels = c("Baseline","V-Breathing Practice", "Pre-stressor1","Post-stressor1","Post-stressor2"))) %>%
mutate(group2 = factor(group2, levels = c("Treatment", "Control"))) %>%
filter(STAI %in% c("Post-stressor1","Post-stressor2")) %>%
filter(group2 == "Treatment") %>%
dplyr::select(id, STAI, STAI_value)
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(df.STAI,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining factors with different levels, coercing to
## character vector#---------------- hist passive & active--------------
hist.active.data = temp2 %>%
filter(group == "Active")
hist(hist.active.data$STAI_value, breaks = 100) 
hist.passive.data = temp2 %>%
filter(group == "Passive")
hist(hist.passive.data$STAI_value, breaks = 100) 
#---------------is active breathing predictor of STAI---------
df.STAI.b= left_join(df.STAI,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active")))## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lmer.1 = lmer(data = df.STAI.b, STAI_value ~ active_perc + (1|id))
summary(fit.lmer.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: STAI_value ~ active_perc + (1 | id)
## Data: df.STAI.b
##
## REML criterion at convergence: 604.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94702 -0.42735 0.05965 0.48070 2.14449
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 68.05 8.249
## Residual 31.78 5.638
## Number of obs: 86, groups: id, 43
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 42.2948 2.6559 41.0000 15.925 <2e-16 ***
## active_perc -0.9726 3.7275 41.0000 -0.261 0.795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## active_perc -0.850 #OR
df.STAI.b= left_join(df.STAI,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,STAI, STAI_value,active_perc) %>%
spread(STAI,STAI_value) %>%
mutate( diff_STAI = `Post-stressor2` - `Post-stressor1`)## Warning: Column `id` joining factors with different levels, coercing to
## character vectorfit.lm = lm(data = df.STAI.b, diff_STAI ~ active_perc)
summary(fit.lm)##
## Call:
## lm(formula = diff_STAI ~ active_perc, data = df.STAI.b)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.3293 -5.0223 0.2004 6.3194 10.2249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.6707 2.1592 -1.700 0.0967 .
## active_perc 0.6846 3.0304 0.226 0.8224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.449 on 41 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.001243, Adjusted R-squared: -0.02312
## F-statistic: 0.05104 on 1 and 41 DF, p-value: 0.8224#----------------------if i bocketize passive and active
#---------------------will i see a signifciant difference?--
fit.lmer = lmer(formula = STAI_value ~ STAI * group + (1 | id), data = temp2)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: STAI_value ~ STAI * group + (1 | id)
## Data: temp2
##
## REML criterion at convergence: 429.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72503 -0.47465 0.00145 0.53128 2.08036
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 68.95 8.303
## Residual 24.25 4.924
## Number of obs: 64, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 40.9804 2.3414 38.7767 17.503 <2e-16
## STAIPost-stressor2 -3.5294 1.6890 30.0000 -2.090 0.0452
## groupPassive 1.2418 3.4198 38.7767 0.363 0.7185
## STAIPost-stressor2:groupPassive 0.1961 2.4670 30.0000 0.079 0.9372
##
## (Intercept) ***
## STAIPost-stressor2 *
## groupPassive
## STAIPost-stressor2:groupPassive
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) STAIPs-2 grpPss
## STAIPst-st2 -0.361
## groupPassiv -0.685 0.247
## STAIPst-2:P 0.247 -0.685 -0.361# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 4 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 41.0 -0.132 2.25 36.5 45.3
## 2 STAIPost-stressor2 -3.53 -0.0265 1.63 -6.75 -0.371
## 3 groupPassive 1.24 0.0685 3.36 -5.38 7.75
## 4 STAIPost-stressor2:groupPassive 0.196 -0.0298 2.44 -4.65 4.84p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = stai, y = stai_value, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.3) +
xlab("") +
ylab("STAI ") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
#—————————————— ——
*CDA.nSCR (P vs A)
#-----------------------------------------------------------
temp = ledalab_results %>%
dplyr::select(id, block, CDA.nSCR, group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040","s048","s064","s014")) %>%
filter( ! id %in% c("s020", "s036", "s094")) %>%
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2"))
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active")## Warning: Column `id` joining factors with different levels, coercing to
## character vectorp = ggplot(data = temp2, mapping = aes(x = block, y = CDA.nSCR, color = group, group = group)) +
# individual data points (jittered horizontally)
geom_point(alpha = 0.8,
position = position_jitter(width = 0.2, height = 0),
size = 2) +
# lines
#geom_smooth(method = "lm", se = F) +
# error bars
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
color = "black",
size = 1) +
# means
stat_summary(fun.y = "mean",
aes(linetype = group),
geom = "line",
shape = 21,
color = "#808080",
size = 1) +
stat_summary(fun.y = "mean",
geom = "point",
shape = 21,
fill = c("#FFFFFF","#FFFFFF","#FFFFFF","#FFFFFF"),
color = "black",
size = 4) +
theme_classic() +
xlab("") +
ylab("SC CDA.nSCR") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 15),
axis.text.y = element_text( size = 15),
axis.title.y = element_text(size = 15),
legend.text=element_text(size=15), legend.title = element_blank()) +
scale_color_manual(values=c("#f4c4a5","grey")) +
ylim(0, 100)## Warning: Ignoring unknown parameters: shapep ## Warning: Removed 3 rows containing non-finite values (stat_summary).## Warning: Removed 3 rows containing non-finite values (stat_summary).
## Warning: Removed 3 rows containing non-finite values (stat_summary).## Warning: Removed 3 rows containing missing values (geom_point).
#ggsave("CDA_nSCR_CHI2021.png", plot = p, width = 6, height = 5, units = "in")Global Mean (active vs passive (S1) then (S2))
# cohen's d
library(lsr)
# t-test
library("ggpubr")
########################### Control vs Treatment for condition stressor 1 ######################
temp = ledalab_results %>%
dplyr::select(id, block, Global.Mean..muS., group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>%
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>%
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2"))
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
df.Global.Mean..muS.= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,block, group, Global.Mean..muS.,active_perc) %>%
spread(group,Global.Mean..muS.) %>%
filter(block == "stressor2")## Warning: Column `id` joining factors with different levels, coercing to
## character vectorcohensD(df.Global.Mean..muS.$Active, df.Global.Mean..muS.$Passive)## [1] 0.7076208t.test(df.Global.Mean..muS.$Active, df.Global.Mean..muS.$Passive, alternative = "two.sided", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: df.Global.Mean..muS.$Active and df.Global.Mean..muS.$Passive
## t = -1.9287, df = 27.167, p-value = 0.06428
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.3822129 0.1965422
## sample estimates:
## mean of x mean of y
## 10.01148 13.10432print("================")## [1] "================"########################### Control vs Treatment for condition stressor 2 ######################
df.Global.Mean..muS.= left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
dplyr::select(id,block, group, Global.Mean..muS.,active_perc) %>%
spread(group,Global.Mean..muS.) %>%
filter(block == "stressor1")## Warning: Column `id` joining factors with different levels, coercing to
## character vectorcohensD(df.Global.Mean..muS.$Active, df.Global.Mean..muS.$Passive)## [1] 0.6213271t.test(df.Global.Mean..muS.$Active, df.Global.Mean..muS.$Passive, alternative = "two.sided", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: df.Global.Mean..muS.$Active and df.Global.Mean..muS.$Passive
## t = -1.688, df = 26.693, p-value = 0.1031
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.2996730 0.6144902
## sample estimates:
## mean of x mean of y
## 10.86212 13.70472SC.Driver (active vs passive (S1) then (S2))
# cohen's d
library(lsr)
# t-test
library("ggpubr")
########################### Control vs Treatment for condition stressor 1 ######################
temp = df.antegral.sc.driver %>%
rename( block = condition) %>%
dplyr::select(id, block, antegral, group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>% # outliers
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>% #bad antegral
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2")) %>%
filter(antegral < 40000)
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active") %>%
spread(group,antegral) %>%
filter(block == "stressor2")## Warning: Column `id` joining character vector and factor, coercing into
## character vectorcohensD(temp2$Active, temp2$Passive)## [1] 0.2654338t.test(temp2$Active, temp2$Passive, alternative = "two.sided", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: temp2$Active and temp2$Passive
## t = -0.73972, df = 28.989, p-value = 0.4654
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -7635.592 3579.424
## sample estimates:
## mean of x mean of y
## 11041.04 13069.12print("================")## [1] "================"########################### Control vs Treatment for condition stressor 2 ######################
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active") %>%
spread(group,antegral) %>%
filter(block == "stressor1")## Warning: Column `id` joining character vector and factor, coercing into
## character vectorcohensD(temp2$Active, temp2$Passive)## [1] 0.2131425t.test(temp2$Active, temp2$Passive, alternative = "two.sided", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: temp2$Active and temp2$Passive
## t = -0.58531, df = 27.877, p-value = 0.5631
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5900.137 3278.097
## sample estimates:
## mean of x mean of y
## 13735.29 15046.31#————————————————— # pasive vs control
temp = df.antegral.sc.driver %>%
rename( block = condition) %>%
dplyr::select(id, block, antegral, group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>% # outliers
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>% #bad antegral
filter(group2 == "Treatment") %>%
filter(block %in% c("stressor1","stressor2")) %>%
filter(antegral < 40000)
#note: just using passive breathing in s2
temp3 = df.breathing_act_pas %>%
rename(block = `condition`) %>%
filter(block == "stressor2") %>%
dplyr::select(id,active_perc)
temp2 = left_join(temp,temp3, by = "id") %>%
mutate(group = if_else(active_perc > .85,"Active",if_else(active_perc < .4,"Passive", "Exc_passive_active"))) %>%
mutate(group = factor(group, levels = c("Active","Passive","Exc_passive_active"))) %>%
filter(! group == "Exc_passive_active") %>%
select(id, block, antegral, group)## Warning: Column `id` joining character vector and factor, coercing into
## character vectortemp4 = df.antegral.sc.driver %>%
rename( block = condition) %>%
dplyr::select(id, block, antegral, group2) %>%
filter(! id %in% problematic_skin_conductance_ids) %>%
filter(! id %in% c("s086", "s040")) %>% # outliers
filter( ! id %in% c("s020", "s036", "s094","s048","s064","s014")) %>% #bad antegral
filter(group2 == "Control") %>%
filter(block %in% c("stressor1","stressor2")) %>%
filter(antegral < 40000) %>%
rename (group = group2)
df.sc.driver = left_join(temp2,temp4)## Joining, by = c("id", "block", "antegral", "group")## Warning: Column `group` joining factors with different levels, coercing to
## character vectorfit.lmer = lmer(formula = antegral ~ block * group + (1 | id), data = df.sc.driver)
summary(fit.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: antegral ~ block * group + (1 | id)
## Data: df.sc.driver
##
## REML criterion at convergence: 1167
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76220 -0.54238 -0.09275 0.47585 1.95098
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 34739287 5894
## Residual 16354831 4044
## Number of obs: 61, groups: id, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13735.29 1787.00 37.57 7.686 3.19e-09 ***
## blockstressor2 -2694.26 1429.81 26.06 -1.884 0.0707 .
## groupPassive 2246.16 2594.32 38.59 0.866 0.3920
## blockstressor2:groupPassive -218.07 2087.06 26.33 -0.104 0.9176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) blcks2 grpPss
## blckstrssr2 -0.400
## groupPassiv -0.689 0.276
## blckstrs2:P 0.274 -0.685 -0.414# bootstrap parameter estimates
boot.lmer = bootMer(fit.lmer,
FUN = fixef,
nsim = 1000)
n = 1 # not running concurrent tests
# compute confidence intervals for all the estimates in one go.
tidy(boot.lmer,conf.int=TRUE,conf.method="perc", conf.level = 1 - .05/n )## # A tibble: 4 x 6
## term statistic bias std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 13735. -101. 1791. 10261. 17214.
## 2 blockstressor2 -2694. 12.0 1431. -5431. 163.
## 3 groupPassive 2246. 125. 2698. -3114. 7500.
## 4 blockstressor2:groupPassive -218. -35.8 2152. -4558. 3615.p1 = fit.lmer %>%
augment() %>%
clean_names() %>%
ggplot(data = ., mapping = aes(x = block, y = antegral
, col = id, group = interaction(group, id))) + facet_wrap(~ group) + geom_point(alpha = 0.2) + geom_line(alpha = 0.5) +
geom_point(aes(y = fitted), color = "red") +
geom_line(aes(y = fitted),
color = "red",
alpha = 0.1) +
xlab("") +
ylab("Sum SC Driver") +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
# ylim(0,100)
p1
#————————— # nSCR and SCR in a more sys? - autorcorlation …