Statistical Rethinking. Chapter 5
5 Multivariate Linear Models
Correlations are not rare in nature but they don’t indicate causal relationships. To distinguish mere association from evidence of causation - much of statistical effort is devoted to multivariate regression. Reasons often given for multivariate models include:
- Statistical “control” for confounds. A confound is a variable that may be correlated with another variable of interest.
- Multiple causation. Even when confounds are absent a phenomenon may arise from multiple causes. When causation is multiple one cause can hide another. Multivariate models can help in such setting.
- Interactions.Even when variables are completely uncorrelated, the importance of each may still depend upone the other.
We’ll focus on two thigs multivariate models can help with:
- revealing spurious correlations
- revealing important correlationg that can be masked by unrevealed correlations with other variables.
5.1. Spurious association
Let’s explore correlation between divorce rate and marriage rate. Another predictor associated with divorce is the median age at marriage. It is also can be a good predictor of divorce rate - higher age at marriage predicts less divorce.
So let’s create a linear regression model:
\[D_i {\sim} Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta_A A_i\\ \alpha {\sim} Normal(10, 10)\\ \beta {\sim} Normal(0, 1)\\ \sigma {\sim} Uniform(0, 10)\]
\(D_i\) is the divorce rate for State i, and \(A_i\) is State’s median age at marriage. We’re going to standardize the predictor here, because it’s a good habit to get into.
library(rethinking)
data(WaffleDivorce)
d <- WaffleDivorce
#standardize predictor
d$MedianAgeMarriage.s <- (d$MedianAgeMarriage - mean(d$MedianAgeMarriage)) / sd(d$MedianAgeMarriage)
#fit model
m5.1 <- map(
alist(
Divorce ~ dnorm(mu, sigma),
mu <- a + bA * MedianAgeMarriage.s,
a ~ dnorm(10, 10),
bA ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)And the following code will compute the shaded confidence region.
precis(m5.1) Mean StdDev 5.5% 94.5%
a 9.69 0.20 9.36 10.01
bA -1.04 0.20 -1.37 -0.72
sigma 1.45 0.14 1.22 1.68#compute percentile interval of mean
MAM.seq <- seq(from = -3, to = 3.5, length.out = 30)
mu <- link(m5.1, data = data.frame(MedianAgeMarriage.s = MAM.seq))[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]mu.PI <- apply(mu, 2, PI)
#plot it all
plot(Divorce ~ MedianAgeMarriage.s, data = d, col = rangi2)
abline(m5.1)only using the first two of 3 regression coefficientsshade(mu.PI, MAM.seq)
As we can see from precis() output each additional standard deviation (sigma?) of delay in marriage (1.04 years) predicts a decrease of about one divorce per thousands abults, with a 89% interval from about -1.4 to -0.7. So it’s reliably negative, even though the magnitude of the difference may vary quite a lot - the upper bound is half the lower bound. The magnitute of difference means that even though the association here is implausibly positive, it could be both much stronger of weaker than the mean.
We can fit a similar regression for another relationship - marriage rate:
#standardize predictor
d$Marriage.s <- (d$Marriage - mean(d$Marriage)) / sd(d$Marriage)
#fit model
m5.2 <- map(
alist(
Divorce ~ dnorm(mu, sigma),
mu <- a + bR * Marriage,
a ~ dnorm(10, 10),
bR ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.2) Mean StdDev 5.5% 94.5%
a 6.16 1.27 4.13 8.20
bR 0.18 0.06 0.08 0.27
sigma 1.67 0.17 1.40 1.94And this shows an increase of 0.18 divorces for every additional standard deviation of marriage rate. So this relationship is not that strong as the previous one.
But comparing just parameter means between different bivariate regressions is not way to decide which predictor is better, so we’ll build a multivariate model with the goal of measuring the partial value of each preditor. The question we want answered is:
What is the predictive value of a variable, one we already know all of the other predictor variables?
In our case this question can be divided into two separate questions: * After we already know marriage rate, what additional value is here in also knowing age at marriage? * After we already know age at marriage, what additional value is there in also knowing marriate rate?
5.1.1. Multivariate notation
The strategy:
- Nominate the predictor variables you want in the linear model of the mean
- For each predictor, make a parameter that will measure its association with the outcome
- Multiply the parameter by the variable and add the term to the linear model
\[D_i {\sim} Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta_RR_i + \beta_aA_i\\ \alpha {\sim} Normal(10, 10)\\ \beta_R {\sim} Normal(0, 1)\\ \beta_A {\sim} Normal(0, 1)\\ \sigma {\sim} Uniform(0, 10)\]
R was chosen here for marriage rate and A for age at marriage.
So what does it mean to assume \(\\mu_i = \alpha + \beta_RR_i + \beta_aA_i\)? It means that the expected outcome for any State with marriage rate \(R_i\) and median age at marriage \(A_i\) is the sume of three independent terms. The first term is a constant, \(\alpha\). The second is the product of marriage rate, \(R_i\), and the coefficient \(\beta_R\), that measures the association between marriage rate and divorce rate. The third is similar, but for association between medin age at marriage instead.
5.1.2. Fitting the model
To fit this model to the divorce data, we just expand our linear model.
m5.3 <- map(
alist(
Divorce ~ dnorm(mu, sigma),
mu <- a + bR*Marriage.s + bA*MedianAgeMarriage.s,
a ~ dnorm(10, 10),
bR ~ dnorm(0, 1),
bA ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.3) Mean StdDev 5.5% 94.5%
a 9.69 0.20 9.36 10.01
bR -0.13 0.28 -0.58 0.31
bA -1.13 0.28 -1.58 -0.69
sigma 1.44 0.14 1.21 1.67The posterior mean for marriage rate, \(B_R\), is now close to zero, with plenty of probability of both sizes of zero. The posterior mean for age at marriage, \(B_A\), has gotten actually slightly farther from zero, but is essentially unchanged.
It will help to visualize these posterir distribution estimates:
rethinking::precis_plot(precis(m5.3))
This is the result, with MAP values shown by the points and the percentile intervals by the solid horizontal lines.
You can interpret these estimates as saying:
Once we know median age at marriage for a State, there is little or no additional predictive power in also knowing the rate of marriage in that State
5.1.3 Plotting multivariate posteriors
There are three useful types of interpretive plots:
- Predictor residuals plots - show the outcome against residual predictor values
- Counterfactual plots - show the impled preictions for imaginary experiments in which the different predictor variables can be changed independently of one another.
- Posterior predicton plots - show model-based predictions agains raw data, or otherwise display the error in prediction.
5.1.3.1. Predictor residual plots
A predictor variable residual is the average prediction error when we use all of the other predictor variables to model a predictor of interest. The benifit of computing these things is that, once plotted against the outcome, we have a bivariate regression of sorts that has already ‘controlled’ for all of the other predictor variables. It just leaves in the variation that is not expected by the model of the mean, \(\mu\), as a function of the other predictors.
In our multivariate model of divorce rate, we have two predictors: 1) marriage rate (Marriage.s) and 2) median age at marriage (MedianAgeMarriage.s). To compute predictor residuals for either, we just use the other predictor to model it. So for marriage rate, this is the model we need:
\[R_i {\sim} Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta A_i\\ \alpha {\sim} Normal(0, 10)\\ \beta {\sim} Normal(0, 1)\\ \sigma {\sim} Uniform(0, 10)\]
As before \(R\) is for marriage rate and \(A\) is for age at marriage.
This code will fit the model:
m5.4 <- map(
alist(
Marriage.s ~ dnorm(mu, sigma),
mu <- a + b*MedianAgeMarriage.s,
a ~ dnorm(0, 10),
b ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)And then we compute the residuals by substracting the observed marriage rate in each State from the predicted rate, base upon usign age at marriage:
#compute expected MAP, for each state
mu <- coef(m5.4)["a"] + coef(m5.4)["b"] * d$MedianAgeMarriage.s
#compute residual for each state
m.resid <- d$Marriage.s - muWhen a residual is positive, that means that the observed rate was in excess of what we’d expect, given the median age at marriage in that State (and below when negative). States with positive redisuals marry fast for their age of marriage and vise versa.
It will help to plot the relationship between these two variables, and show the residuals as well.
plot(Marriage.s ~ MedianAgeMarriage.s, d, col=rangi2)
abline(m5.4)only using the first two of 3 regression coefficients#loop over states
for(i in 1:length(m.resid)) {
x <- d$MedianAgeMarriage.s[i] # x location of line segment
y <- d$Marriage.s[i] # observed endpoint of line segment
#draw the line segment
lines(c(x, x), c(mu[i], y), lwd = 0.5, col = col.alpha("black", 0.7))
}
Redisual marriage rate in each states, after accounting for the linear association with median age at marriage. Each gray line represent a residual, the distance of each observed marriage rate from the expected value, attemping to predict marriage rate with median age at marriage alone. So states that lie above the black regression line have higher rates of marriage than expected, according to age of marriage, and vise versa.
Notice that the residuals are variation in marriage rate that is left over, after taking out the purely linear relationship between two variables.
To use these residuals, let’s put then on a horizontal axis and plot them agains the actual outcome of interest, divorce rate.
PLOT HERE
You can think of this lot as displaying the linear relationship between divorce and marriage rates, having statistically “controlled” for media age of marriage. The vertical dashed line indicates marriage rate that exactly matches the expectation from median age at marriage. So states to the right of the line marry faster than expected; to the left - slower. Average divorce rate on both sides of the line is abot the same, and so the regression line demonstates little relationship between divorce and marriage rates. The slope of this line is -.13, exactly what we found in the multivariate model, m5.3.
PLOT HERE
This plot shows the same kind of calculation, but now for median age at marriage, “controlling” for marriage rate. So states to the right have older-than-expected median age at marriage, than those of the left. Now we find that the average divorce rate on the right is lower than on the left, as indicated by the regression line. States in which people marry older than expected for a given rate of marriage ten to have less divorce. The slope of the regression line here is -1.13, again the same as in the multivariate model, m5.3.
5.1.3.2 Counterfactual plots.
These plots show the implied predictions of the model. We call them counterfactual, because they can be produced for any values of the predictor variable you line, even unobserved or impossible combinations like very high median age or marriage and very high marriage rate.
The simplest use of a counterfactual plot is to see how the predictions change as you change only one predictor at a time. This means holding the values of all predictors constant, except for a single predictor of interest.
Let’s draw a pair of counterfactual plots for the divorce model. Beginning with a plot showing the impact of changes in Marriage.s on predictions:
#prepare new counterfactual data
A.avg <- mean(d$MedianAgeMarriage.s)
R.seq <- seq(from = -3, to = 3, length.out = 30)
pred.data <- data.frame(
Marriage.s=R.seq,
MedianAgeMarriage.s=A.avg
)
#compute counterfactual mean divorce rate (mu)
mu <- link(m5.3, data = pred.data)[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI)
#simulate counterfactual divorce outcomes
R.sim <- rethinking::sim(m5.3, data = pred.data, n = 1e4)[ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]R.PI <- apply(R.sim, 2, PI)
#display predictions, hiding raw data with type = 'n'
plot(Divorce ~ Marriage.s, data = d, type = "n")
mtext("MedianAgeMarriage.s = 0")lines(R.seq, mu.mean)
shade(mu.PI, R.seq)shade(R.PI, R.seq)
The strategy above is to build a new list of data that describe the counterfactual cases we wish to simulate predictions for. The list names pred.data holds these cases. Not that the observed values for MedianAgeMarriage.s are not used. Instead we compute the average value and then use it inside the linear model. So Marriage.S changes across the values of R.seq, while the other predictor is held constant as it mean, A.avg.
Let’s do the same for MedianAgeMarriage.s which are allowed to vary and Marriage.s being held constant:
R.avg <- mean(d$Marriage.s)
A.seq <- seq(from = -3, to = 3.5, length.out = 30)
pred.data2 <- data.frame(
Marriage.s = R.avg,
MedianAgeMarriage.s = A.seq
)
#compute counterfactual mean divorce rate (mu)
mu <- link(m5.3, data = pred.data2)[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI)
#simulate counterfactual divorce outcomes
A.sim <- rethinking::sim(m5.3, data = pred.data2, n = 1e4)[ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]A.PI <- apply(A.sim, 2, PI)
#display predictions, hiding raw data with type = 'n'
plot(Divorce ~ MedianAgeMarriage.s, data = d, type = "n")
mtext("MedianAgeMarriage.s = 0")lines(A.seq, mu.mean)
shade(mu.PI, A.seq)shade(A.PI, A.seq)
These plots have the same splots as the resudual plots, but they don’t display any data, raw or residual, because they are counterfactual. And they also show percentile intervals on the scale of the data, instead of on the redisual scale. As the result, they are dicrect displays of the impact on prediction of a change in each variable.
#5.1.3.3. Posterior prediction plots
In addition to understanting the estimates, it’s importatn to check the model fit agains the observed data.
- Did the model fit correctly? (Compare implied predictions to the raw data)
- how does the model fail? (By inspecting the individual cases where the model maeks poor predictions, you might get an idea of how to improve the model)
Let’s begin by simulation prediction, averaing over the posterior.
#call rethinking::link() without specifying new data, so it uses original data
mu <- link(m5.3)[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]#summarise samples across cases
mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI)
#simulate observations, again no new data, so uses original data
divorce.sim <- rethinking::sim(m5.3, n = 1e4)[ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]divorce.PI <- apply(divorce.sim, 2, PI)The code is similar we’ve seen before but now we use the original observed data.
For multivariate models, there are many different ways to diplay these simulations. The simplest is to just plot predictions against observed.
plot(mu.mean ~ d$Divorce, col = rangi2, ylim = range(mu.PI),
xlab = "Observed divorce", ylab = "Predicted divorce")
abline(a = 0, b = 1, lty = 2)for(i in 1:nrow(d)) {
lines(rep(d$Divorce[i], 2), c(mu.PI[1, i], mu.PI[2, i]), col = rangi2)
}
It’s easy to see from the plot that the model under-predicts for States with very high divorce rate while over-predicts for States with very lower divorce rate.
plot(mu.mean ~ d$Divorce, col = rangi2, ylim = range(mu.PI), xlab = "Observed divorce", ylab = "Predicted divorce")
abline(a = 0, b = 1, lty = 2)for(i in 1:nrow(d)) {
lines(rep(d$Divorce[i], 2), c(mu.PI[1, i], mu.PI[2, i]), col = rangi2)
}
#doesn't work in notebook, only as a part of script or in console (better in console)
#just click on several points and then right-click in console (or ESC) to see the labels
identify(x = d$Divorce, y = mu.mean, labels = as.character(d$Loc), cex = 0.8)integer(0)
This plot makes is hard to see the amount of prediction error. For this reason, we’ll use the residual plot that show the mean prediction error for each row.
#compute residuals
divorce.resid <- d$Divorce - mu.mean
#get ordering by divorce rate
o <- order(divorce.resid)
#make the plot
dotchart(divorce.resid[o], labels = d$Loc[o], xlim = c(-6, 5), cex = 0.6)
abline(v = 0, col = col.alpha("black", 0.2))for(i in 1:nrow(d)) {
j <- o[i] # which state in order
lines(d$Divorce[j] - c(mu.PI[1, j], mu.PI[2, j]), rep(i, 2))
points(d$Divorce[j] - c(divorce.PI[1, j], divorce.PI[2, j]), rep(i, 2), pch = 3, cex=0.6, col = "gray")
}
It’s much easier to see the large model failures, such as ID and ME.
Another common use for these simulations is to constuct novel predictor residual plots. Once you’ve computed the divorce residuals (as done above), you can plot those residuals against new predictor variables. This is a quick way to see if remaining variation in theoutcome is associated with another predictor.
5.2. Masked relationship
The divorce rate example demonstrates that multiple predictor variables are useful for knocking out spurious association. A second reason to use more than one predictor variable is to measure the direct influences of multiple factors on an outcome, when non of those infuences is apparent from bivariate relationships. This problem tends to aries when there’re two predictor variables that are correlated with one another, but one correlated positively and other negatively.
Let’s see an example of milk:
library(rethinking)
data("milk")
d <- milk
glimpse(d)Observations: 29
Variables: 8
$ clade <fct> Strepsirrhine, Strepsirrhine, Strepsirrhine, Strepsirrhine, Strepsirrhine, New Wo...
$ species <fct> Eulemur fulvus, E macaco, E mongoz, E rubriventer, Lemur catta, Alouatta seniculu...
$ kcal.per.g <dbl> 0.49, 0.51, 0.46, 0.48, 0.60, 0.47, 0.56, 0.89, 0.91, 0.92, 0.80, 0.46, 0.71, 0.7...
$ perc.fat <dbl> 16.60, 19.27, 14.11, 14.91, 27.28, 21.22, 29.66, 53.41, 46.08, 50.58, 41.35, 3.93...
$ perc.protein <dbl> 15.42, 16.91, 16.85, 13.18, 19.50, 23.58, 23.46, 15.80, 23.34, 22.33, 20.85, 25.3...
$ perc.lactose <dbl> 67.98, 63.82, 69.04, 71.91, 53.22, 55.20, 46.88, 30.79, 30.58, 27.09, 37.80, 70.7...
$ mass <dbl> 1.95, 2.09, 2.51, 1.62, 2.19, 5.25, 5.37, 2.51, 0.71, 0.68, 0.12, 0.47, 0.32, 0.6...
$ neocortex.perc <dbl> 55.16, NA, NA, NA, NA, 64.54, 64.54, 67.64, NA, 68.85, 58.85, 61.69, 60.32, NA, N...The popular hypothesis has it that primates with larger brains product more energetic milk, so that brains can grow quickly. So the variables we’ll consider now are:
- kcal.per.g - Kilocaliries of energy per gram of milk
- mass - Average female body mass, in kg
- neocortex.perc - the percent of total brain mass that is neocortex mass.
The question here is to what extent energy content of milk, measured here by kilocalories, is related to the percent of the brain mass that is neocortex. Female mass is used here to see the masking that hides relationship among the variable.
But before fitting a model, let’s exclude the missing values from the data.
d <- d[complete.cases(d), ]m5.5 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + bn*neocortex.perc,
a ~ dnorm(0, 100),
bn ~ dnorm(0, 1),
sigma ~ dunif(0, 1)
), data = d
)
precis(m5.5, digits = 3) Mean StdDev 5.5% 94.5%
a 0.353 0.471 -0.399 1.106
bn 0.005 0.007 -0.007 0.016
sigma 0.166 0.028 0.120 0.211Next we can plot the predicted mean and 89% interval:
np.seq <- 0:100
pred.data <- data.frame(neocortex.perc = np.seq)
mu <- link(m5.5, data = pred.data, n = 1e4)[ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI)
plot(kcal.per.g ~ neocortex.perc, data = d, col = rangi2)
lines(np.seq, mu.mean)lines(np.seq, mu.PI[1, ], lty = 2)
lines(np.seq, mu.PI[2, ], lty = 2)
MAP line is weakly positive, but it’s highly imprecise.
Now let’s consider another predictor variable, adult female body mass, mass in our data frame. Let’s use log of mass as a predictor.
d$log.mass <- log(d$mass)
m5.6 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + bm*log.mass,
a ~ dnorm(0, 100),
bm ~ dnorm(0, 1),
sigma ~ dunif(0, 1)
), data = d
)
precis(m5.6, digits = 3) Mean StdDev 5.5% 94.5%
a 0.705 0.049 0.627 0.783
bm -0.032 0.020 -0.064 0.001
sigma 0.157 0.027 0.114 0.200So log-mass is negatively correlated with kilocalories and this influence is stronger than neocortex one.
Now let’s see what happens when we add both predictor variables at the same time to the regression:
m5.7 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + bm*log.mass + bn*neocortex.perc,
a ~ dnorm(0, 100),
bn ~ dnorm(0, 1),
bm ~ dnorm(0, 1),
sigma ~ dunif(0, 1)
), data = d
)
precis(m5.7, digits = 3) Mean StdDev 5.5% 94.5%
a -1.084 0.468 -1.832 -0.337
bn 0.028 0.007 0.016 0.040
bm -0.096 0.022 -0.132 -0.060
sigma 0.115 0.020 0.083 0.146By incorporating both predictor variable in the regression, the estimated association of both with the outcome has increased.
Let’s plot the intervals for the predicted mean kilocalories, for this new model.
mean.log.mass <- mean(log(d$mass))
mean.neocortex.perc = mean(d$neocortex.perc)
np.seq <- 0:100
#neocortex perc
pred.data <- data.frame(
neocortex.perc = np.seq,
log.mass = mean.log.mass
)
mu <- link(m5.7, data = pred.data, n = 1e4)[ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI)
plot(kcal.per.g ~ neocortex.perc, data = d, type = "n", main = "Neocortex Perc vs kcal per gramm ")
lines(np.seq, mu.mean)lines(np.seq, mu.PI[1, ], lty = 2)
lines(np.seq, mu.PI[2, ], lty = 2)
#log mass
np.seq <- -100:100
pred.data <- data.frame(
log.mass = np.seq,
neocortex.perc = mean.neocortex.perc
)
mu <- link(m5.7, data = pred.data, n = 1e4)[ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI)
plot(kcal.per.g ~ log.mass, data = d, type = "n", main = "Log Mass vs kcal per gramm ")
lines(np.seq, mu.mean)lines(np.seq, mu.PI[1, ], lty = 2)
lines(np.seq, mu.PI[2, ], lty = 2)
Since one variable is positively correlated with the outcome and another is negatively correlated, they tend to cancel one another out.
5.3 When adding variables hurts
There’re several good, purely statistical reasons to avoid doing this.
- Multicollinearity - strong correlation between two or more predictor variables. The consequence of it is that the posterior distribution will say that a veray large range of parameter are plausible.
- Post-treatment bias - statistically controlling for consequences of a causal factor.
- Overfitting
5.3.1. Multicollinear legs
N <- 100 # number of individuals
height <- rnorm(N, 10, 2) #sim total height of each
leg_prop <- runif(N, 0.4, 0.5) #leg as proportion of heigh
leg_left <- leg_prop * height + rnorm(N, 0, 0.02) #sim left leg as proportion + error
leg_right <- leg_prop * height + rnorm(N, 0, 0.02) #sim right leg as proportion + error
d <- data.frame(height, leg_left, leg_right)Now let’s analyze these data, prediction the outcome height with both predictions, leg_left and leg_right. Our expectation that beta coefficient is equal to 10 / 4.5 = 2.2 (4.5 = 10 * 0.45 - average proportion of the leg).
m5.8 <- map(
alist(
height ~ dnorm(mu, sigma),
mu <- a + bl*leg_left + br*leg_right,
a ~ dnorm(10, 100),
bl ~ dnorm(2, 10),
br ~ dnorm(2, 10),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.8) Mean StdDev 5.5% 94.5%
a 1.46 0.33 0.93 1.98
bl 2.35 2.16 -1.10 5.81
br -0.45 2.16 -3.90 2.99
sigma 0.62 0.04 0.55 0.70precis_plot(precis(m5.8))
Despite the fact that the precis plot looks odd, model did fit correctly, and the posterior distribution here is the right answer to the question we asked - what is the value of knowing each predictor, after already knowing all the other predictors? The answer to the question is wierd, but perfectly logical. The posterior distributino is the answer to this question, considering every possible combination of the parameters and assigning relative plausibilities to every combination, conditional on this model and these data. It might help to look at the bivariate posterio distribution for bl and br:
post <- extract.samples(m5.8)
plot(bl ~ br, post, col = col.alpha(rangi2, 0.1), pch = 16)
So the posterior distribution for these two parameters is very highly correlated - when br is large, bl is small.
One way to think about this phenomenon is that you have approximated this likelihood:
\[y_i {\sim} Normal(\mu_i, \sigma)\\ m_i = \alpha + \beta_1x_i + \beta_2x_i\]
Here is a single predictor, and it’s used twice. From the computer perspective this likelihood is really:
\[y_i {\sim} Normal(\mu_i, \sigma)\\ m_i = \alpha + (\beta_1 + \beta_2)x_i\]
So neither \(\beta_1\) nor \(\beta_2\) influence \(\mu\), only their sum \(\beta_1 + \beta_2\) influences it. This means that there’re infinite combinations of \(\beta_1\) and \(\beta_2\).
Just to check, let’s plot the sum:
sum_blbr <- post$bl + post$br
dens(sum_blbr, col = rangi2, lwd = 2, xlab = "sum of bl and br")
If you fit a regression with only one of the leg length, you’ll ge approximately the same posterior mean:
m5.9 <- map(
alist(
height ~ dnorm(mu, sigma),
mu <- a + bl * leg_left,
a ~ dnorm(10, 100),
bl ~ dnorm(2, 10),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.9) Mean StdDev 5.5% 94.5%
a 1.46 0.33 0.94 1.99
bl 1.90 0.07 1.78 2.01
sigma 0.62 0.04 0.55 0.70The lesson here is - when two predictor variables are very strongly correlated, including both in a model may lead to confusion.
5.3.2 Multicollinear milk
library(rethinking)
data(milk)
d <- milkIn these data, we have variables perc.fat and perc.lactose that we might use to model the total energy content, kcal.per.g to explore a natural case of multicollinearity.
Let’s start by odeling kcal.per.g as a function of perc.fact and perc.lactose, but in two bivariate regressions:
m5.10 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + bf * perc.fat,
a ~ dnorm(0.6, 10),
bf ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)
m5.11 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + bl * perc.lactose,
a ~ dnorm(0.6, 10),
bl ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.10, digits = 3) Mean StdDev 5.5% 94.5%
a 0.301 0.036 0.244 0.358
bf 0.010 0.001 0.008 0.012
sigma 0.073 0.010 0.058 0.089precis_plot(precis(m5.10))
precis(m5.11, digits = 3) Mean StdDev 5.5% 94.5%
a 1.166 0.043 1.098 1.235
bl -0.011 0.001 -0.012 -0.009
sigma 0.062 0.008 0.049 0.075precis_plot(precis(m5.11))
The posterior means are essentially mirror images of one another.
Let’s fit a multivariate regression:
m5.12 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + bl * perc.lactose + bf * perc.fat,
a ~ dnorm(0.6, 10),
bl ~ dnorm(0, 1),
bf ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.12) Mean StdDev 5.5% 94.5%
a 1.01 0.20 0.69 1.33
bl -0.01 0.00 -0.01 0.00
bf 0.00 0.00 0.00 0.01
sigma 0.06 0.01 0.05 0.07precis_plot(precis(m5.12))
Now the posterior means for both bf and bl are closer to zero. This happens because the variables perc.fact and perc.lactose contain much of the same information. In the case of fat and lactose, there two variables form essentially a single axis of the variation:
pairs(~ kcal.per.g + perc.fat + perc.lactose, data = d, col = rangi2)
Two variables are negatively correlated, and so strongly so that they are nearly redundant. Either helps in predicting kcal.per.g, but neither helps much once you already know the other.
You can compute the correlation between two variables with cor:
cor(d$perc.fat, d$perc.lactose)[1] -0.9416373What can be done about multicollinearity? The best thing to do is be aware of it. You can anticipate this problem by checking the predictor variables against one another in a pairs plot.
5.3.3 Post-treatment bias
It’s routine to worry about mistaken inferences that arise from omitting predictor variables. Such mistakes are often called omitted variable bias. It’s much less routing to worry about mistaken inferences arising from including variables that are consequences of other variables. We’ll call this post-treatment bias.
Plants in a greenhouse example Different in growth under different antifungail soil treatment. Four variables: initial height, final height, treatment, and presence of fungus. If your goal is to make a causal inference about the treatment, you shouldn’t include the presence of fungus, because it is a post-treatment effect.
Let’s simulate some data:
#number of plants
N <- 100
#simulate initial heights
h0 <- rnorm(N, 10, 2)
#assign treatments and simulate fungus and growth
treatment <- rep(0:1, each = N/2)
fungus <- rbinom(N, size = 1, prob = 0.5 - treatment * 0.4)
h1 <- h0 + rnorm(N, 5 - 3 * fungus)
#compose a clean data frame
d <- data.frame(h0 = h0, h1 =h1, treatment = treatment, fungus = fungus)Then fit a model:
m5.13 <- map(
alist(
h1 ~ dnorm(mu, sigma),
mu <- a + bh * h0 + bt * treatment + bf * fungus,
a ~ dnorm(0, 100),
c(bh, bt, bf) ~ dnorm(0, 10),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.13) Mean StdDev 5.5% 94.5%
a 5.49 0.56 4.59 6.39
bh 0.96 0.05 0.89 1.04
bt -0.15 0.20 -0.47 0.18
bf -3.43 0.21 -3.77 -3.09
sigma 0.90 0.06 0.79 1.00precis_plot(precis(m5.13))
The marginal posterior for bt, the effect of treatment, is actually negative here with little effect, while other predictors have important effect. The problem is that fungus is mostly a consequence of treatment. So when we control for fungus, the model is implicitly answering the question: Once we already know whether or not a plant developed fungus, does soil treatment matters?. The answer is ‘no’, because soil treatment has its effects on growth through reducing fungus.
What we actually want to know is the impact of treatment on growth.
m5.14 <- map(
alist(
h1 ~ dnorm(mu, sigma),
mu <- a + bh * h0 + bt * treatment,
a ~ dnorm(0, 100),
c(bh, bt) ~ dnorm(0, 10),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.14) Mean StdDev 5.5% 94.5%
a 1.87 0.97 0.31 3.42
bh 1.12 0.09 0.98 1.27
bt 1.40 0.34 0.86 1.94
sigma 1.69 0.12 1.50 1.88precis_plot(precis(m5.14))
Now the impact of treatment is strong and positive.
5.4. Categorical variables
A common question for statistical methods is to what extent an outcome changes asa result of presence or absence of a category.
5.4.1. Binary categories
Let’s return to Kalahari height data
data("Howell1")
d <- Howell1
str(d)'data.frame': 544 obs. of 4 variables:
$ height: num 152 140 137 157 145 ...
$ weight: num 47.8 36.5 31.9 53 41.3 ...
$ age : num 63 63 65 41 51 35 32 27 19 54 ...
$ male : int 1 0 0 1 0 1 0 1 0 1 ...The male variable is our new predictor, an example of dummy variable.
\[h_i {\sim} Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta_mm_i\\ \alpha {\sim} Normal(178, 100)\\ \beta_m {\sim} Normal(0, 10)\\ \sigma {\sim} Uniform(0, 50)\]
where \(m\) is a dummy variable indicating a male individual. The parameter \(\beta_m\) influences prediction only for those cases where \(m_i = 1\).
Let’s fit it:
m5.15 <- map(
alist(
height ~ dnorm(mu, sigma),
mu <- a + bm * male,
a ~ dnorm(178, 100),
bm ~ dnorm(0, 10),
sigma ~ dunif(0, 50)
), data = d
)
precis(m5.15) Mean StdDev 5.5% 94.5%
a 134.83 1.59 132.28 137.37
bm 7.29 2.28 3.64 10.94
sigma 27.31 0.83 25.99 28.63precis_plot(precis(m5.15))
To interpret these estimates, you have to note that the parameter \(\alpha\) is now the average height among females. The parameter \(\beta_m\) then tells us the average difference between males and females = 7.3 cm.
Let’s sample from the posterior to derive a percentile inverval for average male height.
post <- extract.samples(m5.15)
mu.male <- post$a + post$bm
PI(mu.male) 5% 94%
139.4139 144.7717 5.4.2. Many categories
When we have k categories, we need k - 1 dummy variables.
Get back to primate milk example with 4 types of primates:
data(milk)
d <- milk
unique(d$clade)[1] Strepsirrhine New World Monkey Old World Monkey Ape
Levels: Ape New World Monkey Old World Monkey StrepsirrhineLet’s create dummy variables:
library(dplyr)
d <- d %>%
mutate(clade.NWM = if_else(clade == "New World Monkey", 1, 0),
clade.OWM = if_else(clade == "Old World Monkey", 1, 0),
clade.S = if_else(clade == "Strepsirrhine", 1, 0),
)There’s no need to make another variable for Ape, because it will be the default, “intercept” category.
Let’s specify the model for kcal.per.g:
\[k_i {\sim} Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta_{NWM}NWM_i + \beta_{OWM}OWM_i + \beta_SS_i\\ \alpha {\sim} Normal(0.6, 10)\\ \beta_{NWM}{\sim} Normal(0, 1)\\ \beta_{OWM}{\sim} Normal(0, 1)\\ \beta_{S}{\sim} Normal(0, 1)\\ \sigma {\sim} Uniform(0, 10)\\\]
A linear model like this really defines four different linear models, each corresponsing to a different category.
m5.16 <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a + b.NWM * clade.NWM + b.OWM * clade.OWM + b.S * clade.S,
a ~ dnorm(0.6, 10),
c(b.NWM, b.OWM, b.S) ~ dnorm(0, 1),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.16) Mean StdDev 5.5% 94.5%
a 0.55 0.04 0.49 0.61
b.NWM 0.17 0.05 0.08 0.25
b.OWM 0.24 0.06 0.15 0.34
b.S -0.04 0.06 -0.14 0.06
sigma 0.11 0.02 0.09 0.14precis_plot(precis(m5.16))
To get posterior distribution for the average milk energy in each category, you can again use samples:
#sample posterior
post <- extract.samples(m5.16)
#compute averages for each category
mu.ape <- post$a
mu.NWM <- post$a + post$b.NWM
mu.OWM <- post$a + post$b.OWM
mu.S <- post$a + post$b.S
#summarize using precis
precis(data.frame(mu.ape, mu.NWM, mu.OWM, mu.S)) Mean StdDev |0.89 0.89|
mu.ape 0.55 0.04 0.49 0.61
mu.NWM 0.71 0.04 0.65 0.77
mu.OWM 0.79 0.05 0.71 0.86
mu.S 0.51 0.05 0.43 0.59Suppose we want to know the estimated difference between the two monkey groups. Then just subtract the estimated means to get a difference:
diff.NWM.OWM <- mu.NWM - mu.OWM
quantile(diff.NWM.OWM, probs = c(0.025, 0.5, 0.975)) 2.5% 50% 97.5%
-0.19300099 -0.07438118 0.04406231 5.4.3. Adding regular predictor variables
There’s not obstacle now in including other predictor variables, like perc.fat or log(mass) to the model.
\[\mu_i = \alpha + \beta_{NWM}NWM_i + \beta_{OWM}OWM_i + \beta_SS_i + \beta_FF_i\]
where \(F_i\) is the percent fat for the i-th case in the data and \(\beta_F\) is the slope that measure the influence of \(F\) on predictions about the mean milk energy.
5.4.4. Unique intercepts.
Another way to conceptualize categorical variables is to construct a vector of intercept parameters, one parameter for each category. Then you can create an index variable in your data frame that says which parameter goes with each case.
(d$clade_id <- coerce_index(d$clade)) [1] 4 4 4 4 4 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1This variable just gives the number of unique clade value. Then you tell map to maek a vector of intercepts, one intercept for each unique values in clade_id:
m5.16_alt <- map(
alist(
kcal.per.g ~ dnorm(mu, sigma),
mu <- a[clade_id],
a[clade_id] ~ dnorm(0.6, 10),
sigma ~ dunif(0, 10)
), data = d
)
precis(m5.16_alt, depth = 2) Mean StdDev 5.5% 94.5%
a[1] 0.55 0.04 0.48 0.61
a[2] 0.71 0.04 0.65 0.78
a[3] 0.79 0.05 0.71 0.86
a[4] 0.51 0.05 0.43 0.59
sigma 0.11 0.02 0.09 0.14precis_plot(precis(m5.16_alt, depth = 2))
Now you getn the same averages for each clade but this time directly from the fit model. You need to add depth = 2 to precis call in order to print out vector parameters like these.
5.5. Ordinary least squares and lm
If you ok with flat priors.
5.5.1. Design formulas
\[y_i {\sim} Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta x_i\]
is equivalent to y ~ 1 + x.
Intercepts are optional y ~ x = y = 1 + x
To omit intercept y ~ 0 + x
Categorical variables y ~ 1 + as.factor(x)
Transform variables y ~ 1 + x + I(x^2) + I(x^3)
5.5.3 Building map formulas from lm formulas with rethinking::glimmmer()
data(cars)
glimmer(dist ~ speed, data = cars)alist(
dist ~ dnorm( mu , sigma ),
mu <- Intercept +
b_speed*speed,
Intercept ~ dnorm(0,10),
b_speed ~ dnorm(0,10),
sigma ~ dcauchy(0,2)
)