11 Interpreting Model Uncertainty
So far we have learned how to write down our ideas of the data generation process as a statistical model, automatically fit the model, assess fit, and interpret the \(\beta\) coefficients for what they mean about the role of an explanatory variable. This is useful information both for explanation research goals (e.g. “how big is the effect of \(X\) on \(Y\)?”) and for prediction goals (e.g. “what is the most likely value of \(Y\) given what we know about \(X\)?”).
Stopping our model interpretation here, however, can leave us with an incomplete understanding. Seeing a big number for \(\beta\) might make us think that \(X\) is the most important variable for explaining \(Y\), when in fact this model of the data generation process is critically incomplete and there is still uncertainty in \(Y\) that could be explained better with other information. Similarly, making a prediction for some value of \(Y\) given \(X\) doesn’t tell us how good that guess is - how far off will we be when we make such a prediction?
In order to have a more complete understanding of our statistical model’s performance, we also need to interpret the other parameter that the model fitting procedure gives us: \(\sigma\), in the observation model \(Y \sim N(\mu, \sigma)\). In this chapter we will cover this interpretation step.
11.1 Prediction accuracy
Sometimes in a prediction goal, we want to predict the summary statistic of a large group of data. This is useful to do when our research is about group-level attributes, like overall student graduation rate under different financial aid strategies or the average number of customer service calls to expect at different times of day.
However, other times your goal is to make predictions about individual data points. In this case the estimand is not a summary statistic of a group, such as \(\mu_{Y|X}\), but the value of the target variable for one specific observation, \(y_i\) when \(X=x_i\).
We face some difficulty in doing this when we’re dealing with a probabilistic data generation process. Recall that in this sort of process, we can’t perfectly predict specific data values, only their probability of happening via the parameter of the probability distribution output by the prediction equation. For instance, here is a plot of the distribution of highest_year_of_school_completed values in the GSS dataset for people who have 4 siblings specifically:
Despite having the same value on the explanatory variable, these people still vary in their education attainment. So what do we do if we need to guess one specific value to accomplish our research goal?
If we don’t know someone’s \(y_i\) value, but we do know their \(x_i\) value and we have a best-fitting statistical model that tells us \(P(Y|X)\), we could go with whatever value of \(Y\) is most likely to happen when \(X=x_i\). In the model we’ve been working with the last couple chapters, \(Y \sim \mathcal{N}(\mu, \sigma)\) and \(\mu = \beta_0 + \beta_1 X\), we are assuming the output of the data generation process is a normal probability distribution influenced by \(X\). If this is correct, then the most likely value to happen is the mean of that distribution, \(\mu\). So we can use our model to predict the most likely value of \(Y\) given \(X\) by calculating \(\mu = \beta_0 + \beta_1 X\). We will call our prediction about someone’s \(Y\) value \(\hat{Y}\), or “Y hat,” to distinguish it from the true value of \(Y\) for that person, \(y_i\). Essentially, by making a prediction \(\hat{Y} = \beta_0 + \beta_1 X\), we are treating the data generation process as deterministic.
Here is the predicted \(\hat{Y}\) value for people with \(X=4\), drawn as a red line on the previous histogram:
Even though 13.62 is the most likely value in this conditional probability distribution, if we had guessed this value as the years of education for everyone with 4 siblings, clearly we’d be wrong most of the time1. Written mathematically, \(Y - \hat{Y} \neq 0\). We’d be less wrong on average than if we had guessed any other number, but still wrong because there is a difference between our prediction \(\hat{Y}\) and the true value \(y_i\) for each person. This is known as the prediction error or residual for that person under our current statistical model. It is often denoted as \(\epsilon\), such that:
residual = the difference between what a statistical model predicts for an individual observation and the observation’s true value.
\[\epsilon = Y - \hat{Y} \tag{11.1}\]
How big is our prediction error for each person? We can again plot this, but this time plot the difference between our guess and the true value.
The shape of the histogram didn’t change, but the values on the x-axis did. Now, we are looking at the difference between our guess \(\hat{Y}\) (for which we are using the mean of the conditional probability distribution, \(\mu_{Y|X=4}\)) and each real data point (\(y_i\)). Think of these residual values as how much our guess missed by.
Now let’s plot the residuals for every data point in this data set, not just the people with 4 siblings. We can find these quickly with the function resid(), where the argument is the model object.
TipExercise
Use resid() to calculate the residuals for every data point in the GSS dataset and store them in a new column called residuals.
We can think of each residual as a prediction error. We used a different value for \(\hat{Y}\) for each person depending on what their \(x_i\) value is and what the prediction equation \(\beta_0 + \beta_1 X\) calculated, but we can still compute \(Y - \hat{Y}\) for everyone and thus get a residual value. The distribution of all of them together is called the error distribution.
error distribution = the distribution of residuals from a statistical model.
11.2 Quantifying prediction accuracy
If prediction for any individual observation is our goal, we would like to quantify how good our model is for making predictions regardless of what someone’s value on \(X\) may be. In other words, it would be useful to summarize this error distribution. We can lean on the distribution summary metrics we already know to do this.
Recall when we first learned about spread metrics in Chapter 5. We first saw that trying to quantify spread as simply the sum of all deviations from a central value didn’t work, because some values were negative and positive and they canceled each other out. Because model residuals can also be positive or negative (we can guess too low or too high), this same reality applies:
Instead, if we make sure all deviation values are positive by squaring them, we can calculate the sum of squares. Specifically, when talking about residuals, we can sum up the squared residuals to do the same thing since a residual is a deviation score between a predicted value and a true value. This metric, when talking about model residuals, is called the residual sum of squares (RSS) or sum of squared errors (SSE).
sum of squared errors = sum of squares calculated for an error distribution.
\[SSE = \sum_{i=1}^{n} (Y_i - \hat{Y}_i)^2 \tag{11.2}\]
TipExercise
Calculate the sum of squared errors of the residuals.
Because this is just the sum of squares equation applied to a vector of residuals, the same interpretation applies to this value as we would make with any other sum of squares. For the same set of data, a larger value means more spread in the residuals, which means our model is worse at predicting individual data points. Here’s the sum of squared error calculated for the model predicting education using respondent sex. Which model do you think would make more accurate predictions - the one above using number of siblings, or using sex?
However, a larger \(SSE\) value can also occur just because there are more data points. So we can also calculate the typical amount of squared error per observation, called the mean squared error (MSE). This is similar to calculating the variance of the error distribution, except instead of dividing but \(n-1\), we divide by \(n-k\) where \(k\) is the number of coefficients in the prediction equation of the model2.
mean squared error = variance of the error distribution.
\[MSE = \frac{\sum_{i=1}^{n} (Y_i - \hat{Y}_i)^2}{n-k} \tag{11.3}\]
Finally, because MSE is in units of squared residuals, it is tough to interpret by how much our guesses are typically wrong. So we can take the square root of MSE to get back to the original units of the target variable. This is called the root mean squared error (RMSE), and is like calculating the standard deviation of the error distribution.
root mean squared error = standard deviation of the error distribution.
When we first learned about standard deviation, we learned that the straightforward interpretation of it is the typical amount of deviation from the mean within a data distribution - if someone guessed the mean value for every data point, on average they would be off by plus or minus the standard deviation amount. We can interpret the RMSE of the error distribution the same way - it is the typical amount of prediction error our statistical model makes when predicting individual data points.
Now let’s look back at what our model fitting procedure identified as the \(\sigma\) parameter for the observation model:
Look familiar? This is the same value as the RMSE we just calculated. So we can interpret \(\sigma\) in the observation model as the typical amount of prediction error our statistical model makes when predicting individual data points. The larger the \(\sigma\) value in our model, the less accurate our predictions will be.
11.3 Prediction improvement
The above discussion hopefully makes it clear how we can interpret the fitted estimates from a statistical model when we have prediction goals. The prediction equation coefficients tell you how to generate predictions for new data, while the \(\sigma\) coefficient tells you how much those predictions will be wrong by, on average.
Understanding the \(\sigma\) estimate also helps us with explanation goals too, in that it helps us understand how much uncertainty is in our model of the data generation process. When we were first learning about probability, we saw that if two variables are associated, knowing the value of one variable changes our guesses about the value of another. In Chapter 10 we saw that this change in prediction is a change in the mean of the target’s conditional probability distribution, such that a 1-unit increase in the explanatory variable shifts our expectation about the mean target value in some direction.
But if \(X\) is an explanatory variable of \(Y\), it should also change our predictions about \(Y\) by narrowing the spread of the conditional probability distribution as well. In other words, before we know anything about \(X\), the marginal probability distribution for the \(Y\) values we could draw into our sample is wide. But once we know \(X\), if it contains information about what \(Y\) will be, then we should be able to reduce the options we’d expect for \(Y\) and be more confident that a narrower range of \(Y\) options will happen instead. It’s like if we are in a windowless room, we are pretty uncertain about what the weather is outside. But if we know it is August 15th in California, we can be pretty confident that it’s not snowing outside! Once we knew our date and location, the list of likely options about the weather got narrower, with our guesses more confident and accurate.
Thus the more accurately we can predict individual data points, the more we understand about what influences the target variable, and the less uncertain we will be about what someone’s \(Y\) value will be given their known \(X\) value.
Uncertainty relative to what?
- prediction accuracy as variance explanation / variance reduction
- null model
- variance explained: SSmodel = SStotal - SSerror
11.4 Standardized prediction accuracy
- R2
11.5 Relationship with correlation
11.6 Sources of unexplained uncertainty
Another reason to look at data predictions is to understand how much uncertainty is in the model. Even if the simulated data is a good fit for real data, if the prediction intervals are very wide, that means our explanatory variables don’t contain that much information about the target variable. Taking them into account doesn’t help our predictions of exact data values very much. Alternatively, if the prediction intervals are narrow, that tells us that knowing the explanatory variables changes our expectations for the target variable by a lot. We are confident about what that value will be if we know the explanatory variable values. In other words, we are much closer to understanding all the influences on the target variable and having a deterministic idea of the data generation process.
11.7 Estimands in statistical models
An estimand says precisely what quantity about the population/data-generating process you want to know.
Then the model gives you a way of thinking about the estimand
11.8 connection to correlation
sd when simulating Y|X with given correlation r
\[\sigma_{Y|X}^2 = \sigma_Y^2 - \frac{\text{Cov}(X,Y)^2}{\sigma_X^2}\]
\[Var(Y|X) = Var(Y) - \frac{Cov(X,Y)^2}{Var(X)}\]
\[Var(Y|X) = Var(Y)*(1 - r^2)\]
- Decompose the total variation in Y
Write
\[Y=\hat{Y}+\epsilon\]
where
\[\hat{Y} = \beta_0 + \beta_1X\]
Then
\[Var(Y)=Var(\hat{Y}+\epsilon)\]
Using the variance-of-a-sum formula:
\[Var(Y)=Var(\hat{Y})+Var(\epsilon)+2Cov(\hat{Y},\epsilon)\]
In ordinary least squares, the predicted values and residuals are uncorrelated, so
\[Cov(\hat{Y},\epsilon)=0.\]
Thus:
\[Var(Y)=Var(\hat{Y}) + Var(\epsilon)\]
This is the population version of the familiar regression decomposition:
total variation=explained variation+unexplained variation
- Express the variance of the predictions in terms of r
Because \(\hat{Y} = \beta_0 + \beta_1X\), the intercept doesn’t affect the variance:
\[Var(\hat{Y})=\beta_1^2 Var(X)\]
For simple linear regression,
\[\beta_1 = \frac{Cov(X,Y)}{Var(X)}\]
Therefore,
\[Var(\hat{Y})=\frac{Cov(X,Y)^2}{Var(X)}\]
Now use
\[r= \frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}}\] Squaring:
\[r^2 = \frac{Cov(X,Y)^2}{Var(X)Var(Y)}\]
Therefore,
\[Var(\hat{Y})=r^2 Var(Y)\]
- Substitute into the variance decomposition
We had
\[Var(Y)=Var(\hat{Y})+Var(\epsilon)\]
Substitute
\[Var(\hat{Y})=r^2 Var(Y):\] \[Var(Y)=r^2 Var(Y)+Var(\epsilon)\] Therefore,
\[Var(\epsilon)=Var(Y)−r^2 Var(Y)\] and hence
\[Var(\epsilon)=Var(Y)(1−r^2)\]
Under the homoscedastic regression model,
\[Var(Y∣X)=Var(\epsilon)\],
so:
\[Var(Y∣X)=Var(Y)(1−r^2)\]
This gives a very nice interpretation of r^2:
\[\frac{Var(Y∣X)}{Var(Y)}=1−r^2\]
In words, the proportion of variability in Y that remains within the conditional distribution, after accounting for X, is 1−r^2.
Equivalently,
\[r^2=1-\frac{Var(Y∣X)}{Var(Y)}\] So r^2 isn’t just an abstract “proportion of variance explained.” It tells us how much the predictor reduces the variance of the target distribution relative to its unconditional variance.