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 hierarchal 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
library(afex)
# library(emmeans) # for ANCOVA

# #optional packages:
# library(psych)
library(car) # for ANCOVA
# library(compute.es) # for ANCOVA
# library(msmeans) # outdated package replaced by emmeans

Step 2: Load data

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

d
## # A tibble: 161 × 45
##    voice form_smean form_svar form_rmean form_rvar intense_smean intense_svar
##    <dbl>      <dbl>     <dbl>      <dbl>     <dbl>         <dbl>        <dbl>
##  1     1      1043.    38805.      1259.    68275.          65.5        157. 
##  2     2      1083.    35609.      1278.    54612.          65.5        301. 
##  3     3      1092.    25979.      1334.    55175.          53.6         71.2
##  4     4      1268.    71346.      1298.    74340.          60.7        201. 
##  5     5      1071.    38246.      1256.    67846.          60.1        204. 
##  6     6      1095.    40716.      1278.    76674.          60.1        150. 
##  7     7      1060.    38855.      1310.    88947.          60.5        127. 
##  8     8      1051.    28191.      1223.    52682.          57.3        207. 
##  9     9      1038.    34105.      1256.    68508.          57.9        199. 
## 10    10      1346.    69260.      1314.    55400.          60.3        149. 
## # ℹ 151 more rows
## # ℹ 38 more variables: intense_rmean <dbl>, intense_rvar <dbl>,
## #   pitch_smean <dbl>, pitch_svar <dbl>, pitch_rmean <dbl>, pitch_rvar <dbl>,
## #   pow <dbl>, age <dbl>, sex <chr>, race <chr>, native <chr>, feelpower <dbl>,
## #   plev <dbl>, vsex <dbl>, pitch_rmeanMD <dbl>, pitch_rvarMD <dbl>,
## #   intense_rmeanMD <dbl>, intense_rvarMD <dbl>, formant_rmeanMD <dbl>,
## #   formant_rvarMD <dbl>, pitch_smeanMD <dbl>, pitch_svarMD <dbl>, …
# DT::datatable(d)

#Acoustic variable names
#‘form’ = formant disperson/measure of resonance
#‘intense’ = intensity/loudness
#‘pitch’= pitch
#‘_s’ = script, acoustics recorded while saying negotiation script
#‘_r’ = rainbow, acoustics recorded during rainbow passage
#‘mean’ = average value of acoustic during a clip
#‘var’ = variance for acoustic during a clip
#‘MD’ = mean-deviated or centered
#‘Z’ = standardized version of acoustic

Step 3: Tidy data

filtered_d <- d |> 
  select(c("voice","sex","plev","pitch_rmean","pitch_smean","pitch_rvar","pitch_svar","intense_rmean","intense_smean","intense_rvar","intense_svar","form_rmean","form_smean","form_rvar","form_svar"))

# Males recoded to 1, females recoded to 2
# filtered_d$sex <- ifelse(filtered_d$sex == 'M', "1", "2")
# filtered_d$sex <- as.numeric(filtered_d$sex)

filtered_d
## # A tibble: 161 × 15
##    voice sex    plev pitch_rmean pitch_smean pitch_rvar pitch_svar intense_rmean
##    <dbl> <chr> <dbl>       <dbl>       <dbl>      <dbl>      <dbl>         <dbl>
##  1     1 M         1        109.        116.      2482.       858.          58.7
##  2     2 M         1        133.        146.      1231.      1404.          63.9
##  3     3 M         1        108.        106.       229.       706.          53.0
##  4     4 M         1        108.        102.      1744.       450.          59.8
##  5     5 M         1        115.        122.       946.       478.          57.5
##  6     6 M         1        127.        137.       420.       623.          58.4
##  7     7 M         1        115.        106.      2267.       313.          57.9
##  8     8 M        -1        101.        105.      2898.      1120.          57.3
##  9     9 M         1        162.        161.       697.       444.          58.1
## 10    10 M        -1        118.        132.       185.       640.          58.4
## # ℹ 151 more rows
## # ℹ 7 more variables: intense_smean <dbl>, intense_rvar <dbl>,
## #   intense_svar <dbl>, form_rmean <dbl>, form_smean <dbl>, form_rvar <dbl>,
## #   form_svar <dbl>
any(is.na(filtered_d))
## [1] FALSE

Step 4: Run analysis

Pre-processing

# Prepare data for ANCOVA

# Calculate the means for each condition by subsetting the data
# condition_values <- unique(filtered_d$plev)  # Get unique condition values

# Initialize an empty vector to store the means
# condition_means <- numeric(length(condition_values))

# filtered_d['diff_pitchmean'] = filtered_d['pitch_smean'] - filtered_d['pitch_rmean']
# filtered_d['diff_pitchvar'] = filtered_d['pitch_svar'] - filtered_d['pitch_rvar']
# filtered_d['diff_intensemean'] = filtered_d['intense_smean'] - filtered_d['intense_rmean']
# filtered_d['diff_intensevar'] = filtered_d['intense_svar'] - filtered_d['intense_rvar']
# filtered_d['diff_formmean'] = filtered_d['form_smean'] - filtered_d['form_rmean']
# filtered_d['diff_formvar'] = filtered_d['form_svar'] - filtered_d['form_rvar']
# 
# filtered_d |> 
#   group_by(plev) |> 
#   summarize(
#    mean_pitch = mean(diff_pitchmean),
#    mean_pitchvar = mean(diff_pitchvar),
#    mean_intense = mean(diff_intensemean),
#    mean_intensevar = mean(diff_intensevar),
#    mean_pitchs = mean(diff_formmean),
#    mean_pitchvars = mean(diff_formvar)
#   )

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:

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

# Run the ANCOVA for each of the six hierarchy-based acoustic variables
pitchmean_model <- lm(pitch_smean ~ plev * sex + pitch_rmean, data=filtered_d)
ancova_pitch_result <- Anova(pitchmean_model, type="III")

# Fit an ANCOVA model using aov_ez
# pitchmean_model <- aov_ez(id = "voice", dv = "pitch_smean", data = filtered_d, between = c("plev", "sex"), covariate = c("pitch_rmean"))

pitchvar_model <- lm(pitch_svar ~ plev * sex + pitch_rvar, data=filtered_d)

intensemean_model <- lm(intense_smean ~ plev * sex + intense_rmean, data=filtered_d)

intensevar_model <- lm(intense_svar ~ plev * sex + intense_rvar, data=filtered_d)

formmean_model <- lm(form_smean ~ plev * sex + form_rmean, data=filtered_d)

formvar_model <- lm(form_svar ~ plev * sex + form_rvar, data=filtered_d)

# Display the summary of the ANCOVA for each variable
print(ancova_pitch_result) # summary(pitchmean_model) #anova(pitchmean_model)
## Anova Table (Type III tests)
## 
## Response: pitch_smean
##             Sum Sq  Df  F value  Pr(>F)    
## (Intercept)    379   1   4.4744 0.03600 *  
## plev           590   1   6.9640 0.00916 ** 
## sex            396   1   4.6687 0.03224 *  
## pitch_rmean  43334   1 511.2208 < 2e-16 ***
## plev:sex       222   1   2.6156 0.10783    
## Residuals    13223 156                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(pitchvar_model)
## 
## Call:
## lm(formula = pitch_svar ~ plev * sex + pitch_rvar, data = filtered_d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1133.0  -428.0   -94.1   292.2  3651.7 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1344.79878  110.60029  12.159  < 2e-16 ***
## plev        -127.07528   72.56901  -1.751   0.0819 .  
## sexM        -820.49821  104.40916  -7.858 5.94e-13 ***
## pitch_rvar     0.34360    0.04246   8.092 1.55e-13 ***
## plev:sexM     30.80017  103.06722   0.299   0.7655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 651.1 on 156 degrees of freedom
## Multiple R-squared:  0.5111, Adjusted R-squared:  0.4985 
## F-statistic: 40.77 on 4 and 156 DF,  p-value: < 2.2e-16
summary(intensemean_model)
## 
## Call:
## lm(formula = intense_smean ~ plev * sex + intense_rmean, data = filtered_d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.3607 -1.7424 -0.1007  1.5687  6.8850 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   18.25180    3.37665   5.405 2.38e-07 ***
## plev           0.06325    0.31684   0.200 0.842041    
## sexM          -1.66687    0.46567  -3.580 0.000459 ***
## intense_rmean  0.72368    0.05958  12.146  < 2e-16 ***
## plev:sexM      0.54133    0.45074   1.201 0.231577    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.847 on 156 degrees of freedom
## Multiple R-squared:  0.4987, Adjusted R-squared:  0.4858 
## F-statistic:  38.8 on 4 and 156 DF,  p-value: < 2.2e-16
summary(intensevar_model)
## 
## Call:
## lm(formula = intense_svar ~ plev * sex + intense_rvar, data = filtered_d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -126.73  -22.63   -5.05   20.21  130.57 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  57.56283   11.96221   4.812 3.50e-06 ***
## plev          7.12182    4.31875   1.649    0.101    
## sexM         27.53904    6.19014   4.449 1.63e-05 ***
## intense_rvar  0.64991    0.05874  11.065  < 2e-16 ***
## plev:sexM    -0.99304    6.13391  -0.162    0.872    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.86 on 156 degrees of freedom
## Multiple R-squared:  0.4709, Adjusted R-squared:  0.4573 
## F-statistic: 34.71 on 4 and 156 DF,  p-value: < 2.2e-16
summary(formmean_model)
## 
## Call:
## lm(formula = form_smean ~ plev * sex + form_rmean, data = filtered_d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -112.665  -52.619  -25.835   -0.345  229.149 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 458.6088   194.9988   2.352 0.019929 *  
## plev         -0.7387     9.2891  -0.080 0.936722    
## sexM        -45.2304    13.3128  -3.398 0.000863 ***
## form_rmean    0.5362     0.1500   3.576 0.000466 ***
## plev:sexM     2.0558    13.1922   0.156 0.876362    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83.59 on 156 degrees of freedom
## Multiple R-squared:  0.1537, Adjusted R-squared:  0.132 
## F-statistic: 7.082 on 4 and 156 DF,  p-value: 2.892e-05
summary(formvar_model)
## 
## Call:
## lm(formula = form_svar ~ plev * sex + form_rvar, data = filtered_d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17761  -8627  -4289   2159  59762 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.362e+04  5.726e+03   5.870 2.53e-08 ***
## plev        -1.051e+03  1.517e+03  -0.693    0.489    
## sexM        -1.186e+03  2.403e+03  -0.493    0.622    
## form_rvar    1.542e-01  9.431e-02   1.635    0.104    
## plev:sexM    6.183e+02  2.159e+03   0.286    0.775    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13650 on 156 degrees of freedom
## Multiple R-squared:  0.02184,    Adjusted R-squared:  -0.00324 
## F-statistic: 0.8708 on 4 and 156 DF,  p-value: 0.4829
# Calculate the pitch means for each condition
# for (i in 1:length(condition_values)) {
#   condition_means[i] <- mean(filtered_d$pitch_smean[filtered_d$plev == condition_values[i]])
# }
# 
# Create a 1x2 output table
# output_table <- data.frame(
#   "Condition" = condition_values,
#   "Pitch Mean" = condition_means
# )
# 
# Print the output table
# print(output_table)

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

Step 5: Reflection

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

No. I couldn’t even make Table 4 even after 3 hours.

How difficult was it to reproduce your results?

I tried different ways of computing means, including taking the difference between postmanipulation and baseline, as well as running various ANCOVAs with lm, aov, and aov_ez. Even when reading the supplementary materials on OSF, I could not figure out if the ANCOVAs were performed using the raw values, MD values, or Z-transformed values without resorting to trial and error.

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

Difficulties were detailed above, I didn’t find it easy at all.