Introduction
Before we begin exploring some ideas about variable selection in multiple regression, we need to review
- Why multiple regression may be needed and
- Some of the problems that arise when covariates aren’t independent of one another, i.e., when they are collinear.
It will take us a while to get through the first of these points, so that’s all we’ll cover in this notebook. The next notebook will take up some of the problems that arise when covariates aren’t independent of one another.
To do that we’ll create three toy data sets in which there are 9 covariates, only 3 of which are associated with the response variable. The three data sets will have the same regression coefficients, i.e., the same intercept, the same non-zero value for the 3 covariates that have an association, and 0 for the remaining 6 covariates. They’ll differ in the degree of collinearity among the covariates, ranging from complete independence of all covariates to relativelly modest collinearity to relatively strong collinearity.
To get started, we’ll load all of the libraries we’re going to use in this session and set up sample data sets. Note: I’m going to be generating the covariate data from a multivariate normal, but that’s only for convenience. Linear regression has some important assumptions. (Multi)normality of covariates is at the bottom of the list.
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)
}
Generating the data
OK, with those preliminaries out of the way. Let’s generate data under our three scenarios:
- no association:
ind (rho = 0)
- weak association:
wk (rho = 0.2)
- strong association:
str (rho = 0.8)
We’ll use the same residual standard deviation, resid = 0.2, for all of them. cov_ refers to the covariates. y_ is the corresponding response
## 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)
Now I’ll collect the data into a data frame so that we can look at the covariates and explore how they relate to the response variable.
dat_ind <- data.frame(y_ind, cov_ind, rep("Independent", length(y_ind)))
dat_wk <- data.frame(y_wk, cov_wk, rep("Weak", length(y_wk)))
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_ind) <- column_names
colnames(dat_wk) <- column_names
colnames(dat_str) <- column_names
dat <- rbind(dat_ind, dat_wk, dat_str)
Covariate associations
First, let’s just look at the association of the covariates under the three scenarios.
p <- PairPlot(dat_ind,
meas_vars = column_names[grep("x", column_names)],
title = "Independent")
print(p)
p <- PairPlot(dat_wk,
meas_vars = column_names[grep("x", column_names)],
title = "Weak association")
print(p)
p <- PairPlot(dat_str,
meas_vars = column_names[grep("x", column_names)],
title = "Strong association")
print(p)
When you compare the plot for independent samples with the one for weak associations, you may not see much of a difference, but corrplot shows that the differences are there.
corrplot(cor(dat_ind[,grepl("x", column_names)]))
corrplot(cor(dat_wk[,grepl("x", column_names)]))
corrplot(cor(dat_str[,grepl("x", column_names)]))
Exploring bivariate associations between covariates and the response
Now that we’ve seen how the covariats are related, let’s examine pairwise associations between each of the covariates and the response variable. We first have to change dat from a “wide” format to a “long” format. Then we can use ggplot to produce a color-coded scatterplot with a regression line in which the colors correspond to the scenarios.
for.plot <- melt(dat,
id.vars = c("Scenario", "y"))
p <- ggplot(for.plot, aes(x = value, y = y, color = Scenario, fill = Scenario), group = Scenario) +
geom_point(alpha = 0.3) +
stat_smooth(method = lm) +
xlab("x") +
facet_wrap(~ variable)
print(p)
Remember that when we simulated these data, we set the coefficient on x1 to 1, the coefficient on x2 to 2, and the coefficient on x3 to 1, so the top row of figures looks OK. With either independent covariates or weak associations, there doesn’t seem to be a pairwise association between x4-x9 and y, but what’s going on with strong associations and x4-x9? Look back at the corrplot for strong associations. You’ll see that x4, x6, and x8 are strongly associated with x2 and that x5, x7, and x9 are strongly associated with x1 and x3. The association between x5 and y, for example, arises because x5 is strongly associated with x1 and x3, not because there’s a direct association between x5 and y.
We can check our visual impression that there are “spurious” associations in the strong association scenario by running a simple linear regression for each of the pairwise relationships.
covariate <- character(length(beta))
est <- numeric(length(beta))
lo <- numeric(length(beta))
hi <- numeric(length(beta))
for (i in 1:length(beta)) {
variable_name <- paste("x", i, sep="")
tmp_lm <- lm(y ~ value, data = subset(for.plot, Scenario == "Strong" & variable == variable_name))
covariate[i] <- variable_name
est[i] <- coef(tmp_lm)[2]
tmp_interval <- confint(tmp_lm)
lo[i] <- tmp_interval[2,1]
hi[i] <- tmp_interval[2,2]
}
results <- tibble(covariate = covariate,
estimate = round(est, 3),
lo = round(lo, 3),
hi = round(hi, 3))
print(results)
As you can see, the confidence interval does not include 0 for any of the covariates.
Multiple regression
Let’s compare that to what we get from a multiple regression of y on all of the covariates.
covariates <- paste("x", seq(1, length(beta)), sep = "")
summary(lm(as.formula(paste("y ~ ", paste(covariates, collapse = " + "))), data = subset(dat, Scenario == "Strong")))
Call:
lm(formula = as.formula(paste("y ~ ", paste(covariates, collapse = " + "))),
data = subset(dat, Scenario == "Strong"))
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-16
Here you can see that the data provide very strong evidence for an association between x1, x2, x3 and y. The magnitudes of the estimated coefficients for x1-x3 are close to 1, and the corresponding standard errors are no more than 0.044. In contrast, there is very little evidence for an association between any of the remaining covariates and y. For all of them the magnitude of the estimated association is smaller than its standard error.
Conclusion
We can distinguish between two types of associations, direct associations and indirect association. An indirect association arises when a covariate, for example x9, is associated with a response only because it is associated with one or more other covariates that do have a direct association with the response, in this example x1 and x3. If we look only at x9 and y, we’ll find an association. When we include x1-x9 in a multiple regression, x9 now has only a negligible association. A multivariate regression may allow us to distinguish “real” from “spurious” associations.
Points left as an exercise
I’ve illustrated the phenomena here with only one simulated data set under each scenario. A rigorous analysis would require a lot of fairly complicated math, which you’re welcome to look up in a book on multiple regression. If you don’t want to do that, you might want to set up a small simulation study in which you produce a hundred or a thousand data sets under conditions similar to those above to make sure the pattern seen in this one example is reproducible across a large number of examples. And as long as you’re doing that, you might also want to generate data sets of different sizes to explore how these results depend on sample size relative to the number of predictors.
LS0tCnRpdGxlOiAiV2h5IG11bHRpcGxlIHJlZ3Jlc3Npb24gaXMgbmVlZGVkIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCiMjIEludHJvZHVjdGlvbgoKQmVmb3JlIHdlIGJlZ2luIGV4cGxvcmluZyBzb21lIGlkZWFzIGFib3V0IHZhcmlhYmxlIHNlbGVjdGlvbiBpbiBtdWx0aXBsZSByZWdyZXNzaW9uLCB3ZSBuZWVkIHRvIHJldmlldyAKCiogV2h5IG11bHRpcGxlIHJlZ3Jlc3Npb24gbWF5IGJlIG5lZWRlZCBhbmQgCiogU29tZSBvZiB0aGUgcHJvYmxlbXMgdGhhdCBhcmlzZSB3aGVuIGNvdmFyaWF0ZXMgYXJlbid0IGluZGVwZW5kZW50IG9mIG9uZSBhbm90aGVyLCBpLmUuLCB3aGVuIHRoZXkgYXJlIFtjb2xsaW5lYXJdKGh0dHBzOi8vZW4ud2lraXBlZGlhLm9yZy93aWtpL011bHRpY29sbGluZWFyaXR5KS4gCgpJdCB3aWxsIHRha2UgdXMgYSB3aGlsZSB0byBnZXQgdGhyb3VnaCB0aGUgZmlyc3Qgb2YgdGhlc2UgcG9pbnRzLCBzbyB0aGF0J3MgYWxsIHdlJ2xsIGNvdmVyIGluIHRoaXMgbm90ZWJvb2suIFRoZSBuZXh0IG5vdGVib29rIHdpbGwgdGFrZSB1cCBzb21lIG9mIHRoZSBwcm9ibGVtcyB0aGF0IGFyaXNlIHdoZW4gY292YXJpYXRlcyBhcmVuJ3QgaW5kZXBlbmRlbnQgb2Ygb25lIGFub3RoZXIuCgpUbyBkbyB0aGF0IHdlJ2xsIGNyZWF0ZSB0aHJlZSB0b3kgZGF0YSBzZXRzIGluIHdoaWNoIHRoZXJlIGFyZSA5IGNvdmFyaWF0ZXMsIG9ubHkgMyBvZiB3aGljaCBhcmUgYXNzb2NpYXRlZCB3aXRoIHRoZSByZXNwb25zZSB2YXJpYWJsZS4gVGhlIHRocmVlIGRhdGEgc2V0cyB3aWxsIGhhdmUgdGhlIHNhbWUgcmVncmVzc2lvbiBjb2VmZmljaWVudHMsIGkuZS4sIHRoZSBzYW1lIGludGVyY2VwdCwgdGhlIHNhbWUgbm9uLXplcm8gdmFsdWUgZm9yIHRoZSAzIGNvdmFyaWF0ZXMgdGhhdCBoYXZlIGFuIGFzc29jaWF0aW9uLCBhbmQgMCBmb3IgdGhlIHJlbWFpbmluZyA2IGNvdmFyaWF0ZXMuIFRoZXknbGwgZGlmZmVyIGluIHRoZSBkZWdyZWUgb2YgY29sbGluZWFyaXR5IGFtb25nIHRoZSBjb3ZhcmlhdGVzLCByYW5naW5nIGZyb20gY29tcGxldGUgaW5kZXBlbmRlbmNlIG9mIGFsbCBjb3ZhcmlhdGVzIHRvIHJlbGF0aXZlbGx5IG1vZGVzdCBjb2xsaW5lYXJpdHkgdG8gcmVsYXRpdmVseSBzdHJvbmcgY29sbGluZWFyaXR5LgoKVG8gZ2V0IHN0YXJ0ZWQsIHdlJ2xsIGxvYWQgYWxsIG9mIHRoZSBsaWJyYXJpZXMgd2UncmUgZ29pbmcgdG8gdXNlIGluIHRoaXMgc2Vzc2lvbiBhbmQgc2V0IHVwIHNhbXBsZSBkYXRhIHNldHMuIE5vdGU6IEknbSBnb2luZyB0byBiZSBnZW5lcmF0aW5nIHRoZSBjb3ZhcmlhdGUgZGF0YSBmcm9tIGEgbXVsdGl2YXJpYXRlIG5vcm1hbCwgYnV0IHRoYXQncyBvbmx5IGZvciBjb252ZW5pZW5jZS4gTGluZWFyIHJlZ3Jlc3Npb24gaGFzIHNvbWUgaW1wb3J0YW50IGFzc3VtcHRpb25zLiAoTXVsdGkpbm9ybWFsaXR5IG9mIGNvdmFyaWF0ZXMgaXMgYXQgdGhlIFtib3R0b20gb2YgdGhlIGxpc3RdKGh0dHBzOi8vc3RhdG1vZGVsaW5nLnN0YXQuY29sdW1iaWEuZWR1LzIwMTMvMDgvMDQvMTk0NzAvKS4KYGBge3Igc2V0dXAsIHdhcm5pbmcgPSBGQUxTRSwgbWVzc2FnZSA9IEZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShyZXNoYXBlMikKbGlicmFyeShjb3JycGxvdCkKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGNvd3Bsb3QpCmxpYnJhcnkoV1ZQbG90cykKbGlicmFyeShtdnRub3JtKQpgYGAKYGBge3J9CiMjIEkgYWx3YXlzIGxpa2UgdG8gY2xlYW4gb3V0IG15IHdvcmtzcGFjZSBiZWZvcmUgcnVubmluZyBhIHNjcmlwdCB0byBtYWtlIHN1cmUKIyMgdGhhdCBJJ20gc3RhcnRpbmcgUiBpbiB0aGUgc2FtZSBzdGF0ZS4gVGhpcyBoZWxwcyB0byBlbnN1cmUgdGhhdCBJIGNhbiAKIyMgcmVwcm9kdWNlIG15IHJlc3VsdHMgbGF0ZXIKcm0obGlzdCA9IGxzKCkpCgojIyBpbnRldGNlcHQKIyMKYmV0YTAgPC0gMS4wCiMjIHJlZ3Jlc3Npb24gY29lZmZpY2llbnRzCiMjCmJldGEgPC0gYygxLjAsIC0xLjAsIDEuMCwgMC4wLCAwLjAsIDAuMCwgMC4wLCAwLjAsIDAuMCkKIyMgcGF0dGVybiBvZiBjb3JyZWxhdGlvbiBtYXRyaXgsIGFsbCBub24temVybyBlbnRyaWVzIGFyZSBzZXQgdG8gc2FlbQojIyBjb3JyZWxhdGlvbiwgY292YXJpYW5jZSBtYXRyaXggY2FsZHVsYXRlZCBmcm9tIGluZGl2aWR1YWwgdmFyaWFuY2VzIGFuZCBhIAojIyBzaW5nbGUgYXNzb2NpYXRpb24gcGFyYW1ldGVyIGdvdmVybmluZyB0aGUgbm9uLXplcm8gY29ycmVsYXRpb24gY29lZmZpY2llbnRzCiMjCiMjIE5vdGU6IE5vdCBqdXN0IGFueSBwYXR0ZXJuIHdpbGwgd29yayBoZXJlLiBUaGUgY29ycmVsYXRpb24gbWF0cml4IGFuZAojIyBjb3ZhcmlhbmNlIG1hdHJpeCBnZW5lcmF0ZWQgZnJvbSB0aGlzIHBhdHRlcm4gbXVzdCBiZSBwb3NpdGl2ZSBkZWZpbml0ZS4KIyMgSWYgeW91IGNoYW5nZSB0aGlzIHBhdHRlcm4sIHlvdSBtYXkgZ2V0IGFuIGVycm9yIHdoZW4geW91IHRyeSB0byBnZW5lcmF0ZQojIyBkYXRhIHdpdGggYSBub24temVybyBhc3NvY2lhdGlvbiBwYXJhbWV0ZXIuCiMjClJobyA8LSBtYXRyaXgobnJvdyA9IDksIG5jb2wgPSAsIGJ5cm93ID0gVFJVRSwgCiAgICAgICAgICAgICAgZGF0YSA9IGMoMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEKICAgICAgICAgICAgICAgICAgICAgICApKQojIyB2ZWN0b3Igb2Ygc3RhbmRhcmQgZGV2aWF0aW9ucyBmb3IgY292YXJpYXRlcwojIwpzaWdtYSA8LSByZXAoMSwgOSkKCiMjIGNvbnN0cnVjdCBhIGNvdmFyaWFuY2UgbWF0cml4IGZyb20gdGhlIHBhdHRlcm4sIHN0YW5kYXJkIGRldmlhdGlvbnMsIGFuZAojIyBvbmUgcGFyYW1ldGVyIGluIFstMSwxXSB0aGF0IGdvdmVybnMgdGhlIG1hZ25pdHVkZSBvZiBub24temVybyBjb3JyZWxhdGlvbgojIyBjb2VmZmljaWVudHMKIyMKIyMgUmhvIC0gdGhlIHBhdHRlcm4gb2YgYXNzb2NpYXRpb25zCiMjIHNpZ21hIC0gdGhlIHZlY3RvciBvZiBzdGFuZGFyZCBkZXZpYXRpb25zCiMjIHJobyAtIHRoZSBhc3NvY2lhdGlvbiBwYXJhbWV0ZXIKIyMKY29uc3RydWN0X1NpZ21hIDwtIGZ1bmN0aW9uKFJobywgc2lnbWEsIHJobykgewogICMjIGdldCB0aGUgY29ycmVsYXRpb24gbWF0cmlzCiAgIyMKICBSaG8gPC0gUmhvKnJobwogIGZvciAoaSBpbiAxOm5jb2woUmhvKSkgewogICAgUmhvW2ksaV0gPC0gMS4wCiAgfQogICMjIG5vdGljZSB0aGUgdXNlIG9mIG1hdHJpeCBtdWx0aXBsaWNhdGlvbgogICMjCiAgU2lnbWEgPC0gZGlhZyhzaWdtYSkgJSolIFJobyAlKiUgZGlhZyhzaWdtYSkKICByZXR1cm4oU2lnbWEpCn0KYGBgCgojIyBHZW5lcmF0aW5nIHRoZSBkYXRhCgpPSywgd2l0aCB0aG9zZSBwcmVsaW1pbmFyaWVzIG91dCBvZiB0aGUgd2F5LiBMZXQncyBnZW5lcmF0ZSBkYXRhIHVuZGVyIG91ciB0aHJlZSBzY2VuYXJpb3M6CgoqIG5vIGFzc29jaWF0aW9uOiBgaW5kYCAoYHJobyA9IDBgKQoqIHdlYWsgYXNzb2NpYXRpb246IGB3a2AgKGByaG8gPSAwLjJgKQoqIHN0cm9uZyBhc3NvY2lhdGlvbjogYHN0cmAgKGByaG8gPSAwLjhgKQoKV2UnbGwgdXNlIHRoZSBzYW1lIHJlc2lkdWFsIHN0YW5kYXJkIGRldmlhdGlvbiwgYHJlc2lkID0gMC4yYCwgZm9yIGFsbCBvZiB0aGVtLiBgY292X2AgcmVmZXJzIHRvIHRoZSBjb3ZhcmlhdGVzLiBgeV9gIGlzIHRoZSBjb3JyZXNwb25kaW5nIHJlc3BvbnNlCmBgYHtyfQojIyBzZXQgdGhlIHJhbmRvbSBudW1iZXIgc2VlZCBtYW51YWxseSBzbyB0aGF0IGV2ZXJ5IHJ1biBvZiB0aGUgY29kZSB3aWxsIAojIyBwcm9kdWNlIHRoZSBzYW1lIG51bWJlcnMKIyMKc2V0LnNlZWQoMTIzNCkKCm5fc2FtcCA8LSAxMDAKY292X2luZCA8LSBybXZub3JtKG5fc2FtcCwKICAgICAgICAgICAgICAgICAgIG1lYW4gPSByZXAoMCwgbnJvdyhSaG8pKSwKICAgICAgICAgICAgICAgICAgIHNpZ21hID0gY29uc3RydWN0X1NpZ21hKFJobywgc2lnbWEsIDAuMCkpCmNvdl93ayA8LSBybXZub3JtKG5fc2FtcCwKICAgICAgICAgICAgICAgICAgbWVhbiA9IHJlcCgwLCBucm93KFJobykpLAogICAgICAgICAgICAgICAgICBzaWdtYSA9IGNvbnN0cnVjdF9TaWdtYShSaG8sIHNpZ21hLCAwLjIpKQpjb3Zfc3RyIDwtIHJtdm5vcm0obl9zYW1wLAogICAgICAgICAgICAgICAgICAgbWVhbiA9IHJlcCgwLCBucm93KFJobykpLAogICAgICAgICAgICAgICAgICAgc2lnbWEgPSBjb25zdHJ1Y3RfU2lnbWEoUmhvLCBzaWdtYSwgMC44KSkKCnJlc2lkIDwtIHJlcCgwLjIsIG5fc2FtcCkKeV9pbmQgPC0gcm5vcm0obnJvdyhjb3ZfaW5kKSwgbWVhbiA9IGJldGEwICsgY292X2luZCAlKiUgYmV0YSwgc2QgPSByZXNpZCkKeV93ayA8LSBybm9ybShucm93KGNvdl93ayksIG1lYW4gPSBiZXRhMCArIGNvdl93ayAlKiUgYmV0YSwgc2QgPSByZXNpZCkKeV9zdHIgPC0gcm5vcm0obnJvdyhjb3Zfc3RyKSwgbWVhbiA9IGJldGEwICsgY292X3N0ciAlKiUgYmV0YSwgc2QgPSByZXNpZCkKYGBgCk5vdyBJJ2xsIGNvbGxlY3QgdGhlIGRhdGEgaW50byBhIGBkYXRhIGZyYW1lYCBzbyB0aGF0IHdlIGNhbiBsb29rIGF0IHRoZSBjb3ZhcmlhdGVzIGFuZCBleHBsb3JlIGhvdyB0aGV5IHJlbGF0ZSB0byB0aGUgcmVzcG9uc2UgdmFyaWFibGUuIApgYGB7cn0KZGF0X2luZCA8LSBkYXRhLmZyYW1lKHlfaW5kLCBjb3ZfaW5kLCByZXAoIkluZGVwZW5kZW50IiwgbGVuZ3RoKHlfaW5kKSkpCmRhdF93ayA8LSBkYXRhLmZyYW1lKHlfd2ssIGNvdl93aywgcmVwKCJXZWFrIiwgbGVuZ3RoKHlfd2spKSkKZGF0X3N0ciA8LSBkYXRhLmZyYW1lKHlfc3RyLCBjb3Zfc3RyLCByZXAoIlN0cm9uZyIsIGxlbmd0aCh5X3N0cikpKQoKY29sdW1uX25hbWVzIDwtIGMoInkiLCBwYXN0ZSgieCIsIHNlcSgxLCBsZW5ndGgoYmV0YSkpLCBzZXAgPSAiIiksICJTY2VuYXJpbyIpCmNvbG5hbWVzKGRhdF9pbmQpIDwtIGNvbHVtbl9uYW1lcwpjb2xuYW1lcyhkYXRfd2spIDwtIGNvbHVtbl9uYW1lcwpjb2xuYW1lcyhkYXRfc3RyKSA8LSBjb2x1bW5fbmFtZXMKCmRhdCA8LSByYmluZChkYXRfaW5kLCBkYXRfd2ssIGRhdF9zdHIpCmBgYAoKIyMgQ292YXJpYXRlIGFzc29jaWF0aW9ucwoKRmlyc3QsIGxldCdzIGp1c3QgbG9vayBhdCB0aGUgYXNzb2NpYXRpb24gb2YgdGhlIGNvdmFyaWF0ZXMgdW5kZXIgdGhlIHRocmVlIHNjZW5hcmlvcy4KYGBge3J9CnAgPC0gUGFpclBsb3QoZGF0X2luZCwKICAgICAgICAgICAgICBtZWFzX3ZhcnMgPSBjb2x1bW5fbmFtZXNbZ3JlcCgieCIsIGNvbHVtbl9uYW1lcyldLAogICAgICAgICAgICAgIHRpdGxlID0gIkluZGVwZW5kZW50IikKcHJpbnQocCkKcCA8LSBQYWlyUGxvdChkYXRfd2ssCiAgICAgICAgICAgICAgbWVhc192YXJzID0gY29sdW1uX25hbWVzW2dyZXAoIngiLCBjb2x1bW5fbmFtZXMpXSwKICAgICAgICAgICAgICB0aXRsZSA9ICJXZWFrIGFzc29jaWF0aW9uIikKcHJpbnQocCkKcCA8LSBQYWlyUGxvdChkYXRfc3RyLAogICAgICAgICAgICAgIG1lYXNfdmFycyA9IGNvbHVtbl9uYW1lc1tncmVwKCJ4IiwgY29sdW1uX25hbWVzKV0sCiAgICAgICAgICAgICAgdGl0bGUgPSAiU3Ryb25nIGFzc29jaWF0aW9uIikKcHJpbnQocCkKYGBgCldoZW4geW91IGNvbXBhcmUgdGhlIHBsb3QgZm9yIGluZGVwZW5kZW50IHNhbXBsZXMgd2l0aCB0aGUgb25lIGZvciB3ZWFrIGFzc29jaWF0aW9ucywgeW91IG1heSBub3Qgc2VlIG11Y2ggb2YgYSBkaWZmZXJlbmNlLCBidXQgYGNvcnJwbG90YCBzaG93cyB0aGF0IHRoZSBkaWZmZXJlbmNlcyBhcmUgdGhlcmUuCmBgYHtyfQpjb3JycGxvdChjb3IoZGF0X2luZFssZ3JlcGwoIngiLCBjb2x1bW5fbmFtZXMpXSkpCmNvcnJwbG90KGNvcihkYXRfd2tbLGdyZXBsKCJ4IiwgY29sdW1uX25hbWVzKV0pKQpjb3JycGxvdChjb3IoZGF0X3N0clssZ3JlcGwoIngiLCBjb2x1bW5fbmFtZXMpXSkpCmBgYAojIyBFeHBsb3JpbmcgYml2YXJpYXRlIGFzc29jaWF0aW9ucyBiZXR3ZWVuIGNvdmFyaWF0ZXMgYW5kIHRoZSByZXNwb25zZQoKTm93IHRoYXQgd2UndmUgc2VlbiBob3cgdGhlIGNvdmFyaWF0cyBhcmUgcmVsYXRlZCwgbGV0J3MgZXhhbWluZSBwYWlyd2lzZSBhc3NvY2lhdGlvbnMgYmV0d2VlbiBlYWNoIG9mIHRoZSBjb3ZhcmlhdGVzIGFuZCB0aGUgcmVzcG9uc2UgdmFyaWFibGUuIFdlIGZpcnN0IGhhdmUgdG8gY2hhbmdlIGBkYXRgIGZyb20gYSAid2lkZSIgZm9ybWF0IHRvIGEgImxvbmciIGZvcm1hdC4gW14xXSBUaGVuIHdlIGNhbiB1c2UgYGdncGxvdGAgdG8gcHJvZHVjZSBhIGNvbG9yLWNvZGVkIHNjYXR0ZXJwbG90IHdpdGggYSByZWdyZXNzaW9uIGxpbmUgaW4gd2hpY2ggdGhlIGNvbG9ycyBjb3JyZXNwb25kIHRvIHRoZSBzY2VuYXJpb3MuCmBgYHtyfQpmb3IucGxvdCA8LSBtZWx0KGRhdCwKICAgICAgICAgICAgICAgICBpZC52YXJzID0gYygiU2NlbmFyaW8iLCAieSIpKQoKcCA8LSBnZ3Bsb3QoZm9yLnBsb3QsIGFlcyh4ID0gdmFsdWUsIHkgPSB5LCBjb2xvciA9IFNjZW5hcmlvLCBmaWxsID0gU2NlbmFyaW8pLCBncm91cCA9IFNjZW5hcmlvKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuMykgKwogIHN0YXRfc21vb3RoKG1ldGhvZCA9IGxtKSArCiAgeGxhYigieCIpICsKICBmYWNldF93cmFwKH4gdmFyaWFibGUpCnByaW50KHApCmBgYApSZW1lbWJlciB0aGF0IHdoZW4gd2Ugc2ltdWxhdGVkIHRoZXNlIGRhdGEsIHdlIHNldCB0aGUgY29lZmZpY2llbnQgb24gYHgxYCB0byAxLCB0aGUgY29lZmZpY2llbnQgb24gYHgyYCB0byAyLCBhbmQgdGhlIGNvZWZmaWNpZW50IG9uIGB4M2AgdG8gMSwgc28gdGhlIHRvcCByb3cgb2YgZmlndXJlcyBsb29rcyBPSy4gV2l0aCBlaXRoZXIgaW5kZXBlbmRlbnQgY292YXJpYXRlcyBvciB3ZWFrIGFzc29jaWF0aW9ucywgdGhlcmUgZG9lc24ndCBzZWVtIHRvIGJlIGEgcGFpcndpc2UgYXNzb2NpYXRpb24gYmV0d2VlbiBgeDQteDlgIGFuZCBgeWAsIGJ1dCB3aGF0J3MgZ29pbmcgb24gd2l0aCBzdHJvbmcgYXNzb2NpYXRpb25zIGFuZCBgeDQteDlgPyBMb29rIGJhY2sgYXQgdGhlIGBjb3JycGxvdGAgZm9yIHN0cm9uZyBhc3NvY2lhdGlvbnMuIFlvdSdsbCBzZWUgdGhhdCBgeDRgLCBgeDZgLCBhbmQgYHg4YCBhcmUgc3Ryb25nbHkgYXNzb2NpYXRlZCB3aXRoIGB4MmAgYW5kIHRoYXQgYHg1YCwgYHg3YCwgYW5kIGB4OWAgYXJlIHN0cm9uZ2x5IGFzc29jaWF0ZWQgd2l0aCBgeDFgIGFuZCBgeDNgLiBUaGUgYXNzb2NpYXRpb24gYmV0d2VlbiBgeDVgIGFuZCBgeWAsIGZvciBleGFtcGxlLCBhcmlzZXMgYmVjYXVzZSBgeDVgIGlzIHN0cm9uZ2x5IGFzc29jaWF0ZWQgd2l0aCBgeDFgIGFuZCBgeDNgLCBub3QgYmVjYXVzZSB0aGVyZSdzIGEgZGlyZWN0IGFzc29jaWF0aW9uIGJldHdlZW4gYHg1YCBhbmQgYHlgLiAKCldlIGNhbiBjaGVjayBvdXIgdmlzdWFsIGltcHJlc3Npb24gdGhhdCB0aGVyZSBhcmUgInNwdXJpb3VzIiBhc3NvY2lhdGlvbnMgaW4gdGhlIHN0cm9uZyBhc3NvY2lhdGlvbiBzY2VuYXJpbyBieSBydW5uaW5nIGEgc2ltcGxlIGxpbmVhciByZWdyZXNzaW9uIGZvciBlYWNoIG9mIHRoZSBwYWlyd2lzZSByZWxhdGlvbnNoaXBzLiBbXjJdIApgYGB7cn0KY292YXJpYXRlIDwtIGNoYXJhY3RlcihsZW5ndGgoYmV0YSkpCmVzdCA8LSBudW1lcmljKGxlbmd0aChiZXRhKSkKbG8gPC0gbnVtZXJpYyhsZW5ndGgoYmV0YSkpCmhpIDwtIG51bWVyaWMobGVuZ3RoKGJldGEpKQpmb3IgKGkgaW4gMTpsZW5ndGgoYmV0YSkpIHsKICB2YXJpYWJsZV9uYW1lIDwtIHBhc3RlKCJ4IiwgaSwgc2VwPSIiKQogIHRtcF9sbSA8LSBsbSh5IH4gdmFsdWUsIGRhdGEgPSBzdWJzZXQoZm9yLnBsb3QsIFNjZW5hcmlvID09ICJTdHJvbmciICYgdmFyaWFibGUgPT0gdmFyaWFibGVfbmFtZSkpCiAgY292YXJpYXRlW2ldIDwtIHZhcmlhYmxlX25hbWUKICBlc3RbaV0gPC0gY29lZih0bXBfbG0pWzJdCiAgdG1wX2ludGVydmFsIDwtIGNvbmZpbnQodG1wX2xtKQogIGxvW2ldIDwtIHRtcF9pbnRlcnZhbFsyLDFdCiAgaGlbaV0gPC0gdG1wX2ludGVydmFsWzIsMl0KfQpyZXN1bHRzIDwtIHRpYmJsZShjb3ZhcmlhdGUgPSBjb3ZhcmlhdGUsIAogICAgICAgICAgICAgICAgICBlc3RpbWF0ZSA9IHJvdW5kKGVzdCwgMyksCiAgICAgICAgICAgICAgICAgIGxvID0gcm91bmQobG8sIDMpLAogICAgICAgICAgICAgICAgICBoaSA9IHJvdW5kKGhpLCAzKSkgIAoKcHJpbnQocmVzdWx0cykKYGBgCkFzIHlvdSBjYW4gc2VlLCB0aGUgY29uZmlkZW5jZSBpbnRlcnZhbCBkb2VzIG5vdCBpbmNsdWRlIDAgZm9yIGFueSBvZiB0aGUgY292YXJpYXRlcy4gCgojIyBNdWx0aXBsZSByZWdyZXNzaW9uCgpMZXQncyBjb21wYXJlIHRoYXQgdG8gd2hhdCB3ZSBnZXQgZnJvbSBhIG11bHRpcGxlIHJlZ3Jlc3Npb24gb2YgeSBvbiBhbGwgb2YgdGhlIGNvdmFyaWF0ZXMuCmBgYHtyfQpjb3ZhcmlhdGVzIDwtIHBhc3RlKCJ4Iiwgc2VxKDEsIGxlbmd0aChiZXRhKSksIHNlcCA9ICIiKQpzdW1tYXJ5KGxtKGFzLmZvcm11bGEocGFzdGUoInkgfiAiLCBwYXN0ZShjb3ZhcmlhdGVzLCBjb2xsYXBzZSA9ICIgKyAiKSkpLCBkYXRhID0gc3Vic2V0KGRhdCwgU2NlbmFyaW8gPT0gIlN0cm9uZyIpKSkKYGBgCkhlcmUgeW91IGNhbiBzZWUgdGhhdCB0aGUgZGF0YSBwcm92aWRlIHZlcnkgc3Ryb25nIGV2aWRlbmNlIGZvciBhbiBhc3NvY2lhdGlvbiBiZXR3ZWVuIGB4MWAsIGB4MmAsIGB4M2AgYW5kIGB5YC4gVGhlIG1hZ25pdHVkZXMgb2YgdGhlIGVzdGltYXRlZCBjb2VmZmljaWVudHMgZm9yIGB4MWAtYHgzYCBhcmUgY2xvc2UgdG8gMSwgYW5kIHRoZSBjb3JyZXNwb25kaW5nIHN0YW5kYXJkIGVycm9ycyBhcmUgbm8gbW9yZSB0aGFuIDAuMDQ0LiBJbiBjb250cmFzdCwgdGhlcmUgaXMgdmVyeSBsaXR0bGUgZXZpZGVuY2UgZm9yIGFuIGFzc29jaWF0aW9uIGJldHdlZW4gYW55IG9mIHRoZSByZW1haW5pbmcgY292YXJpYXRlcyBhbmQgYHlgLiBGb3IgYWxsIG9mIHRoZW0gdGhlIG1hZ25pdHVkZSBvZiB0aGUgZXN0aW1hdGVkIGFzc29jaWF0aW9uIGlzIHNtYWxsZXIgdGhhbiBpdHMgc3RhbmRhcmQgZXJyb3IuCgojIyBDb25jbHVzaW9uCgpXZSBjYW4gZGlzdGluZ3Vpc2ggYmV0d2VlbiB0d28gdHlwZXMgb2YgYXNzb2NpYXRpb25zLCBkaXJlY3QgYXNzb2NpYXRpb25zIGFuZCBpbmRpcmVjdCBhc3NvY2lhdGlvbi4gW14zXSBBbiBpbmRpcmVjdCBhc3NvY2lhdGlvbiBhcmlzZXMgd2hlbiBhIGNvdmFyaWF0ZSwgZm9yIGV4YW1wbGUgYHg5YCwgaXMgYXNzb2NpYXRlZCB3aXRoIGEgcmVzcG9uc2Ugb25seSBiZWNhdXNlIGl0IGlzIGFzc29jaWF0ZWQgd2l0aCBvbmUgb3IgbW9yZSBvdGhlciBjb3ZhcmlhdGVzIHRoYXQgZG8gaGF2ZSBhIGRpcmVjdCBhc3NvY2lhdGlvbiB3aXRoIHRoZSByZXNwb25zZSwgaW4gdGhpcyBleGFtcGxlIGB4MWAgYW5kIGB4M2AuIElmIHdlIGxvb2sgb25seSBhdCBgeDlgIGFuZCBgeWAsIHdlJ2xsIGZpbmQgYW4gYXNzb2NpYXRpb24uIFdoZW4gd2UgaW5jbHVkZSBgeDFgLWB4OWAgaW4gYSBtdWx0aXBsZSByZWdyZXNzaW9uLCBgeDlgIG5vdyBoYXMgb25seSBhIG5lZ2xpZ2libGUgYXNzb2NpYXRpb24uIEEgbXVsdGl2YXJpYXRlIHJlZ3Jlc3Npb24gX21heV8gYWxsb3cgdXMgdG8gZGlzdGluZ3Vpc2ggInJlYWwiIGZyb20gInNwdXJpb3VzIiBhc3NvY2lhdGlvbnMuIFteNF0KCiMjIFBvaW50cyBsZWZ0IGFzIGFuIGV4ZXJjaXNlCgpJJ3ZlIGlsbHVzdHJhdGVkIHRoZSBwaGVub21lbmEgaGVyZSB3aXRoIG9ubHkgb25lIHNpbXVsYXRlZCBkYXRhIHNldCB1bmRlciBlYWNoIHNjZW5hcmlvLiBBIHJpZ29yb3VzIGFuYWx5c2lzIHdvdWxkIHJlcXVpcmUgYSBsb3Qgb2YgZmFpcmx5IGNvbXBsaWNhdGVkIG1hdGgsIHdoaWNoIHlvdSdyZSB3ZWxjb21lIHRvIGxvb2sgdXAgaW4gYSBib29rIG9uIG11bHRpcGxlIHJlZ3Jlc3Npb24uIElmIHlvdSBkb24ndCB3YW50IHRvIGRvIHRoYXQsIHlvdSBtaWdodCB3YW50IHRvIHNldCB1cCBhIHNtYWxsIHNpbXVsYXRpb24gc3R1ZHkgaW4gd2hpY2ggeW91IHByb2R1Y2UgYSBodW5kcmVkIG9yIGEgdGhvdXNhbmQgZGF0YSBzZXRzIHVuZGVyIGNvbmRpdGlvbnMgc2ltaWxhciB0byB0aG9zZSBhYm92ZSB0byBtYWtlIHN1cmUgdGhlIHBhdHRlcm4gc2VlbiBpbiB0aGlzIG9uZSBleGFtcGxlIGlzIHJlcHJvZHVjaWJsZSBhY3Jvc3MgYSBsYXJnZSBudW1iZXIgb2YgZXhhbXBsZXMuIEFuZCBhcyBsb25nIGFzIHlvdSdyZSBkb2luZyB0aGF0LCB5b3UgbWlnaHQgYWxzbyB3YW50IHRvIGdlbmVyYXRlIGRhdGEgc2V0cyBvZiBkaWZmZXJlbnQgc2l6ZXMgdG8gZXhwbG9yZSBob3cgdGhlc2UgcmVzdWx0cyBkZXBlbmQgb24gc2FtcGxlIHNpemUgcmVsYXRpdmUgdG8gdGhlIG51bWJlciBvZiBwcmVkaWN0b3JzLgoKW14xXTogSWYgeW91J3JlIG5vdCBmYW1pbGlhciB3aXRoIHRoZSBkaWZmZXJlbmNlIGJldHdlZW4gIndpZGUiIGFuZCAibG9uZyIgZm9ybWF0LCBHb29nbGUgd2lsbCByZXR1cm4gYSB2YXJpZXR5IG9mIGxpbmtzLCBidXQgSGFkbGV5IFdpY2toYW0ncyBbdGlkeSBkYXRhIHBhcGVyXShodHRwczovL3ZpdGEuaGFkLmNvLm56L3BhcGVycy90aWR5LWRhdGEucGRmKSBpcyB0aGUgZGVmaW5pdGl2ZSBzb3VyY2UuClteMl06IFRoZSByZWFzb24gZm9yIHRoZSBzY2FyZSBxdW90ZXMgYXJvdW5kICJzcHVyaW91cyIgd2lsbCBiZWNvbWUgYXBwYXJlbnQgaW4gYSBsYXRlciBwb3N0LgpbXjNdOiBJJ20gbm90IGdvaW5nIHRvIGRlZmluZSAiZGlyZWN0IiBvciAiaW5kaXJlY3QuIiBJJ2xsIGRlYWwgd2l0aCB0aGUgZGlzdGluY3Rpb24gd2hlbiBJIGdldCB0byBkaXNjdXNzaW5nICJzcHVyaW91cyIgY29ycmVsYXRpb25zIGluIGEgbGF0ZXIgcG9zdC4KW140XTogWW91IGNhbiBwcm9iYWJseSBndWVzcyBieSBub3cgdGhhdCBJJ2xsIGJlIHJldHVybmluZyB0byB0aGlzIGxhdGVyLg==