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

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

Step 2: Load data

## Step 2: Load data
datadir <- "/Users/eu/Documents/GitHub/problem_sets/ps3/Group B/Choice 2/data"

# Just Experiment 1
d <- read_csv(file.path(datadir, "S1_voice_level_Final.csv"))

str(d)
## spc_tbl_ [161 × 45] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ voice          : num [1:161] 1 2 3 4 5 6 7 8 9 10 ...
##  $ form_smean     : num [1:161] 1043 1083 1092 1268 1071 ...
##  $ form_svar      : num [1:161] 38805 35609 25979 71346 38246 ...
##  $ form_rmean     : num [1:161] 1259 1278 1334 1298 1256 ...
##  $ form_rvar      : num [1:161] 68275 54612 55175 74340 67846 ...
##  $ intense_smean  : num [1:161] 65.5 65.5 53.6 60.7 60.1 ...
##  $ intense_svar   : num [1:161] 157.3 301 71.2 201.2 203.7 ...
##  $ intense_rmean  : num [1:161] 58.7 63.9 53 59.8 57.5 ...
##  $ intense_rvar   : num [1:161] 174.1 271.4 87.4 198.8 202.4 ...
##  $ pitch_smean    : num [1:161] 116 146 106 102 122 ...
##  $ pitch_svar     : num [1:161] 858 1404 706 450 478 ...
##  $ pitch_rmean    : num [1:161] 109 133 108 108 115 ...
##  $ pitch_rvar     : num [1:161] 2482 1231 229 1744 946 ...
##  $ pow            : num [1:161] 1 3 3 3 4 1 2 1 1 1 ...
##  $ age            : num [1:161] 19 21 18 22 20 21 18 18 23 19 ...
##  $ sex            : chr [1:161] "M" "M" "M" "M" ...
##  $ race           : chr [1:161] "X" "C" "A" "C" ...
##  $ native         : chr [1:161] "Y" "Y" "Y" "Y" ...
##  $ feelpower      : num [1:161] 6 5 3 6 6 4 6 2 5 3 ...
##  $ plev           : num [1:161] 1 1 1 1 1 1 1 -1 1 -1 ...
##  $ vsex           : num [1:161] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ pitch_rmeanMD  : num [1:161] -40.7 -16.4 -41.3 -42 -34.5 ...
##  $ pitch_rvarMD   : num [1:161] 729.28 -521.37 -1523.73 -7.96 -806.19 ...
##  $ intense_rmeanMD: num [1:161] 1.2404 6.4302 -4.4406 2.2973 0.0353 ...
##  $ intense_rvarMD : num [1:161] -8.66 88.6 -95.39 16.04 19.64 ...
##  $ formant_rmeanMD: num [1:161] -33.4 -14.7 41.2 5.2 -36.8 ...
##  $ formant_rvarMD : num [1:161] 4145 -9519 -8955 10210 3716 ...
##  $ pitch_smeanMD  : num [1:161] -41.1 -11.4 -51.6 -55 -34.6 ...
##  $ pitch_svarMD   : num [1:161] -680 -133 -831 -1088 -1060 ...
##  $ intense_smeanMD: num [1:161] 6.47 6.46 -5.42 1.69 1.05 ...
##  $ intense_svarMD : num [1:161] -32.8 110.8 -119 11 13.6 ...
##  $ formant_smeanMD: num [1:161] -86.4 -46.2 -37 138.3 -57.8 ...
##  $ formant_svarMD : num [1:161] -4107 -7303 -16933 28433 -4667 ...
##  $ Zpitch_rmean   : num [1:161] -0.242 1.353 -0.279 -0.322 0.166 ...
##  $ Zpitch_rvar    : num [1:161] 0.619 -0.195 -0.847 0.139 -0.38 ...
##  $ Zform_rmean    : num [1:161] -0.602 -0.186 1.056 0.256 -0.678 ...
##  $ Zform_rvar     : num [1:161] -0.112 -1.125 -1.083 0.338 -0.143 ...
##  $ Zintense_rmean : num [1:161] 0.0481 1.3553 -1.3828 0.3143 -0.2554 ...
##  $ Zintense_rvar  : num [1:161] -0.0243 1.6525 -1.5195 0.4016 0.4636 ...
##  $ Zpitch_smean   : num [1:161] -0.101 1.683 -0.734 -0.939 0.287 ...
##  $ Zpitch_svar    : num [1:161] -0.229 0.439 -0.415 -0.728 -0.694 ...
##  $ Zform_smean    : num [1:161] -0.684 -0.228 -0.124 1.865 -0.36 ...
##  $ Zform_svar     : num [1:161] -0.28 -0.486 -1.104 1.81 -0.316 ...
##  $ Zintense_smean : num [1:161] 1.593 1.59 -1.307 0.427 0.272 ...
##  $ Zintense_svar  : num [1:161] -0.6972 1.6798 -2.1231 0.028 0.0698 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   voice = col_double(),
##   ..   form_smean = col_double(),
##   ..   form_svar = col_double(),
##   ..   form_rmean = col_double(),
##   ..   form_rvar = col_double(),
##   ..   intense_smean = col_double(),
##   ..   intense_svar = col_double(),
##   ..   intense_rmean = col_double(),
##   ..   intense_rvar = col_double(),
##   ..   pitch_smean = col_double(),
##   ..   pitch_svar = col_double(),
##   ..   pitch_rmean = col_double(),
##   ..   pitch_rvar = col_double(),
##   ..   pow = col_double(),
##   ..   age = col_double(),
##   ..   sex = col_character(),
##   ..   race = col_character(),
##   ..   native = col_character(),
##   ..   feelpower = col_double(),
##   ..   plev = col_double(),
##   ..   vsex = col_double(),
##   ..   pitch_rmeanMD = col_double(),
##   ..   pitch_rvarMD = col_double(),
##   ..   intense_rmeanMD = col_double(),
##   ..   intense_rvarMD = col_double(),
##   ..   formant_rmeanMD = col_double(),
##   ..   formant_rvarMD = col_double(),
##   ..   pitch_smeanMD = col_double(),
##   ..   pitch_svarMD = col_double(),
##   ..   intense_smeanMD = col_double(),
##   ..   intense_svarMD = col_double(),
##   ..   formant_smeanMD = col_double(),
##   ..   formant_svarMD = col_double(),
##   ..   Zpitch_rmean = col_double(),
##   ..   Zpitch_rvar = col_double(),
##   ..   Zform_rmean = col_double(),
##   ..   Zform_rvar = col_double(),
##   ..   Zintense_rmean = col_double(),
##   ..   Zintense_rvar = col_double(),
##   ..   Zpitch_smean = col_double(),
##   ..   Zpitch_svar = col_double(),
##   ..   Zform_smean = col_double(),
##   ..   Zform_svar = col_double(),
##   ..   Zintense_smean = col_double(),
##   ..   Zintense_svar = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>
names(d)
##  [1] "voice"           "form_smean"      "form_svar"       "form_rmean"     
##  [5] "form_rvar"       "intense_smean"   "intense_svar"    "intense_rmean"  
##  [9] "intense_rvar"    "pitch_smean"     "pitch_svar"      "pitch_rmean"    
## [13] "pitch_rvar"      "pow"             "age"             "sex"            
## [17] "race"            "native"          "feelpower"       "plev"           
## [21] "vsex"            "pitch_rmeanMD"   "pitch_rvarMD"    "intense_rmeanMD"
## [25] "intense_rvarMD"  "formant_rmeanMD" "formant_rvarMD"  "pitch_smeanMD"  
## [29] "pitch_svarMD"    "intense_smeanMD" "intense_svarMD"  "formant_smeanMD"
## [33] "formant_svarMD"  "Zpitch_rmean"    "Zpitch_rvar"     "Zform_rmean"    
## [37] "Zform_rvar"      "Zintense_rmean"  "Zintense_rvar"   "Zpitch_smean"   
## [41] "Zpitch_svar"     "Zform_smean"     "Zform_svar"      "Zintense_smean" 
## [45] "Zintense_svar"
head(d)
## # A tibble: 6 × 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. 
## # ℹ 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>,
## #   intense_smeanMD <dbl>, intense_svarMD <dbl>, formant_smeanMD <dbl>, …

Step 3: Tidy data

d <- d %>%
  mutate(
    condition = ifelse(plev == 1, "high", "low"),
    speaker_sex = ifelse(vsex == 1, "female", "male")
  )

Step 4: Run analysis

Pre-processing

d <- d %>%
  mutate(
    condition = factor(plev, levels = c(-1, 1), labels = c("low", "high")),
    sex       = factor(vsex, levels = c(-1, 1), labels = c("male", "female"))
  )

pairs <- tibble::tribble(
  ~metric,            ~post,           ~base,
  "pitch_mean",       "pitch_smean",   "pitch_rmean",
  "pitch_var",        "pitch_svar",    "pitch_rvar",
  "loudness_mean",    "intense_smean", "intense_rmean",
  "loudness_var",     "intense_svar",  "intense_rvar",
  "resonance_mean",   "form_smean",    "form_rmean",
  "resonance_var",    "form_svar",     "form_rvar"
)

table(d$condition, useNA = "ifany")
## 
##  low high 
##   79   82
table(d$sex, useNA = "ifany")
## 
##   male female 
##     80     81

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:

get_adj_means <- function(post, base, metric, dat) {
  di <- dat %>%
    dplyr::select(dplyr::all_of(c(post, base, "condition", "sex"))) %>%
    dplyr::filter(!is.na(.data[[post]]), !is.na(.data[[base]]))

  # if nothing to fit, return NAs
  if (nrow(di) < 2) {
    return(tibble::tibble(
      metric    = metric,
      condition = factor(c("high","low"), levels = c("high","low")),
      adj_mean  = NA_real_
    ))
  }

  di$condition <- droplevels(di$condition)
  di$sex       <- droplevels(di$sex)

  include_cond <- nlevels(di$condition) >= 2
  include_sex  <- nlevels(di$sex)       >= 2

  rhs <- c(base)
  if (include_cond) rhs <- c("condition", rhs)
  if (include_sex)  rhs <- c(rhs, "sex")

  fit <- tryCatch(
    lm(reformulate(rhs, response = post), data = di),
    error = function(e) NULL
  )

  base_mean <- mean(di[[base]], na.rm = TRUE)

  # prediction grid
  newdat <- list()
  if (include_cond) newdat$condition <- levels(di$condition) else newdat$condition <- "high"
  if (include_sex)  newdat$sex       <- levels(di$sex)       else newdat$sex       <- "female"
  newdat <- expand.grid(newdat, KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE)
  newdat[[base]] <- base_mean

  preds <- if (!is.null(fit)) as.numeric(predict(fit, newdata = newdat)) else NA_real_

  out <- tibble::tibble(
    metric    = metric,
    condition = factor(newdat$condition, levels = c("high","low")),
    adj_mean  = preds
  )

  if (!include_cond) {
    out <- tibble::tibble(
      metric    = metric,
      condition = factor(c("high","low"), levels = c("high","low")),
      adj_mean  = mean(out$adj_mean, na.rm = TRUE)
    )
  }

  out
}

pairs <- tibble::tribble(
  ~metric,            ~post,           ~base,
  "pitch_mean",       "pitch_smean",   "pitch_rmean",
  "pitch_var",        "pitch_svar",    "pitch_rvar",
  "loudness_mean",    "intense_smean", "intense_rmean",
  "loudness_var",     "intense_svar",  "intense_rvar",
  "resonance_mean",   "form_smean",    "form_rmean",
  "resonance_var",    "form_svar",     "form_rvar"
)

adj_means <- purrr::map_dfr(
  seq_len(nrow(pairs)),
  ~ get_adj_means(
      post   = pairs$post[.x],
      base   = pairs$base[.x],
      metric = pairs$metric[.x],
      dat    = d
    )
)

adj_table <- adj_means %>%
  dplyr::group_by(metric, condition) %>%
  dplyr::summarise(adj_mean = mean(adj_mean, na.rm = TRUE), .groups = "drop") %>%
  tidyr::pivot_wider(names_from = condition, values_from = adj_mean) %>%
  dplyr::rename(
    `High-rank condition` = high,
    `Low-rank condition`  = low
  )
knitr::kable(
  adj_table %>% mutate(across(c(`High-rank condition`,`Low-rank condition`), ~round(.x, 2))),
  caption = "Results From Experiment 1: Adjusted Means for the Hierarchy-Based Acoustic Cues and Effect Sizes for Between-Condition Differences in These Means"
)
Results From Experiment 1: Adjusted Means for the Hierarchy-Based Acoustic Cues and Effect Sizes for Between-Condition Differences in These Means
metric High-rank condition Low-rank condition
loudness_mean 59.34 58.68
loudness_var 196.72 183.46
pitch_mean 158.57 155.49
pitch_var 1425.42 1648.93
resonance_mean 1129.41 1128.85
resonance_var 42176.00 43668.90

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

ancova_out <- purrr::pmap_dfr(
  list(pairs$metric, pairs$post, pairs$base),
  function(metric, post, base) {
    fit <- lm(reformulate(c("condition", "sex", base), response = post), data = d)
    a   <- anova(fit)
    tibble::tibble(
      metric   = metric,
      F        = unname(a["condition", "F value"]),
      p_value  = unname(a["condition", "Pr(>F)"]),
      df1      = unname(a["condition", "Df"]),
      df2      = fit$df.residual,
      N        = stats::nobs(fit)
    )
  }
)

knitr::kable(ancova_out, digits = 4, caption = "ANCOVA results")
ANCOVA results
metric F p_value df1 df2 N
pitch_mean 2.4306 0.1210 1 157 161
pitch_var 7.2255 0.0080 1 157 161
loudness_mean 2.4140 0.1223 1 157 161
loudness_var 7.7377 0.0061 1 157 161
resonance_mean 0.0823 0.7745 1 157 161
resonance_var 0.7172 0.3983 1 157 161

Step 5: Reflection

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

Same variables reached significance & conditions effects were in same direction as reported, but my F values for pitch and loudness variability were larger, and the pitch mean effect was slightly weaker than originally reported. List of discrepancies here: - Pitch mean: F(1, 157) = 2.43, p = .12 (smaller than reported F = 4.48, p < .05) - Pitch variability: F(1, 157) = 7.23, p = .008 (larger than reported F = 4.73, p < .05) - Loudness mean: F(1, 157) = 2.41, p = .12 (not significant, consistent with report) - Loudness variability: F(1, 157) = 7.74, p = .006 (larger than reported F = 4.66, p < .05) - Resonance mean: F(1, 157) = 0.08, p = .77 (similar to F < 1 reported) - Resonance variability: F(1, 157) = 0.72, p = .40 (similar to F < 1 reported)

How difficult was it to reproduce your results?

This one was more difficult than the one in Group A because I had to troubleshoot several data-wrangling issues before the models could run successfully.

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

Incomplete documentation and inconsistent naming conventions made this more difficult than necessary. For example, the paper didn’t clearly specify which baseline vars corresponded to each post-manipulation measure. Also, the preprocessing details such as how the hierarchy and sex variables were coded should have been included. It just led to unnecessary model and factor-level errors. The code book for all data sets document helped clarify the meaning of each variable, but a more organized format for this would also help.