What is multiple regression doing?
If you came to this page directly, please take a moment to look at this blog post for the comment. I’ll wait until you’re back.
OK. Now you’ve read the blogpost, and you know why we’re here. I’m going to illustrate how it is that multiple regression separates the effects that are “really there” from those that are only there because of statistical associations. We’ll use exactly the same data aw we did last time, but we’ll only analyse the strong association scenario.{^1] Everything in this block of R code is just setting up the data set.
library(tidyverse)
library(reshape2)
library(corrplot)
library(ggplot2)
library(cowplot)
library(WVPlots)
library(mvtnorm)## I always like to clean out my workspace before running a script to make sure
## that I'm starting R in the same state. This helps to ensure that I can
## reproduce my results later
rm(list = ls())
## intetcept
##
beta0 <- 1.0
## regression coefficients
##
beta <- c(1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
## pattern of correlation matrix, all non-zero entries are set to saem
## correlation, covariance matrix caldulated from individual variances and a
## single association parameter governing the non-zero correlation coefficients
##
## Note: Not just any pattern will work here. The correlation matrix and
## covariance matrix generated from this pattern must be positive definite.
## If you change this pattern, you may get an error when you try to generate
## data with a non-zero association parameter.
##
Rho <- matrix(nrow = 9, ncol = , byrow = TRUE,
data = c(1,0,1,0,1,0,1,0,1,
0,1,0,1,0,1,0,1,0,
1,0,1,0,1,0,1,0,1,
0,1,0,1,0,1,0,1,0,
1,0,1,0,1,0,1,0,1,
0,1,0,1,0,1,0,1,0,
1,0,1,0,1,0,1,0,1,
0,1,0,1,0,1,0,1,0,
1,0,1,0,1,0,1,0,1
))
## vector of standard deviations for covariates
##
sigma <- rep(1, 9)
## construct a covariance matrix from the pattern, standard deviations, and
## one parameter in [-1,1] that governs the magnitude of non-zero correlation
## coefficients
##
## Rho - the pattern of associations
## sigma - the vector of standard deviations
## rho - the association parameter
##
construct_Sigma <- function(Rho, sigma, rho) {
## get the correlation matris
##
Rho <- Rho*rho
for (i in 1:ncol(Rho)) {
Rho[i,i] <- 1.0
}
## notice the use of matrix multiplication
##
Sigma <- diag(sigma) %*% Rho %*% diag(sigma)
return(Sigma)
}
## set the random number seed manually so that every run of the code will
## produce the same numbers
##
set.seed(1234)
n_samp <- 100
cov_ind <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.0))
cov_wk <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.2))
cov_str <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.8))
resid <- rep(0.2, n_samp)
y_ind <- rnorm(nrow(cov_ind), mean = beta0 + cov_ind %*% beta, sd = resid)
y_wk <- rnorm(nrow(cov_wk), mean = beta0 + cov_wk %*% beta, sd = resid)
y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_str <- data.frame(y_str, cov_str, rep("Strong", length(y_str)))
column_names <- c("y", paste("x", seq(1, length(beta)), sep = ""), "Scenario")
colnames(dat_str) <- column_namesFirst, let’s double check the multiple regression analysis to make sure that the data we’re using really are the same.
covariates <- paste("x", seq(1, length(beta)), sep = "")
summary(lm(as.formula(paste("y ~ ", paste(covariates, collapse = " + "))), data = dat_str))
Call:
lm(formula = as.formula(paste("y ~ ", paste(covariates, collapse = " + "))),
data = dat_str)
Residuals:
Min 1Q Median 3Q Max
-0.47277 -0.13864 0.01532 0.13314 0.50245
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.009106 0.021042 47.957 <2e-16 ***
x1 0.981123 0.044866 21.868 <2e-16 ***
x2 -1.000192 0.037447 -26.710 <2e-16 ***
x3 0.975234 0.039012 24.998 <2e-16 ***
x4 0.022840 0.045962 0.497 0.620
x5 0.036022 0.043129 0.835 0.406
x6 -0.012031 0.041184 -0.292 0.771
x7 0.002427 0.040711 0.060 0.953
x8 0.002452 0.042054 0.058 0.954
x9 -0.006829 0.037163 -0.184 0.855
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2003 on 90 degrees of freedom
Multiple R-squared: 0.9913, Adjusted R-squared: 0.9905
F-statistic: 1144 on 9 and 90 DF, p-value: < 2.2e-16If you compare those results with the ones we obtained before,1 you’ll see that they’re identical to 6 decimal places. That makes it pretty likely that they’re the same data sets.
We’ll focus this example on x1, x3, and x9, where we know that the “real” associations are between x1, x3, and y and that the “apparent” association between x9 and y arises only because x9 is strongly associated with x1 and x3. Let’s start by looking at the bivariate association between x1 and y, including a plot of the residuals.
lm_str <- lm(y ~ x1, data = dat_str)
residual <- residuals(lm_str)
for.plot <- data.frame(x = dat_str$x1, y = residual)
p <- ggplot(for.plot, aes(x = x, y = y)) +
geom_point(fill = "tomato", color = "tomato") +
geom_hline(yintercept = 0.0) +
xlab("Observed (x1)") +
ylab("Residual")
print(p)
As you can see, the residuals seem pretty randomly distributed, which is what we hope to see when we do a residual plot. But let’s try this. Let’s try regressing the residuals on x3. Think of this as looking at the influence of x3 on y once we’ve removed the influence of x1.
res_x3 <- lm(residual ~ x3, data = dat_str)
print(summary(res_x3))
Call:
lm(formula = residual ~ x3, data = dat_str)
Residuals:
Min 1Q Median 3Q Max
-2.8942 -0.8188 0.1138 0.8285 2.3859
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.349e-05 1.127e-01 -0.001 0.9995
x3 3.419e-01 1.182e-01 2.893 0.0047 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.127 on 98 degrees of freedom
Multiple R-squared: 0.0787, Adjusted R-squared: 0.0693
F-statistic: 8.371 on 1 and 98 DF, p-value: 0.004697Look at that. We have good evidence that the residuals from regressing y on x1 show a relationship with x3 - the larger x3, the larger the residual. We can see this both by coloring the residuals above and by plotting the residuals against x3.
for.plot <- data.frame(x = dat_str$x1, y = residual, x3 = dat_str$x3)
p <- ggplot(for.plot, aes(x = x, y = y, fill = x3, color = x3)) +
geom_point() +
scale_fill_gradient2() +
scale_color_gradient2() +
geom_hline(yintercept = 0.0) +
xlab("Observed (x1)") +
ylab("Residual")
print(p)
p <- ggplot(for.plot, aes(x = x3, y = residual)) +
geom_point(color = "tomato", fill = "tomato") +
geom_smooth(method = "lm") +
xlab("Observed (x3)") +
ylab("Residual (from x1)")
print(p)
Now let’s see what happens if we do the same thing with x9.
res_x9 <- lm(residual ~ x9, data = dat_str)
print(summary(res_x9))
Call:
lm(formula = residual ~ x9, data = dat_str)
Residuals:
Min 1Q Median 3Q Max
-3.5085 -0.8446 0.0678 0.7904 2.7547
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01417 0.11669 -0.121 0.904
x9 0.15777 0.11294 1.397 0.166
Residual standard error: 1.163 on 98 degrees of freedom
Multiple R-squared: 0.01952, Adjusted R-squared: 0.00952
F-statistic: 1.952 on 1 and 98 DF, p-value: 0.1656for.plot <- data.frame(x = dat_str$x1, y = residual, x9 = dat_str$x9)
p <- ggplot(for.plot, aes(x = x, y = y, fill = x9, color = x9)) +
geom_point() +
scale_fill_gradient2() +
scale_color_gradient2() +
geom_hline(yintercept = 0.0) +
xlab("Observed (x1)") +
ylab("Residual")
print(p)
p <- ggplot(for.plot, aes(x = x9, y = residual)) +
geom_point(color = "tomato", fill = "tomato") +
geom_smooth(method = "lm") +
xlab("Observed (x9)") +
ylab("Residual (from x1)")
print(p)
Here we see that there is a weak positive association between x9 and y, but it isn’t convincing. Once we’ve removed the influence of x1 the remaining association between x9 and y presumably arises because x3 and x9 are associated. Let’s try this again after removing the influence of both x1 and x3.
lm_x1_x3 <- lm(residual ~ x1 +x3, data = dat_str)
print(summary(lm_x1_x3))
Call:
lm(formula = residual ~ x1 + x3, data = dat_str)
Residuals:
Min 1Q Median 3Q Max
-2.8345 -0.9209 0.1885 0.7596 2.1640
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01948 0.10439 0.187 0.852
x1 -0.73398 0.17580 -4.175 6.50e-05 ***
x3 0.95163 0.18244 5.216 1.04e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.043 on 97 degrees of freedom
Multiple R-squared: 0.219, Adjusted R-squared: 0.2029
F-statistic: 13.6 on 2 and 97 DF, p-value: 6.199e-06residual <- residuals(lm_x1_x3)
res_x9 <- lm(residual ~ x9, data = dat_str)
print(summary(res_x9))
Call:
lm(formula = residual ~ x9, data = dat_str)
Residuals:
Min 1Q Median 3Q Max
-2.6948 -0.9243 0.1956 0.7409 2.2469
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.008041 0.103726 -0.078 0.938
x9 0.089548 0.100388 0.892 0.375
Residual standard error: 1.033 on 98 degrees of freedom
Multiple R-squared: 0.008054, Adjusted R-squared: -0.002068
F-statistic: 0.7957 on 1 and 98 DF, p-value: 0.3746for.plot <- data.frame(x = dat_str$x1, y = residual, x9 = dat_str$x9)
p <- ggplot(for.plot, aes(x = x9, y = residual)) +
geom_point(color = "tomato", fill = "tomato") +
geom_smooth(method = "lm") +
xlab("Observed (x9)") +
ylab("Residual (from x1)")
print(p)
Notice that the coefficient on x9 is now only about 0.09 rather than 0.16. There’s still a bit of an association, but it’s weak and poorly supported.
So that’s what lm() is doing. It is “statistically controlling” for the effects of other covariates when estimating the effect of each one individually. It doesn’t fit the regression sequentially, as I’ve done it here. It fits all of the coefficients simultaneously. As a result, it doesn’t make any difference what order you put variables into the model statement. You’ll get the same result.2
Scroll all the way to the bottom of the page, just before the Conclusion at http://darwin.eeb.uconn.edu/pages/variable-selection/multiple-regression-basics.nb.html.↩
If you don’t believe me, try manually changing the order of the variables in this multiple regression and see what you get.↩
LS0tCnRpdGxlOiAiV2hhdCBpcyBtdWx0aXBsZSByZWdyZXNzaW9uIGRvaW5nPyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKSWYgeW91IGNhbWUgdG8gdGhpcyBwYWdlIGRpcmVjdGx5LCBwbGVhc2UgdGFrZSBhIG1vbWVudCB0byBsb29rIGF0IFt0aGlzIGJsb2cgcG9zdF0oIGh0dHA6Ly9kYXJ3aW4uZWViLnVjb25uLmVkdS91bmNvbW1vbi1ncm91bmQvYmxvZy8yMDE5LzA4LzE0L3doYXQtaXMtbXVsdGlwbGUtcmVncmVzc2lvbi1kb2luZy8pIGZvciB0aGUgY29tbWVudC4gSSdsbCB3YWl0IHVudGlsIHlvdSdyZSBiYWNrLgoKT0suIE5vdyB5b3UndmUgcmVhZCB0aGUgYmxvZ3Bvc3QsIGFuZCB5b3Uga25vdyB3aHkgd2UncmUgaGVyZS4gSSdtIGdvaW5nIHRvIGlsbHVzdHJhdGUgaG93IGl0IGlzIHRoYXQgbXVsdGlwbGUgcmVncmVzc2lvbiBzZXBhcmF0ZXMgdGhlIGVmZmVjdHMgdGhhdCBhcmUgInJlYWxseSB0aGVyZSIgZnJvbSB0aG9zZSB0aGF0IGFyZSBvbmx5IHRoZXJlIGJlY2F1c2Ugb2Ygc3RhdGlzdGljYWwgYXNzb2NpYXRpb25zLiBXZSdsbCB1c2UgZXhhY3RseSB0aGUgc2FtZSBkYXRhIGF3IHdlIGRpZCBsYXN0IHRpbWUsIGJ1dCB3ZSdsbCBvbmx5IGFuYWx5c2UgdGhlIHN0cm9uZyBhc3NvY2lhdGlvbiBzY2VuYXJpby57XjFdIEV2ZXJ5dGhpbmcgaW4gdGhpcyBibG9jayBvZiBSIGNvZGUgaXMganVzdCBzZXR0aW5nIHVwIHRoZSBkYXRhIHNldC4KCmBgYHtyIHNldHVwLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkoY29ycnBsb3QpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShjb3dwbG90KQpsaWJyYXJ5KFdWUGxvdHMpCmxpYnJhcnkobXZ0bm9ybSkKYGBgCmBgYHtyfQojIyBJIGFsd2F5cyBsaWtlIHRvIGNsZWFuIG91dCBteSB3b3Jrc3BhY2UgYmVmb3JlIHJ1bm5pbmcgYSBzY3JpcHQgdG8gbWFrZSBzdXJlCiMjIHRoYXQgSSdtIHN0YXJ0aW5nIFIgaW4gdGhlIHNhbWUgc3RhdGUuIFRoaXMgaGVscHMgdG8gZW5zdXJlIHRoYXQgSSBjYW4gCiMjIHJlcHJvZHVjZSBteSByZXN1bHRzIGxhdGVyCnJtKGxpc3QgPSBscygpKQoKIyMgaW50ZXRjZXB0CiMjCmJldGEwIDwtIDEuMAojIyByZWdyZXNzaW9uIGNvZWZmaWNpZW50cwojIwpiZXRhIDwtIGMoMS4wLCAtMS4wLCAxLjAsIDAuMCwgMC4wLCAwLjAsIDAuMCwgMC4wLCAwLjApCiMjIHBhdHRlcm4gb2YgY29ycmVsYXRpb24gbWF0cml4LCBhbGwgbm9uLXplcm8gZW50cmllcyBhcmUgc2V0IHRvIHNhZW0KIyMgY29ycmVsYXRpb24sIGNvdmFyaWFuY2UgbWF0cml4IGNhbGR1bGF0ZWQgZnJvbSBpbmRpdmlkdWFsIHZhcmlhbmNlcyBhbmQgYSAKIyMgc2luZ2xlIGFzc29jaWF0aW9uIHBhcmFtZXRlciBnb3Zlcm5pbmcgdGhlIG5vbi16ZXJvIGNvcnJlbGF0aW9uIGNvZWZmaWNpZW50cwojIwojIyBOb3RlOiBOb3QganVzdCBhbnkgcGF0dGVybiB3aWxsIHdvcmsgaGVyZS4gVGhlIGNvcnJlbGF0aW9uIG1hdHJpeCBhbmQKIyMgY292YXJpYW5jZSBtYXRyaXggZ2VuZXJhdGVkIGZyb20gdGhpcyBwYXR0ZXJuIG11c3QgYmUgcG9zaXRpdmUgZGVmaW5pdGUuCiMjIElmIHlvdSBjaGFuZ2UgdGhpcyBwYXR0ZXJuLCB5b3UgbWF5IGdldCBhbiBlcnJvciB3aGVuIHlvdSB0cnkgdG8gZ2VuZXJhdGUKIyMgZGF0YSB3aXRoIGEgbm9uLXplcm8gYXNzb2NpYXRpb24gcGFyYW1ldGVyLgojIwpSaG8gPC0gbWF0cml4KG5yb3cgPSA5LCBuY29sID0gLCBieXJvdyA9IFRSVUUsIAogICAgICAgICAgICAgIGRhdGEgPSBjKDEsMCwxLDAsMSwwLDEsMCwxLAogICAgICAgICAgICAgICAgICAgICAgIDAsMSwwLDEsMCwxLDAsMSwwLAogICAgICAgICAgICAgICAgICAgICAgIDEsMCwxLDAsMSwwLDEsMCwxLAogICAgICAgICAgICAgICAgICAgICAgIDAsMSwwLDEsMCwxLDAsMSwwLAogICAgICAgICAgICAgICAgICAgICAgIDEsMCwxLDAsMSwwLDEsMCwxLAogICAgICAgICAgICAgICAgICAgICAgIDAsMSwwLDEsMCwxLDAsMSwwLAogICAgICAgICAgICAgICAgICAgICAgIDEsMCwxLDAsMSwwLDEsMCwxLAogICAgICAgICAgICAgICAgICAgICAgIDAsMSwwLDEsMCwxLDAsMSwwLAogICAgICAgICAgICAgICAgICAgICAgIDEsMCwxLDAsMSwwLDEsMCwxCiAgICAgICAgICAgICAgICAgICAgICAgKSkKIyMgdmVjdG9yIG9mIHN0YW5kYXJkIGRldmlhdGlvbnMgZm9yIGNvdmFyaWF0ZXMKIyMKc2lnbWEgPC0gcmVwKDEsIDkpCgojIyBjb25zdHJ1Y3QgYSBjb3ZhcmlhbmNlIG1hdHJpeCBmcm9tIHRoZSBwYXR0ZXJuLCBzdGFuZGFyZCBkZXZpYXRpb25zLCBhbmQKIyMgb25lIHBhcmFtZXRlciBpbiBbLTEsMV0gdGhhdCBnb3Zlcm5zIHRoZSBtYWduaXR1ZGUgb2Ygbm9uLXplcm8gY29ycmVsYXRpb24KIyMgY29lZmZpY2llbnRzCiMjCiMjIFJobyAtIHRoZSBwYXR0ZXJuIG9mIGFzc29jaWF0aW9ucwojIyBzaWdtYSAtIHRoZSB2ZWN0b3Igb2Ygc3RhbmRhcmQgZGV2aWF0aW9ucwojIyByaG8gLSB0aGUgYXNzb2NpYXRpb24gcGFyYW1ldGVyCiMjCmNvbnN0cnVjdF9TaWdtYSA8LSBmdW5jdGlvbihSaG8sIHNpZ21hLCByaG8pIHsKICAjIyBnZXQgdGhlIGNvcnJlbGF0aW9uIG1hdHJpcwogICMjCiAgUmhvIDwtIFJobypyaG8KICBmb3IgKGkgaW4gMTpuY29sKFJobykpIHsKICAgIFJob1tpLGldIDwtIDEuMAogIH0KICAjIyBub3RpY2UgdGhlIHVzZSBvZiBtYXRyaXggbXVsdGlwbGljYXRpb24KICAjIwogIFNpZ21hIDwtIGRpYWcoc2lnbWEpICUqJSBSaG8gJSolIGRpYWcoc2lnbWEpCiAgcmV0dXJuKFNpZ21hKQp9CgoKIyMgc2V0IHRoZSByYW5kb20gbnVtYmVyIHNlZWQgbWFudWFsbHkgc28gdGhhdCBldmVyeSBydW4gb2YgdGhlIGNvZGUgd2lsbCAKIyMgcHJvZHVjZSB0aGUgc2FtZSBudW1iZXJzCiMjCnNldC5zZWVkKDEyMzQpCgpuX3NhbXAgPC0gMTAwCmNvdl9pbmQgPC0gcm12bm9ybShuX3NhbXAsCiAgICAgICAgICAgICAgICAgICBtZWFuID0gcmVwKDAsIG5yb3coUmhvKSksCiAgICAgICAgICAgICAgICAgICBzaWdtYSA9IGNvbnN0cnVjdF9TaWdtYShSaG8sIHNpZ21hLCAwLjApKQpjb3Zfd2sgPC0gcm12bm9ybShuX3NhbXAsCiAgICAgICAgICAgICAgICAgIG1lYW4gPSByZXAoMCwgbnJvdyhSaG8pKSwKICAgICAgICAgICAgICAgICAgc2lnbWEgPSBjb25zdHJ1Y3RfU2lnbWEoUmhvLCBzaWdtYSwgMC4yKSkKY292X3N0ciA8LSBybXZub3JtKG5fc2FtcCwKICAgICAgICAgICAgICAgICAgIG1lYW4gPSByZXAoMCwgbnJvdyhSaG8pKSwKICAgICAgICAgICAgICAgICAgIHNpZ21hID0gY29uc3RydWN0X1NpZ21hKFJobywgc2lnbWEsIDAuOCkpCgpyZXNpZCA8LSByZXAoMC4yLCBuX3NhbXApCnlfaW5kIDwtIHJub3JtKG5yb3coY292X2luZCksIG1lYW4gPSBiZXRhMCArIGNvdl9pbmQgJSolIGJldGEsIHNkID0gcmVzaWQpCnlfd2sgPC0gcm5vcm0obnJvdyhjb3Zfd2spLCBtZWFuID0gYmV0YTAgKyBjb3Zfd2sgJSolIGJldGEsIHNkID0gcmVzaWQpCnlfc3RyIDwtIHJub3JtKG5yb3coY292X3N0ciksIG1lYW4gPSBiZXRhMCArIGNvdl9zdHIgJSolIGJldGEsIHNkID0gcmVzaWQpCgpkYXRfc3RyIDwtIGRhdGEuZnJhbWUoeV9zdHIsIGNvdl9zdHIsIHJlcCgiU3Ryb25nIiwgbGVuZ3RoKHlfc3RyKSkpCgpjb2x1bW5fbmFtZXMgPC0gYygieSIsIHBhc3RlKCJ4Iiwgc2VxKDEsIGxlbmd0aChiZXRhKSksIHNlcCA9ICIiKSwgIlNjZW5hcmlvIikKY29sbmFtZXMoZGF0X3N0cikgPC0gY29sdW1uX25hbWVzCmBgYApGaXJzdCwgbGV0J3MgZG91YmxlIGNoZWNrIHRoZSBtdWx0aXBsZSByZWdyZXNzaW9uIGFuYWx5c2lzIHRvIG1ha2Ugc3VyZSB0aGF0IHRoZSBkYXRhIHdlJ3JlIHVzaW5nIHJlYWxseSBhcmUgdGhlIHNhbWUuCmBgYHtyfQpjb3ZhcmlhdGVzIDwtIHBhc3RlKCJ4Iiwgc2VxKDEsIGxlbmd0aChiZXRhKSksIHNlcCA9ICIiKQpzdW1tYXJ5KGxtKGFzLmZvcm11bGEocGFzdGUoInkgfiAiLCBwYXN0ZShjb3ZhcmlhdGVzLCBjb2xsYXBzZSA9ICIgKyAiKSkpLCBkYXRhID0gZGF0X3N0cikpCmBgYApJZiB5b3UgY29tcGFyZSB0aG9zZSByZXN1bHRzIHdpdGggdGhlIG9uZXMgd2Ugb2J0YWluZWQgYmVmb3JlLFteMl0geW91J2xsIHNlZSB0aGF0IHRoZXkncmUgaWRlbnRpY2FsIHRvIDYgZGVjaW1hbCBwbGFjZXMuIFRoYXQgbWFrZXMgaXQgcHJldHR5IGxpa2VseSB0aGF0IHRoZXkncmUgdGhlIHNhbWUgZGF0YSBzZXRzLgoKV2UnbGwgZm9jdXMgdGhpcyBleGFtcGxlIG9uIGB4MWAsIGB4M2AsIGFuZCBgeDlgLCB3aGVyZSB3ZSBrbm93IHRoYXQgdGhlICJyZWFsIiBhc3NvY2lhdGlvbnMgYXJlIGJldHdlZW4gYHgxYCwgYHgzYCwgYW5kIGB5YCBhbmQgdGhhdCB0aGUgImFwcGFyZW50IiBhc3NvY2lhdGlvbiBiZXR3ZWVuIGB4OWAgYW5kIGB5YCBhcmlzZXMgb25seSBiZWNhdXNlIGB4OWAgaXMgc3Ryb25nbHkgYXNzb2NpYXRlZCB3aXRoIGB4MWAgYW5kIGB4M2AuIExldCdzIHN0YXJ0IGJ5IGxvb2tpbmcgYXQgdGhlIGJpdmFyaWF0ZSBhc3NvY2lhdGlvbiBiZXR3ZWVuIGB4MWAgYW5kIGB5YCwgaW5jbHVkaW5nIGEgcGxvdCBvZiB0aGUgcmVzaWR1YWxzLgpgYGB7cn0KbG1fc3RyIDwtIGxtKHkgfiB4MSwgZGF0YSA9IGRhdF9zdHIpCnJlc2lkdWFsIDwtIHJlc2lkdWFscyhsbV9zdHIpCmZvci5wbG90IDwtIGRhdGEuZnJhbWUoeCA9IGRhdF9zdHIkeDEsIHkgPSByZXNpZHVhbCkKcCA8LSBnZ3Bsb3QoZm9yLnBsb3QsIGFlcyh4ID0geCwgeSA9IHkpKSArCiAgZ2VvbV9wb2ludChmaWxsID0gInRvbWF0byIsIGNvbG9yID0gInRvbWF0byIpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAwLjApICsKICB4bGFiKCJPYnNlcnZlZCAoeDEpIikgKwogIHlsYWIoIlJlc2lkdWFsIikKcHJpbnQocCkKYGBgCkFzIHlvdSBjYW4gc2VlLCB0aGUgcmVzaWR1YWxzIHNlZW0gcHJldHR5IHJhbmRvbWx5IGRpc3RyaWJ1dGVkLCB3aGljaCBpcyB3aGF0IHdlIGhvcGUgdG8gc2VlIHdoZW4gd2UgZG8gYSByZXNpZHVhbCBwbG90LiBCdXQgbGV0J3MgdHJ5IHRoaXMuIExldCdzIHRyeSByZWdyZXNzaW5nIHRoZSByZXNpZHVhbHMgb24gYHgzYC4gVGhpbmsgb2YgdGhpcyBhcyBsb29raW5nIGF0IHRoZSBpbmZsdWVuY2Ugb2YgYHgzYCBvbiBgeWAgb25jZSB3ZSd2ZSByZW1vdmVkIHRoZSBpbmZsdWVuY2Ugb2YgYHgxYC4KYGBge3J9CnJlc194MyA8LSBsbShyZXNpZHVhbCB+IHgzLCBkYXRhID0gZGF0X3N0cikKcHJpbnQoc3VtbWFyeShyZXNfeDMpKQpgYGAKTG9vayBhdCB0aGF0LiBXZSBoYXZlIGdvb2QgZXZpZGVuY2UgdGhhdCB0aGUgcmVzaWR1YWxzIGZyb20gcmVncmVzc2luZyBgeWAgb24gYHgxYCBzaG93IGEgcmVsYXRpb25zaGlwIHdpdGggYHgzYCAtIHRoZSBsYXJnZXIgYHgzYCwgdGhlIGxhcmdlciB0aGUgcmVzaWR1YWwuIFdlIGNhbiBzZWUgdGhpcyBib3RoIGJ5IGNvbG9yaW5nIHRoZSByZXNpZHVhbHMgYWJvdmUgYW5kIGJ5IHBsb3R0aW5nIHRoZSByZXNpZHVhbHMgYWdhaW5zdCBgeDNgLgpgYGB7cn0KZm9yLnBsb3QgPC0gZGF0YS5mcmFtZSh4ID0gZGF0X3N0ciR4MSwgeSA9IHJlc2lkdWFsLCB4MyA9IGRhdF9zdHIkeDMpCnAgPC0gZ2dwbG90KGZvci5wbG90LCBhZXMoeCA9IHgsIHkgPSB5LCBmaWxsID0geDMsIGNvbG9yID0geDMpKSArCiAgZ2VvbV9wb2ludCgpICsKICBzY2FsZV9maWxsX2dyYWRpZW50MigpICsKICBzY2FsZV9jb2xvcl9ncmFkaWVudDIoKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMC4wKSArCiAgeGxhYigiT2JzZXJ2ZWQgKHgxKSIpICsKICB5bGFiKCJSZXNpZHVhbCIpCnByaW50KHApCnAgPC0gZ2dwbG90KGZvci5wbG90LCBhZXMoeCA9IHgzLCB5ID0gcmVzaWR1YWwpKSArCiAgZ2VvbV9wb2ludChjb2xvciA9ICJ0b21hdG8iLCBmaWxsID0gInRvbWF0byIpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iKSArCiAgeGxhYigiT2JzZXJ2ZWQgKHgzKSIpICsKICB5bGFiKCJSZXNpZHVhbCAoZnJvbSB4MSkiKQpwcmludChwKQpgYGAKTm93IGxldCdzIHNlZSB3aGF0IGhhcHBlbnMgaWYgd2UgZG8gdGhlIHNhbWUgdGhpbmcgd2l0aCBgeDlgLgpgYGB7cn0KcmVzX3g5IDwtIGxtKHJlc2lkdWFsIH4geDksIGRhdGEgPSBkYXRfc3RyKQpwcmludChzdW1tYXJ5KHJlc194OSkpCgpmb3IucGxvdCA8LSBkYXRhLmZyYW1lKHggPSBkYXRfc3RyJHgxLCB5ID0gcmVzaWR1YWwsIHg5ID0gZGF0X3N0ciR4OSkKcCA8LSBnZ3Bsb3QoZm9yLnBsb3QsIGFlcyh4ID0geCwgeSA9IHksIGZpbGwgPSB4OSwgY29sb3IgPSB4OSkpICsKICBnZW9tX3BvaW50KCkgKwogIHNjYWxlX2ZpbGxfZ3JhZGllbnQyKCkgKwogIHNjYWxlX2NvbG9yX2dyYWRpZW50MigpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAwLjApICsKICB4bGFiKCJPYnNlcnZlZCAoeDEpIikgKwogIHlsYWIoIlJlc2lkdWFsIikKcHJpbnQocCkKcCA8LSBnZ3Bsb3QoZm9yLnBsb3QsIGFlcyh4ID0geDksIHkgPSByZXNpZHVhbCkpICsKICBnZW9tX3BvaW50KGNvbG9yID0gInRvbWF0byIsIGZpbGwgPSAidG9tYXRvIikgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpICsKICB4bGFiKCJPYnNlcnZlZCAoeDkpIikgKwogIHlsYWIoIlJlc2lkdWFsIChmcm9tIHgxKSIpCnByaW50KHApCmBgYApIZXJlIHdlIHNlZSB0aGF0IHRoZXJlIGlzIGEgd2VhayBwb3NpdGl2ZSBhc3NvY2lhdGlvbiBiZXR3ZWVuIGB4OWAgYW5kIGB5YCwgYnV0IGl0IGlzbid0IGNvbnZpbmNpbmcuIE9uY2Ugd2UndmUgcmVtb3ZlZCB0aGUgaW5mbHVlbmNlIG9mIGB4MWAgdGhlIHJlbWFpbmluZyBhc3NvY2lhdGlvbiBiZXR3ZWVuIGB4OWAgYW5kIGB5YCBwcmVzdW1hYmx5IGFyaXNlcyBiZWNhdXNlIGB4M2AgYW5kIGB4OWAgYXJlIGFzc29jaWF0ZWQuIExldCdzIHRyeSB0aGlzIGFnYWluIGFmdGVyIHJlbW92aW5nIHRoZSBpbmZsdWVuY2Ugb2YgXypib3RoKl8gYHgxYCBhbmQgYHgzYC4KYGBge3J9CmxtX3gxX3gzIDwtIGxtKHJlc2lkdWFsIH4geDEgK3gzLCBkYXRhID0gZGF0X3N0cikKcHJpbnQoc3VtbWFyeShsbV94MV94MykpCgpyZXNpZHVhbCA8LSByZXNpZHVhbHMobG1feDFfeDMpCgpyZXNfeDkgPC0gbG0ocmVzaWR1YWwgfiB4OSwgZGF0YSA9IGRhdF9zdHIpCnByaW50KHN1bW1hcnkocmVzX3g5KSkKCmZvci5wbG90IDwtIGRhdGEuZnJhbWUoeCA9IGRhdF9zdHIkeDEsIHkgPSByZXNpZHVhbCwgeDkgPSBkYXRfc3RyJHg5KQpwIDwtIGdncGxvdChmb3IucGxvdCwgYWVzKHggPSB4OSwgeSA9IHJlc2lkdWFsKSkgKwogIGdlb21fcG9pbnQoY29sb3IgPSAidG9tYXRvIiwgZmlsbCA9ICJ0b21hdG8iKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIikgKwogIHhsYWIoIk9ic2VydmVkICh4OSkiKSArCiAgeWxhYigiUmVzaWR1YWwgKGZyb20geDEpIikKcHJpbnQocCkKYGBgCk5vdGljZSB0aGF0IHRoZSBjb2VmZmljaWVudCBvbiBgeDlgIGlzIG5vdyBvbmx5IGFib3V0IDAuMDkgcmF0aGVyIHRoYW4gMC4xNi4gVGhlcmUncyBzdGlsbCBhIGJpdCBvZiBhbiBhc3NvY2lhdGlvbiwgYnV0IGl0J3Mgd2VhayBhbmQgcG9vcmx5IHN1cHBvcnRlZC4gCgpTbyB0aGF0J3Mgd2hhdCBgbG0oKWAgaXMgZG9pbmcuIEl0IGlzICJzdGF0aXN0aWNhbGx5IGNvbnRyb2xsaW5nIiBmb3IgdGhlIGVmZmVjdHMgb2Ygb3RoZXIgY292YXJpYXRlcyB3aGVuIGVzdGltYXRpbmcgdGhlIGVmZmVjdCBvZiBlYWNoIG9uZSBpbmRpdmlkdWFsbHkuIEl0IGRvZXNuJ3QgZml0IHRoZSByZWdyZXNzaW9uIHNlcXVlbnRpYWxseSwgYXMgSSd2ZSBkb25lIGl0IGhlcmUuIEl0IGZpdHMgYWxsIG9mIHRoZSBjb2VmZmljaWVudHMgc2ltdWx0YW5lb3VzbHkuIEFzIGEgcmVzdWx0LCBpdCBkb2Vzbid0IG1ha2UgYW55IGRpZmZlcmVuY2Ugd2hhdCBvcmRlciB5b3UgcHV0IHZhcmlhYmxlcyBpbnRvIHRoZSBtb2RlbCBzdGF0ZW1lbnQuIFlvdSdsbCBnZXQgdGhlIHNhbWUgcmVzdWx0LlteM10KClteMV06IFRvIG1ha2Ugc3VyZSB0aGF0IHRoZSBkYXRhIGFyZSB0aGUgc2FtZSwgd2UgaGF2ZSB0byBnZW5lcmF0ZSBhbGwgdGhyZWUgZGF0YSBzZXRzLgpbXjJdOiBTY3JvbGwgYWxsIHRoZSB3YXkgdG8gdGhlIGJvdHRvbSBvZiB0aGUgcGFnZSwganVzdCBiZWZvcmUgdGhlIENvbmNsdXNpb24gYXQgW2h0dHA6Ly9kYXJ3aW4uZWViLnVjb25uLmVkdS9wYWdlcy92YXJpYWJsZS1zZWxlY3Rpb24vbXVsdGlwbGUtcmVncmVzc2lvbi1iYXNpY3MubmIuaHRtbF0oaHR0cDovL2Rhcndpbi5lZWIudWNvbm4uZWR1L3BhZ2VzL3ZhcmlhYmxlLXNlbGVjdGlvbi9tdWx0aXBsZS1yZWdyZXNzaW9uLWJhc2ljcy5uYi5odG1sKS4KW14zXTogSWYgeW91IGRvbid0IGJlbGlldmUgbWUsIHRyeSBtYW51YWxseSBjaGFuZ2luZyB0aGUgb3JkZXIgb2YgdGhlIHZhcmlhYmxlcyBpbiB0aGlzIG11bHRpcGxlIHJlZ3Jlc3Npb24gYW5kIHNlZSB3aGF0IHlvdSBnZXQu