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.
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.
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).
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(emmeans)
# #optional packages:
#library(psych)
library(car) # for ANCOVA
#library(compute.es) # for ANCOVA
#library(lsmeans) # for ANCOVA
library(effectsize)# 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)df_clean <- d %>% select(voice, plev, vsex,
starts_with("pitch"),
starts_with("intense"),
starts_with("form"))
df_clean## # A tibble: 161 × 27
## voice plev vsex pitch_smean pitch_svar pitch_rmean pitch_rvar pitch_rmeanMD
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 -1 116. 858. 109. 2482. -40.7
## 2 2 1 -1 146. 1404. 133. 1231. -16.4
## 3 3 1 -1 106. 706. 108. 229. -41.3
## 4 4 1 -1 102. 450. 108. 1744. -42.0
## 5 5 1 -1 122. 478. 115. 946. -34.5
## 6 6 1 -1 137. 623. 127. 420. -22.7
## 7 7 1 -1 106. 313. 115. 2267. -34.9
## 8 8 -1 -1 105. 1120. 101. 2898. -48.5
## 9 9 1 -1 161. 444. 162. 697. 12.2
## 10 10 -1 -1 132. 640. 118. 185. -31.2
## # ℹ 151 more rows
## # ℹ 19 more variables: pitch_rvarMD <dbl>, pitch_smeanMD <dbl>,
## # pitch_svarMD <dbl>, intense_smean <dbl>, intense_svar <dbl>,
## # intense_rmean <dbl>, intense_rvar <dbl>, intense_rmeanMD <dbl>,
## # intense_rvarMD <dbl>, intense_smeanMD <dbl>, intense_svarMD <dbl>,
## # form_smean <dbl>, form_svar <dbl>, form_rmean <dbl>, form_rvar <dbl>,
## # formant_rmeanMD <dbl>, formant_rvarMD <dbl>, formant_smeanMD <dbl>, …table4 <- list()
ancova_df <- list()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:
## FOR PITCH MEAN
key1 <- "pitch mean"
row1 <- list()
m1 <- lm(pitch_smean ~ plev * vsex + pitch_rmean,data = df_clean)
ancova1 <- car::Anova(m1, type = "III")
adj_means1 <- emmeans::emmeans(m1, ~plev,data=df_clean)
ancova_df[[key1]] <- ancova1
es1 <- effectsize::eta_squared(ancova1, partial = F)
row1[["effect size"]] <- round(es1$Eta2[1], 2)
es1_partial <- effectsize::eta_squared(ancova1, partial = T)
row1[["partial effect size"]] <- round(es1_partial$Eta2_partial[1], 2)
high_adj_mean1 <- summary(adj_means1)$emmean[2]
low_adj_mean1 <- summary(adj_means1)$emmean[1]
row1[["high-rank condition"]] <- high_adj_mean1
row1[["low-rank condition"]] <- low_adj_mean1
table4[[key1]] <- row1
# FOR PITCH VARIABILITY
key2 <- "pitch variability"
row2 <- list()
m2 <- lm(pitch_svar ~ plev * vsex + pitch_rvar,data = df_clean)
ancova2 <- car::Anova(m2, type = "III")
adj_means2 <- emmeans::emmeans(m2, ~plev,data=df_clean)
ancova_df[[key2]] <- ancova2
es2 <- effectsize::eta_squared(ancova2, partial = F)
row2[["effect size"]] <-round(es2$Eta2[1], 2)
es2_partial <- effectsize::eta_squared(ancova2, partial = T)
row2[["partial effect size"]] <-round(es2_partial$Eta2_partial[1], 2)
high_adj_mean2 <- summary(adj_means2)$emmean[2]
low_adj_mean2 <- summary(adj_means2)$emmean[1]
row2[["high-rank condition"]] <- high_adj_mean2
row2[["low-rank condition"]] <- low_adj_mean2
table4[[key2]] <- row2
# FOR LOUDNESS/intensity
key3 <- "loudness"
row3 <- list()
m3 <- lm(intense_smean ~ plev * vsex + intense_rmean,data = df_clean)
ancova3 <- car::Anova(m3, type = "III")
adj_means3 <- emmeans::emmeans(m3, ~plev,data=df_clean)
ancova_df[[key3]] <- ancova3
es3 <- effectsize::eta_squared(ancova3, partial = F)
row3[["effect size"]] <-round(es3$Eta2[1], 2)
es3_partial <- effectsize::eta_squared(ancova3, partial = T)
row3[["partial effect size"]] <-round(es3_partial$Eta2_partial[1], 2)
high_adj_mean3 <- summary(adj_means3)$emmean[2]
low_adj_mean3 <- summary(adj_means3)$emmean[1]
row3[["high-rank condition"]] <- high_adj_mean3
row3[["low-rank condition"]] <- low_adj_mean3
table4[[key3]] <- row3
# FOR LOUDNESS/intensity variability
key4 <- "loudness variability"
row4 <- list()
m4 <- lm(intense_svar ~ plev * vsex + intense_rvar,data = df_clean)
ancova4 <- car::Anova(m4, type = "III")
adj_means4 <- emmeans::emmeans(m4, ~plev,data=df_clean)
ancova_df[[key4]] <- ancova4
es4 <- effectsize::eta_squared(ancova4, partial = F)
row4[["effect size"]] <-round(es4$Eta2[1], 2)
es4_partial <- effectsize::eta_squared(ancova4, partial = T)
row4[["partial effect size"]] <-round(es4_partial$Eta2_partial[1], 2)
high_adj_mean4 <- summary(adj_means4)$emmean[2]
low_adj_mean4 <- summary(adj_means4)$emmean[1]
row4[["high-rank condition"]] <- high_adj_mean4
row4[["low-rank condition"]] <- low_adj_mean4
table4[[key4]] <- row4
# FOR RESONANCE/FORM
key5 <- "resonance"
row5 <- list()
m5 <- lm(form_smean ~ plev * vsex + form_rmean,data = df_clean)
ancova5 <- car::Anova(m5, type = "III")
adj_means5 <- emmeans::emmeans(m5, ~plev,data=df_clean)
ancova_df[[key5]] <- ancova5
es5 <- effectsize::eta_squared(ancova5, partial = F)
row5[["effect size"]] <-round(es5$Eta2[1], 2)
es5_partial <- effectsize::eta_squared(ancova5, partial = T)
row5[["partial effect size"]] <-round(es5_partial$Eta2_partial[1], 2)
high_adj_mean5 <- summary(adj_means5)$emmean[2]
low_adj_mean5 <- summary(adj_means5)$emmean[1]
row5[["high-rank condition"]] <- high_adj_mean5
row5[["low-rank condition"]] <- low_adj_mean5
table4[[key5]] <- row5
# FOR RESONANCE/FORM variability
key6 <- "resonance variability"
row6 <- list()
m6 <- lm(form_svar ~ plev * vsex + form_rvar,data = df_clean)
ancova6 <- car::Anova(m6, type = "III")
adj_means6 <- emmeans::emmeans(m6, ~plev,data=df_clean)
ancova_df[[key6]] <- ancova6
es6 <- effectsize::eta_squared(ancova6, partial = F)
row6[["effect size"]] <-round(es6$Eta2[1], 2)
es6_partial <- effectsize::eta_squared(ancova6, partial = T)
row6[["partial effect size"]] <-round(es6_partial$Eta2_partial[1], 2)
high_adj_mean6 <- summary(adj_means6)$emmean[2]
low_adj_mean6 <- summary(adj_means6)$emmean[1]
row6[["high-rank condition"]] <- high_adj_mean6
row6[["low-rank condition"]] <- low_adj_mean6
table4[[key6]] <- row6# printing the results now:
table4[key1]## $`pitch mean`
## $`pitch mean`$`effect size`
## [1] 0.01
##
## $`pitch mean`$`partial effect size`
## [1] 0.03
##
## $`pitch mean`$`high-rank condition`
## [1] 158.6098
##
## $`pitch mean`$`low-rank condition`
## [1] 155.5227table4[key2]## $`pitch variability`
## $`pitch variability`$`effect size`
## [1] 0.02
##
## $`pitch variability`$`partial effect size`
## [1] 0.03
##
## $`pitch variability`$`high-rank condition`
## [1] 1425.016
##
## $`pitch variability`$`low-rank condition`
## [1] 1648.367table4[key3]## $loudness
## $loudness$`effect size`
## [1] 0.01
##
## $loudness$`partial effect size`
## [1] 0.01
##
## $loudness$`high-rank condition`
## [1] 59.33567
##
## $loudness$`low-rank condition`
## [1] 58.66784table4[key4]## $`loudness variability`
## $`loudness variability`$`effect size`
## [1] 0.02
##
## $`loudness variability`$`partial effect size`
## [1] 0.03
##
## $`loudness variability`$`high-rank condition`
## [1] 196.7301
##
## $`loudness variability`$`low-rank condition`
## [1] 183.4795table4[key5]## $resonance
## $resonance$`effect size`
## [1] 0
##
## $resonance$`partial effect size`
## [1] 0
##
## $resonance$`high-rank condition`
## [1] 1129.384
##
## $resonance$`low-rank condition`
## [1] 1128.806table4[key6]## $`resonance variability`
## $`resonance variability`$`effect size`
## [1] 0
##
## $`resonance variability`$`partial effect size`
## [1] 0
##
## $`resonance variability`$`high-rank condition`
## [1] 42170.78
##
## $`resonance variability`$`low-rank condition`
## [1] 43654.54print("all of the results matched the table")## [1] "all of the results matched the table"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
print_output <- function(key,rownum) {# key is a string variable, row num is an int that signifies row in table 4
print(sprintf("results for %s : ",key))
currancova <- ancova_df[[key]]
F <- currancova[3,3]
p_val <- currancova[3,4]
cat(paste0(sprintf("\nF(%d, %d) = %.2f p = %.3f\n", currancova[3,2], currancova[6,2], F, p_val)))
}
# for pitch mean
print_output(key1,1)## [1] "results for pitch mean : "
##
## F(1, 156) = 4.67 p = 0.032# for pitch var
print_output(key2,2)## [1] "results for pitch variability : "
##
## F(1, 156) = 61.76 p = 0.000# for loudness mean
print_output(key3,3)## [1] "results for loudness : "
##
## F(1, 156) = 12.81 p = 0.000# for loudness var
print_output(key4,4)## [1] "results for loudness variability : "
##
## F(1, 156) = 19.79 p = 0.000# for resonance mean
print_output(key5,5)## [1] "results for resonance : "
##
## F(1, 156) = 11.54 p = 0.001# for resonance var
print_output(key6,6)## [1] "results for resonance variability : "
##
## F(1, 156) = 0.24 p = 0.622Were 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 most of it, but for part 2 some of the results were incorrect, and I coudn’t figure out why even though I spent up to 3 hours on it.
How difficult was it to reproduce your results?
The code was okay, just took some googling of different statistical R functions (like in group A), but trying to debug it/figure out why part 2 was wrong was hard.
What aspects made it difficult? What aspects made it easy?
Debugging was difficult because I’m unsure about where I went wrong.Also this dataset was a bit confusing/hard to understand at first since the variables were strangely named and it was so big.