For this exercise, please try to reproduce the results from Experiment 1 of the associated paper (Ko, Sadler & Galinsky, 2015). The PDF of the paper is included in the same folder as this Rmd file.

Methods summary:

A sense of power has often been tied to how we perceive each other’s voice. Social hierarchy is embedded into the structure of society and provides a metric by which others relate to one another. In 1956, the Brunswik Lens Model was introduced to examine how vocal cues might influence hierarchy. In “The Sound of Power: Conveying and Detecting Hierarchical Rank Through Voice,” Ko and colleagues investigated how manipulation of hierarchical rank within a situation might impact vocal acoustic cues. Using the Brunswik Model, six acoustic metrics were utilized (pitch mean & variability, loudness mean & variability, and resonance mean & variability) to isolate a potential contribution between individuals of different hierarchal rank. In the first experiment, Ko, Sadler & Galinsky examined the vocal acoustic cues of individuals before and after being assigned a hierarchal rank in a sample of 161 subjects (80 male). Each of the six hierarchy acoustic cues were analyzed with a 2 (high vs. low rank condition) x 2 (male vs. female) analysis of covariance, controlling for the baseline of the respective acoustic cue.


Target outcomes:

Below is the specific result you will attempt to reproduce (quoted directly from the results section of Experiment 1):

The impact of hierarchical rank on speakers’ acoustic cues. Each of the six hierarchy-based (i.e., postmanipulation) acoustic variables was submitted to a 2 (condition: high rank, low rank) × 2 (speaker’s sex: female, male) between-subjects analysis of covariance, controlling for the corresponding baseline acoustic variable. Table 4 presents the adjusted means by condition. Condition had a significant effect on pitch, pitch variability, and loudness variability. Speakers’ voices in the high-rank condition had higher pitch, F(1, 156) = 4.48, p < .05; were more variable in loudness, F(1, 156) = 4.66, p < .05; and were more monotone (i.e., less variable in pitch), F(1, 156) = 4.73, p < .05, compared with speakers’ voices in the low-rank condition (all other Fs < 1; see the Supplemental Material for additional analyses of covariance involving pitch and loudness). (from Ko et al., 2015, p. 6; emphasis added)

The adjusted means for these analyses are reported in Table 4 (Table4_AdjustedMeans.png, included in the same folder as this Rmd file).


Step 1: Load packages

library(tidyverse) # for data munging
library(knitr) # for kable table formating
library(haven) # import and export 'SPSS', 'Stata' and 'SAS' Files
library(readxl) # import excel files

#optional packages:
library(psych)
library(car) # for ANCOVA
library(compute.es) # for ANCOVA
library(lsmeans) # for ANCOVA

Step 2: Load data

# Just Experiment 1
d <-read_csv("data/S1_voice_level_Final.csv")
# DT::datatable(d)

Step 3: Tidy data

# Remove columns not needed for reproducibility analyses 
d <- d %>% 
  select(plev, vsex, pitch_smean, pitch_rmean, 
         pitch_svar, pitch_rvar, intense_smean,
         intense_rmean, intense_svar, intense_rvar, 
         form_smean, form_rmean, form_svar, 
         form_rvar)

Step 4: Run analysis

Pre-processing

Inferential statistics

The impact of hierarchical rank on speakers’ acoustic cues. Each of the six hierarchy-based (i.e., postmanipulation) acoustic variables was submitted to a 2 (condition: high rank, low rank) × 2 (speaker’s sex: female, male) between-subjects analysis of covariance, controlling for the corresponding baseline acoustic variable. […] Condition had a significant effect on pitch, pitch variability, and loudness variability. Speakers’ voices in the high-rank condition had higher pitch, F(1, 156) = 4.48, p < .05; were more variable in loudness, F(1, 156) = 4.66, p < .05; and were more monotone (i.e., less variable in pitch), F(1, 156) = 4.73, p < .05, compared with speakers’ voices in the low-rank condition (all other Fs < 1; see the Supplemental Material for additional analyses of covariance involving pitch and loudness).

# reproduce the above results here

model_pitch_mean <- lm(pitch_smean ~ plev*vsex + pitch_rmean, data = d)
car::Anova(model_pitch_mean, type = 3)
## Anova Table (Type III tests)
## 
## Response: pitch_smean
##             Sum Sq  Df  F value  Pr(>F)    
## (Intercept)    359   1   4.2359 0.04124 *  
## plev           380   1   4.4837 0.03581 *  
## vsex           396   1   4.6687 0.03224 *  
## pitch_rmean  43334   1 511.2208 < 2e-16 ***
## plev:vsex      222   1   2.6156 0.10783    
## Residuals    13223 156                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F(1, 156) = 4.48, p < .05

model_pitch_var <- lm(pitch_svar ~ plev*vsex + pitch_rvar, data = d)
car::Anova(model_pitch_var, type = 3)
## Anova Table (Type III tests)
## 
## Response: pitch_svar
##               Sum Sq  Df  F value    Pr(>F)    
## (Intercept) 45439310   1 107.1878 < 2.2e-16 ***
## plev         2003882   1   4.7270    0.0312 *  
## vsex        26179681   1  61.7558 5.939e-13 ***
## pitch_rvar  27760551   1  65.4850 1.548e-13 ***
## plev:vsex      37857   1   0.0893    0.7655    
## Residuals   66131879 156                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F(1, 156) = 4.73, p < .05

model_loud_mean <- lm(intense_smean ~ plev*vsex + intense_rmean, data = d)
car::Anova(model_loud_mean, type = 3)
## Anova Table (Type III tests)
## 
## Response: intense_smean
##                Sum Sq  Df  F value    Pr(>F)    
## (Intercept)    208.90   1  25.7761  1.08e-06 ***
## plev            17.93   1   2.2118 0.1389771    
## vsex           103.84   1  12.8131 0.0004591 ***
## intense_rmean 1195.54   1 147.5165 < 2.2e-16 ***
## plev:vsex       11.69   1   1.4424 0.2315772    
## Residuals     1264.30 156                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F(1, 156) > 1 (different)

model_loud_var <- lm(intense_svar ~ plev*vsex + intense_rvar, data = d)
car::Anova(model_loud_var, type = 3)
## Anova Table (Type III tests)
## 
## Response: intense_svar
##              Sum Sq  Df  F value    Pr(>F)    
## (Intercept)   61767   1  40.8930 1.777e-09 ***
## plev           7042   1   4.6620   0.03236 *  
## vsex          29895   1  19.7923 1.634e-05 ***
## intense_rvar 184927   1 122.4323 < 2.2e-16 ***
## plev:vsex        40   1   0.0262   0.87160    
## Residuals    235630 156                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F(1, 156) = 4.66, p < .05

model_resonance_mean <- lm(form_smean ~ plev*vsex + form_rmean, data = d)
car::Anova(model_resonance_mean, type = 3)
## Anova Table (Type III tests)
## 
## Response: form_smean
##              Sum Sq  Df F value    Pr(>F)    
## (Intercept)   35300   1  5.0525 0.0259916 *  
## plev             13   1  0.0019 0.9651132    
## vsex          80646   1 11.5431 0.0008630 ***
## form_rmean    89319   1 12.7845 0.0004656 ***
## plev:vsex       170   1  0.0243 0.8763623    
## Residuals   1089894 156                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F(1, 156) < 1

model_resonance_var <- lm(form_svar ~ plev*vsex + form_rvar, data = d)
car::Anova(model_resonance_var, type = 3)
## Anova Table (Type III tests)
## 
## Response: form_svar
##                 Sum Sq  Df F value    Pr(>F)    
## (Intercept) 5.3682e+09   1 28.8227 2.826e-07 ***
## plev        8.7633e+07   1  0.4705    0.4938    
## vsex        4.5357e+07   1  0.2435    0.6224    
## form_rvar   4.9801e+08   1  2.6739    0.1040    
## plev:vsex   1.5279e+07   1  0.0820    0.7749    
## Residuals   2.9055e+10 156                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F(1, 156) < 1

Descriptive statistics

In the paper, the adjusted means by condition are reported (see Table 4, or Table4_AdjustedMeans.png, included in the same folder as this Rmd file). Reproduce these values below:

emmeans(model_pitch_mean, ~ plev) # -1 = 156, 1 = 159 
##  plev emmean   SE  df lower.CL upper.CL
##    -1    156 1.04 156      153      158
##     1    159 1.02 156      157      161
## 
## Results are averaged over the levels of: vsex 
## Confidence level used: 0.95
es_model_pitch_mean <- effectsize::eta_squared(model_pitch_mean)
print(es_model_pitch_mean, digits = 6) # eta squared = .02 (.03)
## # Effect Size for ANOVA (Type I)
## 
## Parameter   | Eta2 (partial) |               95% CI
## ---------------------------------------------------
## plev        |       0.015497 | [0.000000, 1.000000]
## vsex        |       0.949046 | [0.938005, 1.000000]
## pitch_rmean |       0.765403 | [0.716788, 1.000000]
## plev:vsex   |       0.016490 | [0.000000, 1.000000]
## 
## - One-sided CIs: upper bound fixed at [1.000000].
emmeans(model_pitch_var, ~ plev) # -1 = 1648, 1 = 1425
##  plev emmean   SE  df lower.CL upper.CL
##    -1   1648 73.3 156     1504     1793
##     1   1425 71.9 156     1283     1567
## 
## Results are averaged over the levels of: vsex 
## Confidence level used: 0.95
es_model_pitch_var <- effectsize::eta_squared(model_pitch_var)
print(es_model_pitch_var, digits = 6) # eta squared = .04 (.03)
## # Effect Size for ANOVA (Type I)
## 
## Parameter  | Eta2 (partial) |               95% CI
## --------------------------------------------------
## plev       |       0.044022 | [0.006510, 1.000000]
## vsex       |       0.364350 | [0.269042, 1.000000]
## pitch_rvar |       0.298469 | [0.204051, 1.000000]
## plev:vsex  |       0.000572 | [0.000000, 1.000000]
## 
## - One-sided CIs: upper bound fixed at [1.000000].
emmeans(model_loud_mean, ~ plev) # -1 = 58.7, 1 = 59.3
##  plev emmean    SE  df lower.CL upper.CL
##    -1   58.7 0.321 156     58.0     59.3
##     1   59.3 0.314 156     58.7     60.0
## 
## Results are averaged over the levels of: vsex 
## Confidence level used: 0.95
es_model_loud_mean <- effectsize::eta_squared(model_loud_mean)
print(es_model_loud_mean, digits = 6) # eta squared = .02 (.01)
## # Effect Size for ANOVA (Type I)
## 
## Parameter     | Eta2 (partial) |               95% CI
## -----------------------------------------------------
## plev          |       0.015281 | [0.000000, 1.000000]
## vsex          |       0.000709 | [0.000000, 1.000000]
## intense_rmean |       0.492203 | [0.403709, 1.000000]
## plev:vsex     |       0.009161 | [0.000000, 1.000000]
## 
## - One-sided CIs: upper bound fixed at [1.000000].
emmeans(model_loud_var, ~ plev) # -1 = 183, 1 = 197 
##  plev emmean   SE  df lower.CL upper.CL
##    -1    183 4.38 156      175      192
##     1    197 4.30 156      188      205
## 
## Results are averaged over the levels of: vsex 
## Confidence level used: 0.95
es_model_loud_var <- effectsize::eta_squared(model_loud_var)
print(es_model_loud_var, digits = 6) # eta squared = .05 (.03)
## # Effect Size for ANOVA (Type I)
## 
## Parameter    | Eta2 (partial) |               95% CI
## ----------------------------------------------------
## plev         |       0.046978 | [0.007739, 1.000000]
## vsex         |       0.052667 | [0.010275, 1.000000]
## intense_rvar |       0.439772 | [0.347291, 1.000000]
## plev:vsex    |       0.000168 | [0.000000, 1.000000]
## 
## - One-sided CIs: upper bound fixed at [1.000000].
emmeans(model_resonance_mean, ~ plev) # -1 = 1129, 1 = 1129 
##  plev emmean   SE  df lower.CL upper.CL
##    -1   1129 9.42 156     1110     1147
##     1   1129 9.24 156     1111     1148
## 
## Results are averaged over the levels of: vsex 
## Confidence level used: 0.95
es_model_resonance_mean <- effectsize::eta_squared(model_resonance_mean)
print(es_model_resonance_mean, digits = 6) # eta squared = .00
## # Effect Size for ANOVA (Type I)
## 
## Parameter  | Eta2 (partial) |               95% CI
## --------------------------------------------------
## plev       |       0.000524 | [0.000000, 1.000000]
## vsex       |       0.090168 | [0.031183, 1.000000]
## form_rmean |       0.075616 | [0.022352, 1.000000]
## plev:vsex  |       0.000156 | [0.000000, 1.000000]
## 
## - One-sided CIs: upper bound fixed at [1.000000].
emmeans(model_resonance_var, ~ plev) # -1 = 43655, 1 = 42171 
##  plev emmean   SE  df lower.CL upper.CL
##    -1  43655 1540 156    40611    46698
##     1  42171 1510 156    39187    45155
## 
## Results are averaged over the levels of: vsex 
## Confidence level used: 0.95
es_model_resonance_var <- effectsize::eta_squared(model_resonance_var) 
print(es_model_resonance_var, digits = 6) # eta squared = .00
## # Effect Size for ANOVA (Type I)
## 
## Parameter | Eta2 (partial) |               95% CI
## -------------------------------------------------
## plev      |       0.004550 | [0.000000, 1.000000]
## vsex      |       0.000446 | [0.000000, 1.000000]
## form_rvar |       0.016509 | [0.000000, 1.000000]
## plev:vsex |       0.000526 | [0.000000, 1.000000]
## 
## - One-sided CIs: upper bound fixed at [1.000000].

Step 5: Reflection

Were you able to reproduce the results you attempted to reproduce? If not, what part(s) were you unable to reproduce?

I was able to reproduce all adjusted means and all F and p values in the inferential statistics section, apart from loudness mean model, where F was greater than 1, despite the paper reporting that all models apart from the pitch mean, pitch variability, and loudness variability models had an F value of less than 1. I also was unable to reproduce the majority of the eta squared effect sizes for the models; only those for the resonance mean and resonance variability models were reproduced.

How difficult was it to reproduce your results?

It was somewhat difficult to reproduce the results due to some gaps in my conceptual knowledge. I am still unsure of why I was unable to reproduce the majority of the effect sizes; it’s possible that I made a coding error that I am not currently able to identify or that the authors made a calculation or coding error.

What aspects made it difficult? What aspects made it easy?

The codebook the authors provided helped a lot with building the models in the inferential statistics section. On the other hand, I had to do some separate research to conceptually grasp what adjusted means were. I understood them to come out of ANCOVA models, so I ran the models before the descriptive statistics, which felt less intuitive given the original sequence of this .rmd. I also had to grasp the difference between the types of sums of squares an ANOVA vs. an ANCOVA requires for interpretation. ChatGPT provided me with coding assistance for the Type III ANOVAs as well as for the adjusted means.