Supplementary Analyses

The emodiversity construct was put forth by Quoidbach and colleagues (2014) describing the variety and relative abundance of individuals’ discrete emotion experiences. The original study used two samples and asked individuals to provide a single occasion rating of a number of discrete emotions. Results indicated that global, positive, and negative emodiversity were related to a variety of psychological and physical health items after accounting for mean levels of positive and negative emotions. Building on this initial work, Benson and colleagues (2017) published a methodological investigation of emodiversity using 30 days of diary data. Most recently, Brown & Coyne (2017) published a critique of the intial emodiversity paper, including critiques involving hierarchical regression and suppression. The following are a set of supplementary analyses utilizing the same As U Live data (Sturgeon, Zautra, & Okun, 2014) presented in the paper by Benson and colleagues (2017). Individuals reading these supplementary analyses are also encouraged to read the original articles on emodiversity.

References

Quoidbach, J., Gruber, J., Mikolajczak, M., Kogan, A., Kotsou, I., & Norton, M. I. (2014). Emodiversity and the emotional ecosystem. Journal of experimental psychology: General, 143(6), 2057. DOI: doi:10.1037/a0038025

Benson, L., Ram, N., Almeida, D. M., Zautra, A. J., & Ong, A. D. (2017). Fusing Biodiversity Metrics into Investigations of Daily Life: Illustrations and Recommendations With Emodiversity. The Journals of Gerontology: Series B. DOI: https://doi.org/10.1093/geronb/gbx025

Brown, N. J. L., & Coyne, J. C. (2017). Emodiversity: Robust Predictor of Outcomes or Statistical Artifact? Journal of Experimental Psychology. General. DOI: 10.1037/xge0000330

Summary of Supplementary Analyses

First, emodiversity calculated using the diary emotion data was examined (utilizing just the Gini coefficient):

  • model 1 hierarchical series investigating (a) mean positive emotion, (b) mean positive emotion and positive emodiversity, and (c) mean positive emotion, positive emodiversity, and the interaction. Results suggested no suppression was occuring. Model comparisons suggested that model1a with just mean positive emotion provided the best fit to these data.
  • model 2 hierarchical series investigating (a) mean negative emotion, (b) mean negative emotion and negative emodiversity, and (c) mean negative emotion, negative emodiversity, and the interaction. Results suggested suppression was occuring in terms of the magnitude of the beta coefficient increasing for mean negative emotion and the sign flipping for negative emodiversity. Model comparisons suggested that model2c with all three terms provided the best fit to these data.
  • model 3 hierarchical series investigating (a) mean positive and mean negative emotion, (b) mean positive emotion, mean negative emotion, positive emodiversity, and negative emodiversity, and (c) mean positive emotion, mean negative emotion, positive emodiversity, and negative emodiversity, and the interactions. Results suggested suppression was occuring in terms of the magnitude of the beta coefficient increasing for mean negative emotion and the sign flipping for negative emodiversity. Model comparisons suggested that model3c including all terms provided the best fit to these data.

Second, emodiversity calculated using the single occasion emotion data was examined (utilizing both the Gini coefficient and Shannon’s entropy):

  • model 4 hierarchical series utilizing the Gini coefficient, investigating (a) mean positive emotion, (b) mean positive emotion and positive emodiversity, and (c) mean positive emotion, positive emodiversity, and the interaction. Results suggested no suppression was occuring. Model comparisons suggested that model4b with mean positive emotion and positive emodiversity provided the best fit to these data.
  • model 4 hierarchical series utilizing Shannon’s entropy, investigating (d) mean positive emotion, (e) mean positive emotion and positive emodiversity, and (f) mean positive emotion, positive emodiversity, and the interaction. Results suggested no suppression was occuring. Model comparisons suggested that model4e with mean positive emotion and positive emodiversity provided the best fit to these data.
  • model 5 hierarchical series utilizing the Gini coefficient, investigating (a) mean negative emotion, (b) mean negative emotion and negative emodiversity, and (c) mean negative emotion, negative emodiversity, and the interaction. Results suggested suppression was occuring in terms of the magnitude of the beta coefficient increasing for mean negative emotion and the sign flipping for negative emodiversity. Model comparisons suggested that model5a with only mean negative emotion, provided the best fit to these data.
  • model 5 hierarchical series utilizing Shannon’s entropy, investigating (d) mean negative emotion, (e) mean negative emotion and negative emodiversity, and (f) mean negative emotion, negative emodiversity, and the interaction. Results suggested suppression was occuring in terms of the magnitude of the beta coefficient increasing for mean negative emotion and the sign flipping for negative emodiversity. Model comparisons suggested that model5d with only mean negative emotion, provided the best fit to these data.
  • model 6 hierarchical series utilizing the Gini coefficient, investigating (a) mean positive and mean negative emotion, (b) mean positive emotion, mean negative emotion, positive emodiversity, and negative emodiversity, and (c) mean positive emotion, mean negative emotion, positive emodiversity, and negative emodiversity, and the interactions. Results suggested suppression was occuring in terms of the magnitude of the beta coefficient increasing for mean negative emotion and the sign flipping for negative emodiversity. Model comparisons suggested that model6a including only mean positive and mean negative emotion, provided the best fit to these data.
  • model 6 hierarchical series utilizing Shannon’s entropy, investigating (d) mean positive and mean negative emotion, (e) mean positive emotion, mean negative emotion, positive emodiversity, and negative emodiversity, and (f) mean positive emotion, mean negative emotion, positive emodiversity, and negative emodiversity, and the interactions. Results suggested suppression was occuring in terms of the magnitude of the beta coefficient increasing for mean negative emotion and the sign flipping for negative emodiversity. Model comparisons suggested that model6d including only mean positive and mean negative emotion, provided the best fit to these data.

Set-up

Read-in Data

zautra_full <- read.csv("6_EmoDiversity_Emotion and Biomarker & Followup Data_Zautra_2016-12-09.csv",header=TRUE)

Load packages:

library(psych)
library(ggplot2)
library(lavaan)
library(foreign)
library(car)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(ineq)
library(reshape)

Subset of individuals who have 6 or more days of data (for diary analyses)

zautra <- zautra_full[which(zautra_full$days >= 6),] 

Subsets for the single occasion data analyses

Reduce dataset so same number of observations used in all models (for model fit comparison purposes later on) - this is for single occasion data using just the positive emotion items:

zautra2 <- zautra_full[which(!is.na(zautra_full$panaspos_c_mean_4wk) & 
                                 !is.na(zautra_full$panaspos_c_gini_4wk)),]

Reduce dataset so same number of observations used in all models (for model fit comparison purposes later on) - this is for single occasion data using just the negative emotion items:

zautra3 <- zautra_full[which(!is.na(zautra_full$panasneg_c_mean_4wk) & 
                                 !is.na(zautra_full$panasneg_c_gini_4wk)),]

Reduce dataset so same number of observations used in all models (for model fit comparison purposes later on) - this is for single occasion data using both positive and negative emotion items:

zautra4 <- zautra_full[which(!is.na(zautra_full$panaspos_c_mean_4wk) & 
                                 !is.na(zautra_full$panaspos_c_gini_4wk) &
                                 !is.na(zautra_full$panasneg_c_mean_4wk) & 
                                 !is.na(zautra_full$panasneg_c_gini_4wk)),]

Regression Models

Please see Benson, Ram, Almeida, Zautra, & Ong (2017) for further details on the individual emotion items used to calculated positive and negative emodiversity, mean positive and negative emotion, and physical health.

Center Variables

# Center diary data variables
zautra$pos_b_gini.c <- scale(zautra$pos_b_gini,center=TRUE, scale=FALSE)
zautra$neg_b_gini.c <- scale(zautra$neg_b_gini,center=TRUE, scale=FALSE)

zautra$pos_c_mean.c <- scale(zautra$pos_c_mean,center=TRUE, scale=FALSE)
zautra$neg_c_mean.c <- scale(zautra$neg_c_mean,center=TRUE, scale=FALSE)

# Center single occasion data variables - participant reflecting over past 4 weeks
zautra2$panaspos_c_gini_4wk.c <- scale(zautra2$panaspos_c_gini_4wk,center=TRUE, scale=FALSE)
zautra4$panaspos_c_gini_4wk.c <- scale(zautra4$panaspos_c_gini_4wk,center=TRUE, scale=FALSE)
zautra3$panasneg_c_gini_4wk.c <- scale(zautra3$panasneg_c_gini_4wk,center=TRUE, scale=FALSE)
zautra4$panasneg_c_gini_4wk.c <- scale(zautra4$panasneg_c_gini_4wk,center=TRUE, scale=FALSE)

zautra2$panaspos_c_shannon_4wk.c <- scale(zautra2$panaspos_c_shannon_4wk,center=TRUE, scale=FALSE)
zautra4$panaspos_c_shannon_4wk.c <- scale(zautra4$panaspos_c_shannon_4wk,center=TRUE, scale=FALSE)
zautra3$panasneg_c_shannon_4wk.c <- scale(zautra3$panasneg_c_shannon_4wk,center=TRUE, scale=FALSE)
zautra4$panasneg_c_shannon_4wk.c <- scale(zautra4$panasneg_c_shannon_4wk,center=TRUE, scale=FALSE)

zautra2$panaspos_c_mean_4wk.c <- scale(zautra2$panaspos_c_mean_4wk,center=TRUE, scale=FALSE)
zautra4$panaspos_c_mean_4wk.c <- scale(zautra4$panaspos_c_mean_4wk,center=TRUE, scale=FALSE)
zautra3$panasneg_c_mean_4wk.c <- scale(zautra3$panasneg_c_mean_4wk,center=TRUE, scale=FALSE)
zautra4$panasneg_c_mean_4wk.c <- scale(zautra4$panasneg_c_mean_4wk,center=TRUE, scale=FALSE)

Mean Positive Emotion and Positive Emodiversity Predicting Physical Health - Diary Data

Descritpives:

describe(zautra[,c("physhealth","pos_c_mean","pos_b_gini")])
##            vars   n  mean    sd median trimmed   mad   min    max range
## physhealth    1 141 77.74 21.33  84.44   81.02 16.40 12.06 100.00 87.94
## pos_c_mean    2 175  2.09  0.77   2.10    2.09  0.69  0.20   3.96  3.76
## pos_b_gini    3 175  0.92  0.11   0.97    0.95  0.04  0.39   1.00  0.61
##             skew kurtosis   se
## physhealth -1.23     0.78 1.80
## pos_c_mean -0.04    -0.08 0.06
## pos_b_gini -2.37     5.97 0.01

Plot bivariate associations:

pairs.panels(zautra[,c("physhealth","pos_c_mean","pos_b_gini")], 
             lm = TRUE, hist = "gray")

Fit regression model:

mod1a <- lm(physhealth ~ pos_c_mean.c,
           data = zautra)
summary(mod1a)
## 
## Call:
## lm(formula = physhealth ~ pos_c_mean.c, data = zautra)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.543 -11.753   4.911  13.625  35.336 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    77.655      1.647   47.16  < 2e-16 ***
## pos_c_mean.c   11.302      2.153    5.25 5.56e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.55 on 139 degrees of freedom
##   (34 observations deleted due to missingness)
## Multiple R-squared:  0.1655, Adjusted R-squared:  0.1595 
## F-statistic: 27.57 on 1 and 139 DF,  p-value: 5.557e-07
mod1a_adjrsq <- summary(mod1a)$adj.r.squared
mod1a_beta1 <- summary(mod1a)$coefficients[2]

Fit regression model:

mod1b <- lm(physhealth ~ pos_c_mean.c + pos_b_gini.c,
           data = zautra)
summary(mod1b)
## 
## Call:
## lm(formula = physhealth ~ pos_c_mean.c + pos_b_gini.c, data = zautra)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -59.01 -11.62   4.58  12.80  44.28 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    77.641      1.644   47.22   <2e-16 ***
## pos_c_mean.c    8.619      3.101    2.78   0.0062 ** 
## pos_b_gini.c   25.322     21.098    1.20   0.2321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.52 on 138 degrees of freedom
##   (34 observations deleted due to missingness)
## Multiple R-squared:  0.1741, Adjusted R-squared:  0.1622 
## F-statistic: 14.55 on 2 and 138 DF,  p-value: 1.85e-06
mod1b_adjrsq <- summary(mod1b)$adj.r.squared
mod1b_beta1 <- summary(mod1b)$coefficients[2]
mod1b_beta2 <- summary(mod1b)$coefficients[3]

Fit regression model:

mod1c <- lm(physhealth ~ pos_c_mean.c + pos_b_gini.c + I(pos_c_mean.c*pos_b_gini.c),
           data = zautra)
summary(mod1c)
## 
## Call:
## lm(formula = physhealth ~ pos_c_mean.c + pos_b_gini.c + I(pos_c_mean.c * 
##     pos_b_gini.c), data = zautra)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.718 -10.598   5.938  12.918  43.723 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      79.086      2.369  33.384  < 2e-16 ***
## pos_c_mean.c                     10.030      3.521   2.848  0.00507 ** 
## pos_b_gini.c                     -2.527     39.037  -0.065  0.94848    
## I(pos_c_mean.c * pos_b_gini.c)  -22.882     26.975  -0.848  0.39777    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.54 on 137 degrees of freedom
##   (34 observations deleted due to missingness)
## Multiple R-squared:  0.1784, Adjusted R-squared:  0.1605 
## F-statistic: 9.919 on 3 and 137 DF,  p-value: 5.82e-06
mod1c_adjrsq <- summary(mod1c)$adj.r.squared
mod1c_beta1 <- summary(mod1c)$coefficients[2]
mod1c_beta2 <- summary(mod1c)$coefficients[3]

Model fit comparisons:

anova(mod1a,mod1b); fit1 <- anova(mod1a,mod1b)[6][[2,1]]
## Analysis of Variance Table
## 
## Model 1: physhealth ~ pos_c_mean.c
## Model 2: physhealth ~ pos_c_mean.c + pos_b_gini.c
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1    139 53146                           
## 2    138 52597  1    549.03 1.4405 0.2321
anova(mod1b,mod1c); fit2 <- anova(mod1b,mod1c)[6][[2,1]]
## Analysis of Variance Table
## 
## Model 1: physhealth ~ pos_c_mean.c + pos_b_gini.c
## Model 2: physhealth ~ pos_c_mean.c + pos_b_gini.c + I(pos_c_mean.c * pos_b_gini.c)
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1    138 52597                           
## 2    137 52322  1     274.8 0.7195 0.3978

Summary:

The adjusted \(R^2\) for

  • mod1a = 0.1595
  • mod1b = 0.1622
  • mod1c = 0.1605
  • Indicating an increase of 0.0027 from mod1a to mod1b
  • And an additional increase of -0.0017 from mod1b to mod1c in the overall variance explained.

Model fit comparisons:

  • The change in fit from mod1a to mod1b was 0.232
  • The change in fit from mod1b to mod1c was 0.398

The beta coefficient for mean positive emotion was

  • 11.3021 for mod1a
  • 8.6193 for mod1b
  • 10.03 for mod1c
  • In this case, the beta coefficient did not increase from mod1a to mod1b

The beta coefficient for positive emodiversity was

  • 25.3217 for mod1b
  • -2.5272 for mod1c
  • The association between positive emodiversity and physical health remained the same (in terms of sign direction) when comparing the bivariate correlation to the regression model results

Conclusions:

  • Given there was no sign flip for positive emodiversity and no increase in the beta coefficient value for mean positive emotion, it does not seem that statistical suppression is occurring here. The model fit information also suggests that the mod1a (just mean positive emotion) provides the most best/ most parsimonious fit to these data.

Mean Negative Emotion and Negative Emodiversity Predicting Physical Health - Diary Data

Descritpives:

describe(zautra[,c("physhealth","neg_c_mean","neg_b_gini")])
##            vars   n  mean    sd median trimmed   mad   min    max range
## physhealth    1 141 77.74 21.33  84.44   81.02 16.40 12.06 100.00 87.94
## neg_c_mean    2 173  0.38  0.38   0.26    0.31  0.21  0.01   2.51  2.50
## neg_b_gini    3 173  0.49  0.18   0.48    0.48  0.19  0.12   0.95  0.83
##             skew kurtosis   se
## physhealth -1.23     0.78 1.80
## neg_c_mean  2.69     9.66 0.03
## neg_b_gini  0.13    -0.50 0.01

Plot bivariate associations:

pairs.panels(zautra[,c("physhealth","neg_c_mean","neg_b_gini")], 
             lm = TRUE, hist = "gray")