Analysis for Gatzke-Kopp & Ram: Developmental Dynamics of Autonomic Function in Childhood
Overview
Little is known about the relative stability versus developmental plasticity of autonomic function, and the extent of coordination between autonomic branches, across childhood. This script works through analysis of phsychophysiological data collected from 339 children annually from kindergarten to 2nd grade. Autonomic measures of sympathetic (cardiac pre-ejection period; PEP, electrodermal activity; EDA) and parasympathetic (respiratory sinus arrhythmia; RSA) function were obtianed each year as children completed an emotion-induction task wherein they watched a series of vidio clips. After developing the data descriptives, a series of multilevel models are used to identify the portion of variance in psychophysiological function across the 31 task epochs, repeated across 3 visits, that could be attributable to the individual (trait), visit (potential developmental change), or epoch (reactivity to task).
Outline
The script covers …
A. Data Preparation and Selectivity Tests
B. Descriptives and Tests of mean level differences
C. Multilevel models for variance decompositions of each measure
D. Multilevel models for variance decompositions of bivariate couplings
E. Multilevel model examining situational differences
Prelim - Loading libraries used in this script.
Loading libraries
library(plyr)
library(psych)
library(MASS)
library(nlme)
library(effects)## Loading required package: carData## lattice theme set by effectsTheme()
## See ?effectsTheme for details.library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserReading in Data
Reading in the data
############################
####### Reading in the Data
############################
#Input a comma delimited data files:
#epoch-level data
epochdata <- read.csv(file="lionking_ALL_epochlevel_2017_0115.csv", header=TRUE)
#visit-level data
visitdata <- read.csv(file="lionking_visitlevel_2017_0329.csv", header=TRUE)
#person-level data
persondata <- read.csv(file="lionking_personlevel_2017_0103.csv", header=TRUE)
#id data
iddata <- read.csv(file="T1sample.csv", header=TRUE)
iddata$insample <- 1
names(iddata) <- c("id","insample")This step read in a set of datafiles
epochdata (epoch-level repeated measures)
visitdata (visit-level repeated measures)
persondata (person-level variables)
iddata (master list of ids with designation for inclusion in sample)
A. Data Preparation and Selectivity Tests
Calculating visit-level summary data imean,isd, ietc …
#Computing Individual parameters
istats.visit <- ddply(epochdata, c("id","visit"), summarize,
imean.rsa=mean(rsa, na.rm=TRUE),
isd.rsa=sd(rsa, na.rm=TRUE),
iskew.rsa=skew(rsa, na.rm=TRUE),
ikurt.rsa=kurtosi(rsa, na.rm=TRUE),
icount.rsa=sum(!is.na(rsa)),
imean.pep=mean(pep, na.rm=TRUE),
isd.pep=sd(pep, na.rm=TRUE),
iskew.pep=skew(pep, na.rm=TRUE),
ikurt.pep=kurtosi(pep, na.rm=TRUE),
icount.pep=sum(!is.na(pep)),
imean.ns_scr=mean(ns_scr, na.rm=TRUE),
isd.ns_scr=sd(ns_scr, na.rm=TRUE),
iskew.ns_scr=skew(ns_scr, na.rm=TRUE),
ikurt.ns_scr=kurtosi(ns_scr, na.rm=TRUE),
icount.ns_scr=sum(!is.na(ns_scr)))
#Calculating the resting baseline scores
baselineepochdata <- epochdata[which(epochdata$epoch <= 1),]
istats.resting <- ddply(baselineepochdata, c("id","visit"), summarize,
resting.rsa=mean(rsa,na.rm=TRUE),
resting.pep=mean(pep,na.rm=TRUE),
resting.ns_scr=mean(ns_scr,na.rm=TRUE))Merging visit-level and visit-level physio summary variables.
istats.visit <- merge(istats.resting,istats.visit,by=c("id","visit"), all=TRUE)
visitdataplus <- merge(visitdata,istats.visit,by=c("id","visit"), all=TRUE)
#adding in visit dummy variables
visitdataplus$visit0 <- ifelse(visitdataplus$visit==0,1,0)
visitdataplus$visit1 <- ifelse(visitdataplus$visit==1,1,0)
visitdataplus$visit2 <- ifelse(visitdataplus$visit==2,1,0)Making person-level summaries
#making person-level means
istats.person <- ddply(visitdataplus, c("id"), summarize,
ipersonmean.rsa=mean(imean.rsa, na.rm=TRUE),
ipersonsd.rsa=sd(imean.rsa, na.rm=TRUE),
ipersonmean.pep=mean(imean.pep, na.rm=TRUE),
ipersonsd.pep=sd(imean.pep, na.rm=TRUE),
ipersonmean.ns_scr=mean(imean.ns_scr, na.rm=TRUE),
ipersonsd.ns_scr=sd(imean.ns_scr, na.rm=TRUE),
ipersonmean.emo=mean(emo, na.rm=TRUE))
#datavisit2 <- visitdataplus[which(visitdataplus$visit ==2),]
#emo2data <- datavisit2[,c("id","emo","cp")]
#names(emo2data) <- c("id","emo2","cp2")
#emo2data$emo2_c <- emo2data$emo2 - mean(emo2data$)
#istats.person <- merge(istats.person,emo2data,by="id",all=TRUE)
#reshaping visit data to wide
visitvarnames <- (names(visitdata))[3:12]
visitwide <- reshape(data=visitdata,
timevar=c("visit"),
idvar=c("id"),
v.names=visitvarnames,
direction="wide", sep="_")
names(visitwide)## [1] "id" "visitdate_0" "visitage_0" "emo_0" "pro_0"
## [6] "soc_0" "ext_0" "es_0" "cp_0" "hy_0"
## [11] "pp_0" "visitdate_1" "visitage_1" "emo_1" "pro_1"
## [16] "soc_1" "ext_1" "es_1" "cp_1" "hy_1"
## [21] "pp_1" "visitdate_2" "visitage_2" "emo_2" "pro_2"
## [26] "soc_2" "ext_2" "es_2" "cp_2" "hy_2"
## [31] "pp_2"#combining visit-level wide data
visitdataplus1 <- merge(visitdataplus,visitwide,by="id",all=TRUE)
visitdataplus2 <- merge(visitdataplus1,istats.person,by="id",all=TRUE)
visitdataplus3 <- merge(visitdataplus2,iddata,by.x ="id") #, by.xall=TRUE)
#checking all persons are in-sample persons
table(visitdataplus3$insample, useNA="always")##
## 1 <NA>
## 1017 0#making wide visit data
visitvarnames <- names(visitdataplus3)[c(3:15,16,21,26)]
visitdataplus3 <- visitdataplus3[order(visitdataplus3$id,visitdataplus3$visit),]
visitplus3wide <- reshape(data=visitdataplus3[,c(1:15,16,21,26)],
timevar=c("visit"),
idvar=c("id"),
v.names=visitvarnames,
direction="wide", sep="_")
#indiciating all in sample
visitplus3wide$insample <- 1
#making missing trackers
#visitplus3wide$present_emo <- rowSums(!is.na(visitplus3wide[,c("emo_0","emo_1","emo_2")]))
visitplus3wide$present_resting.rsa <- rowSums(!is.na(visitplus3wide[,c("resting.rsa_0","resting.rsa_1","resting.rsa_2")]))
visitplus3wide$present_imean.rsa <- rowSums(!is.na(visitplus3wide[,c("imean.rsa_0","imean.rsa_1","imean.rsa_2")]))
visitplus3wide$present_visit0 <- ifelse(#!is.na(visitplus3wide$emo_0) == TRUE |
!is.na(visitplus3wide$resting.rsa_0) == TRUE |
!is.na(visitplus3wide$imean.rsa_0) == TRUE, 1, 0)
visitplus3wide$present_visit1 <- ifelse(#!is.na(visitplus3wide$emo_1) == TRUE |
!is.na(visitplus3wide$resting.rsa_1) == TRUE |
!is.na(visitplus3wide$imean.rsa_1) == TRUE, 1, 0)
visitplus3wide$present_visit2 <- ifelse(#!is.na(visitplus3wide$emo_2) == TRUE |
!is.na(visitplus3wide$resting.rsa_2) == TRUE |
!is.na(visitplus3wide$imean.rsa_2) == TRUE, 1, 0)
visitplus3wide$present_all <- visitplus3wide$present_visit0 + visitplus3wide$present_visit1 + visitplus3wide$present_visit2Tables about data contribution
#tables about data contribution
table(visitplus3wide$present_all)##
## 0 1 2 3
## 5 54 64 216216/339## [1] 0.6371681(216+64)/339## [1] 0.8259587#table(visitplus3wide$present_emo)
table(visitplus3wide$present_resting.rsa)##
## 0 1 2 3
## 10 54 90 185#THIS IS THE SAMPLE THAT COMES OUT IN LATER ANALYSIS
#54 + 90 + 185 = 329
table(visitplus3wide$present_imean.rsa)##
## 0 1 2 3
## 5 54 64 216table(visitplus3wide$present_visit0, useNA = "always")##
## 0 1 <NA>
## 26 313 0table(visitplus3wide$present_visit1, useNA = "always")##
## 0 1 <NA>
## 69 270 0table(visitplus3wide$present_visit2, useNA = "always")##
## 0 1 <NA>
## 92 247 0#table(visitplus3wide$present_emo,visitplus3wide$present_resting.rsa)
table(visitplus3wide$present_resting.rsa, visitplus3wide$present_imean.rsa)##
## 0 1 2 3
## 0 5 5 0 0
## 1 0 49 3 2
## 2 0 0 61 29
## 3 0 0 0 185table(visitplus3wide$present_visit0, visitplus3wide$present_visit1)##
## 0 1
## 0 6 20
## 1 63 250table(visitplus3wide$present_visit0, visitplus3wide$present_visit2)##
## 0 1
## 0 10 16
## 1 82 231table(visitplus3wide$present_visit1, visitplus3wide$present_visit2)##
## 0 1
## 0 53 16
## 1 39 231table(visitplus3wide$present_visit0,visitplus3wide$present_visit1, visitplus3wide$present_visit2)## , , = 0
##
##
## 0 1
## 0 5 5
## 1 48 34
##
## , , = 1
##
##
## 0 1
## 0 1 15
## 1 15 216#combining data with demographics
persondataplus <- merge(visitplus3wide,persondata,by="id", all=TRUE)
table(persondataplus$insample, useNA = "always")##
## 1 <NA>
## 339 15#keeping to in sample participants
persondata339 <- persondataplus[which(persondataplus$insample==1),]
table(persondata339$present_all)##
## 0 1 2 3
## 5 54 64 216#216+64+54+5 = 339
describe(persondata339$present_all)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 339 2.45 0.81 3 2.58 0 0 3 3 -1.15 0.01 0.04Checking demographic differences in data provision.
#sex differences in visit presence
chisq.test(table(persondata339$present_visit0,persondata339$sex))## Warning in chisq.test(table(persondata339$present_visit0,
## persondata339$sex)): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(persondata339$present_visit0, persondata339$sex)
## X-squared = 0.084908, df = 1, p-value = 0.7708chisq.test(table(persondata339$present_visit1,persondata339$sex))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(persondata339$present_visit1, persondata339$sex)
## X-squared = 9.8591, df = 1, p-value = 0.00169chisq.test(table(persondata339$present_visit2,persondata339$sex))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(persondata339$present_visit2, persondata339$sex)
## X-squared = 3.5667, df = 1, p-value = 0.05895table(persondata339$present_visit0)##
## 0 1
## 26 313table(persondata339$present_visit1)##
## 0 1
## 69 270table(persondata339$present_visit2)##
## 0 1
## 92 247table(persondata339$sex)##
## 0 1
## 116 205#sex differences in visits
allcheck1 <- lm(present_all ~ sex,
data=persondata339,
na.action=na.exclude)
summary(allcheck1)##
## Call:
## lm(formula = present_all ~ sex, data = persondata339, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6049 -0.3621 0.3951 0.3951 0.6379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.36207 0.06892 34.271 < 2e-16 ***
## sex 0.24281 0.08625 2.815 0.00518 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7423 on 319 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.02424, Adjusted R-squared: 0.02118
## F-statistic: 7.926 on 1 and 319 DF, p-value: 0.005177#ethnic differences in visit presence
chisq.test(table(persondata339$present_visit0,persondata339$ethniccode))## Warning in chisq.test(table(persondata339$present_visit0,
## persondata339$ethniccode)): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: table(persondata339$present_visit0, persondata339$ethniccode)
## X-squared = 1.9573, df = 3, p-value = 0.5813chisq.test(table(persondata339$present_visit1,persondata339$ethniccode))## Warning in chisq.test(table(persondata339$present_visit1,
## persondata339$ethniccode)): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: table(persondata339$present_visit1, persondata339$ethniccode)
## X-squared = 2.5854, df = 3, p-value = 0.4601chisq.test(table(persondata339$present_visit2,persondata339$ethniccode))## Warning in chisq.test(table(persondata339$present_visit2,
## persondata339$ethniccode)): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: table(persondata339$present_visit2, persondata339$ethniccode)
## X-squared = 1.2415, df = 3, p-value = 0.7431#group differences in visit presence
chisq.test(table(persondata339$present_visit0,persondata339$group))## Warning in chisq.test(table(persondata339$present_visit0,
## persondata339$group)): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(persondata339$present_visit0, persondata339$group)
## X-squared = 1.4067, df = 1, p-value = 0.2356chisq.test(table(persondata339$present_visit1,persondata339$group))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(persondata339$present_visit1, persondata339$group)
## X-squared = 0.088582, df = 1, p-value = 0.766chisq.test(table(persondata339$present_visit2,persondata339$group))##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(persondata339$present_visit2, persondata339$group)
## X-squared = 0.69889, df = 1, p-value = 0.4032# #sex differences in emo ratings provision
# emocheck1 <- lm(present_emo ~ sex,
# data=persondata339,
# na.action=na.exclude)
# summary(emocheck1)
#
# #ethnic differences in emo ratings provision
# emocheck2 <- lm(present_emo ~ factor(ethniccode),
# data=persondata339,
# na.action=na.exclude)
# summary(emocheck2)
#
# #emo differences in emo ratings provision
# emocheck3 <- lm(present_emo ~ emo_0,
# data=persondata339,
# na.action=na.exclude)
# summary(emocheck3)
#
# #group differences in emo ratings provision
# emocheck4 <- lm(present_emo ~ group,
# data=persondata339,
# na.action=na.exclude)
# summary(emocheck4)
#sex differences in rsa ratings provision
rsacheck1 <- lm(present_resting.rsa ~ sex,
data=persondata339,
na.action=na.exclude)
summary(rsacheck1)##
## Call:
## lm(formula = present_resting.rsa ~ sex, data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4829 -0.4829 0.5171 0.5171 0.7759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.22414 0.07377 30.151 < 2e-16 ***
## sex 0.25879 0.09231 2.804 0.00536 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7945 on 319 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.02405, Adjusted R-squared: 0.02099
## F-statistic: 7.86 on 1 and 319 DF, p-value: 0.005364#ethnic differences in rsa ratings provision
rsacheck2 <- lm(present_resting.rsa ~ factor(ethniccode),
data=persondata339,
na.action=na.exclude)
summary(rsacheck2)##
## Call:
## lm(formula = present_resting.rsa ~ factor(ethniccode), data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4009 -0.4009 0.5991 0.5991 0.6897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.40090 0.05412 44.365 <2e-16 ***
## factor(ethniccode)2 -0.09056 0.15921 -0.569 0.570
## factor(ethniccode)3 -0.01284 0.11240 -0.114 0.909
## factor(ethniccode)4 -0.06757 0.46867 -0.144 0.885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8063 on 317 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.001067, Adjusted R-squared: -0.008387
## F-statistic: 0.1129 on 3 and 317 DF, p-value: 0.9525#emo differences in rsa ratings provision
rsacheck3 <- lm(present_resting.rsa ~ emo_0,
data=persondata339,
na.action=na.exclude)
summary(rsacheck3)##
## Call:
## lm(formula = present_resting.rsa ~ emo_0, data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3987 -0.4197 0.5383 0.6153 0.7063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.23766 0.15745 14.212 <2e-16 ***
## emo_0 0.04201 0.03975 1.057 0.292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.792 on 283 degrees of freedom
## (54 observations deleted due to missingness)
## Multiple R-squared: 0.003931, Adjusted R-squared: 0.0004109
## F-statistic: 1.117 on 1 and 283 DF, p-value: 0.2915rsacheck3 <- lm(present_resting.rsa ~ ext_0,
data=persondata339,
na.action=na.exclude)
summary(rsacheck3)##
## Call:
## lm(formula = present_resting.rsa ~ ext_0, data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4119 -0.4281 0.5233 0.6044 0.7584
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.53346 0.10818 23.419 <2e-16 ***
## ext_0 -0.05675 0.04040 -1.405 0.161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7908 on 283 degrees of freedom
## (54 observations deleted due to missingness)
## Multiple R-squared: 0.006923, Adjusted R-squared: 0.003414
## F-statistic: 1.973 on 1 and 283 DF, p-value: 0.1612rsacheck3 <- lm(present_resting.rsa ~ hy_0,
data=persondata339,
na.action=na.exclude)
summary(rsacheck3)##
## Call:
## lm(formula = present_resting.rsa ~ hy_0, data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4069 -0.3998 0.5895 0.6073 0.6251
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.410488 0.080200 30.056 <2e-16 ***
## hy_0 -0.003555 0.014310 -0.248 0.804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.794 on 282 degrees of freedom
## (55 observations deleted due to missingness)
## Multiple R-squared: 0.0002188, Adjusted R-squared: -0.003327
## F-statistic: 0.0617 on 1 and 282 DF, p-value: 0.804#group differences in rsa ratings provision
rsacheck4 <- lm(present_resting.rsa ~ group,
data=persondata339,
na.action=na.exclude)
summary(rsacheck4)##
## Call:
## lm(formula = present_resting.rsa ~ group, data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.301 -0.301 0.472 0.699 0.699
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.52800 0.07124 35.48 <2e-16 ***
## group -0.22698 0.09117 -2.49 0.0133 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7965 on 319 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.01906, Adjusted R-squared: 0.01598
## F-statistic: 6.198 on 1 and 319 DF, p-value: 0.0133#sex differences in rsa ratings provision
rsacheck1 <- lm(present_imean.rsa ~ sex,
data=persondata339,
na.action=na.exclude)
summary(rsacheck1)##
## Call:
## lm(formula = present_imean.rsa ~ sex, data = persondata339, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6049 -0.3621 0.3951 0.3951 0.6379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.36207 0.06892 34.271 < 2e-16 ***
## sex 0.24281 0.08625 2.815 0.00518 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7423 on 319 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.02424, Adjusted R-squared: 0.02118
## F-statistic: 7.926 on 1 and 319 DF, p-value: 0.005177#ethnic differences in rsa ratings provision
rsacheck2 <- lm(present_imean.rsa ~ factor(ethniccode),
data=persondata339,
na.action=na.exclude)
summary(rsacheck2)##
## Call:
## lm(formula = present_imean.rsa ~ factor(ethniccode), data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5451 -0.5451 0.4550 0.4550 0.6207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.54505 0.05047 50.427 <2e-16 ***
## factor(ethniccode)2 -0.16573 0.14848 -1.116 0.265
## factor(ethniccode)3 -0.06743 0.10482 -0.643 0.520
## factor(ethniccode)4 0.12162 0.43709 0.278 0.781
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.752 on 317 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.004971, Adjusted R-squared: -0.004446
## F-statistic: 0.5279 on 3 and 317 DF, p-value: 0.6634# #emo differences in rsa ratings provision
# rsacheck3 <- lm(present_imean.rsa ~ emo_0,
# data=persondata339,
# na.action=na.exclude)
# summary(rsacheck3)
#group differences in rsa ratings provision
rsacheck4 <- lm(present_imean.rsa ~ group,
data=persondata339,
na.action=na.exclude)
summary(rsacheck4)##
## Call:
## lm(formula = present_imean.rsa ~ group, data = persondata339,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5760 -0.4796 0.4240 0.5204 0.5204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.57600 0.06708 38.400 <2e-16 ***
## group -0.09641 0.08585 -1.123 0.262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.75 on 319 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.003938, Adjusted R-squared: 0.0008152
## F-statistic: 1.261 on 1 and 319 DF, p-value: 0.2623Merging visit-level data into long file, creating deviation scores and centering
#epochdataplus <- merge(epochdata,visitdataplus3,by=c("id","visit"),all=TRUE)
#table(epochdataplus$insample, useNA = "always")
#nots <- epochdataplus[which(is.na(epochdataplus$insample) == TRUE),]
epochdataplus <- merge(epochdata,visitdataplus3,by.x=c("id","visit"))
table(epochdataplus$insample, useNA = "always")##
## 1 <NA>
## 32544 0#creating epoch-level deviation scores
epochdataplus$rsa_d <- epochdataplus$rsa - epochdataplus$imean.rsa
epochdataplus$ns_scr_d <- epochdataplus$ns_scr - epochdataplus$imean.ns_scr
epochdataplus$pep_d <- epochdataplus$pep - epochdataplus$imean.pep
#creating visit-level deviation scores
epochdataplus$imean.rsa_d <- epochdataplus$imean.rsa - epochdataplus$ipersonmean.rsa
epochdataplus$imean.ns_scr_d <- epochdataplus$imean.ns_scr - epochdataplus$ipersonmean.ns_scr
epochdataplus$imean.pep_d <- epochdataplus$imean.pep - epochdataplus$ipersonmean.pep
#creating visit-level deviation
#epochdataplus$emo_v <- epochdataplus$emo - epochdataplus$ipersonmean.emo
#centering imeans
epochdataplus$ipersonmean.rsa_c <- epochdataplus$ipersonmean.rsa - mean(istats.person$ipersonmean.rsa, na.rm=TRUE)
epochdataplus$ipersonmean.pep_c <- epochdataplus$ipersonmean.pep - mean(istats.person$ipersonmean.pep, na.rm=TRUE)
epochdataplus$ipersonmean.ns_scr_c <- epochdataplus$ipersonmean.ns_scr - mean(istats.person$ipersonmean.ns_scr, na.rm=TRUE)
# epochdataplus$ipersonmean.emo_c <- epochdataplus$ipersonmean.emo - mean(istats.person$ipersonmean.emo, na.rm=TRUE)
# epochdataplus$emo2_c <- epochdataplus$emo2 - mean(istats.person$emo2, na.rm=TRUE)
# epochdataplus$cp2_c <- epochdataplus$cp2 - mean(istats.person$cp2, na.rm=TRUE)
# describe(epochdataplus)
#need an additional epoch time variable to plot all at once
epochdataplus$epochtime <- (34*epochdataplus$visit)+epochdataplus$epoch
#reordering for easy viewing
epochdataplus <- epochdataplus[order(epochdataplus$id,epochdataplus$visit,epochdataplus$epochtime),]
#checking values of epochtime variable
unique(epochdataplus$epochtime)## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
## [24] 23 24 25 26 27 28 29 30 31 34 35 36 37 38 39 40 41 42 43 44 45 46 47
## [47] 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 68 69 70 71 72
## [70] 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
## [93] 96 97 98 99Plots of the data structure
For Figure 1 in paper - A plot that explains the data structure.
#creating new segement variables
epochdataplus$bigsegment <- ifelse(epochdataplus$segment == "initialbase" |
epochdataplus$segment == "postlion" |
epochdataplus$segment == "prelion" |
epochdataplus$segment =="fixation", "baseline", as.character(epochdataplus$segment))
epochdataplus$bigsegment <- ifelse(epochdataplus$epoch == 10, "fear", epochdataplus$bigsegment)
epochdataplus$bigsegment <- ifelse(epochdataplus$epoch == 16, "sad", epochdataplus$bigsegment)
epochdataplus$bigsegment <- ifelse(epochdataplus$epoch == 23, "happy", epochdataplus$bigsegment)
epochdataplus$bigsegment <- ifelse(epochdataplus$epoch == 29, "anger", epochdataplus$bigsegment)
epochdataplus$bigsegmentf <- factor(epochdataplus$bigsegment, levels= c("baseline","fear","sad","happy","anger"))
#checking variable type
str(epochdataplus$bigsegmentf)## Factor w/ 5 levels "baseline","fear",..: 1 1 1 1 2 2 2 2 2 2 ...#another bigsegment coding for emotions (for later analysis)
epochdataplus$bigsegmentemo <- epochdataplus$bigsegment
epochdataplus$bigsegmentemo <- ifelse(epochdataplus$epoch == 11, "sad", epochdataplus$bigsegmentemo)
epochdataplus$bigsegmentemo <- ifelse(epochdataplus$epoch == 17, "happy", epochdataplus$bigsegmentemo)
epochdataplus$bigsegmentemo <- ifelse(epochdataplus$epoch == 24, "anger", epochdataplus$bigsegmentemo)
epochdataplus$bigsegmentemo <- factor(epochdataplus$bigsegmentemo)
levels(epochdataplus$bigsegmentemo)## [1] "anger" "baseline" "fear" "happy" "sad"#selecting one individual for plotting
dataindiv <- subset(epochdataplus, id ==101176) #113730)
#setting colors for background
cols <- c("baseline" ="#bdbdbd", "fear"="#edf8e9", "sad"="#c7e9c0","happy"="#f2f0f7","anger"="#dadaeb")
#making plot of individual data
ggplot(data = dataindiv, aes(x=epochtime), legend=FALSE) +
#background conditions
geom_rect(mapping=aes(xmin=epochtime, xmax=epochtime+1, ymin=3, ymax=10, fill=bigsegment), alpha=0.8) +
#visit breaks
geom_vline(aes(xintercept = 33), color="gray") +
geom_vline(aes(xintercept = 67), color="gray") +
#visit-level means
geom_line(aes(x=epochtime,y = imean.rsa), lty=1, size=1,colour="black") +
geom_point(aes(x=epochtime+.5, y=ipersonmean.rsa), color="black", shape=15,size=2) +
geom_line(aes(x=epochtime+.5, y=ipersonmean.rsa), color="black",size=1.5) +
#epoch-level
geom_step(data = subset(dataindiv, visit == 0), aes(x=epochtime,y = rsa), lty=1, size=.5,colour="black") +
geom_step(data = subset(dataindiv, visit == 1), aes(x=epochtime,y = rsa), lty=1, size=.5,colour="black") +
geom_step(data = subset(dataindiv, visit == 2), aes(x=epochtime,y = rsa), lty=1, size=.5,colour="black") +
geom_segment(data = subset(dataindiv, epoch==31), aes(x=epochtime, xend=epochtime+1,y = rsa, yend=rsa), lty=1, size=.5,colour="black") +
xlab("EpochTime") +
ylab("RSA") + #ylim(1,7) +
guides(fill=FALSE) +
scale_x_continuous(breaks=0.5+c(0:31,34:65,68:99)) +
scale_fill_manual(values = cols) +
theme_classic() +
theme(axis.title.y = element_text(size = rel(1.5), angle = 90),
axis.text.y = element_text(size = rel(1.2)),
axis.text.x = element_blank())
Saving figure as a jpg file
ggsave(filename = "Figure1.png",
width = 7.5, height = 3, dpi = 300, units = "in", device='png')Our objective is also to look at the coupling between two variables. within person, across visits, and across persons.
So - A few sample plots to get an idea of how those data structures look.
ggplot(data = subset(epochdataplus, id == 101176), aes(x=epochtime), legend=FALSE, group=id) +
geom_rect(mapping=aes(xmin=epochtime-.5, xmax=epochtime+.5, ymin=-3, ymax=3, fill=visit), alpha=0.8) +
#geom_line(aes(x=epochtime,y = ipersonmean.rsa, group=id), lty=1, size=1,colour="black") +
#geom_line(aes(x=epochtime,y = imean.rsa_d, group=id), lty=1, size=1.5,colour="blue") +
geom_line(aes(x=epochtime,y = scale(rsa_d), group=id), lty=1, size=.5,colour="red") +
geom_line(aes(x=epochtime,y = scale(pep_d), group=id), lty=1, size=.5,colour="blue") +
#geom_point(aes(x=epochtime,y = rsa, group=id), size=2,colour="red") +
xlab("EpochTime") +
ylab("Physio Level") #+ #ylim(1,7) +
#scale_x_continuous(breaks=seq(0,100,by=10))
ggplot(data = subset(epochdataplus, id == 101176), aes(x=epochtime), legend=FALSE, group=id) +
geom_rect(mapping=aes(xmin=epochtime-.5, xmax=epochtime+.5, ymin=-3, ymax=3, fill=visit), alpha=0.8) +
#geom_line(aes(x=epochtime,y = ipersonmean.rsa, group=id), lty=1, size=1,colour="black") +
#geom_line(aes(x=epochtime,y = imean.rsa_d, group=id), lty=1, size=1.5,colour="blue") +
geom_line(aes(x=epochtime,y = scale(imean.rsa_d), group=id), lty=1, size=.5,colour="red") +
geom_line(aes(x=epochtime,y = scale(imean.pep_d), group=id), lty=1, size=.5,colour="blue") +
#geom_point(aes(x=epochtime,y = rsa, group=id), size=2,colour="red") +
xlab("EpochTime") +
ylab("Physio Level") #+ #ylim(1,7) +
#scale_x_continuous(breaks=seq(0,100,by=10))
ggplot(data = subset(epochdataplus, id <= 106000), aes(x=ipersonmean.pep_c,y = ipersonmean.rsa_c), legend=FALSE) +
geom_point(aes(x=ipersonmean.pep_c,y = ipersonmean.rsa_c, group=id), size=2,colour="black") +
geom_smooth(method="lm", color="red") +
xlab("PEP (sample centered)") +
ylab("RSA (sample centered)") #+ #ylim(1,7) +
#scale_x_continuous(breaks=seq(0,100,by=10)) B. Descriptives and Tests of mean level differences
Descriptives
#Getting Ns
length(unique(visitdataplus2[which(visitdataplus3$visit ==0 & is.na(visitdataplus3$resting.rsa) == FALSE), ]$id))## [1] 286length(unique(visitdataplus2[which(visitdataplus3$visit ==1 & is.na(visitdataplus3$resting.rsa) == FALSE), ]$id))## [1] 262length(unique(visitdataplus2[which(visitdataplus3$visit ==2 & is.na(visitdataplus3$resting.rsa) == FALSE), ]$id))## [1] 241#physio variables
describeBy(visitdataplus3$resting.rsa,group=visitdataplus3$visit, mat=TRUE, digits=3)## item group1 vars n mean sd median trimmed mad min max
## X11 1 0 1 286 7.283 1.285 7.290 7.289 1.215 2.487 10.362
## X12 2 1 1 262 7.406 1.306 7.432 7.445 1.235 3.626 10.266
## X13 3 2 1 241 7.497 1.275 7.522 7.529 1.319 4.059 10.259
## range skew kurtosis se
## X11 7.875 -0.275 0.618 0.076
## X12 6.639 -0.270 0.008 0.081
## X13 6.199 -0.216 -0.317 0.082describeBy(visitdataplus3$resting.pep,group=visitdataplus3$visit, mat=TRUE, digits=3)## item group1 vars n mean sd median trimmed mad min
## X11 1 0 1 272 103.685 14.797 100.925 102.911 13.232 57.90
## X12 2 1 1 254 104.795 15.483 103.200 104.944 16.568 68.75
## X13 3 2 1 232 102.961 15.105 100.000 102.688 16.309 69.00
## max range skew kurtosis se
## X11 165.75 107.85 0.554 0.550 0.897
## X12 136.00 67.25 -0.018 -0.804 0.972
## X13 142.00 73.00 0.193 -0.821 0.992describeBy(visitdataplus3$resting.ns_scr,group=visitdataplus3$visit, mat=TRUE, digits=3)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 0 1 254 4.774 3.736 4.0 4.395 4.448 0 19.0 19.0
## X12 2 1 1 225 2.649 2.970 1.5 2.144 2.224 0 12.5 12.5
## X13 3 2 1 218 1.550 2.313 0.5 1.048 0.741 0 13.0 13.0
## skew kurtosis se
## X11 0.778 0.090 0.234
## X12 1.321 1.016 0.198
## X13 2.263 5.764 0.157# #emotion regulation variable
# describeBy(visitdataplus3$emo,group=visitdataplus3$visit, mat=TRUE, digits=3)Testing Mean Differences across Visits
#Setting up for anova
visitdataplus3$factorvisit <- factor(visitdataplus3$visit)
#mean diff tests and post-hoc for RSA
lme_restingrsa = lme(resting.rsa ~ factorvisit,
random = ~1|id,
correlation=corCompSymm(form=~1|id),
data=visitdataplus3,
na.action=na.exclude)
#summary(lme_restingrsa)
anova(lme_restingrsa)## numDF denDF F-value p-value
## (Intercept) 1 458 14679.568 <.0001
## factorvisit 2 458 3.624 0.0274summary(glht(lme_restingrsa, linfct=mcp(factorvisit = "Tukey")), test = adjusted(type = "bonferroni"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = resting.rsa ~ factorvisit, data = visitdataplus3,
## random = ~1 | id, correlation = corCompSymm(form = ~1 | id),
## na.action = na.exclude)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 0.11654 0.07949 1.466 0.4278
## 2 - 0 == 0 0.21859 0.08137 2.686 0.0217 *
## 2 - 1 == 0 0.10205 0.08138 1.254 0.6296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)#mean diff tests and post-hoc for PEP
lme_restingpep = lme(resting.pep ~ factorvisit,
random = ~1|id,
correlation=corCompSymm(form=~1|id),
data=visitdataplus3,
na.action=na.exclude)
#summary(lme_restingpep)
anova(lme_restingpep)## numDF denDF F-value p-value
## (Intercept) 1 434 20678.856 <.0001
## factorvisit 2 434 2.286 0.1029summary(glht(lme_restingpep, linfct=mcp(factorvisit = "Tukey")), test = adjusted(type = "bonferroni"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = resting.pep ~ factorvisit, data = visitdataplus3,
## random = ~1 | id, correlation = corCompSymm(form = ~1 | id),
## na.action = na.exclude)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 1.0972 0.9616 1.141 0.7616
## 2 - 0 == 0 -0.9973 0.9832 -1.014 0.9312
## 2 - 1 == 0 -2.0945 0.9810 -2.135 0.0983 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)#mean diff tests and post-hoc for NS_SCR
lme_restingns_scr = lme(resting.ns_scr ~ factorvisit,
random = ~1|id,
correlation=corCompSymm(form=~1|id),
data=visitdataplus3,
na.action=na.exclude)
#summary(lme_restingns_scr)
anova(lme_restingns_scr)## numDF denDF F-value p-value
## (Intercept) 1 375 517.4040 <.0001
## factorvisit 2 375 81.7549 <.0001summary(glht(lme_restingns_scr, linfct=mcp(factorvisit = "Tukey")), test = adjusted(type = "bonferroni"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lme.formula(fixed = resting.ns_scr ~ factorvisit, data = visitdataplus3,
## random = ~1 | id, correlation = corCompSymm(form = ~1 | id),
## na.action = na.exclude)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 - 0 == 0 -2.1309 0.2543 -8.379 < 2e-16 ***
## 2 - 0 == 0 -3.2098 0.2572 -12.478 < 2e-16 ***
## 2 - 1 == 0 -1.0790 0.2622 -4.114 0.000117 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)# #mean diff tests and post-hoc for EMO
# lme_emo = lme(emo ~ factorvisit,
# random = ~1|id,
# correlation=corCompSymm(form=~1|id),
# data=visitdataplus3,
# na.action=na.exclude)
# #summary(lme_restingns_scr)
# anova(lme_emo)
# summary(glht(lme_emo, linfct=mcp(factorvisit = "Tukey")), test = adjusted(type = "bonferroni"))
#Checking mean differences across emotion conditions
#setting reference category for emotion variable = baseline
epochdataplus$bigsegmentemo <- relevel(epochdataplus$bigsegmentemo, ref='baseline')
str(epochdataplus$bigsegmentemo)## Factor w/ 5 levels "baseline","anger",..: 1 1 1 1 3 3 3 3 3 3 ...#testing big segment differences
model_emodiffs <- lme(rsa ~ 1 + bigsegmentemo,
random = ~ 1 |id,
data = epochdataplus, #[which(epochdataplus$visit == 0),],
na.action = na.exclude)
summary(model_emodiffs)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 64317.76 64374.58 -32151.88
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.071096 0.8574493
##
## Fixed effects: rsa ~ 1 + bigsegmentemo
## Value Std.Error DF t-value p-value
## (Intercept) 7.310861 0.06001833 24455 121.81047 0.0000
## bigsegmentemoanger -0.090479 0.01785835 24455 -5.06648 0.0000
## bigsegmentemofear -0.036414 0.01716791 24455 -2.12105 0.0339
## bigsegmentemohappy 0.032133 0.01716695 24455 1.87182 0.0612
## bigsegmentemosad -0.005436 0.01796456 24455 -0.30258 0.7622
## Correlation:
## (Intr) bgsgmntmn bgsgmntmf bgsgmntmh
## bigsegmentemoanger -0.142
## bigsegmentemofear -0.148 0.521
## bigsegmentemohappy -0.148 0.522 0.543
## bigsegmentemosad -0.141 0.498 0.518 0.518
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.54735060 -0.56607333 0.05485858 0.61708607 6.17207777
##
## Number of Observations: 24793
## Number of Groups: 334model_emodiffs <- lme(ns_scr ~ 1 + bigsegmentemo,
random = ~ 1 |id,
data = epochdataplus, #[which(epochdataplus$visit == 0),],
na.action = na.exclude)
summary(model_emodiffs)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 107346.3 107402.5 -53666.13
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.438184 2.471874
##
## Fixed effects: ns_scr ~ 1 + bigsegmentemo
## Value Std.Error DF t-value p-value
## (Intercept) 4.411359 0.08842162 22538 49.89005 0
## bigsegmentemoanger -1.396161 0.05346847 22538 -26.11185 0
## bigsegmentemofear -1.805515 0.05137444 22538 -35.14422 0
## bigsegmentemohappy -2.432590 0.05133630 22538 -47.38538 0
## bigsegmentemosad -2.462837 0.05445161 22538 -45.22983 0
## Correlation:
## (Intr) bgsgmntmn bgsgmntmf bgsgmntmh
## bigsegmentemoanger -0.303
## bigsegmentemofear -0.315 0.526
## bigsegmentemohappy -0.315 0.527 0.548
## bigsegmentemosad -0.297 0.496 0.517 0.517
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.5846879 -0.6170616 -0.1431053 0.4819266 10.3683011
##
## Number of Observations: 22870
## Number of Groups: 328Correlations
#reshaping visit-level data to wide for easy calculations of person-level correlations
restinglong <- visitdataplus3[,c("id","visit","resting.rsa","resting.pep","resting.ns_scr")]
restinglong <- restinglong[order(restinglong$id,restinglong$visit),]
restingwide <- reshape(data=restinglong,
timevar=c("visit"),
idvar=c("id"),
v.names=c("resting.rsa","resting.pep","resting.ns_scr"),
direction="wide", sep="_")
#resting RSA correlations
print(corr.test(restingwide[,c("resting.rsa_0","resting.rsa_1","resting.rsa_2")], use = "pairwise",method="pearson",adjust="holm", alpha=.05,ci=TRUE), short=FALSE)## Call:corr.test(x = restingwide[, c("resting.rsa_0", "resting.rsa_1",
## "resting.rsa_2")], use = "pairwise", method = "pearson",
## adjust = "holm", alpha = 0.05, ci = TRUE)
## Correlation matrix
## resting.rsa_0 resting.rsa_1 resting.rsa_2
## resting.rsa_0 1.00 0.59 0.47
## resting.rsa_1 0.59 1.00 0.50
## resting.rsa_2 0.47 0.50 1.00
## Sample Size
## resting.rsa_0 resting.rsa_1 resting.rsa_2
## resting.rsa_0 286 220 207
## resting.rsa_1 220 262 218
## resting.rsa_2 207 218 241
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## resting.rsa_0 resting.rsa_1 resting.rsa_2
## resting.rsa_0 0 0 0
## resting.rsa_1 0 0 0
## resting.rsa_2 0 0 0
##
## To see confidence intervals of the correlations, print with the short=FALSE option
##
## Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
## raw.lower raw.r raw.upper raw.p lower.adj upper.adj
## rs._0-rs._1 0.50 0.59 0.67 0 0.47 0.69
## rs._0-rs._2 0.36 0.47 0.57 0 0.36 0.57
## rs._1-rs._2 0.39 0.50 0.59 0 0.37 0.60#resting PEP correlations
print(corr.test(restingwide[,c("resting.pep_0","resting.pep_1","resting.pep_2")], use = "pairwise",method="pearson",adjust="holm", alpha=.05,ci=TRUE), short=FALSE)## Call:corr.test(x = restingwide[, c("resting.pep_0", "resting.pep_1",
## "resting.pep_2")], use = "pairwise", method = "pearson",
## adjust = "holm", alpha = 0.05, ci = TRUE)
## Correlation matrix
## resting.pep_0 resting.pep_1 resting.pep_2
## resting.pep_0 1.00 0.54 0.40
## resting.pep_1 0.54 1.00 0.64
## resting.pep_2 0.40 0.64 1.00
## Sample Size
## resting.pep_0 resting.pep_1 resting.pep_2
## resting.pep_0 272 204 194
## resting.pep_1 204 254 205
## resting.pep_2 194 205 232
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## resting.pep_0 resting.pep_1 resting.pep_2
## resting.pep_0 0 0 0
## resting.pep_1 0 0 0
## resting.pep_2 0 0 0
##
## To see confidence intervals of the correlations, print with the short=FALSE option
##
## Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
## raw.lower raw.r raw.upper raw.p lower.adj upper.adj
## rs._0-rs._1 0.43 0.54 0.63 0 0.42 0.64
## rs._0-rs._2 0.27 0.40 0.51 0 0.27 0.51
## rs._1-rs._2 0.55 0.64 0.71 0 0.53 0.73#resting NS_SCR correlations
print(corr.test(restingwide[,c("resting.ns_scr_0","resting.ns_scr_1","resting.ns_scr_2")], use = "pairwise",method="pearson",adjust="holm", alpha=.05,ci=TRUE), short=FALSE)## Call:corr.test(x = restingwide[, c("resting.ns_scr_0", "resting.ns_scr_1",
## "resting.ns_scr_2")], use = "pairwise", method = "pearson",
## adjust = "holm", alpha = 0.05, ci = TRUE)
## Correlation matrix
## resting.ns_scr_0 resting.ns_scr_1 resting.ns_scr_2
## resting.ns_scr_0 1.00 0.27 0.29
## resting.ns_scr_1 0.27 1.00 0.26
## resting.ns_scr_2 0.29 0.26 1.00
## Sample Size
## resting.ns_scr_0 resting.ns_scr_1 resting.ns_scr_2
## resting.ns_scr_0 254 172 163
## resting.ns_scr_1 172 225 170
## resting.ns_scr_2 163 170 218
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## resting.ns_scr_0 resting.ns_scr_1 resting.ns_scr_2
## resting.ns_scr_0 0 0 0
## resting.ns_scr_1 0 0 0
## resting.ns_scr_2 0 0 0
##
## To see confidence intervals of the correlations, print with the short=FALSE option
##
## Confidence intervals based upon normal theory. To get bootstrapped values, try cor.ci
## raw.lower raw.r raw.upper raw.p lower.adj upper.adj
## r.__0-r.__1 0.13 0.27 0.40 0 0.11 0.42
## r.__0-r.__2 0.14 0.29 0.42 0 0.11 0.45
## r.__1-r.__2 0.11 0.26 0.39 0 0.11 0.39# #EMO correlations
# print(corr.test(visitwide[,c("emo_0","emo_1","emo_2")], use = "pairwise",method="pearson",adjust="holm", alpha=.05,ci=TRUE), short=FALSE)C. Multilevel models for variance decompositions of each measure
RSA:Calculating some variance breakdowns for epoch, visit, person
model0 <- lme(rsa ~ 1,
random = ~ 1|id/visit,
data = epochdataplus,
na.action = na.exclude)
summary(model0)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 47155.55 47188.03 -23573.78
##
## Random effects:
## Formula: ~1 | id
## (Intercept)
## StdDev: 0.9138973
##
## Formula: ~1 | visit %in% id
## (Intercept) Residual
## StdDev: 0.8083329 0.5799978
##
## Fixed effects: rsa ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 7.290404 0.05826116 23963 125.1332 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.42777164 -0.56098303 0.01005193 0.56404670 9.22775377
##
## Number of Observations: 24793
## Number of Groups:
## id visit %in% id
## 334 830VarCorr(model0)## Variance StdDev
## id = pdLogChol(1)
## (Intercept) 0.8352082 0.9138973
## visit = pdLogChol(1)
## (Intercept) 0.6534021 0.8083329
## Residual 0.3363974 0.5799978totalvar <- as.numeric(VarCorr(model0)[2]) + as.numeric(VarCorr(model0)[4]) + as.numeric(VarCorr(model0)[5])
totalvar## [1] 1.825008#percentages
round(as.numeric(VarCorr(model0)[2])/totalvar,2)## [1] 0.46round(as.numeric(VarCorr(model0)[4])/totalvar,2)## [1] 0.36round(as.numeric(VarCorr(model0)[5])/totalvar,2)## [1] 0.18#calculating ICC
totalvar <- (0.9138973)^2 + (0.8083329)^2 + (0.5799978)^2
(0.9138973)^2/totalvar + (0.8083329)^2/totalvar + (0.5799978)^2/totalvar## [1] 1RSA: 46% personlevel, 36% visitlevel, 18% epochlevel
pct_rsa <- c(46, 36, 18)
pie(pct_rsa, main="RSA", col=rainbow(length(pct_rsa)),
labels=c("Persons 46%","Visits 36%","Epochs 18%"))
PEP: Calculating some variance breakdowns for epoch, visit, person
model0 <- lme(pep ~ 1,
random = ~ 1|id/visit,
data = epochdataplus,
na.action = na.exclude)
summary(model0)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 139027 139059.3 -69509.52
##
## Random effects:
## Formula: ~1 | id
## (Intercept)
## StdDev: 11.57705
##
## Formula: ~1 | visit %in% id
## (Intercept) Residual
## StdDev: 9.791796 4.097968
##
## Fixed effects: pep ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 104.3605 0.7353749 22954 141.9147 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.092861e+01 -4.506675e-01 -3.197161e-04 4.564533e-01 1.792600e+01
##
## Number of Observations: 23760
## Number of Groups:
## id visit %in% id
## 329 806VarCorr(model0)## Variance StdDev
## id = pdLogChol(1)
## (Intercept) 134.02818 11.577054
## visit = pdLogChol(1)
## (Intercept) 95.87928 9.791796
## Residual 16.79334 4.097968totalvar <- as.numeric(VarCorr(model0)[2]) + as.numeric(VarCorr(model0)[4]) + as.numeric(VarCorr(model0)[5])
totalvar## [1] 246.7008#percentages
round(as.numeric(VarCorr(model0)[2])/totalvar,2)## [1] 0.54round(as.numeric(VarCorr(model0)[4])/totalvar,2)## [1] 0.39round(as.numeric(VarCorr(model0)[5])/totalvar,2)## [1] 0.07#calculating ICC
totalvar <- (11.577054)^2 + (9.791796)^2 + (4.097968)^2
(11.577054)^2/totalvar + (9.791796)^2/totalvar + (4.097968)^2/totalvar## [1] 1PEP: 54% personlevel, 39% visitlevel, 7% epochlevel
pct_pep <- c(54, 39, 7)
pie(pct_pep, main="PEP", col=rainbow(length(pct_pep)),
labels=c("Persons 54%","Visits 39%","Epochs 7%"))
EDA: Calculating some variance breakdowns for epoch, visit, person
model0 <- lme(ns_scr ~ 1,
random = ~ 1|id/visit,
data = epochdataplus,
na.action = na.exclude)
summary(model0)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 106520.5 106552.7 -53256.27
##
## Random effects:
## Formula: ~1 | id
## (Intercept)
## StdDev: 0.9298363
##
## Formula: ~1 | visit %in% id
## (Intercept) Residual
## StdDev: 1.511979 2.368507
##
## Fixed effects: ns_scr ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 2.694398 0.07744505 22053 34.79109 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.7176994 -0.5761958 -0.1501375 0.3873147 9.9916404
##
## Number of Observations: 22870
## Number of Groups:
## id visit %in% id
## 328 817VarCorr(model0)## Variance StdDev
## id = pdLogChol(1)
## (Intercept) 0.8645955 0.9298363
## visit = pdLogChol(1)
## (Intercept) 2.2860813 1.5119793
## Residual 5.6098246 2.3685068totalvar <- as.numeric(VarCorr(model0)[2]) + as.numeric(VarCorr(model0)[4]) + as.numeric(VarCorr(model0)[5])
totalvar## [1] 8.760501#percentages
round(as.numeric(VarCorr(model0)[2])/totalvar,2)## [1] 0.1round(as.numeric(VarCorr(model0)[4])/totalvar,2)## [1] 0.26round(as.numeric(VarCorr(model0)[5])/totalvar,2)## [1] 0.64#calculating ICC
totalvar <- (0.9298363)^2 + (1.5119793)^2 + (2.3685068)^2
(0.9298363)^2/totalvar + (1.5119793)^2/totalvar + (2.3685068)^2/totalvar## [1] 1EDA: 10% personlevel, 26% visitlevel, 64% epochlevel
pct_eda <- c(10, 26, 64)
pie(pct_eda, main="NS_SCR", col=rainbow(length(pct_eda)),
labels=c("Persons 10%","Visits 26%","Epochs 64%"))
D. Multilevel models for variance decompositions of bivariate couplings
Looking at breakdown of variance in coupling for person and visit RSA ~ PEP Coupling
#coupling model
model1 <- lme(rsa ~ 1 + pep_d,
random = ~ 1 + pep_d |id/visit,
data = epochdataplus,
na.action = na.exclude)
summary(model1)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 43920.29 43992.77 -21951.15
##
## Random effects:
## Formula: ~1 + pep_d | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.91232060 (Intr)
## pep_d 0.01594729 0.118
##
## Formula: ~1 + pep_d | visit %in% id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.81318858 (Intr)
## pep_d 0.02540923 -0.176
## Residual 0.56890008
##
## Fixed effects: rsa ~ 1 + pep_d
## Value Std.Error DF t-value p-value
## (Intercept) 7.290434 0.05884653 22412 123.88892 0
## pep_d 0.014944 0.00178459 22412 8.37392 0
## Correlation:
## (Intr)
## pep_d 0.008
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.393993270 -0.560522452 0.008195362 0.558202079 9.296836814
##
## Number of Observations: 23218
## Number of Groups:
## id visit %in% id
## 329 805VarCorr(model1)## Variance StdDev Corr
## id = pdLogChol(1 + pep_d)
## (Intercept) 0.8323288838 0.91232060 (Intr)
## pep_d 0.0002543160 0.01594729 0.118
## visit = pdLogChol(1 + pep_d)
## (Intercept) 0.6612756711 0.81318858 (Intr)
## pep_d 0.0006456291 0.02540923 -0.176
## Residual 0.3236473061 0.56890008totalvar <- as.numeric(VarCorr(model1)[3]) + as.numeric(VarCorr(model1)[6])
totalvar## [1] 0.0008999451#percentages
round(as.numeric(VarCorr(model1)[3])/totalvar,2)## [1] 0.28round(as.numeric(VarCorr(model1)[6])/totalvar,2)## [1] 0.72couplingvar <- (0.01594729)^2 + (0.02540923)^2
(0.01594729)^2/couplingvar + (0.02540923)^2/couplingvar## [1] 1RSA ~ PEP: 28% personlevel, 72% visitlevel
pct_rsapep <- c(28, 72)
pie(pct_rsapep, main="RSA ~ PEP", col=rainbow(length(pct_rsapep)),
labels=c("Persons 28%","Visits 72%"))
Looking at breakdown of variance in coupling for person and visit RSA ~ EDA Coupling
#coupling model
model1 <- lme(rsa ~ 1 + ns_scr_d,
random = ~ 1 + ns_scr_d |id/visit,
data = epochdataplus,
na.action = na.exclude)
summary(model1)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 42103.83 42175.89 -21042.92
##
## Random effects:
## Formula: ~1 + ns_scr_d | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.91906910 (Intr)
## ns_scr_d 0.01361602 -0.183
##
## Formula: ~1 + ns_scr_d | visit %in% id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.81461819 (Intr)
## ns_scr_d 0.03462067 -0.369
## Residual 0.57259106
##
## Fixed effects: rsa ~ 1 + ns_scr_d
## Value Std.Error DF t-value p-value
## (Intercept) 7.286400 0.05944467 21371 122.57449 0
## ns_scr_d -0.012041 0.00239039 21371 -5.03707 0
## Correlation:
## (Intr)
## ns_scr_d -0.147
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.800086874 -0.559780948 0.009054321 0.561311196 8.966249173
##
## Number of Observations: 22168
## Number of Groups:
## id visit %in% id
## 328 796VarCorr(model1)## Variance StdDev Corr
## id = pdLogChol(1 + ns_scr_d)
## (Intercept) 0.844688018 0.91906910 (Intr)
## ns_scr_d 0.000185396 0.01361602 -0.183
## visit = pdLogChol(1 + ns_scr_d)
## (Intercept) 0.663602791 0.81461819 (Intr)
## ns_scr_d 0.001198590 0.03462067 -0.369
## Residual 0.327860517 0.57259106totalvar <- as.numeric(VarCorr(model1)[3]) + as.numeric(VarCorr(model1)[6])
totalvar## [1] 0.001383986#percentages
round(as.numeric(VarCorr(model1)[3])/totalvar,2)## [1] 0.13round(as.numeric(VarCorr(model1)[6])/totalvar,2)## [1] 0.87couplingvar <- (0.01361602)^2 + (0.03462067)^2
(0.01361602)^2/couplingvar## [1] 0.1339579(0.03462067)^2/couplingvar## [1] 0.8660421RSA ~ EDA: 13% person, 87% visit
pct_rsaeda <- c(13, 87)
pie(pct_rsaeda, main="RSA ~ NS_SCR", col=rainbow(length(pct_rsaeda)),
labels=c("Persons 13%","Visits 87%"))
Looking at breakdown of variance in coupling for person and visit PEP ~ EDA Coupling
#coupling model
model1 <- lme(pep ~ 1 + ns_scr_d,
random = ~ 1 + ns_scr_d |id/visit,
data = epochdataplus,
na.action = na.exclude)
summary(model1)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 123696.8 123768.4 -61839.39
##
## Random effects:
## Formula: ~1 + ns_scr_d | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 11.5753372 (Intr)
## ns_scr_d 0.1042482 0.116
##
## Formula: ~1 + ns_scr_d | visit %in% id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 9.9589284 (Intr)
## ns_scr_d 0.4926138 -0.122
## Residual 3.9726131
##
## Fixed effects: pep ~ 1 + ns_scr_d
## Value Std.Error DF t-value p-value
## (Intercept) 104.42139 0.7487536 20422 139.46029 0.0000
## ns_scr_d -0.06654 0.0245596 20422 -2.70937 0.0067
## Correlation:
## (Intr)
## ns_scr_d -0.019
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -10.154923043 -0.448122427 -0.002462409 0.454518396 18.320224554
##
## Number of Observations: 21189
## Number of Groups:
## id visit %in% id
## 323 766VarCorr(model1)## Variance StdDev Corr
## id = pdLogChol(1 + ns_scr_d)
## (Intercept) 133.98843061 11.5753372 (Intr)
## ns_scr_d 0.01086768 0.1042482 0.116
## visit = pdLogChol(1 + ns_scr_d)
## (Intercept) 99.18025538 9.9589284 (Intr)
## ns_scr_d 0.24266838 0.4926138 -0.122
## Residual 15.78165448 3.9726131totalvar <- as.numeric(VarCorr(model1)[3]) + as.numeric(VarCorr(model1)[6])
totalvar## [1] 0.2535361#percentages
round(as.numeric(VarCorr(model1)[3])/totalvar,2)## [1] 0.04round(as.numeric(VarCorr(model1)[6])/totalvar,2)## [1] 0.96couplingvar <- (0.1042482)^2 + (0.4926138)^2
(0.1042482)^2/couplingvar## [1] 0.04286447(0.4926138)^2/couplingvar## [1] 0.9571355PEP ~ EDA: 3% person, 97% visit
pct_pepeda <- c(4, 96)
pie(pct_pepeda, main="PEP ~ NS_SCR", col=rainbow(length(pct_pepeda)),
labels=c("Persons 4%","Visits 96%"))
EMO: Looking at breakdown of variance in EMO for person and visit This one removed from paper
# #coupling model
# model1 <- lme(emo ~ 1 ,
# random = ~ 1 |id,
# data = visitdataplus3,
# na.action = na.exclude)
# summary(model1)
# VarCorr(model1)
# totalvar <- as.numeric(VarCorr(model1)[1]) + as.numeric(VarCorr(model1)[2])
# totalvar
#
# #percentages
# round(as.numeric(VarCorr(model1)[1])/totalvar,2)
# #54% person
# round(as.numeric(VarCorr(model1)[2])/totalvar,2)
# #46% visitComputing within-person correlations
# icorr.visit <- ddply(epochdataplus, c("id","visit"), summarize,
# icorr.nsscr_rsa=cor(x=rsa_d,y=ns_scr_d, #correlation
# use="pairwise.complete.obs", method="pearson"))
# describe(icorr.visit[which(icorr.visit$visit==0),])
# visitdataplus3 <- merge(visitdataplus2,icorr.visit,by=c("id","visit"),all=TRUE)
#
# ggplot(data=subset(visitdataplus3,visit==0), aes(x=emo, y=icorr.nsscr_rsa, group=id, color="gray"), legend=FALSE) +
# geom_point(aes(color=factor(id))) + guides(color=FALSE) +
# #geom_line(aes(color=factor(id))) +
# #geom_smooth(method=lm, se=FALSE, fullrange=FALSE, lty=1, size=.5, color="gray40") +
# geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="red") +
# ylab("Coupling") + scale_x_continuous(breaks=seq(1,6,by=1)) +
# xlab("Emotion Regulation (Kinder)") + #ylim(1,7) +
# ggtitle("Predicted Coupling NS_SCR ~ EDA")
#
# ggplot(data=icorr.visit, aes(x=visit,y=icorr.nsscr_rsa, group=id, color="gray"), legend=FALSE) +
# geom_point(aes(color=factor(id))) + guides(color=FALSE) +
# geom_line(aes(color=factor(id))) +
# #geom_smooth(method=lm, se=FALSE, fullrange=FALSE, lty=1, size=.5, color="gray40") +
# geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="red") +
# xlab("Visit") + #scale_x_continuous(breaks=seq(0,1,by=1)) +
# ylab("iCORR") + #ylim(1,7) +
# ggtitle("Coupling overtime")E. Multilevel model examining situational differences in coupling
Setting up emotion variable
#Coding the kind of emotion = approach or avoid
epochdataplus$approach <- ifelse(epochdataplus$epoch>=17,1,0)
epochdataplus$avoid <- ifelse(epochdataplus$epoch>=17,0,1)Setting up multilevel model to examine situational predictor We check avoid (vs approach emotions) and visit
#NEW MODEL
model_situation <- lme(rsa_d ~ 0 + ns_scr_d + avoid + visit +
ns_scr_d:avoid + ns_scr_d:visit, # + ns_scr_d:visit:avoid,
random = ~ 0 + ns_scr_d |id,
data = epochdataplus, #[which(epochdataplus$visit == 0),],
na.action = na.exclude)
summary(model_situation)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 37923.48 37979.53 -18954.74
##
## Random effects:
## Formula: ~0 + ns_scr_d | id
## ns_scr_d Residual
## StdDev: 0.02784436 0.5659631
##
## Fixed effects: rsa_d ~ 0 + ns_scr_d + avoid + visit + ns_scr_d:avoid + ns_scr_d:visit
## Value Std.Error DF t-value p-value
## ns_scr_d -0.014694132 0.003248845 21836 -4.522879 0.0000
## avoid 0.002909575 0.006259582 21836 0.464819 0.6421
## visit -0.000359274 0.003655251 21836 -0.098290 0.9217
## ns_scr_d:avoid 0.008808164 0.003325794 21836 2.648439 0.0081
## ns_scr_d:visit -0.001498397 0.002216718 21836 -0.675953 0.4991
## Correlation:
## ns_sc_ avoid visit ns_scr_d:vd
## avoid 0.001
## visit -0.014 -0.538
## ns_scr_d:avoid -0.534 0.015 0.017
## ns_scr_d:visit -0.441 0.017 -0.009 0.029
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.810817587 -0.573473709 0.005484974 0.568662965 9.157402972
##
## Number of Observations: 22168
## Number of Groups: 328#extra info on results
VarCorr(model_situation)## id = pdLogChol(0 + ns_scr_d)
## Variance StdDev
## ns_scr_d 0.0007753084 0.02784436
## Residual 0.3203142617 0.56596313intervals(model_situation)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## ns_scr_d -0.021062104 -0.0146941322 -0.008326161
## avoid -0.009359661 0.0029095748 0.015178811
## visit -0.007523831 -0.0003592743 0.006805282
## ns_scr_d:avoid 0.002289366 0.0088081637 0.015326962
## ns_scr_d:visit -0.005843326 -0.0014983971 0.002846531
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd(ns_scr_d) 0.02341928 0.02784436 0.03310556
##
## Within-group standard error:
## lower est. upper
## 0.5606841 0.5659631 0.5712918Note the removal of the intercept (outcome variable has been person and visit centered).
Parallel model for RSA~PEP
model_situation.pep <- lme(rsa_d ~ 0 + pep_d + avoid + visit +
pep_d:avoid + pep_d:visit, # + ns_scr_d:visit:avoid,
random = ~ 0 + pep_d |id,
data = epochdataplus, #[which(epochdataplus$visit == 0),],
na.action = na.exclude)
summary(model_situation.pep)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 39584.97 39641.34 -19785.49
##
## Random effects:
## Formula: ~0 + pep_d | id
## pep_d Residual
## StdDev: 0.02226021 0.5633206
##
## Fixed effects: rsa_d ~ 0 + pep_d + avoid + visit + pep_d:avoid + pep_d:visit
## Value Std.Error DF t-value p-value
## pep_d 0.012543200 0.002253487 22885 5.566130 0.0000
## avoid 0.007245408 0.006154891 22885 1.177179 0.2391
## visit -0.002500105 0.003545220 22885 -0.705204 0.4807
## pep_d:avoid 0.000295086 0.001917337 22885 0.153904 0.8777
## pep_d:visit 0.001014887 0.001299636 22885 0.780901 0.4349
## Correlation:
## pep_d avoid visit pp_d:vd
## avoid 0.002
## visit -0.012 -0.557
## pep_d:avoid -0.460 0.009 0.017
## pep_d:visit -0.471 0.007 -0.004 0.001
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.50887070 -0.57190930 0.00456113 0.56342391 9.56710159
##
## Number of Observations: 23218
## Number of Groups: 329#extra info on results
VarCorr(model_situation)## id = pdLogChol(0 + ns_scr_d)
## Variance StdDev
## ns_scr_d 0.0007753084 0.02784436
## Residual 0.3203142617 0.56596313Follow-ups with specfic emotions (different reference categories)
#Checking for particular emotion conditions
#setting reference category for emotion variable = baseline
epochdataplus$bigsegmentemo <- relevel(epochdataplus$bigsegmentemo, ref='baseline')
str(epochdataplus$bigsegmentemo)## Factor w/ 5 levels "baseline","anger",..: 1 1 1 1 3 3 3 3 3 3 ...#trying other emotion distinction
model_situation.emo <- lme(rsa_d ~ 0 + ns_scr_d + bigsegmentemo + visit +
ns_scr_d:bigsegmentemo + ns_scr_d:visit,
# + ns_scr_d:visit:bigsegmentemo,
random = ~ 0 + ns_scr_d |id,
data = epochdataplus, #[which(epochdataplus$visit == 0),],
na.action = na.exclude)
#Higher order interaction trimmed
summary(model_situation.emo)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 37869 37981.08 -18920.5
##
## Random effects:
## Formula: ~0 + ns_scr_d | id
## ns_scr_d Residual
## StdDev: 0.02748625 0.564527
##
## Fixed effects: rsa_d ~ 0 + ns_scr_d + bigsegmentemo + visit + ns_scr_d:bigsegmentemo + ns_scr_d:visit
## Value Std.Error DF t-value
## ns_scr_d -0.00941642 0.003409489 21829 -2.761828
## bigsegmentemobaseline 0.03417340 0.010847658 21829 3.150302
## bigsegmentemoanger -0.05780222 0.009824361 21829 -5.883561
## bigsegmentemofear -0.02375888 0.009145081 21829 -2.597995
## bigsegmentemohappy 0.05264876 0.009993085 21829 5.268520
## bigsegmentemosad 0.01672138 0.011108549 21829 1.505271
## visit -0.00119378 0.004662715 21829 -0.256027
## ns_scr_d:bigsegmentemoanger -0.01724256 0.005500497 21829 -3.134728
## ns_scr_d:bigsegmentemofear -0.00253543 0.005216562 21829 -0.486035
## ns_scr_d:bigsegmentemohappy 0.00203644 0.005640710 21829 0.361026
## ns_scr_d:bigsegmentemosad 0.01109739 0.006074090 21829 1.827005
## ns_scr_d:visit 0.00028206 0.002219086 21829 0.127105
## p-value
## ns_scr_d 0.0058
## bigsegmentemobaseline 0.0016
## bigsegmentemoanger 0.0000
## bigsegmentemofear 0.0094
## bigsegmentemohappy 0.0000
## bigsegmentemosad 0.1323
## visit 0.7979
## ns_scr_d:bigsegmentemoanger 0.0017
## ns_scr_d:bigsegmentemofear 0.6269
## ns_scr_d:bigsegmentemohappy 0.7181
## ns_scr_d:bigsegmentemosad 0.0677
## ns_scr_d:visit 0.8989
## Correlation:
## ns_sc_ bgsgmntmb bgsgmntmn bgsgmntmf bgsgmntmh
## bigsegmentemobaseline -0.342
## bigsegmentemoanger -0.034 0.195
## bigsegmentemofear -0.038 0.209 0.204
## bigsegmentemohappy -0.046 0.202 0.195 0.210
## bigsegmentemosad -0.038 0.179 0.172 0.187 0.181
## visit 0.089 -0.452 -0.435 -0.469 -0.451
## ns_scr_d:bigsegmentemoanger -0.364 0.214 -0.066 0.045 0.036
## ns_scr_d:bigsegmentemofear -0.394 0.233 0.048 0.119 0.046
## ns_scr_d:bigsegmentemohappy -0.368 0.221 0.047 0.052 0.391
## ns_scr_d:bigsegmentemosad -0.340 0.206 0.044 0.049 0.047
## ns_scr_d:visit -0.383 0.022 -0.031 -0.025 0.014
## bgsgmntms visit ns_scr_d:bgsgmntmn
## bigsegmentemobaseline
## bigsegmentemoanger
## bigsegmentemofear
## bigsegmentemohappy
## bigsegmentemosad
## visit -0.400
## ns_scr_d:bigsegmentemoanger 0.031 -0.082
## ns_scr_d:bigsegmentemofear 0.039 -0.100 0.262
## ns_scr_d:bigsegmentemohappy 0.042 -0.104 0.236
## ns_scr_d:bigsegmentemosad 0.399 -0.099 0.215
## ns_scr_d:visit 0.011 0.003 -0.066
## ns_scr_d:bgsgmntmf ns_scr_d:bgsgmntmh
## bigsegmentemobaseline
## bigsegmentemoanger
## bigsegmentemofear
## bigsegmentemohappy
## bigsegmentemosad
## visit
## ns_scr_d:bigsegmentemoanger
## ns_scr_d:bigsegmentemofear
## ns_scr_d:bigsegmentemohappy 0.251
## ns_scr_d:bigsegmentemosad 0.228 0.211
## ns_scr_d:visit -0.041 -0.010
## ns_scr_d:bgsgmntms
## bigsegmentemobaseline
## bigsegmentemoanger
## bigsegmentemofear
## bigsegmentemohappy
## bigsegmentemosad
## visit
## ns_scr_d:bigsegmentemoanger
## ns_scr_d:bigsegmentemofear
## ns_scr_d:bigsegmentemohappy
## ns_scr_d:bigsegmentemosad
## ns_scr_d:visit 0.001
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.757422266 -0.572517213 0.003633263 0.570220541 9.194947131
##
## Number of Observations: 22168
## Number of Groups: 328anova(model_situation.emo)## numDF denDF F-value p-value
## ns_scr_d 1 21829 22.452979 <.0001
## bigsegmentemo 5 21829 22.716654 <.0001
## visit 1 21829 0.069266 0.7924
## ns_scr_d:bigsegmentemo 4 21829 4.237849 0.0020
## ns_scr_d:visit 1 21829 0.016156 0.8989#setting reference category for emotion variable = fear
epochdataplus$bigsegmentemofear <- relevel(epochdataplus$bigsegmentemo, ref='fear')
str(epochdataplus$bigsegmentemofear)## Factor w/ 5 levels "fear","baseline",..: 2 2 2 2 1 1 1 1 1 1 ...#trying other emotion distinction
model_situation.emofear <- lme(rsa_d ~ 0 + ns_scr_d + bigsegmentemofear + visit +
ns_scr_d:bigsegmentemofear + ns_scr_d:visit,
# + ns_scr_d:visit:bigsegmentemofear,
random = ~ 0 + ns_scr_d |id,
data = epochdataplus, #[which(epochdataplus$visit == 0),],
na.action = na.exclude)
#Higher order interaction trimmed
summary(model_situation.emofear)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 37869 37981.08 -18920.5
##
## Random effects:
## Formula: ~0 + ns_scr_d | id
## ns_scr_d Residual
## StdDev: 0.02748623 0.564527
##
## Fixed effects: rsa_d ~ 0 + ns_scr_d + bigsegmentemofear + visit + ns_scr_d:bigsegmentemofear + ns_scr_d:visit
## Value Std.Error DF t-value
## ns_scr_d -0.01195185 0.004983822 21829 -2.398130
## bigsegmentemofearfear -0.02375888 0.009145081 21829 -2.597995
## bigsegmentemofearbaseline 0.03417340 0.010847658 21829 3.150302
## bigsegmentemofearanger -0.05780222 0.009824361 21829 -5.883561
## bigsegmentemofearhappy 0.05264876 0.009993085 21829 5.268520
## bigsegmentemofearsad 0.01672138 0.011108549 21829 1.505271
## visit -0.00119378 0.004662715 21829 -0.256027
## ns_scr_d:bigsegmentemofearbaseline 0.00253543 0.005216562 21829 0.486035
## ns_scr_d:bigsegmentemofearanger -0.01470713 0.006516160 21829 -2.257024
## ns_scr_d:bigsegmentemofearhappy 0.00457187 0.006651021 21829 0.687394
## ns_scr_d:bigsegmentemofearsad 0.01363282 0.007048694 21829 1.934092
## ns_scr_d:visit 0.00028206 0.002219086 21829 0.127104
## p-value
## ns_scr_d 0.0165
## bigsegmentemofearfear 0.0094
## bigsegmentemofearbaseline 0.0016
## bigsegmentemofearanger 0.0000
## bigsegmentemofearhappy 0.0000
## bigsegmentemofearsad 0.1323
## visit 0.7979
## ns_scr_d:bigsegmentemofearbaseline 0.6269
## ns_scr_d:bigsegmentemofearanger 0.0240
## ns_scr_d:bigsegmentemofearhappy 0.4918
## ns_scr_d:bigsegmentemofearsad 0.0531
## ns_scr_d:visit 0.8989
## Correlation:
## ns_sc_ bgsgmntmfrf bgsgmntmfrb
## bigsegmentemofearfear 0.098
## bigsegmentemofearbaseline 0.010 0.209
## bigsegmentemofearanger 0.026 0.204 0.195
## bigsegmentemofearhappy 0.017 0.210 0.202
## bigsegmentemofearsad 0.015 0.187 0.179
## visit -0.044 -0.469 -0.452
## ns_scr_d:bigsegmentemofearbaseline -0.777 -0.119 -0.233
## ns_scr_d:bigsegmentemofearanger -0.602 -0.057 -0.006
## ns_scr_d:bigsegmentemofearhappy -0.600 -0.048 0.005
## ns_scr_d:bigsegmentemofearsad -0.570 -0.045 0.004
## ns_scr_d:visit -0.305 -0.025 0.022
## bgsgmntmfrn bgsgmntmfrh bgsgmntmfrs
## bigsegmentemofearfear
## bigsegmentemofearbaseline
## bigsegmentemofearanger
## bigsegmentemofearhappy 0.195
## bigsegmentemofearsad 0.172 0.181
## visit -0.435 -0.451 -0.400
## ns_scr_d:bigsegmentemofearbaseline -0.048 -0.046 -0.039
## ns_scr_d:bigsegmentemofearanger -0.094 -0.006 -0.005
## ns_scr_d:bigsegmentemofearhappy 0.003 0.296 0.005
## ns_scr_d:bigsegmentemofearsad 0.003 0.006 0.315
## ns_scr_d:visit -0.031 0.014 0.011
## visit ns_scr_d:bgsgmntmfrb
## bigsegmentemofearfear
## bigsegmentemofearbaseline
## bigsegmentemofearanger
## bigsegmentemofearhappy
## bigsegmentemofearsad
## visit
## ns_scr_d:bigsegmentemofearbaseline 0.100
## ns_scr_d:bigsegmentemofearanger 0.011 0.580
## ns_scr_d:bigsegmentemofearhappy -0.010 0.571
## ns_scr_d:bigsegmentemofearsad -0.012 0.544
## ns_scr_d:visit 0.003 0.041
## ns_scr_d:bgsgmntmfrn
## bigsegmentemofearfear
## bigsegmentemofearbaseline
## bigsegmentemofearanger
## bigsegmentemofearhappy
## bigsegmentemofearsad
## visit
## ns_scr_d:bigsegmentemofearbaseline
## ns_scr_d:bigsegmentemofearanger
## ns_scr_d:bigsegmentemofearhappy 0.453
## ns_scr_d:bigsegmentemofearsad 0.428
## ns_scr_d:visit -0.022
## ns_scr_d:bgsgmntmfrh
## bigsegmentemofearfear
## bigsegmentemofearbaseline
## bigsegmentemofearanger
## bigsegmentemofearhappy
## bigsegmentemofearsad
## visit
## ns_scr_d:bigsegmentemofearbaseline
## ns_scr_d:bigsegmentemofearanger
## ns_scr_d:bigsegmentemofearhappy
## ns_scr_d:bigsegmentemofearsad 0.423
## ns_scr_d:visit 0.024
## ns_scr_d:bgsgmntmfrs
## bigsegmentemofearfear
## bigsegmentemofearbaseline
## bigsegmentemofearanger
## bigsegmentemofearhappy
## bigsegmentemofearsad
## visit
## ns_scr_d:bigsegmentemofearbaseline
## ns_scr_d:bigsegmentemofearanger
## ns_scr_d:bigsegmentemofearhappy
## ns_scr_d:bigsegmentemofearsad
## ns_scr_d:visit 0.031
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.757422106 -0.572517200 0.003633265 0.570220539 9.194947298
##
## Number of Observations: 22168
## Number of Groups: 328anova(model_situation.emofear)## numDF denDF F-value p-value
## ns_scr_d 1 21829 22.452992 <.0001
## bigsegmentemofear 5 21829 22.716655 <.0001
## visit 1 21829 0.069266 0.7924
## ns_scr_d:bigsegmentemofear 4 21829 4.237850 0.0020
## ns_scr_d:visit 1 21829 0.016156 0.8989#setting reference category for emotion variable = anger
epochdataplus$bigsegmentemoanger <- relevel(epochdataplus$bigsegmentemo, ref='anger')
str(epochdataplus$bigsegmentemoanger)## Factor w/ 5 levels "anger","baseline",..: 2 2 2 2 3 3 3 3 3 3 ...#trying other emotion distinction
model_situation.emoanger <- lme(rsa_d ~ 0 + ns_scr_d + bigsegmentemoanger + visit +
ns_scr_d:bigsegmentemoanger + ns_scr_d:visit,
random = ~ 0 + ns_scr_d |id,
data = epochdataplus,
na.action = na.exclude)
summary(model_situation.emoanger)## Linear mixed-effects model fit by REML
## Data: epochdataplus
## AIC BIC logLik
## 37869 37981.08 -18920.5
##
## Random effects:
## Formula: ~0 + ns_scr_d | id
## ns_scr_d Residual
## StdDev: 0.02748626 0.564527
##
## Fixed effects: rsa_d ~ 0 + ns_scr_d + bigsegmentemoanger + visit + ns_scr_d:bigsegmentemoanger + ns_scr_d:visit
## Value Std.Error DF
## ns_scr_d -0.02665898 0.005312232 21829
## bigsegmentemoangeranger -0.05780222 0.009824361 21829
## bigsegmentemoangerbaseline 0.03417340 0.010847658 21829
## bigsegmentemoangerfear -0.02375888 0.009145081 21829
## bigsegmentemoangerhappy 0.05264876 0.009993085 21829
## bigsegmentemoangersad 0.01672138 0.011108549 21829
## visit -0.00119378 0.004662715 21829
## ns_scr_d:bigsegmentemoangerbaseline 0.01724256 0.005500497 21829
## ns_scr_d:bigsegmentemoangerfear 0.01470713 0.006516160 21829
## ns_scr_d:bigsegmentemoangerhappy 0.01927901 0.006885656 21829
## ns_scr_d:bigsegmentemoangersad 0.02833995 0.007265835 21829
## ns_scr_d:visit 0.00028206 0.002219086 21829
## t-value p-value
## ns_scr_d -5.018415 0.0000
## bigsegmentemoangeranger -5.883561 0.0000
## bigsegmentemoangerbaseline 3.150302 0.0016
## bigsegmentemoangerfear -2.597995 0.0094
## bigsegmentemoangerhappy 5.268520 0.0000
## bigsegmentemoangersad 1.505271 0.1323
## visit -0.256027 0.7979
## ns_scr_d:bigsegmentemoangerbaseline 3.134728 0.0017
## ns_scr_d:bigsegmentemoangerfear 2.257024 0.0240
## ns_scr_d:bigsegmentemoangerhappy 2.799879 0.0051
## ns_scr_d:bigsegmentemoangersad 3.900440 0.0001
## ns_scr_d:visit 0.127105 0.8989
## Correlation:
## ns_sc_ bgsgmntmngrn bgsgmntmngrb
## bigsegmentemoangeranger -0.090
## bigsegmentemoangerbaseline 0.002 0.195
## bigsegmentemoangerfear 0.022 0.204 0.209
## bigsegmentemoangerhappy 0.008 0.195 0.202
## bigsegmentemoangersad 0.008 0.172 0.179
## visit -0.028 -0.435 -0.452
## ns_scr_d:bigsegmentemoangerbaseline -0.802 0.066 -0.214
## ns_scr_d:bigsegmentemoangerfear -0.662 0.094 0.006
## ns_scr_d:bigsegmentemoangerhappy -0.633 0.092 0.010
## ns_scr_d:bigsegmentemoangersad -0.603 0.087 0.010
## ns_scr_d:visit -0.314 -0.031 0.022
## bgsgmntmngrf bgsgmntmngrh bgsgmntmngrs
## bigsegmentemoangeranger
## bigsegmentemoangerbaseline
## bigsegmentemoangerfear
## bigsegmentemoangerhappy 0.210
## bigsegmentemoangersad 0.187 0.181
## visit -0.469 -0.451 -0.400
## ns_scr_d:bigsegmentemoangerbaseline -0.045 -0.036 -0.031
## ns_scr_d:bigsegmentemoangerfear 0.057 0.006 0.005
## ns_scr_d:bigsegmentemoangerhappy 0.007 0.291 0.010
## ns_scr_d:bigsegmentemoangersad 0.007 0.011 0.310
## ns_scr_d:visit -0.025 0.014 0.011
## visit ns_scr_d:bgsgmntmngrb
## bigsegmentemoangeranger
## bigsegmentemoangerbaseline
## bigsegmentemoangerfear
## bigsegmentemoangerhappy
## bigsegmentemoangersad
## visit
## ns_scr_d:bigsegmentemoangerbaseline 0.082
## ns_scr_d:bigsegmentemoangerfear -0.011 0.635
## ns_scr_d:bigsegmentemoangerhappy -0.020 0.605
## ns_scr_d:bigsegmentemoangersad -0.021 0.577
## ns_scr_d:visit 0.003 0.066
## ns_scr_d:bgsgmntmngrf
## bigsegmentemoangeranger
## bigsegmentemoangerbaseline
## bigsegmentemoangerfear
## bigsegmentemoangerhappy
## bigsegmentemoangersad
## visit
## ns_scr_d:bigsegmentemoangerbaseline
## ns_scr_d:bigsegmentemoangerfear
## ns_scr_d:bigsegmentemoangerhappy 0.509
## ns_scr_d:bigsegmentemoangersad 0.481
## ns_scr_d:visit 0.022
## ns_scr_d:bgsgmntmngrh
## bigsegmentemoangeranger
## bigsegmentemoangerbaseline
## bigsegmentemoangerfear
## bigsegmentemoangerhappy
## bigsegmentemoangersad
## visit
## ns_scr_d:bigsegmentemoangerbaseline
## ns_scr_d:bigsegmentemoangerfear
## ns_scr_d:bigsegmentemoangerhappy
## ns_scr_d:bigsegmentemoangersad 0.459
## ns_scr_d:visit 0.044
## ns_scr_d:bgsgmntmngrs
## bigsegmentemoangeranger
## bigsegmentemoangerbaseline
## bigsegmentemoangerfear
## bigsegmentemoangerhappy
## bigsegmentemoangersad
## visit
## ns_scr_d:bigsegmentemoangerbaseline
## ns_scr_d:bigsegmentemoangerfear
## ns_scr_d:bigsegmentemoangerhappy
## ns_scr_d:bigsegmentemoangersad
## ns_scr_d:visit 0.050
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.757422316 -0.572517216 0.003633262 0.570220542 9.194947078
##
## Number of Observations: 22168
## Number of Groups: 328anova(model_situation.emofear)## numDF denDF F-value p-value
## ns_scr_d 1 21829 22.452992 <.0001
## bigsegmentemofear 5 21829 22.716655 <.0001
## visit 1 21829 0.069266 0.7924
## ns_scr_d:bigsegmentemofear 4 21829 4.237850 0.0020
## ns_scr_d:visit 1 21829 0.016156 0.8989Follow-up for examination of “habituation”.
#Checking habituation
#merging to get participation variable
epochdataplus.habit <- merge(epochdataplus,persondata339,by="id")
mean(epochdataplus.habit$present_all)## [1] 2.448378model_situation.habit <- lme(rsa_d ~ 0 + ns_scr_d + avoid + present_all +
ns_scr_d:avoid + ns_scr_d:present_all, #+ ns_scr_d:present_all:avoid,
random = ~ 0 + ns_scr_d |id,
data = epochdataplus.habit,
na.action = na.exclude)
summary(model_situation.habit)## Linear mixed-effects model fit by REML
## Data: epochdataplus.habit
## AIC BIC logLik
## 37924.27 37980.31 -18955.13
##
## Random effects:
## Formula: ~0 + ns_scr_d | id
## ns_scr_d Residual
## StdDev: 0.02793115 0.5659574
##
## Fixed effects: rsa_d ~ 0 + ns_scr_d + avoid + present_all + ns_scr_d:avoid + ns_scr_d:present_all
## Value Std.Error DF t-value p-value
## ns_scr_d -0.017840623 0.009568159 21837 -1.8645826 0.0623
## avoid 0.002682667 0.007491404 21837 0.3580994 0.7203
## present_all -0.000019660 0.001928358 327 -0.0101954 0.9919
## ns_scr_d:avoid 0.008884500 0.003324741 21837 2.6722378 0.0075
## ns_scr_d:present_all 0.000829876 0.003484960 21837 0.2381307 0.8118
## Correlation:
## ns_sc_ avoid prsnt_ ns_s_:
## avoid 0.006
## present_all -0.009 -0.710
## ns_scr_d:avoid -0.187 0.008 0.017
## ns_scr_d:present_all -0.952 0.000 0.001 0.011
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.810180982 -0.574301243 0.004946749 0.567966293 9.150727265
##
## Number of Observations: 22168
## Number of Groups: 328