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.

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

**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.

**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.

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

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

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

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)),]
```

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 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)
```

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.

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")
```