Fitting Linear Models

Agenda

Today, we’re going to start with trying to understand and fit some linear regressions.

Setup

~2 Minute Setup

Create your R Notebook for today and double check that your workspace is clear from last time.

Load the important packages for today:

library("ggplot2")
library("scales")
library("lsa2017")
library("tidyverse")
library("broom")

iy <- iy_ah %>%
        filter(plt_vclass == "iy")

ah <- iy_ah %>%
        filter(plt_vclass == "ah")

Understanding Linear Regression

I find the best way to understand linear regressions is as recipes for drawing a graph. Let’s start off with an example relevant to my own life: the relationship between temperatures in Celcius and temperatures in Farentheit. There are two important numeric constants to convert back and forth between them:

1. Freezing = 0°C = 32°F
2. A change of 1°C = 1.8°F

This gives us all of the necessary components for drawing a graph of the relationship between Farentheit and Celcius. First, we know that a line plotting this relationship must cross 0 on the x-axis at 32 on the y-axis, because freezing = 0°C = 32°F.

Next, we know that for 1 every degree we move to the right, we must move up 1.8.

The line, representing this relationship, thus has the formula:

\[ y = 32 + (1.8\times x) \]

And looks like

This line has the following properties:

1. Its Intercept, or where the line crosses 0 on the x-axis, is 32.
2. Its Slope, is 1.8 (that is, a change of 1 on the x-axis corresponds to a change in 1.8 on the y-axis)

And, in fact, all lines, including the ones we fit using linear regression, are described in terms of their intercept (where the line crosses 0 on the x-axis) and slope.

But, let’s say you didn’t know the formula to convert back and forth between Celcius and Farentheit. Maybe you weren’t even sure that there was such a formula. And all you had were a bunch of kind of crappy thermometers that gave you both Celcius and Farentheit readings. And maybe the thermometers are kind of hard to read, or the lighting conditions are dark, or you only have limited experience in reading thermometers. In that case, the relationship between Celcius and Farentheit might look a lot more like this:

In this case, you might want to fit a linear regression. A linear regression assumes that the relationship between two variables can be captured using some line, and it tries to estimate the Intercept and Slope of this assumed line.

The way it does this is by trying to find parameters for a line that minimize the “error” or the “residuals” between the observed data, and the line. For example, here are examples of bad lines. The first assumes that all temperatures Celcius are basically 50°F. The other assumes that all temperatures in Farentheit are basically 3\(\times\)C.

Here is a plot of the best-fitting line, that has the smallest average residuals.

There are actually many different ways to figure out what the best line is (remembering that by “line” we really mean the pairing of Intercept and Slope). For example, one really simplistic approach¹ would be just to draw a bunch of different lines, and see which one worked best.

However, there have been some long established methods for solving this problem, at is, for finding the one and truly best Intercept and Slope for a given set of data points.

Idiom

This isn’t R idiom, but stats idiom. Here are some translations for how people describe solving this problem.

• analytic, closed form = “We solved it with complicated math.”
• iterative, stochastic = “We made a bunch of guesses.”

For the sake of finishing the illustration, here is the solution that a linear regression found for this sample data:

Call:
lm(formula = f_observe ~ c_observe, data = temp_df)

Coefficients:
(Intercept)    c_observe  
     41.414        1.344  

The model has estimated that the Intercept is \(\approx\) 41 and the slope is \(\approx\) 1.3 . That is to say, based on our “sample”, it thinks 0°C = 41°F, and that every 1°C increase corresponds to a 1.3°F increase.

We are lucky to know that this is, strictly speaking, wrong, but we won’t usually, or typically, be able to in such a situation. However, the estimate of both the intercept and slope have the right sign (that is, it hasn’t guessed the intercept is -10, or the slope goes the other way), and they’re just about the right magnitude (that is, it hasn’t guessed that the slope is something like 10, or 100).

Now, let’s consider how the linear model is describing the value of this highlighted dot.

First, you take it’s value in Celcius, and multiply it by 1.3 (which is 22.55)

Then,you add 40 to it (which is 63.96)

Then, you add the residual error to give us the observed data point.

The generic formula for finding the value for any data point is

\[ F_i = 41.4 + (1.3\times C_i) + \epsilon_i. \] In this formula the \(i\) is convention for referring to some generic, \(i^{th}\) data point and \(\epsilon\) is conventional for referring to the residual errors.

Some additional conventions:

• Referring to the outcome variable as \(y\).
• Referring to the predictors as \(x\).
• Referring to the remaining parameters as \(\beta\).

So you may sometimes see people referring to \(\beta_0\), or the “betas” in a model. The fully generic formula for linear regression is sometimes given as

\[ y = \beta_0 + \beta x + \epsilon \]

Idiom

• We call the Intercept and Slopes the parameters or coefficients of the model.
• For every data point, the value you get by only using the slopes and intercepts are called the fitted values, and they all lie along the lines you draw.
• Almost every data point’s actual value is different from it’s fitted value, and these differences are called the residuals.

Thinking about a realer case

Here is some actual data you might work with. Below is a graph of the relationship between vowel duration and normalized F1 for the LOT vowel (i.e. "ah" in the iy_ah data frame). Each point is one vowel token. I predict that as the duration gets longer, the vowels should also get lower (i.e. have a higher F1) for a bunch of different reasons.

~5 Minute Discussion

This real data on duration and F1 is different from the example Celcius and Farentheit in a bunch of different (and important ways). In small groups, discuss what they might be.

Fitting a linear model

In R, the syntax for fitting a simple linear model is very straight forward. The function is called lm(), for linear model.

lm(outcome ~ predictor, data)

• outcome
  • This is the variable whose value that you’re trying to predict, e.g. vowel F1.
• predictor
  • This is the variable you’re usuing the predict the outcome value
• data
  • the data frame that contains the data we’re modelling with.

We can fit a linear regression for the ah data set like so:

ah_model <- lm(F1_n ~ dur, data = ah)

~2 minute activity

Fit the ah_model using the code above. Then look at the summary like so:

summary(ah_model)

Call:
lm(formula = F1_n ~ dur, data = ah)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8049 -0.3773  0.0133  0.3927  3.3113 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.60431    0.01083   55.82   <2e-16 ***
dur          3.50071    0.07421   47.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6302 on 21117 degrees of freedom
Multiple R-squared:  0.09533,   Adjusted R-squared:  0.09529 
F-statistic:  2225 on 1 and 21117 DF,  p-value: < 2.2e-16

Look over these numbers, and (without looking down too much further!), discuss in small groups what you understand here, and how to interpret the outcome.

Visualizing this model.

Before we even get to those “p-values” and what they even mean, let’s talk about the coefficients.

• (Intercept) 0.60431
  • This is the intercept of the line, meaning when dur = 0, F1_n is predicted to be about 0.6.
• dur 3.50071
  • This is the slope, meaning when dur increases by 1 second, we predict F1_n to increase by about 3.5

Here’s how that looks visualized on our data set.

Thinking about the model

Some things not so great about these coefficients:

1. They’re overly precise.
   • At the best of times, linguistics is a one or two decimal place field.
2. The Intercept makes no sense.
   • Not only are there no vowels with 0 duration in this data set, there can’t be a vowel with 0 duration.
3. The slope makes very little sense.
   • It’s possible for a vowel to be 1 second long, but it’s not likely. The longest vowels in this data set (0.8s) are probably forced alignment errors. The vast majority of vowels in this data are <0.25s.

The good news is that these problems are fixable by “rescaling” or “standardizing” our predictor.

Fixing the intercept.

The simplest thing to do to fix the intercept would be to try to figure out where the line crosses the x-axis at some other value, instead of 0. Why not try to estimate where it crosses x at the median duration (0.12) instead?

The bad news that we can’t easilly make a linear model try to estimate this value. All the math surrounding fitting linear models is based on trying to estimate where the line crosses the x-axis at 0.

The good news is that we can trick the model into estimating this number, by just subtracting out the median duration from duration!

ah <- ah %>%
        mutate(dur_center = dur - median(dur))

Doing this is not only “allowed”, it is encouraged. Simply adding or subtracting a constant value from a continuous predictor won’t affect it’s statistical significance. But it does make the intercept of the model more interpretable.

ah_model_c <- lm(F1_n ~ dur_center, data = ah)
summary(ah_model_c)

Call:
lm(formula = F1_n ~ dur_center, data = ah)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8049 -0.3773  0.0133  0.3927  3.3113 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.024392   0.004453  230.03   <2e-16 ***
dur_center  3.500713   0.074211   47.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6302 on 21117 degrees of freedom
Multiple R-squared:  0.09533,   Adjusted R-squared:  0.09529 
F-statistic:  2225 on 1 and 21117 DF,  p-value: < 2.2e-16

Fixing the slope

The next issue to deal with is the slope parameter. The interpretation of the slope parameter as it is is that for every 1 second vowel duration increases, normalize F1 increases 3.5. While this may be accurate, it’s not exactly useful, because the full range of the data, including the extreme outliers doesn’t span a whole 1 second.

summary(ah$dur)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0500  0.0900  0.1200  0.1337  0.1610  0.8700 

Again, what we would like the model to do is tell us not how much normalized F1 will change when the duration changes by 1 second, but maybe something more like 50ms. And again, we can’t directly make the model do this, but we can trick it into doing so, but dividing duration by some value.

What is standard in these circumstances is to divide the predictor variable by its standard deviation, a measure of how spread out the data is from its center point. Gelman & Hill (2006) recommend actually dividing by the standard deviation times 2.

ah <- ah %>%
          mutate(dur_z = dur_center/(sd(dur_center)*2))

ah_model_z <- lm(F1_n ~ dur_z, data = ah)
summary(ah_model_z)

Call:
lm(formula = F1_n ~ dur_z, data = ah)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8049 -0.3773  0.0133  0.3927  3.3113 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.024392   0.004453  230.03   <2e-16 ***
dur_z       0.409123   0.008673   47.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6302 on 21117 degrees of freedom
Multiple R-squared:  0.09533,   Adjusted R-squared:  0.09529 
F-statistic:  2225 on 1 and 21117 DF,  p-value: < 2.2e-16

Now, the slope is describing how much normalized F1 increases every time the vowel duration increases 2 standard deviations, \(\approx\) 100ms

Describing what we just did

Table XX displays the estimated parameters predicting normalized F1 of /ɑ/. Duration was centered at the median and scaled by two times its standard deviation (117) per the recommendation of Gelman & Hill (2006).

term	estimate	std.error	t_value	p.value
(Intercept)	1.02	0.004	230	0
Scaled Duration	0.41	0.009	47	0

Pause for Questions/Debugging

Evaluating our model

So, is this model any good? This question can be approached from a few different angles.

1. Is it possible that any of our model parameters are 0?
2. How great is our uncertainty about the specific model parameters? That is, the model has estimated an intercept of 1.02. How much wiggle room is there?
3. Do all of our model parameters have the right sign? That is, are they predicting an effect in the “right” direction?
4. Are our model parameters a reasonable size? Is it possible that any of them are being dramatically overestimated?

summary(ah_model_z)

Call:
lm(formula = F1_n ~ dur_z, data = ah)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8049 -0.3773  0.0133  0.3927  3.3113 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.024392   0.004453  230.03   <2e-16 ***
dur_z       0.409123   0.008673   47.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6302 on 21117 degrees of freedom
Multiple R-squared:  0.09533,   Adjusted R-squared:  0.09529 
F-statistic:  2225 on 1 and 21117 DF,  p-value: < 2.2e-16

Some of these questions fall under the rubric of “null hypothesis testing.” When people say their results are “significant”, or p<0.05 this is what they’re talking about. In its current use, people often try to divide their predictor variables into those which “have an effect” and those that “don’t have an effect” and they often use p-values to do this.

The summary of the linear model contains the p-values in the final column called Pr(>|t|). Here my best shot at a definition of p-values:

P Values

If we assume that the slope of the line is actually 0 there is no relationship present in the data, the probability we would observe a slope as extreme as, or more extreme than this, given the size of the sample and the variance of the data, is the p-value.

Typically, we take p < 0.05 as the cut off for saying it would be highly improbable to observe a slope of this size if the slope were actually 0 there was no relationship in the data, so the slope is probably not 0, we reject the null hypothesis that there is no relationship in the data.

A problem with p-values is that even if you’ve read that paragraph, you probably still don’t actually understand what p-values are. And when you do really understand what p-values are, you’ll probably forget what they really are. In fact, I tweeted out what I was planning to say about p-values, and based on the feedback I had to make some changes (as evidenced by the strike throughs).

My hypothesis is that truly understanding what p-values are requires a kind of counterfactual thinking and conditional probability that people are actually very bad at, even people who work with p-values all the time.

The fact of the matter is that you must include p-values in your writeups of your analyses because society demands it. However I think it is important to be judicious. If you have a p-value greater than 0.05, your evaluation of the effect should be downgraded. If you have a p-value less than 0.05, it’s also important not to stop there and conclude that the effect is “real”.

With that in mind it’s crucial to also ask whether the estimated effect is in the expected direction, and whether the effect size is sensible. Addressing these two questions requires having some background to draw one (previously reported effect sizes in the lit), a good theoretical reason for fitting this particular model, and perhaps some contextual info.

If you don’t have any reason to hypothesize an effect in any particular direction (could be positive, could be negative), then you might not be ready, as far as your inferences go, to use statistical modelling yet.

For our particular model, I had some pretty good reasons to assume a positive effect, which is what we got. Is 0.41 a “big” or a “small” effect? Compared to what?

Multiple Regression

We often want to figure out how multiple possible predictors influence the data we’re investigating. Such models are called “multiple regression”. For the ah regression, we’re next going to add the word context in which the vowel appears:

table(ah$context)

coextensive       final     initial    internal 
          2         239         304       20574 

We’re going to focus on just word final and word internal contexts.

ah_2 <- ah %>% filter(context %in% c("internal", "final"))

Adding context to the model is as easy as + context

ah_mod_mult <- lm(F1_n ~ dur_z + context, data = ah_2)
summary(ah_mod_mult)

Call:
lm(formula = F1_n ~ dur_z + context, data = ah_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8182 -0.3768  0.0128  0.3916  3.3103 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.805895   0.040797  19.754  < 2e-16 ***
dur_z           0.413973   0.008755  47.283  < 2e-16 ***
contextinternal 0.219005   0.040959   5.347 9.04e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6283 on 20810 degrees of freedom
Multiple R-squared:  0.09723,   Adjusted R-squared:  0.09714 
F-statistic:  1121 on 2 and 20810 DF,  p-value: < 2.2e-16

Here’s how to interpret these results. First, when you add a categorical predictor to a model, one of the categories is chosen as a reference level.² The (Intercept) is now the estimated F1_n when all of the predictors are 0 or the reference level. For this model, that means the predicted value when dur_z==0 and context=="final".

We can, for now, interpret the slope dur_z as being the slope for word final ah.

The effect of context=internal can be interpeted as how much offset the intercept of the line for the internal context is:

And the slope of the word internal context is the same as the word final context, as far as this model is concerned.

Evaluating this model.

~5 Minute Discussion

What do we think of this model? Does context have an effect?

summary(ah_mod_mult)

Call:
lm(formula = F1_n ~ dur_z + context, data = ah_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8182 -0.3768  0.0128  0.3916  3.3103 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.805895   0.040797  19.754  < 2e-16 ***
dur_z           0.413973   0.008755  47.283  < 2e-16 ***
contextinternal 0.219005   0.040959   5.347 9.04e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6283 on 20810 degrees of freedom
Multiple R-squared:  0.09723,   Adjusted R-squared:  0.09714 
F-statistic:  1121 on 2 and 20810 DF,  p-value: < 2.2e-16

Interactions

The fact that the two slopes of word-final and word-internal ah were the same above was an assumption of the model. To put it another way, it didn’t find that the two slopes were different because it didn’t look. But it is possible that the slopes of these two contexts are different, and that we should explore this possibility. This is done by entering context into the model like so:

ah_mod_interact <- lm(F1_n ~ dur_z * context, data = ah_2)
summary(ah_mod_interact)

Call:
lm(formula = F1_n ~ dur_z * context, data = ah_2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8278 -0.3762  0.0130  0.3917  3.3104 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            0.84728    0.04475  18.934  < 2e-16 ***
dur_z                  0.31189    0.04622   6.748 1.54e-11 ***
contextinternal        0.17721    0.04497   3.940 8.16e-05 ***
dur_z:contextinternal  0.10589    0.04707   2.249   0.0245 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6282 on 20809 degrees of freedom
Multiple R-squared:  0.09745,   Adjusted R-squared:  0.09732 
F-statistic: 748.9 on 3 and 20809 DF,  p-value: < 2.2e-16

Here’s how to interpret this model. First, the (Intercept) is the estimated F1_n for ah when dur_z==0 and context=="final".

The dur_z slope is the estimated slope for context=="final".

The effect of context=="final" is how offset the intercept for the final context is

And the interaction term is how different the slope of the final context is. That is, the dashed blue line shows the final context slope if it were the same as the internal context, the solid blue line is the estimated slope of the final context, and the difference is what is represented in the interaction term.

Reading an even Realer Example

Here is an example regression table from one of my own papers (Fruehwald 2016). I was looking at the F1 nucleus of /ay/ using the following predictors:

• Speaker’s year of birth
• Voicing of the following segment
• The frequency of the word

Discussion

How should we read these results?

Next Week - Mixed Effects Models and Model Selection

---

1. Actually, this “simplistic approach” is more like the sophisticated approach these days, because there are complicated maths behind making good guesses which are given semantically opaque names like “Simulated Annealing” and “Hamiltonian Monte Carlo”.↩

2. There are actually other approaches to encoding categorical predictors that don’t involve choosing a reference level that we’re not going to get to.↩