This vignette uses the voicequality package to replicate two of the figures in this paper:

Setup

First, load the tidyverse packages and the data from the voicequality package.

library("tidyverse")

data("acoustics", package = "voicequality")
data("egg", package = "voicequality")

For the Gujarati part of the UCLA Voice Quality Project, Sameer Khan recorded ten speakers producing a list of words. There is one Gujarati-language audio file per speaker for each target word. Each file has several repetitions of a sentence that includes the word. One target vowel was extracted from each repetition, and acoustic and EGG measurements for all of the vowel tokens across all of the speakers and words were collected in the acoustics and egg tables.

Subset the Gujarati data from the acoustics table and join it with the egg table.

gujarati <- filter(acoustics, language == "Gujarati") %>%
  left_join(egg)

List all of the columns in the data. To read the documentation for these columns, use help("acoustics", package = "voicequality") and help("egg", package = "voicequality").

names(gujarati)
#>  [1] "speaker_id"             "speaker_number"        
#>  [3] "speaker_sex"            "language"              
#>  [5] "language_variety"       "filename"              
#>  [7] "segment_start"          "segment_end"           
#>  [9] "duration"               "timepoint"             
#> [11] "textgrid_label"         "segment_type"          
#> [13] "vowel_quality"          "phonation"             
#> [15] "vowel_nasality"         "tone_contour"          
#> [17] "number_of_tone_targets" "tone_category"         
#> [19] "language_specific_tone" "consonant_aspiration"  
#> [21] "consonant_nasality"     "strF0"                 
#> [23] "sF0"                    "pF0"                   
#> [25] "shrF0"                  "sF1"                   
#> [27] "sB1"                    "sF2"                   
#> [29] "sB2"                    "sF3"                   
#> [31] "sB3"                    "sF4"                   
#> [33] "sB4"                    "pF1"                   
#> [35] "pF2"                    "pF3"                   
#> [37] "pF4"                    "H1u"                   
#> [39] "H1c"                    "H2u"                   
#> [41] "H2c"                    "H4u"                   
#> [43] "H4c"                    "A1u"                   
#> [45] "A1c"                    "A2u"                   
#> [47] "A2c"                    "A3u"                   
#> [49] "A3c"                    "H1H2u"                 
#> [51] "H1H2c"                  "H2H4u"                 
#> [53] "H2H4c"                  "H1A1u"                 
#> [55] "H1A1c"                  "H1A2u"                 
#> [57] "H1A2c"                  "H1A3u"                 
#> [59] "H1A3c"                  "Energy"                
#> [61] "HNR05"                  "HNR15"                 
#> [63] "HNR25"                  "HNR35"                 
#> [65] "CPP"                    "SHR"                   
#> [67] "tone_misc"              "CQ"                    
#> [69] "CQ_H"                   "CQ_HT"                 
#> [71] "CQ_PM"                  "min_Vel"               
#> [73] "min_Vel_Time"           "peak_Vel"              
#> [75] "peak_Vel_Time"          "SQ2_SQ1"               
#> [77] "SQ4_SQ3"                "ratio"

Show a few of the columns.

select(gujarati, speaker_id, filename, segment_start, timepoint, CQ_H, H1H2c) %>%
  as_tibble()
#> # A tibble: 29,898 x 6
#>    speaker_id  filename                 segment_start timepoint  CQ_H H1H2c
#>    <chr>       <chr>                            <dbl>     <int> <dbl> <dbl>
#>  1 Gujarati_M1 2008-11-19-1300-01_Audio          441.         1 0.47  10.9 
#>  2 Gujarati_M1 2008-11-19-1300-01_Audio          441.         2 0.459  7.34
#>  3 Gujarati_M1 2008-11-19-1300-01_Audio          441.         3 0.47   5.12
#>  4 Gujarati_M1 2008-11-19-1300-01_Audio          441.         4 0.481  4.54
#>  5 Gujarati_M1 2008-11-19-1300-01_Audio          441.         5 0.48   4.08
#>  6 Gujarati_M1 2008-11-19-1300-01_Audio          441.         6 0.477  3.64
#>  7 Gujarati_M1 2008-11-19-1300-01_Audio          441.         7 0.473  3.37
#>  8 Gujarati_M1 2008-11-19-1300-01_Audio          441.         8 0.473  3.27
#>  9 Gujarati_M1 2008-11-19-1300-01_Audio          441.         9 0.466  3.38
#> 10 Gujarati_M1 2008-11-19-1300-01_Audio         2538.         1 0.503  9.76
#> # ... with 29,888 more rows

Each vowel token can be uniquely identified by its combination of filename and segment_start values. Each vowel token was divided into nine equal windows of time. This table has a row of measurements for each timepoint, so it includes nine rows per vowel token. Each row of measurements can be uniquely identified by its combination of filename, segment_start, and timepoint values.

Each Gujarati word has a file number assigned to it. The last two digits in each filename give the file number for the word that was recorded in that file.

Extract the file numbers from each filename into a separate column.

gujarati <- extract(gujarati,
  filename, "file_number", "^.*(\\d{2})_Audio$", remove = FALSE
)

Khan (2012) used a subset of the Gujarati data in the voicequality package. The appendix has a list of the words that were used for the paper. The file numbers for those words in the gujarati_word_list table are:

khan_2012_stimuli <- tibble(
  file_number = c(
    "01", "03", "06", "08", "09", "11", "12", "13", "30", "31", "32", "33",
    "34", "35", "38", "39", "40", "41", "46", "47", "48", "49", "56", "57",
    "59", "60"
  ),
  vowel_group = c(
    "a", "a", "a", "a", "a", "a", "a", "a", "e", "e", "e", "e",
    "e", "e", "e", "e", "e", "e", "e", "e", "e", "e", "o", "o",
    "o", "o"
  )
)

Join the complete Gujarati data with the list of file numbers to get the subset used in Khan (2012).

khan_2012 <- inner_join(gujarati, khan_2012_stimuli)

Reorder the levels of the phonation variable so that they will be in the correct order for the graph.

khan_2012 <- mutate(khan_2012,
    phonation = fct_relevel(phonation, c("Modal", "Breathy"))
)

Timecourse of contact quotient measure (Figure 10)

The left panel of Figure 10 in Khan (2012) compares the timecourses of the contact quotient for modal versus breathy vowels.

The voicequality data has measurements at nine timepoints within each vowel token. This figure plots the mean and standard error of per-speaker averages for the two phonation types at each timepoint. To recreate the figure, first calculate the average contact quotient for each phonation type, for each speaker. Speakers are uniquely identified by the speaker_id column.

mean_cq <- khan_2012 %>%
  group_by(speaker_id, timepoint, phonation) %>%
  summarize(CQ_H = mean(CQ_H, na.rm = TRUE)) %>%
  ungroup()

mean_cq
#> # A tibble: 180 x 4
#>    speaker_id  timepoint phonation  CQ_H
#>    <chr>           <int> <fct>     <dbl>
#>  1 Gujarati_F2         1 Modal     0.525
#>  2 Gujarati_F2         1 Breathy   0.428
#>  3 Gujarati_F2         2 Modal     0.532
#>  4 Gujarati_F2         2 Breathy   0.416
#>  5 Gujarati_F2         3 Modal     0.540
#>  6 Gujarati_F2         3 Breathy   0.400
#>  7 Gujarati_F2         4 Modal     0.547
#>  8 Gujarati_F2         4 Breathy   0.414
#>  9 Gujarati_F2         5 Modal     0.553
#> 10 Gujarati_F2         5 Breathy   0.420
#> # ... with 170 more rows

Next, plot the means and standard errors. ggplot2 can calculate these itself before plotting by using the summary stat with the appropriate functions.

ggplot(mean_cq, aes(timepoint, CQ_H, shape = phonation, linetype = phonation)) +
  geom_point(stat = "summary", fun.y = mean) +
  geom_line(stat = "summary", fun.y = mean) +
  geom_errorbar(stat = "summary", fun.data = mean_se, width = 0.1) +
  scale_linetype_manual(values = c("Modal" = "solid", "Breathy" = "longdash")) +
  scale_x_continuous(breaks = seq(1, 9)) +
  theme_minimal() +
  theme(legend.position = "bottom")

The results look very close to the left panel of the published Figure 10, except that all of the values are shifted up by about 0.05.

Average spectral tilt measures (Figure 4)

Figure 4 in Khan (2012) shows the average values of four spectral tilt measures for modal versus breathy vowels in each of the three vowel groups.

To recreate this figure, first calculate the average spectral tilt for each phonetic category for each speaker. There are six phonetic categories, with one category for each combination of phonation (Modal or Breathy) and vowel_group (a, e, or o). Speakers are uniquely identified by the speaker_id column. Note that this calculation averages over all nine timepoints in each vowel token.

mean_spectral_tilt <- khan_2012 %>%
  group_by(speaker_id, vowel_group, phonation) %>%
  summarize(
    H1H2c = mean(H1H2c, na.rm = TRUE),
    H2H4c = mean(H2H4c, na.rm = TRUE),
    H1A1c = mean(H1A1c, na.rm = TRUE),
    H1A3c = mean(H1A3c, na.rm = TRUE)
  )

mean_spectral_tilt
#> # A tibble: 60 x 7
#> # Groups:   speaker_id, vowel_group [30]
#>    speaker_id  vowel_group phonation H1H2c   H2H4c H1A1c H1A3c
#>    <chr>       <chr>       <fct>     <dbl>   <dbl> <dbl> <dbl>
#>  1 Gujarati_F2 a           Modal      2.29   6.84   22.8 15.5 
#>  2 Gujarati_F2 a           Breathy    6.77  10.6    31.9 22.3 
#>  3 Gujarati_F2 e           Modal      9.61  -5.13   11.7 13.9 
#>  4 Gujarati_F2 e           Breathy    9.13  -2.38   17.3 17.5 
#>  5 Gujarati_F2 o           Modal      8.34  -2.91   12.9  7.58
#>  6 Gujarati_F2 o           Breathy   11.0    1.26   21.4 17.0 
#>  7 Gujarati_F4 a           Modal      6.82   9.93   28.7 24.6 
#>  8 Gujarati_F4 a           Breathy    7.90  11.2    30.4 25.7 
#>  9 Gujarati_F4 e           Modal      9.86  -2.69   13.2 20.2 
#> 10 Gujarati_F4 e           Breathy   10.5    0.475  21.3 23.0 
#> # ... with 50 more rows

One way to make a single four-panel plot is to gather() the variables to be plotted, then plot each one in a separate facet.

mean_spectral_tilt <- mean_spectral_tilt %>%
  gather(measure, measurement, H1H2c, H2H4c, H1A1c, H1A3c) %>%
  mutate(measure = fct_relevel(measure, c("H1H2c", "H2H4c", "H1A1c", "H1A3c")))

ggplot(
  data = mean_spectral_tilt,
  mapping = aes(
    x = vowel_group,
    y = measurement,
    fill = phonation
    )
  ) +
  geom_bar(
    stat = "summary",
    fun.y = mean,
    position = position_dodge(0.9)
  ) +
  geom_errorbar(
    stat = "summary",
    fun.data = mean_se,
    position = position_dodge(0.9),
    width = 0.1
  ) +
  facet_wrap(~measure, scales = "free") +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    legend.position = "bottom"
  )
#> Warning: Removed 4 rows containing non-finite values (stat_summary).

#> Warning: Removed 4 rows containing non-finite values (stat_summary).

The resulting values look similar to the published Figure 4.