### Data Preparation included cleaning, and transformations, see cleaning .rmd for details
#### Load Relevant Libraries and Functions
#### Import data
#### Data exclusion / filtering - all participants who did did nothave five measures of stress and hr are in this data - kept those who only had stroop/math as their tp1 and 2 tasks dropped others, droppped NAs at hradn stress data aswellneeded all 5
#### Prepare data for analysis - create columns etc.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(knitr)
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(readr)
library(gtsummary)
library(png)
#load wide and long data sets

psych_251_report_cleaned_wide <- read_csv("/Users/krise/R Mac working folder/sommerfeldt2019/psych.251.report.cleaned.wide.csv")
## Warning: Missing column names filled in: 'X1' [1]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
psych_251_report_cleaned_long <- read_csv("/Users/krise/R Mac working folder/sommerfeldt2019/psych.251.report.cleaned.long.csv")
## Warning: Missing column names filled in: 'X1' [1]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
wDF <- psych_251_report_cleaned_wide
f.Ldf <- psych_251_report_cleaned_long

Introduction

This study aimed to establish and validate the value of a particular construct described as coherence. Coherence is a measure of accordance between self-reported stress and heart rate. Evidence suggests coherence is associated with healthy mind-body functioning and adaptivity. The authors hypothesize that the within-participants association between self-reported stress and heart rate is positively related to psychological and physical well-being and negatively related to denial coping.

I chose to reproduce the analyses for this study due to the rich, multi-modal, and dimensional data available and my interest in adding the methods to my toolkit. I’m particularly interested in interventions that use this measure as a means of enhancement of meta-awareness, health, and well-being.

Key analysis

Using linear mixed-effects models, I will use the coherence model below to investigate the relationship of coherence with measures of health and well-being. Fundamental analyses include replicating the stress-heart rate model to predict well-being across various measures (psychological well-being, depression, anxiety, interleukin-6, and denial coping).

I will focus on the interaction effect in this model, representing the degree to which within-participants associations between self-reported stress and heart rate are related to the well-being indicator (PWB, depression, anxiety, IL-6, and CRP) or denial coping.

lmer(heartRate stressClusterMeanCentered ∼ * wellbeingCentered + StressClusterMeanCentered * ageCentered(1 + stressClusterMeanCentered | subject) + (1 + d | subject) + () = 1 | family data dfLong )

Justification for choice of study

I chose to reproduce this study for three main reasons. One is because I get to use a rich and robust data set that would be impossible to replicate due to the expansive and longitudinal nature. Additionally, this data set’s multi-modal nature may stimulate future inquiry and publication that spans across domains. The second and third reason is my interest in gaining a better grasp and practice using the statistical and operational methods in this study for possible future use or exploratory analyses.

Anticipated challenges

I foresee some challenges related to data analysis and interpretation based on adequate cleaning and wrangling. I do have access to the raw and processed data, which should reduce that challenge, but because this requires cleaning, tidying, and transforming many different variables, it introduces chances of error. It may be challenging to understand and justify the methods of the models fully.

Methods

"Participants completed a standardized laboratory-based experimental stress-induction paradigm designed to measure cardiovascular reactivity and recovery from stress (Crowley et al., 2011) detailed documentation of the study protocol is publicly available at http://www.midus.wisc.edu/midus2/project4/). The data were collected at the University of California, Los Angeles; Georgetown University; and the University of Wisconsin and processed at the Columbia University Medical Center in the laboratory of Richard Sloan.

The stress-induction paradigm involved a resting baseline (11 min); two cognitive-psychological stressor tasks (6 min each; counterbalanced across participants); a seated, resting period after each task (recovery period; 6 min each); and an orthostatic challenge, which involved moving from a seated to a standing position and remaining standing (6 min). The orthostatic phase of the task was not included in the analyses because changes in heart rate during this phase are confounded with physical movement. Thus, we examined five phases of interest: baseline, first stressor task, first recovery, second stressor task, and second recovery.

Participants’ heart rate was measured using electrocardiograph electrodes placed on the left and right shoulders and in the left lower quadrant. Heart rate was measured continuously over every phase of the task. Heart rate was calculated as an average of all valid interbeat intervals and converted from interbeat-interval units (milliseconds) to beats-per-minute units. The average of a 5-min epoch was analyzed for each of the five phases of the task. Each epoch was scored for quality, and only epochs containing a full 5 min of good signal quality, without any designated invalid intervals of data that had to be omitted, were included in the analysis. We chose to examine the average heart rate for each phase of the task because the precise timing of each subjective report was not recorded on the physiological time series, and subjective reports did not necessarily occur during the peak physiological response.

Participants were informed at the beginning of the session that, periodically, they would be asked for a verbal stress rating on a scale from 1 (not stressed at all) to 10 (extremely stressed). The experimenter prompted each participant to verbally report his or her level of stress approximately 20 to 30 s before the end of each phase of the task. Thus, a total of six self-reports of stress were collected during the session, near the end of each phase: baseline, during each stressor task, during the recovery period following each stressor task, and after the orthostatic challenge. The first five self-reports of stress were used, excluding the ortho-static time point.".

One can estimate a linear mixed-effects model (LMEM) to examine whether the (statistical) effect of one of the Level 1 variables (e.g., subjective stress) on the other Level 1 variable (e.g., heart rate) is moderated by the individual-differences variable. If, for example, the effect of subjective stress on heart rate is stronger for participants high in psychological well-being, then

Age was included as a covariate because of the broad age range of the sample, extending from early to late adulthood and because older participants had lower stress–heart rate coherence, b = −0.008, F(1, 843.0) = 7.754, p = .005. Gender was not associated with stress–heart rate coherence, b = 0.051, F(1, 850.0) = 0.560, p = .455, and so was not included as a covariate in the analyses. We fitted a separate model for each of the five well-being indicators of interest and denial coping (six total tests). The Anova() function in the car package (Version 3.0.0; Fox & Weisberg, 2011) provided estimates of F, error df (via Kenward-Roger approximation), and p. Multiple comparisons of the six different tests were corrected using the Holm-Bonferroni method the within-participants association is positively related to psychological well-being

Stroop color-word task. Participants completed a modified Stroop color-word task (Stroop, 1935). One of four color-name words was presented in a font color that was either congruent or incongruent with the word itself. The colored color-name stimulus appeared on screen, and participants pressed one of four keys on a keypad cor- responding to the color of the letters in the word, not the color name. The rate of stimuli was modified according to participant performance to roughly standardize the degree of stressfulness. This standardization was set so that

For the LMEM approach, I will regress heart rate on self-reported stress (centered around each participant’s own mean), the well-being indicator under consideration (mean centered; e.g., PWB), and their interaction, adjusting for age, the interaction between self-reported stress and age, and nonindependence due to participants and families (Brauer & Curtin, 2018). Our model thus includes six fixed effects: self-reported stress (Level 1), the well-being indicator of interest (Level 2), their interaction, age (Level 2), the interaction of self- reported stress and age, and the intercept. The model includes a by-participant random intercept, a by-participant random slope for stress, and a by-family random intercept. The two by-participant random effects were allowed to correlate.

This model was represented in R as follows: lmer(heartRate stressClusterMeanCentered ∼ * wellbeingCentered + StressClusterMeanCentered * ageCentered(1 + stressClusterMeanCentered | subject) + (1 + d | subject) + () = 1 | family data dfLong )

Description of the steps required to reproduce the results

Retrieve, clean, analyze, interpret, and visualize key data, models and tables.

Differences from original study

For this reproduction, the final N for analyses was 967. This population count is lower than the original paper (n=1065) reported and is likely due to different exclusion parameters resulting in a difference of 98 participants. Unlike the authors, I removed participants who did not have five self-reported stress and heart rate measures, as opposed to just stress. Additionally, I did not change scores on the PWB measure to NA if they equaled 98, which was not reported in the paper but shown in the authors’ coding script. I likely eliminated more participants in cleaning and transformations, which may have impacted variables differently.

Results

In general, the results support the critical hypotheses in the paper. Stress-heart rate coherence is significantly associated with measures of health and well-being. However, although trending the same way, this was not significant for C-reactive protein. Below I report estimates, standard errors, degrees of freedom, F and P values. See figures at bottom of document.

The statistical effect of stress on heart rate was found to be moderated by PWB, b = 0.04, SE =.011, F(1, 532) = 13.47, p < .0001; participants with higher stress–heart rate coherence also reported higher psychological well-being. The opposite was true for depressive symptoms, b = −0.23, SE=.05, F(1, 500) = 20.63, p < .0001, and trait anxiety, b = −0.18, SE=.04 F(1, 488) = 18.45, p < .0001; individuals with higher stress–heart rate coherence reported fewer depressive symptoms and had lower trait anxiety.

For physical well-being, the statistical effect of stress on heart rate was found to be significantly moderated by IL-6 and CRP; participants with higher stress–heart rate coherence also had lower IL-6, b = −0.131, SE=.037, F(1, 485) = 12.62, p < .0001, and lower CRP, b = −0.147, SE=.077 F(1, 522) = 3.64, p = .057. (trending-almost significant)

I also investigated whether stress–heart rate coherence was associated with the use of denial as a coping strategy. The statistical effect of stress on heart rate was found to be moderated by denial; higher stress–heart rate coherence was associated with less tendency toward the use of denial as a coping strategy, b = −0.069, SE=.018 F(1, 536) = 13.79, p < .0001.

#basic summary stats of data
f.Ldf %>% tbl_summary(statistic = list(all_continuous() ~ "{mean} ({sd})")) #Mean and SDs fo all variables, those that were clustered, mean centered, and log transformed for analyses
Characteristic N = 4,8351
X1 2,418 (1,396)
M2ID 14,766 (2,658)
timepoint
1 967 (20%)
2 967 (20%)
3 967 (20%)
4 967 (20%)
5 967 (20%)
stress_value 3.11 (2.12)
hr_value 75 (11)
pwb2.x 233 (35)
Unknown 100
COPE_denial.x 6.09 (2.22)
Unknown 95
IL6.x 2.99 (2.87)
Unknown 25
CRP.x 2.88 (4.21)
Unknown 55
P4_CESD.x 9 (8)
Unknown 35
P4_STAItrait.x 34 (9)
Unknown 40
M2FAMNUM.x 110,503 (8,314)
Unknown 895
P4_age.x 56 (11)
stress_CMC 0.00 (1.63)
hr_CMC 0.00 (2.69)
pwb2_CMC 0 (0%)
Unknown 100
P4_CESD_CMC 0 (0%)
Unknown 35
COPE_d_CMC 0 (0%)
Unknown 95
IL6_CMC 0 (0%)
Unknown 25
CRP_CMC 0 (0%)
Unknown 55
STAI_CMC 0 (0%)
Unknown 40
P4_age_MC 0 (11)
pwb2_MC 0 (35)
Unknown 100
P4_CESD_MC 0 (8)
Unknown 35
P4_STAItrait_MC 0 (9)
Unknown 40
COPE_denial_MC 0.00 (2.22)
Unknown 95
pwb2_MC_d10 0.0 (3.5)
Unknown 100
P4_CESD_MC_d10 0.00 (0.82)
Unknown 35
P4_STAItrait_MC_d10 0.00 (0.90)
Unknown 40
IL6_MC 0.00 (2.87)
Unknown 25
CRP_MC 0.00 (4.21)
Unknown 55
IL6_T 1.16 (1.05)
Unknown 25
CRP_T 0.17 (0.51)
Unknown 55
IL6_T_MC 0.00 (1.05)
Unknown 25
CRP_T_MC 0.00 (0.51)
Unknown 55

1 Statistics presented: Mean (SD); n (%)

Data preparation

Data preparation included a thorough investigation of coding books, variables, cleaning, wrangling, and transforming (see cleaning .rmd for details). The variables in this data set are only related to key variables. Participants were excluded if they did not complete Stroop or MATH task for stress tasks or did not have five heart rates and self-reported stress measurements.

Key analysis - Origiinal test and results

The statistical effect of stress on heart rate was found to be moderated by PWB, b = 0.050, F(1, 822.8) = 26.70, p < .0001; participants with higher stress–heart rate coherence also reported higher psychological well- being. The opposite was true for depressive symptoms, b = −0.249, F(1, 783.7) = 36.77, p < .0001, and trait anxiety, b = −0.211, F(1, 769.4) = 32.49, p < .0001; individuals with higher stress–heart rate coherence reported fewer depressive symptoms and had lower trait anxiety.

For physical well- being, the statistical effect of stress on heart rate was found to be significantly moderated by IL-6 and CRP; participants with higher stress–heart rate coherence also had lower IL-6, b = −0.145, F(1, 762.3) = 22.20, p < .0001, and lower CRP, b = −0.175, F(1, 827.2) = 7.16, p = .008.

I also investigated whether stress–heart rate coherence was associated with use of denial as a coping strategy.The statistical effect of stress on heart rate was found to be moderated by denial; higher stress–heart rate coherence was associated with less tendency toward the use of denial as a coping strategy, b = −0.069, F(1, 853.3) = 20.69, p < .0001

Side-by-side graph with original graph is ideal here Original author did not graph but listed the above numbers in a table. See below for tables and original or put table here.

Exploratory analyses

None to report.

Discussion

Primarily, the results indicate a significant relationship between heart rate and the interaction of stress and measures of health and well being. Specifically, individuals with higher stress-heart rate coherence also report higher psychological well-being, fewer depressive symptoms, lower trait anxiety, IL-6, CRP and are less likely to use denial as a coping strategy. I had similar results as the original authors for all the variables except CRP. Participants with higher stress-heart-rate coherence had lower CRP than those with lower stress-heart rate coherence, but this was non-significant but trending in the same direction as the original paper.

Summary of Reproduction Attempt

This reproduction was successful for two reasons. 1) I was able to retrieve, prepare, and analyze the data set using a similar methodology, and 2) this resulted in similar results regarding the main hypotheses. I reproduced the two central figures of the paper along with the models that reflect the main hypotheses that stress-heart rate coherence is associated with psychophysiological health and well-being. While the figures differ slightly, the overall relationships are significantly trending in the same direction, and models produce similar estimates, F, and P values. Although we did differ on degrees of freedom listed in the models, this is likely due to the population size and difference in cleaning/elimination methods that may have resulted in different population sizes. The summary statistics were similar despite a smaller sample. The difference in populations may explain the difference in results related to CRP.

Commentary

Overall, the original authors did an excellent job providing insight into their methods and approach and made this mostly available in an online repository. Access to materials and methods much supported this reproductive endeavor and is a critical part of open science. I only had one objection to a scoring method with the psychological well-being variable that I am not sure I understood, but this did not result in meaningful differences in the model results.

While the authors provided their code, there were no data to work from, so I had the opportunity to practice retrieving and digging through a massive data set with 10k+ observers and variables. This process helped me learn much more about the study and appreciate providing clean data sets when choosing open science methods. It also highlighted the original scientists’ incredible efforts and rigor that created and shared the original data because they had to anonymize track code and share these data sets. Their efforts have catalyzed hundreds of papers and highlights how open science can create new knowledge.

This reproduction taught me the value of using code, names, comments, and packages together optimally. This integrative consideration makes it easy to understand thinking processes and reproduce work from others and me! I learned that I need more practice wrangling data. I need to practice understanding why you want data to look a certain way and how to make it look that way for proper or most effective analyses. I now have some workflows, processes, and packages that will support the scientific journey.

TESTS

Stress-heart rate coherence association models

Psychological Well-being

# Run the Lmer tests
lmerM.pwb <- lmerTest::lmer(hr_value ~ stress_CMC * pwb2_MC_d10 + P4_age_MC*stress_CMC + (1 + stress_CMC|M2ID) + (1|M2FAMNUM.x), data = f.Ldf)
a.pwb <- anova(lmerM.pwb, type=3, test="F")
a.pwb
## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## stress_CMC             3347.5  3347.5     1 534.69 594.8275 < 2.2e-16 ***
## pwb2_MC_d10               0.7     0.7     1 784.39   0.1193 0.7298802    
## P4_age_MC               100.1   100.1     1 687.25  17.7861 2.803e-05 ***
## stress_CMC:pwb2_MC_d10   75.9    75.9     1 532.16  13.4785 0.0002657 ***
## stress_CMC:P4_age_MC     76.0    76.0     1 534.19  13.4961 0.0002631 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#pwb.sum <- summary(lmerM.pwb)


tab_model(lmerM.pwb)
  hr_value
Predictors Estimates CI p
(Intercept) 74.74 73.97 – 75.51 <0.001
stress_CMC 0.94 0.87 – 1.02 <0.001
pwb2_MC_d10 0.04 -0.19 – 0.27 0.730
P4_age_MC -0.15 -0.22 – -0.08 <0.001
stress_CMC * pwb2_MC_d10 0.04 0.02 – 0.06 <0.001
stress_CMC * P4_age_MC -0.01 -0.02 – -0.01 <0.001
Random Effects
σ2 5.63
τ00 M2ID 92.34
τ00 M2FAMNUM.x 22.97
τ11 M2ID.stress_CMC 0.56
ρ01 M2ID 0.18
ICC 0.95
N M2ID 788
N M2FAMNUM.x 686
Observations 3940
Marginal R2 / Conditional R2 0.041 / 0.956

Depression

#run the tests
lmerM.dep <- lmerTest::lmer(hr_value ~ stress_CMC *P4_CESD_MC_d10 + P4_age_MC*stress_CMC + (1 + stress_CMC|M2ID) + (1|M2FAMNUM.x), data = f.Ldf)
a.dep<- anova(lmerM.dep, type=3, test="F")
a.dep
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## stress_CMC                3333.2  3333.2     1 524.06 589.9670 < 2.2e-16 ***
## P4_CESD_MC_d10               3.0     3.0     1 778.93   0.5386 0.4632335    
## P4_age_MC                   95.4    95.4     1 671.52  16.8817  4.47e-05 ***
## stress_CMC:P4_CESD_MC_d10  116.5   116.5     1 500.89  20.6267  7.00e-06 ***
## stress_CMC:P4_age_MC        77.6    77.6     1 519.08  13.7335 0.0002333 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(lmerM.dep)

tab_model(lmerM.dep)
  hr_value
Predictors Estimates CI p
(Intercept) 74.75 73.97 – 75.52 <0.001
stress_CMC 0.94 0.86 – 1.01 <0.001
P4_CESD_MC_d10 0.38 -0.63 – 1.38 0.463
P4_age_MC -0.15 -0.22 – -0.08 <0.001
stress_CMC *
P4_CESD_MC_d10 		-0.23 -0.33 – -0.13 <0.001
stress_CMC * P4_age_MC -0.01 -0.02 – -0.01 <0.001
Random Effects
σ2 5.65
τ00 M2ID 91.23
τ00 M2FAMNUM.x 23.75
τ11 M2ID.stress_CMC 0.54
ρ01 M2ID 0.20
ICC 0.95
N M2ID 785
N M2FAMNUM.x 684
Observations 3925
Marginal R2 / Conditional R2 0.042 / 0.956

Anxiety

# Run the tests

lmerM.anx <- lmerTest::lmer(hr_value ~ stress_CMC *P4_STAItrait_MC_d10 + P4_age_MC*stress_CMC + (1 + stress_CMC|M2ID) + (1|M2FAMNUM.x), data = f.Ldf)
a.anx<- anova(lmerM.anx, type=3, test="F")
a.anx
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)
## stress_CMC                     3341.2  3341.2     1 527.46 592.5379 < 2.2e-16
## P4_STAItrait_MC_d10               0.2     0.2     1 771.68   0.0437 0.8344189
## P4_age_MC                       105.2   105.2     1 681.00  18.6648 1.791e-05
## stress_CMC:P4_STAItrait_MC_d10  104.1   104.1     1 488.85  18.4588 2.096e-05
## stress_CMC:P4_age_MC             72.9    72.9     1 523.17  12.9324 0.0003535
##                                   
## stress_CMC                     ***
## P4_STAItrait_MC_d10               
## P4_age_MC                      ***
## stress_CMC:P4_STAItrait_MC_d10 ***
## stress_CMC:P4_age_MC           ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(lmerM.anx)

tab_model(lmerM.anx)
  hr_value
Predictors Estimates CI p
(Intercept) 74.71 73.94 – 75.49 <0.001
stress_CMC 0.94 0.86 – 1.02 <0.001
P4_STAItrait_MC_d10 0.09 -0.78 – 0.97 0.834
P4_age_MC -0.16 -0.23 – -0.09 <0.001
stress_CMC *
P4_STAItrait_MC_d10 		-0.19 -0.27 – -0.10 <0.001
stress_CMC * P4_age_MC -0.01 -0.02 – -0.01 <0.001
Random Effects
σ2 5.64
τ00 M2ID 92.29
τ00 M2FAMNUM.x 22.48
τ11 M2ID.stress_CMC 0.55
ρ01 M2ID 0.19
ICC 0.95
N M2ID 785
N M2FAMNUM.x 683
Observations 3925
Marginal R2 / Conditional R2 0.043 / 0.956

IL6

# Run the tests
lmerM.IL6<- lmerTest::lmer(hr_value ~ stress_CMC *IL6_T_MC + P4_age_MC*stress_CMC + (1 + stress_CMC|M2ID) + (1|M2FAMNUM.x), data = f.Ldf)
a.IL6<- anova(lmerM.IL6, type=3, test="F")
a.IL6
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## stress_CMC           3240.2  3240.2     1 522.37 575.6749 < 2.2e-16 ***
## IL6_T_MC               69.8    69.8     1 778.52  12.3937 0.0004558 ***
## P4_age_MC             133.1   133.1     1 679.16  23.6530 1.435e-06 ***
## stress_CMC:IL6_T_MC    71.1    71.1     1 485.74  12.6262 0.0004175 ***
## stress_CMC:P4_age_MC   32.6    32.6     1 519.31   5.7972 0.0163995 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(lmerM.IL6)

tab_model(lmerM.IL6)
  hr_value
Predictors Estimates CI p
(Intercept) 74.90 74.13 – 75.68 <0.001
stress_CMC 0.92 0.85 – 1.00 <0.001
IL6_T_MC 1.36 0.60 – 2.11 <0.001
P4_age_MC -0.17 -0.24 – -0.10 <0.001
stress_CMC * IL6_T_MC -0.13 -0.20 – -0.06 <0.001
stress_CMC * P4_age_MC -0.01 -0.02 – -0.00 0.016
Random Effects
σ2 5.63
τ00 M2ID 93.44
τ00 M2FAMNUM.x 20.08
τ11 M2ID.stress_CMC 0.53
ρ01 M2ID 0.22
ICC 0.95
N M2ID 785
N M2FAMNUM.x 683
Observations 3925
Marginal R2 / Conditional R2 0.055 / 0.956

CRP

#Run the tests
lmerM.CRP<- lmerTest::lmer(hr_value ~ stress_CMC *CRP_T_MC + P4_age_MC*stress_CMC + (1 + stress_CMC|M2ID) + (1|M2FAMNUM.x), data = f.Ldf)
a.CRP <- anova(lmerM.CRP, type=3, test="F")
a.CRP
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## stress_CMC           3262.7  3262.7     1 522.59 577.814 < 2.2e-16 ***
## CRP_T_MC              130.5   130.5     1 770.10  23.113 1.837e-06 ***
## P4_age_MC             102.4   102.4     1 680.50  18.139 2.343e-05 ***
## stress_CMC:CRP_T_MC    20.5    20.5     1 522.06   3.638  0.057023 .  
## stress_CMC:P4_age_MC   53.2    53.2     1 520.04   9.422  0.002256 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(lmerM.CRP)

tab_model(lmerM.CRP)
  hr_value
Predictors Estimates CI p
(Intercept) 74.89 74.12 – 75.66 <0.001
stress_CMC 0.93 0.86 – 1.01 <0.001
CRP_T_MC 3.71 2.20 – 5.23 <0.001
P4_age_MC -0.15 -0.22 – -0.08 <0.001
stress_CMC * CRP_T_MC -0.15 -0.30 – 0.00 0.056
stress_CMC * P4_age_MC -0.01 -0.02 – -0.00 0.002
Random Effects
σ2 5.65
τ00 M2ID 86.17
τ00 M2FAMNUM.x 25.93
τ11 M2ID.stress_CMC 0.55
ρ01 M2ID 0.21
ICC 0.95
N M2ID 782
N M2FAMNUM.x 681
Observations 3910
Marginal R2 / Conditional R2 0.067 / 0.956

Denial

#Run the tests
lmerM.den<- lmerTest::lmer(hr_value ~ stress_CMC *COPE_denial_MC+ P4_age_MC*stress_CMC + (1 + stress_CMC|M2ID) + (1|M2FAMNUM.x), data = f.Ldf)

a.den<- anova(lmerM.den, type=3, test="F")
a.den
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## stress_CMC                3270.9  3270.9     1 535.84 584.1222 < 2.2e-16 ***
## COPE_denial_MC               0.6     0.6     1 783.91   0.1059 0.7449404    
## P4_age_MC                   99.6    99.6     1 680.81  17.7864 2.806e-05 ***
## stress_CMC:COPE_denial_MC   77.3    77.3     1 536.14  13.7959 0.0002251 ***
## stress_CMC:P4_age_MC        41.3    41.3     1 533.08   7.3787 0.0068143 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(lmerM.den)

tab_model(lmerM.den)
  hr_value
Predictors Estimates CI p
(Intercept) 74.71 73.94 – 75.49 <0.001
stress_CMC 0.92 0.85 – 0.99 <0.001
COPE_denial_MC -0.06 -0.43 – 0.31 0.745
P4_age_MC -0.15 -0.22 – -0.08 <0.001
stress_CMC *
COPE_denial_MC 		-0.07 -0.10 – -0.03 <0.001
stress_CMC * P4_age_MC -0.01 -0.02 – -0.00 0.007
Random Effects
σ2 5.60
τ00 M2ID 93.74
τ00 M2FAMNUM.x 21.32
τ11 M2ID.stress_CMC 0.52
ρ01 M2ID 0.17
ICC 0.95
N M2ID 787
N M2FAMNUM.x 686
Observations 3935
Marginal R2 / Conditional R2 0.041 / 0.956

Summary of the mixed lmer models testing the main hypotheses, that stress-heart rate coherence is associated with improved measures of psychophysiological health and well-being

PWB.regression.tbl  <- tbl_regression(lmerM.pwb, label = list(stress_CMC ~ "Stress", pwb2_MC_d10 ~ "Psychological Well-Being", P4_age_MC ~ "Age", stress_CMC ~ "Stress")) %>%
   italicize_labels() %>%
   as_kable_extra()
## Registered S3 method overwritten by 'broom.mixed':
##   method      from 
##   tidy.gamlss broom
DEP.regression.tbl  <- tbl_regression(lmerM.dep, label = list(stress_CMC ~ "Stress", P4_CESD_MC_d10 ~ "Depression", P4_age_MC ~ "Age", stress_CMC ~ "Stress")) %>%
   italicize_labels() %>%
   as_kable_extra() 
  
ANX.regression.tbl  <- tbl_regression(lmerM.anx,label = list(stress_CMC ~ "Stress", P4_STAItrait_MC_d10 ~ "Anxiety", P4_age_MC ~ "Age", stress_CMC ~ "Stress")) %>%
    italicize_labels() %>%
   as_kable_extra() 

IL6.regression.tbl <- tbl_regression(lmerM.IL6, label = list(stress_CMC ~ "Stress", IL6_T_MC ~ "Interleukin-6", P4_age_MC ~ "Age", stress_CMC ~ "Stress")) %>%
    italicize_labels() %>%
   as_kable_extra() 

CRP.regression.tbl <- tbl_regression(lmerM.CRP, list(stress_CMC ~ "Stress", CRP_T_MC ~ "C-Reactive Protein", P4_age_MC ~ "Age", stress_CMC ~ "Stress")) %>%
    italicize_labels() %>%
   as_kable_extra()

DEN.regression.tbl <- tbl_regression(lmerM.den, list(stress_CMC ~ "Stress", COPE_denial_MC ~ "Denial", P4_age_MC ~ "Age", stress_CMC ~ "Stress")) %>%
    italicize_labels() %>%
   as_kable_extra()

Table2 - Summary of key variable models

PWB.regression.tbl
Characteristic Beta 95% CI p-value
Stress 0.94 0.87, 1.0 <0.001
Psychological Well-Being 0.04 -0.19, 0.27 0.7
Age -0.15 -0.22, -0.08 <0.001
Stress * Psychological Well-Being 0.04 0.02, 0.06 <0.001
Stress * Age -0.01 -0.02, -0.01 <0.001
1 CI = Confidence Interval 
DEP.regression.tbl
Characteristic Beta 95% CI p-value
Stress 0.94 0.86, 1.0 <0.001
Depression 0.38 -0.63, 1.4 0.5
Age -0.15 -0.22, -0.08 <0.001
Stress * Depression -0.23 -0.33, -0.13 <0.001
Stress * Age -0.01 -0.02, -0.01 <0.001
1 CI = Confidence Interval 
ANX.regression.tbl
Characteristic Beta 95% CI p-value
Stress 0.94 0.86, 1.0 <0.001
Anxiety 0.09 -0.78, 1.0 0.8
Age -0.16 -0.23, -0.09 <0.001
Stress * Anxiety -0.19 -0.27, -0.10 <0.001
Stress * Age -0.01 -0.02, -0.01 <0.001
1 CI = Confidence Interval 
IL6.regression.tbl
Characteristic Beta 95% CI p-value
Stress 0.92 0.85, 1.0 <0.001
Interleukin-6 1.4 0.60, 2.1 <0.001
Age -0.17 -0.24, -0.10 <0.001
Stress * Interleukin-6 -0.13 -0.20, -0.06 <0.001
Stress * Age -0.01 -0.02, 0.00 0.016
1 CI = Confidence Interval 
CRP.regression.tbl
Characteristic Beta 95% CI p-value
Stress 0.93 0.86, 1.0 <0.001
C-Reactive Protein 3.7 2.2, 5.2 <0.001
Age -0.15 -0.22, -0.08 <0.001
Stress * C-Reactive Protein -0.15 -0.30, 0.00 0.057
Stress * Age -0.01 -0.02, 0.00 0.002
1 CI = Confidence Interval 
DEN.regression.tbl
Characteristic Beta 95% CI p-value
Stress 0.92 0.85, 1.0 <0.001
Denial -0.06 -0.43, 0.31 0.7
Age -0.15 -0.22, -0.08 <0.001
Stress * Denial -0.07 -0.10, -0.03 <0.001
Stress * Age -0.01 -0.02, 0.00 0.007
1 CI = Confidence Interval

FIGURE 1: Stress and heart rate by phase histograms

ylimits = c(0, 610)
stressHist=ggplot()+
geom_histogram(data=f.Ldf, aes(stress_value), fill="green", binwidth=1, color="black") +
facet_wrap(~timepoint, ncol=5) +
labs(x="Self-reported stress", y="Number of subjects") +
ylim(ylimits)+
scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9,10) )+
theme_bw(base_size=18)
stressHist

xlimits = c(30,150)
ylimits = c(0, 240)
hrHist=ggplot()+
geom_histogram(data=f.Ldf, aes(hr_value), fill="red", binwidth=6, color="black") +
facet_wrap(~timepoint, ncol=5) +
labs(x="Heart rate", y="Number of subjects") +
ylim(ylimits)+
xlim(xlimits)+
geom_vline(xintercept=75, size=.5, color="blue")+
theme_bw(base_size=18)
hrHist
## Warning: Removed 10 rows containing missing values (geom_bar).