LING 3300: Linear Regression, continued
Thomas Kettig
October 29, 2026
Our goal in fitting a regression model will be to obtain three major pieces of information:
- the best-fit line (y = mx + b) – we’ve done this before.
- how much variation in vowel F2 can be accounted for by vowel F1 (R^2) – we’ve also talked about this.
- whether the effects of m (slope) and b (intercept) are significantly different from 0 – we haven’t done this yet.
First step: let’s load in the data and required packages. This is the McMurray & Jongman dataset of fricative and vowel measures that we used for the midterm.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.2
## ✔ purrr 1.2.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsmcm <- read_csv("/Users/Thomas/Library/CloudStorage/OneDrive-YorkUniversity/LING 3300/Datasets/McMurrayJongmanFricativeAcoustics.csv")## Rows: 2880 Columns: 29
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): Talker, Fricative, Vowel
## dbl (26): ID, Rep, F0, F1, F2, F3, F4, F5, maxpf, rms_f, rms_v, dur_f, dur_v...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Linear regression with one continuous predictor
In this example, we’ll predict a speaker’s mean F2 value from their mean F1. Because each person’s physiology is different, their F1 and F2 values won’t be comparable to each other. So rather than throwing all the raw data together, we have to use the speaker means due to the independence of observations assumption.
spkr_means <- mcm %>%
group_by(Talker) %>%
summarise(mean_f1 = mean(F1, na.rm = TRUE),
mean_f2 = mean(F2, na.rm = TRUE))Note that if we wanted to get the by-talker and by-vowel mean F1 and F2 (which is for sure more interesting than just knowing the mean F1/F2 from all the vowels uttered by a particular speaker in a sample) we would add ‘Vowel’ in the group_by() function. But for now let’s just lump all the vowels together and see whether, overall, there’s a relationship between the by-speaker mean F1 and F2.
ggplot(spkr_means, aes(x = mean_f1, y = mean_f2)) +
geom_point() +
geom_smooth(method = "lm") +
theme_classic() +
xlab("speaker mean F1 (Hz)") +
ylab("speaker mean F2 (Hz)")## `geom_smooth()` using formula = 'y ~ x'
In a scatterplot that shows this relationship along the best-fit line, we’ll put F2 on the y-axis because that is the variable we are predicting from F1. We can see from the line that we should expect a positive effect of F1 on F2 (positive slope).
We can also see that if we were to extend the the intercept all the way to the left until it reaches the y-axis (where x = 0), the y-intercept would be well above 0.
ggplot(spkr_means, aes(x = mean_f1, y = mean_f2)) +
geom_point() +
geom_smooth(method = "lm", fullrange = TRUE) +
scale_x_continuous(limits = c(0, NA), expand = c(0, 0)) +
theme_classic() +
xlab("speaker mean F1 (Hz)") +
ylab("speaker mean F2 (Hz)")## `geom_smooth()` using formula = 'y ~ x'
However, this isn’t all that meaningful because F1 and F2 values are always well above 0 anyway, and we’re never interested in knowing what happens when F1 equals 0, since that’s physically impossible.
To make the y-intercept a little more interpretable, people will often center a continuous predictor variable on 0. We can use the scale() function to do this. scale() takes three arguments: the vector to operate on, whether centering should be applied, and whether scaling should be applied. If scaling is set to TRUE, then it gives a z-score; otherwise, if it’s set to false, it just subtracts the mean from each point. For now, we’ll set scale to FALSE and just center our x variable, F1:
spkr_means$centered_f1 <- scale(spkr_means$mean_f1, center = TRUE, scale = FALSE)
ggplot(spkr_means, aes(x = centered_f1, y = mean_f2)) +
geom_point() +
geom_smooth(method = "lm") +
theme_classic() +
xlab("centered speaker mean F1 (Hz)") +
ylab("speaker mean F2 (Hz)")## `geom_smooth()` using formula = 'y ~ x'
From this visualization, we can now see that when the predictor (F1 here) is 0, we get the mean of the F2 dimension (around 1778 Hz). But otherwise, the points and the line look identical to what we plotted at first.
Notice what’s happening: instead of the y-intercept being all the way off the chart and sort of meaningless, now the very center of our plot is the y-intercept! So the intercept is just the mean of the dependent (predicted) variable when the predictor is set to 0.
# Don't believe me?
mean(spkr_means$mean_f2)## [1] 1778.078Now let’s actually get the information from the linear regression model:
fit_centered <- lm(mean_f2 ~ centered_f1, spkr_means)
summary(fit_centered)##
## Call:
## lm(formula = mean_f2 ~ centered_f1, data = spkr_means)
##
## Residuals:
## Min 1Q Median 3Q Max
## -254.61 -79.79 -28.12 94.60 235.74
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1778.0778 25.8224 68.858 <2e-16 ***
## centered_f1 1.2393 0.9602 1.291 0.213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 115.5 on 18 degrees of freedom
## Multiple R-squared: 0.08471, Adjusted R-squared: 0.03386
## F-statistic: 1.666 on 1 and 18 DF, p-value: 0.2131From this output, the key info is in the Coefficients table. You might also be interested in the Adjusted R-squared value. The adjusted R-squared value ranges from 0 to 1. In this case, it’s 0.03, which is ridiculously low: it means that F1 accounts for approximately 3% of the variation in F2. (My hunch is that vowel quality does a much better job predicting F2 than F1… the phoneticians should have the intuition that F1 varies independently of F2).
This leads us to the next critical component of the model: is F1 a significant predictor of F2? That is, are they significantly related? Note that the Coefficients table has two lines: intercept and centered_f1, in other words, the b and the m in our equation y = mx + b. The value of “m” or “b” is in the Estimate column, so here it would by y = 1.24x + 1778.08. In other words, to get a speaker’s F2 from their F1, the best linear prediction is to multiply the F1 value by 1.24 and add 1778.08 Hz to it. This is cool; unfortunately, the model fit as shown in the adjusted R-squared value is awful, so we won’t be all that accurate.
Next to the Estimate column in the table, we get our standard error, our t-value, and our p-value. The p-value for our predictor (F1) is p = 0.21. If our alpha is set at 0.05 as is standard, then p = 0.21 is not less than 0.05, so F1 is not significantly related to F2. In other words, the F1 estimate of 1.24 is not significantly different from 0, and we can’t trust that this is a meaningful relationship.
The intercept is significantly different from 0, with p < 0.05, but this is not that interesting. F2 will always be considerably higher than 0. Most studies will usually omit discussion of the intercept for this reason. That said, it’s cool to see that when you scale your predictors, the intercept is interpretable as the mean of your dependent variable.
Linear regression with one categorical predictor: Treatment/dummy coding
In this example, we’ll use a subset of the dataset with just /s/ and /ʃ/ fricative tokens. We’re interested in predicting M1 (spectral center of gravity, or COG) from the fricative identity (either /s/ or /ʃ/).
Let’s start out with our original raw dataset, and make a new subset with just the voiceless sibilants /s/ and /ʃ/.
voiceless_sibilants <- subset(mcm, Fricative %in% c("s", "sh"))Now we need to get independent means of the two sibilant COGs for each speaker – again, because each speaker has a different physiology and their raw values can’t be directly compared to each other. The column in the raw dataset that contains each token’s mean COG is M1.
vcl_sib <- voiceless_sibilants %>%
group_by(Talker, Fricative) %>%
summarise(meanCOG = mean(M1, na.rm = T))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Talker and Fricative.
## ℹ Output is grouped by Talker.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Talker, Fricative))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.Next, we’ll convert the categorical predictor to a factor and set up the contrast matrix for the linear regression using treatment (dummy) coding. We can then check to see we set it up properly by looking at the contrasts.
The key R code here when dealing with categorical variables is to convert the categorical variable to a factor using the factor() function.
vcl_sib$Fricative <- factor(vcl_sib$Fricative, levels = c("s", "sh"))
# contrasts(vcl_sib$Fricative) <- contr.treatment(2)
contrasts(vcl_sib$Fricative)## sh
## s 0
## sh 1If you want to use treatment coding like we do here, the contr.treatment() line that I’ve commented out above is optional because R will automatically assume treatment coding as its default.
Before we dive into the linear regression, we should visualize our data, though.
ggplot(vcl_sib) +
geom_boxplot(aes(x = Fricative, y = meanCOG, color = Fricative)) +
theme(legend.position = "none")
It looks like there’s definitely a major effect of fricative identity: /ʃ/ has a way lower COG than /s/. So already we expect this effect to come out as significant in our regression model.
We can also use pirateplot() to see this with a really simple line of code. Remember also that we can easily get it to show us the mean and 95% confidence interval by setting inf.method to ‘ci’. Eyeballing the CI will be even more useful than seeing the distribution of the data in a boxplot!
library(yarrr)## Loading required package: jpeg## Loading required package: BayesFactor## Loading required package: coda## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## ************
## Welcome to BayesFactor 0.9.12-4.7. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************## Loading required package: circlize## ========================================
## circlize version 0.4.17
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================##
## Attaching package: 'yarrr'## The following object is masked from 'package:ggplot2':
##
## diamondspirateplot(meanCOG ~ Fricative, vcl_sib,
inf.method = 'ci')
Indeed, we see that the mean of /s/ is well outside the CI of /ʃ/ and vice-versa. So when we run our regression, we will almost certainly find that this difference is significant.
fit_dummy <- lm(meanCOG ~ Fricative, vcl_sib)
summary(fit_dummy)##
## Call:
## lm(formula = meanCOG ~ Fricative, data = vcl_sib)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2772.59 -391.80 -3.11 361.55 1517.24
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6860.8 174.4 39.328 < 2e-16 ***
## Fricativesh -2091.8 246.7 -8.479 2.71e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 780.2 on 38 degrees of freedom
## Multiple R-squared: 0.6542, Adjusted R-squared: 0.6451
## F-statistic: 71.89 on 1 and 38 DF, p-value: 2.711e-10First, let’s check out the intercept - it’s the mean of means of /s/ COGs. Because s comes before sh alphabetically, R has by default set s as the control condition – that is, /s/ is considered to be the x = 0 condition, so its mean will by definition be the y-intercept.
just_s <- subset(vcl_sib, Fricative == "s")
mean(just_s$meanCOG)## [1] 6860.75We also see that the slope of -2091.8 is significantly different from zero. We know this because the p value is really low. 2.71e-10, or 0.000000000271, is far less than 0.05. To save space we’d usually describe this as p < 0.001. This means that we’ve found a significant effect of /s/ vs. /ʃ/ fricative identity on the by-speaker mean COG.
In addition, we can look at the Adjusted R-squared to figure out just how much knowing /s/ vs. /ʃ/ identity helps us in predicting the fricative COG. It turns out that 64.5% of the variation in COG between /s/ and /ʃ/ can be explained simply from knowing which of them we’re talking about. That’s quite a bit!
Practice: Linear regression with categorical variables
Calculate the by-speaker means of M1 and MaxPF. Using these means, plot the relationship between these two measures of fricative frequency (both in Hz). Are they significantly correlated? How strongly?
Run a simple linear regression predicting a talker’s MaxPF from their M1. Is MaxPF significantly predicted by M1? (Note that we did something very similar for the midterm – can you identify the slight difference with the way we’re conceiving of the relationship here?)
What is the multiple R-squared (unadjusted) from the linear regression model in Q2? Take its square root. Compare it to the Pearson’s R found for Q1. Feels good, huh?
Now calculate the by-speaker, by-fricative means of M1 and MaxPF. Using these means, plot the relationship between these two measures of fricative frequency, with the colour of the best-fit regression line differing by fricative. Do some of the fricatives appear to have a stronger correlation between the two measures than others?
Run a simple linear regression predicting a talker’s MaxPF from their M1, but separately for each fricative. For which fricative does a 1 Hz increase in M1 result in the largest predicted increase in MaxPF (that is, which slope is the most highly positive)? For which fricative does M1 best predict MaxPF?
Practice: Linear regression with a categorical predictor
- Try to answer the questions: Do sibilants have an overall higher intensity than non-sibilants? and Do voiced fricatives have an overall higher intensity than voiceless fricatives? using linear regressions, with the predictors treatment coded. You’ll need to construct two models, one regarding sibilants vs. non-sibilants and one for the voiced vs. voiceless question.