Do we do better when covariates are less associated?
Introduction
Unlike previous notebooks, there isn’t anything new here. What you’ll see instead is analysis of two data sets generated with less association among the covariates than the data sets you’ve seen in earlier notebooks. The analysis will use rstanarm, horseshoe priors, and projection prediction variable selection.
NOTE: Data the data not scaled so that regression coefficients are more comparable across the two data sets. You may want to return to some of the earlier notebooks where the data were scaled and re-run the analyses to see how the results change.
library(tidyverse)
library(reshape2)
library(ggplot2)
library(cowplot)
library(mvtnorm)
library(corrplot)
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_str <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.8))
resid <- rep(2.0, n_samp)
y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_1 <- data.frame(y_str, cov_str, rep("Strong", length(y_str)))
cov_str <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.8))
y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_2 <- 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_1) <- column_names
colnames(dat_2) <- column_namesCorrelation plots with rho = 0.8
corrplot(cor(dat_1[, -c(1,11)]))
corrplot(cor(dat_2[, -c(1,11)]))
corr_08_1 <- max(abs(cor(dat_1[, -c(1,11)]) - diag(1, 9)))
corr_08_2 <- max(abs(cor(dat_2[, -c(1,11)]) - diag(1, 9)))Correlation plots with rho = 0.2
n_samp <- 100
cov_str <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.2))
resid <- rep(2.0, n_samp)
y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_1 <- data.frame(y_str, cov_str, rep("Strong", length(y_str)))
cov_str <- rmvnorm(n_samp,
mean = rep(0, nrow(Rho)),
sigma = construct_Sigma(Rho, sigma, 0.2))
y_str <- rnorm(nrow(cov_str), mean = beta0 + cov_str %*% beta, sd = resid)
dat_2 <- 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_1) <- column_names
colnames(dat_2) <- column_names
corrplot(cor(dat_1[, -c(1,11)]))
corrplot(cor(dat_2[, -c(1,11)]))
corr_02_1 <- max(abs(cor(dat_1[, -c(1,11)]) - diag(1, 9)))
corr_02_2 <- max(abs(cor(dat_2[, -c(1,11)]) - diag(1, 9)))The correlation plots make it clear that the magnitude of associations among covariates is much greater with rho = 0.8 than it is with rho = 0.2, but let’s look at the numerical results.
The maximum correlation in data set 1 (rho = 0.8): 0.834
The maximum correlation in data set 2 (rho = 0.8): 0.861
The maximum correlation in data set 1 (rho = 0.2): 0.365
The maximum correlation in data set 2 (rho = 0.2): 0.328
In everything that follows, we’ll be working with the data sets generated when rho = 0.2.
Horseshoe priors
We’ll be working with the data set generated with relatively weak associations among covariates for the rest of this exercise, and we’ll start with horsehoe priors.1
library(rstanarm)
options(mc.cores = parallel::detectCores())
set.seed(1234)
## n is the number of observations
## D is the number of covariates
## p0 is the expected number of important covariates
##
n <- nrow(dat_1)
D <- ncol(dat_1[,2:10])
p0 <- 3
tau0 <- p0/(D - p0) * 1/sqrt(n)
prior_coeff <- hs(global_scale = tau0, slab_scale = 1)
fit_1 <- stan_glm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9,
data = dat_1,
prior = prior_coeff,
refresh = 0,
adapt_delta = 0.999)
## Note: I use the same prior scale here because both data sets have the same number of observations
## and the same expected number of important covariates
##
fit_2 <- stan_glm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9,
data = dat_2,
prior = prior_coeff,
refresh = 0,
adapt_delta = 0.999)
predict_1 <- posterior_predict(fit_1)
predict_2 <- posterior_predict(fit_2)
for.plot <- data.frame(Observed = dat_1$y, Predicted = apply(predict_1, 2, mean))
p <- ggplot(for.plot, aes(x = Observed, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0) +
geom_smooth(method = "lm") +
ggtitle("Data set 1")
print(p)
for.plot <- data.frame(Observed = dat_2$y, Predicted = apply(predict_2, 2, mean))
p <- ggplot(for.plot, aes(x = Observed, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0) +
geom_smooth(method = "lm") +
ggtitle("Data set 2")
print(p)
p <- plot(fit_1) + ggtitle("Estimates from data set 1")
print(p)
p <- plot(fit_2) + ggtitle("Estimates from data set 2")
print(p)
summary(fit_1, digits = 3)
Model Info:
function: stan_glm
family: gaussian [identity]
formula: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
algorithm: sampling
priors: see help('prior_summary')
sample: 4000 (posterior sample size)
observations: 100
predictors: 10
Estimates:
mean sd 2.5% 25% 50% 75% 97.5%
(Intercept) 1.108 0.206 0.709 0.970 1.111 1.244 1.514
x1 1.149 0.198 0.762 1.013 1.148 1.284 1.535
x2 -0.892 0.217 -1.309 -1.039 -0.892 -0.746 -0.470
x3 1.073 0.202 0.662 0.941 1.074 1.210 1.459
x4 -0.150 0.177 -0.550 -0.261 -0.112 -0.012 0.120
x5 -0.003 0.147 -0.327 -0.069 -0.001 0.065 0.307
x6 -0.001 0.138 -0.307 -0.061 -0.001 0.056 0.315
x7 0.000 0.120 -0.264 -0.056 0.000 0.055 0.262
x8 -0.008 0.147 -0.339 -0.075 0.000 0.063 0.306
x9 -0.018 0.141 -0.327 -0.083 -0.006 0.047 0.278
sigma 1.981 0.152 1.711 1.878 1.971 2.074 2.316
mean_PPD 1.014 0.281 0.444 0.829 1.011 1.199 1.563
log-posterior -270.264 4.657 -280.473 -273.091 -269.894 -267.013 -262.167
Diagnostics:
mcse Rhat n_eff
(Intercept) 0.003 1.001 3728
x1 0.003 1.000 4292
x2 0.003 1.001 4166
x3 0.003 1.000 4498
x4 0.003 1.000 2570
x5 0.002 1.000 4963
x6 0.002 1.000 4470
x7 0.002 1.001 4923
x8 0.002 0.999 4766
x9 0.002 0.999 4788
sigma 0.002 0.999 4160
mean_PPD 0.005 1.002 3799
log-posterior 0.139 1.002 1127
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(fit_2, digits = 3)
Model Info:
function: stan_glm
family: gaussian [identity]
formula: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
algorithm: sampling
priors: see help('prior_summary')
sample: 4000 (posterior sample size)
observations: 100
predictors: 10
Estimates:
mean sd 2.5% 25% 50% 75% 97.5%
(Intercept) 0.902 0.197 0.508 0.770 0.898 1.038 1.284
x1 0.637 0.221 0.185 0.490 0.637 0.786 1.066
x2 -0.726 0.205 -1.108 -0.867 -0.733 -0.590 -0.312
x3 0.705 0.263 0.133 0.536 0.715 0.886 1.184
x4 -0.097 0.157 -0.461 -0.188 -0.062 0.002 0.166
x5 0.211 0.206 -0.097 0.042 0.187 0.349 0.659
x6 -0.054 0.144 -0.389 -0.123 -0.025 0.023 0.205
x7 0.045 0.143 -0.229 -0.030 0.020 0.117 0.378
x8 -0.065 0.141 -0.391 -0.138 -0.036 0.014 0.191
x9 -0.002 0.131 -0.276 -0.060 0.000 0.058 0.273
sigma 1.942 0.149 1.681 1.837 1.932 2.034 2.261
mean_PPD 1.025 0.276 0.482 0.842 1.027 1.212 1.569
log-posterior -266.824 4.755 -276.855 -269.812 -266.545 -263.574 -258.373
Diagnostics:
mcse Rhat n_eff
(Intercept) 0.003 0.999 4897
x1 0.003 1.000 4039
x2 0.003 1.000 4121
x3 0.005 1.001 3248
x4 0.002 1.001 4304
x5 0.004 1.000 3225
x6 0.002 1.001 4771
x7 0.002 1.000 4911
x8 0.002 1.000 4458
x9 0.002 1.000 4932
sigma 0.002 1.000 3627
mean_PPD 0.004 1.000 4479
log-posterior 0.129 1.001 1357
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).new_data <- data.frame(x1 = 4.0, x2 = 4.0, x3 = 4.0,
x4 = 4.0, x5 = 4.0, x6 = 4.0,
x7 = 4.0, x8 = 4.0, x9 = 4.0)
predict_1 <- posterior_predict(fit_1, new_data)
predict_2 <- posterior_predict(fit_2, new_data)
summarize_posterior <- function(x, credible = 0.95, digits = 3) {
lo_p <- (1.0 - credible)/2.0
hi_p <- credible + lo_p
ci <- quantile(x, c(lo_p, hi_p))
cat(round(mean(x), 3), " (", round(ci[1], 3), ",", round(ci[2], 3), ")\n", sep = "")
}
summarize_posterior(predict_1)5.687 (0.533,10.652)cat(" True answer: ", beta0 + as.matrix(new_data) %*% beta, "\n",
sep = "") True answer: 5summarize_posterior(predict_2)3.533 (-1.668,8.765)cat(" True answer: ", beta0 + as.matrix(new_data) %*% beta, "\n",
sep = "") True answer: 5That looks pretty encouraging. In both cases, x1, x2, and x3 are picked out as being important, and the estimated coefficients for each are relatively close to the values used to generate the data. The evidence for an association involving any of the other covariates is pretty weak.
Projection prediction
Now let’s try projection prediction to see how many variables we’d include in a model and which ones they are if we wanted explicitly to simplify the model rather than directly interpreting coefficients of the full model.
library(rstanarm)
library(projpred)
library(bayesplot)
options(mc.cores = parallel::detectCores())
projection_prediction <- function(dat) {
n <- nrow(dat)
D <- ncol(dat[,2:10])
p0 <- 3
tau0 <- p0/(D - p0) * 1/sqrt(n)
prior_coeff <- hs(global_scale = tau0, slab_scale = 1)
fit <- stan_glm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9, data = dat,
prior = prior_coeff, refresh = 0)
vs <- varsel(fit, method = "forward")
cvs <- cv_varsel(fit, method = "forward", verbose = FALSE)
proj_size <- suggest_size(cvs)
proj <- project(vs, nv = proj_size, ns = 2000)
mcmc_intervals(as.matrix(proj),
pars = c("(Intercept)", names(vs$vind[1:proj_size]), "sigma"))
pred <- proj_linpred(vs, xnew = dat[, 2:10], ynew = dat$y,
nv = proj_size, integrated = TRUE)
for.plot <- data.frame(Observed = dat$y,
Predicted = pred$pred)
p <- ggplot(for.plot, aes(x = Observed, y = Predicted)) +
geom_point() +
geom_smooth(method = "lm") +
geom_abline(slope = 1, intercept = 0) +
labs(x = "Observed", y = "Predicted")
print(p)
return(list(proj = proj, vs = vs, cvs = cvs, proj_size = proj_size))
}
summarize_posterior <- function(x, credible = 0.95, digits = 3) {
lo_p <- (1.0 - credible)/2.0
hi_p <- credible + lo_p
ci <- quantile(x, c(lo_p, hi_p))
results <- sprintf("% 5.3f (% 5.3f, % 5.3f)\n", mean(x), ci[1], ci[2])
cat(results, sep = "")
}
summarize_results <- function(x, credible = 0.95, digits = 3) {
vars <- colnames(as.matrix(x))
for (i in 1:length(vars)) {
label <- sprintf("%11s: ", vars[i])
cat(label, sep = "")
summarize_posterior(as.matrix(x)[,i])
}
}
proj_1 <- projection_prediction(dat_1)
proj_2 <- projection_prediction(dat_2)
print(varsel_plot(proj_1$cvs, stats = c("elpd", "rmse"), deltas = TRUE))
print(varsel_plot(proj_2$cvs, stats = c("elpd", "rmse"), deltas = TRUE))
summarize_results(proj_1$proj)(Intercept): 1.087 ( 0.711, 1.466)
x3: 1.081 ( 0.709, 1.479)
x1: 1.159 ( 0.770, 1.529)
x2: -0.937 (-1.363, -0.496)
sigma: 2.018 ( 1.729, 2.342)summarize_results(proj_2$proj)(Intercept): 0.920 ( 0.534, 1.304)
x1: 0.682 ( 0.244, 1.113)
x2: -0.794 (-1.145, -0.403)
x3: 0.788 ( 0.282, 1.242)
sigma: 1.986 ( 1.716, 2.321)Good news! Using projection prediction to select the “important” covariates identified the same three covariates our examination of the full model suggested to us and those are the variables we know ore important.
Bottom line: If associations among your covariates are modest, if you have a reasonably large amount of data, and if your covariates include all of those that have a real relationship with the response variable, then you have a good chance of detecting the relationships and of estimating their magnitude. But what happens if the “true” predictors are among the covariates included in your study?
What if the “true” predictors aren’t in the data?
In our case the true predictors are x1, x2, and x3. We know that they influence influence y and that no other covariates do, because we generated the data that way. What happens if we fit a model that excludes x1, x2, and x3 and includes only covariates that we know don’t have a “real” relationship to y?
n <- nrow(dat_1)
D <- ncol(dat_1[,5:10])
p0 <- 3
tau0 <- p0/(D - p0) * 1/sqrt(n)
prior_coeff <- hs(global_scale = tau0, slab_scale = 1)
fit_1 <- stan_glm(y ~ x4 + x5 + x6 + x7 + x8 + x9,
data = dat_1,
prior = prior_coeff,
refresh = 0,
adapt_delta = 0.999)
## Note: I use the same prior scale here because both data sets have the same number of observations
## and the same expected number of important covariates
##
fit_2 <- stan_glm(y ~ x4 + x5 + x6 + x7 + x8 + x9,
data = dat_2,
prior = prior_coeff,
refresh = 0,
adapt_delta = 0.999)
predict_1 <- posterior_predict(fit_1)
predict_2 <- posterior_predict(fit_2)
for.plot <- data.frame(Observed = dat_1$y, Predicted = apply(predict_1, 2, mean))
p <- ggplot(for.plot, aes(x = Observed, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0) +
geom_smooth(method = "lm") +
ggtitle("Data set 1")
print(p)
for.plot <- data.frame(Observed = dat_2$y, Predicted = apply(predict_2, 2, mean))
p <- ggplot(for.plot, aes(x = Observed, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0) +
geom_smooth(method = "lm") +
ggtitle("Data set 2")
print(p)
p <- plot(fit_1) + ggtitle("Estimates from data set 1")
print(p)
p <- plot(fit_2) + ggtitle("Estimates from data set 2")
print(p)
summary(fit_1, digits = 3)
Model Info:
function: stan_glm
family: gaussian [identity]
formula: y ~ x4 + x5 + x6 + x7 + x8 + x9
algorithm: sampling
priors: see help('prior_summary')
sample: 4000 (posterior sample size)
observations: 100
predictors: 7
Estimates:
mean sd 2.5% 25% 50% 75% 97.5%
(Intercept) 1.054 0.280 0.501 0.872 1.060 1.239 1.607
x4 -0.305 0.288 -0.906 -0.513 -0.266 -0.037 0.074
x5 0.037 0.170 -0.304 -0.029 0.009 0.101 0.448
x6 0.040 0.172 -0.286 -0.031 0.010 0.100 0.476
x7 0.030 0.152 -0.280 -0.033 0.007 0.088 0.386
x8 -0.198 0.259 -0.822 -0.359 -0.107 -0.003 0.128
x9 0.194 0.244 -0.115 0.003 0.120 0.339 0.794
sigma 2.833 0.216 2.449 2.682 2.821 2.970 3.298
mean_PPD 1.020 0.393 0.219 0.764 1.022 1.283 1.786
log-posterior -287.483 4.124 -296.477 -290.003 -287.099 -284.641 -280.352
Diagnostics:
mcse Rhat n_eff
(Intercept) 0.005 0.999 3389
x4 0.006 1.000 2153
x5 0.003 1.000 3722
x6 0.003 1.000 3958
x7 0.002 1.000 4489
x8 0.005 1.000 2390
x9 0.005 1.001 2438
sigma 0.003 0.999 4213
mean_PPD 0.007 1.000 3247
log-posterior 0.130 1.003 1001
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(fit_2, digits = 3)
Model Info:
function: stan_glm
family: gaussian [identity]
formula: y ~ x4 + x5 + x6 + x7 + x8 + x9
algorithm: sampling
priors: see help('prior_summary')
sample: 4000 (posterior sample size)
observations: 100
predictors: 7
Estimates:
mean sd 2.5% 25% 50% 75% 97.5%
(Intercept) 0.961 0.240 0.493 0.800 0.957 1.119 1.452
x4 -0.174 0.199 -0.628 -0.308 -0.124 -0.013 0.096
x5 0.409 0.271 -0.026 0.200 0.417 0.596 0.944
x6 -0.225 0.232 -0.746 -0.380 -0.178 -0.024 0.081
x7 0.102 0.175 -0.171 -0.006 0.051 0.193 0.522
x8 -0.057 0.152 -0.439 -0.128 -0.023 0.019 0.216
x9 0.049 0.146 -0.220 -0.023 0.017 0.115 0.393
sigma 2.278 0.167 1.982 2.158 2.265 2.385 2.634
mean_PPD 1.024 0.327 0.395 0.803 1.024 1.242 1.676
log-posterior -265.191 4.167 -274.183 -267.880 -264.860 -262.205 -258.010
Diagnostics:
mcse Rhat n_eff
(Intercept) 0.003 1.000 5107
x4 0.004 1.000 3129
x5 0.006 1.002 2074
x6 0.005 1.001 2398
x7 0.003 1.001 3115
x8 0.002 0.999 4387
x9 0.002 0.999 4632
sigma 0.003 1.000 4136
mean_PPD 0.005 1.000 4396
log-posterior 0.127 1.002 1074
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).new_data <- data.frame(x4 = 4.0, x5 = 4.0, x6 = 4.0,
x7 = 4.0, x8 = 4.0, x9 = 4.0)
predict_1 <- posterior_predict(fit_1, new_data)
predict_2 <- posterior_predict(fit_2, new_data)
summarize_posterior <- function(x, credible = 0.95, digits = 3) {
lo_p <- (1.0 - credible)/2.0
hi_p <- credible + lo_p
ci <- quantile(x, c(lo_p, hi_p))
cat(round(mean(x), 3), " (", round(ci[1], 3), ",", round(ci[2], 3), ")\n", sep = "")
}
summarize_posterior(predict_1)0.214 (-6.214,6.717)summarize_posterior(predict_2)1.371 (-4.324,7.054)The results from analysis of data set 1 suggest (weakly) that x4 and x8 have a negative association with y, and that x9 has a positive associations. That’s right in the sense that (a) x4 and x8 are positively associated with x2, which we know has a negative influence on y and (b) x9 is positively associated with both x1 and x3, which we know have positive influences on y. In data set 2, x4 and x6 have a weakly supported negative association with y, while x5 has more strongly supported positive association with y. Again, those associations are in the direction we’d expect. But the associations are weak. Our evidence for them is shaky.
What happens if we try the projection predictions method for variable selection?
projection_prediction <- function(dat) {
n <- nrow(dat)
D <- ncol(dat) - 1
p0 <- 3
tau0 <- p0/(D - p0) * 1/sqrt(n)
prior_coeff <- hs(global_scale = tau0, slab_scale = 1)
fit <- stan_glm(y ~ x4 + x5 + x6 + x7 + x8 + x9, data = dat,
prior = prior_coeff, refresh = 0)
vs <- varsel(fit, method = "forward")
cvs <- cv_varsel(fit, method = "forward", verbose = FALSE)
proj_size <- suggest_size(cvs)
proj <- project(vs, nv = proj_size, ns = 2000)
mcmc_intervals(as.matrix(proj))
pred <- proj_linpred(vs, xnew = dat[, -1], ynew = dat$y,
nv = proj_size, integrated = TRUE)
for.plot <- data.frame(Observed = dat$y,
Predicted = pred$pred)
p <- ggplot(for.plot, aes(x = Observed, y = Predicted)) +
geom_point() +
geom_smooth(method = "lm") +
geom_abline(slope = 1, intercept = 0) +
labs(x = "Observed", y = "Predicted")
print(p)
return(list(proj = proj, vs = vs, cvs = cvs, proj_size = proj_size))
}
proj_1 <- projection_prediction(dat_1[, -c(2,3,4,11)])
proj_2 <- projection_prediction(dat_2[, -c(2,3,4,11)])
print(varsel_plot(proj_1$cvs, stats = c("elpd", "rmse"), deltas = TRUE))
print(varsel_plot(proj_2$cvs, stats = c("elpd", "rmse"), deltas = TRUE))
summarize_results(proj_1$proj)(Intercept): 1.011 (0.456,1.542)
sigma: 2.902 (2.539,3.341)summarize_results(proj_2$proj)(Intercept): 1.007 (0.564,1.452)
x5: 0.482 (0.016,0.963)
sigma: 2.327 (2.032,2.695)I said that the data set 1 suggested some relationships weakly, and the projection predictions show that the relationships are so weak that we don’t have good evidence for including any of them. There’s reasonable evidence for including one covariate, x5, in the second data set, which means we don’t have any indication of any negative associations with y.[^2]
Conclusions
If associations among our covariates are relatively weak and we have a reasonable amount of data, it’s tempting to think that associations we detect in our data are likely to be real. It’s certainly the case that we’re less likely to be misled when associations are wake, but try this simple thought experiment to realize why you should be cautious.
Imagine that you’re interested in how leaf mass per area (LMA) varies across environmental gradients. Suppose you’re working in a mountainous region of the North Temperate Zone and that you include elevation as a covariate without including any covariates related to temperature (e.g., mean annual temperature, January minimum temperature, number of frost-free days). I suspect you;ll find a positive association between LMA and elevation, but I seriously doubt that it’s elevation per se that’s important. Rather, I suspect that the relationship, assuming that it exists, would reflect the tight association between different aspects of temperature and elevation. Even if you included a temperature covariate, e.g., frost-free days, and you now found an association between LMA and frost-free days while the elevation association was now weak or non-existent, I’d suggest that you think very carefully before concluding that frost-free days really matters. It might be one of several temperature-related covariates that matter, or it might simply be tightly associated with something else that really matters. In the absence of a controlled experiment or a very careful causal analysis,{^3} all we can really say from the data alone is that there’s an association.
In short, I urge you to be very cautious in interpreting the results of a multiple regression. If the associations you identify have been found repeatedly in other data sets and there are good, principled reasons to believe that those associations arose through the process you are trying to attribute them to, I’d say you’re on reasonably solid ground. The less often the associations you see have been seen before and the less they might have been predicted from first principles, the shakier is the ground on which you stand. But being cautious in interpreting the meaning of associations you find, doesn’t mean that those associations are unimportant. Associations that appear important in your analysis are worth further study, and that further study may lead both to finding them repeatedly and developing an understanding of mechanisms and processes that could produce them.
{^2]: I havern’t tried running the models that exclude x1, x2, and x3 with the strong associations (0.8 in the call to rmvnorm()), but I encourage you to try it on your own. I suspect you’ll find some reasonably well supported associations if you do.
You should find it pretty easy to try the Lasso using the least squares regression in
lm()if you’re so inclined.↩
LS0tCnRpdGxlOiAiRG8gd2UgZG8gYmV0dGVyIHdoZW4gY292YXJpYXRlcyBhcmUgbGVzcyBhc3NvY2lhdGVkPyIKb3V0cHV0OiAKICBodG1sX25vdGVib29rOgogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6IAogICAgICBjb2xsYXBzZWQ6IGZhbHNlCi0tLQojIyBJbnRyb2R1Y3Rpb24KClVubGlrZSBwcmV2aW91cyBub3RlYm9va3MsIHRoZXJlIGlzbid0IGFueXRoaW5nIG5ldyBoZXJlLiBXaGF0IHlvdSdsbCBzZWUgaW5zdGVhZCBpcyBhbmFseXNpcyBvZiB0d28gZGF0YSBzZXRzIGdlbmVyYXRlZCB3aXRoIGxlc3MgYXNzb2NpYXRpb24gYW1vbmcgdGhlIGNvdmFyaWF0ZXMgdGhhbiB0aGUgZGF0YSBzZXRzIHlvdSd2ZSBzZWVuIGluIGVhcmxpZXIgbm90ZWJvb2tzLiBUaGUgYW5hbHlzaXMgd2lsbCB1c2UgYHJzdGFuYXJtYCwgaG9yc2VzaG9lIHByaW9ycywgYW5kIHByb2plY3Rpb24gcHJlZGljdGlvbiB2YXJpYWJsZSBzZWxlY3Rpb24uCgpOT1RFOiBEYXRhIHRoZSBkYXRhIG5vdCBzY2FsZWQgc28gdGhhdCByZWdyZXNzaW9uIGNvZWZmaWNpZW50cyBhcmUgX21vcmVfIGNvbXBhcmFibGUgYWNyb3NzIHRoZSB0d28gZGF0YSBzZXRzLiBZb3UgbWF5IHdhbnQgdG8gcmV0dXJuIHRvIHNvbWUgb2YgdGhlIGVhcmxpZXIgbm90ZWJvb2tzIHdoZXJlIHRoZSBkYXRhIHdlcmUgc2NhbGVkIGFuZCByZS1ydW4gdGhlIGFuYWx5c2VzIHRvIHNlZSBob3cgdGhlIHJlc3VsdHMgY2hhbmdlLgoKYGBge3Igc2V0dXAsIHdhcm5pbmcgPSBGQUxTRSwgbWVzc2FnZSA9IEZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShyZXNoYXBlMikKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGNvd3Bsb3QpCmxpYnJhcnkobXZ0bm9ybSkKbGlicmFyeShjb3JycGxvdCkKCnJtKGxpc3QgPSBscygpKQpgYGAKCmBgYHtyfQojIyBpbnRldGNlcHQKIyMKYmV0YTAgPC0gMS4wCiMjIHJlZ3Jlc3Npb24gY29lZmZpY2llbnRzCiMjCmJldGEgPC0gYygxLjAsIC0xLjAsIDEuMCwgMC4wLCAwLjAsIDAuMCwgMC4wLCAwLjAsIDAuMCkKIyMgcGF0dGVybiBvZiBjb3JyZWxhdGlvbiBtYXRyaXgsIGFsbCBub24temVybyBlbnRyaWVzIGFyZSBzZXQgdG8gc2FlbQojIyBjb3JyZWxhdGlvbiwgY292YXJpYW5jZSBtYXRyaXggY2FsZHVsYXRlZCBmcm9tIGluZGl2aWR1YWwgdmFyaWFuY2VzIGFuZCBhIAojIyBzaW5nbGUgYXNzb2NpYXRpb24gcGFyYW1ldGVyIGdvdmVybmluZyB0aGUgbm9uLXplcm8gY29ycmVsYXRpb24gY29lZmZpY2llbnRzCiMjCiMjIE5vdGU6IE5vdCBqdXN0IGFueSBwYXR0ZXJuIHdpbGwgd29yayBoZXJlLiBUaGUgY29ycmVsYXRpb24gbWF0cml4IGFuZAojIyBjb3ZhcmlhbmNlIG1hdHJpeCBnZW5lcmF0ZWQgZnJvbSB0aGlzIHBhdHRlcm4gbXVzdCBiZSBwb3NpdGl2ZSBkZWZpbml0ZS4KIyMgSWYgeW91IGNoYW5nZSB0aGlzIHBhdHRlcm4sIHlvdSBtYXkgZ2V0IGFuIGVycm9yIHdoZW4geW91IHRyeSB0byBnZW5lcmF0ZQojIyBkYXRhIHdpdGggYSBub24temVybyBhc3NvY2lhdGlvbiBwYXJhbWV0ZXIuCiMjClJobyA8LSBtYXRyaXgobnJvdyA9IDksIG5jb2wgPSAsIGJ5cm93ID0gVFJVRSwgCiAgICAgICAgICAgICAgZGF0YSA9IGMoMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEsCiAgICAgICAgICAgICAgICAgICAgICAgMCwxLDAsMSwwLDEsMCwxLDAsCiAgICAgICAgICAgICAgICAgICAgICAgMSwwLDEsMCwxLDAsMSwwLDEKICAgICAgICAgICAgICAgICAgICAgICApKQojIyB2ZWN0b3Igb2Ygc3RhbmRhcmQgZGV2aWF0aW9ucyBmb3IgY292YXJpYXRlcwojIwpzaWdtYSA8LSByZXAoMSwgOSkKCiMjIGNvbnN0cnVjdCBhIGNvdmFyaWFuY2UgbWF0cml4IGZyb20gdGhlIHBhdHRlcm4sIHN0YW5kYXJkIGRldmlhdGlvbnMsIGFuZAojIyBvbmUgcGFyYW1ldGVyIGluIFstMSwxXSB0aGF0IGdvdmVybnMgdGhlIG1hZ25pdHVkZSBvZiBub24temVybyBjb3JyZWxhdGlvbgojIyBjb2VmZmljaWVudHMKIyMKIyMgUmhvIC0gdGhlIHBhdHRlcm4gb2YgYXNzb2NpYXRpb25zCiMjIHNpZ21hIC0gdGhlIHZlY3RvciBvZiBzdGFuZGFyZCBkZXZpYXRpb25zCiMjIHJobyAtIHRoZSBhc3NvY2lhdGlvbiBwYXJhbWV0ZXIKIyMKY29uc3RydWN0X1NpZ21hIDwtIGZ1bmN0aW9uKFJobywgc2lnbWEsIHJobykgewogICMjIGdldCB0aGUgY29ycmVsYXRpb24gbWF0cmlzCiAgIyMKICBSaG8gPC0gUmhvKnJobwogIGZvciAoaSBpbiAxOm5jb2woUmhvKSkgewogICAgUmhvW2ksaV0gPC0gMS4wCiAgfQogICMjIG5vdGljZSB0aGUgdXNlIG9mIG1hdHJpeCBtdWx0aXBsaWNhdGlvbgogICMjCiAgU2lnbWEgPC0gZGlhZyhzaWdtYSkgJSolIFJobyAlKiUgZGlhZyhzaWdtYSkKICByZXR1cm4oU2lnbWEpCn0KCiMjIHNldCB0aGUgcmFuZG9tIG51bWJlciBzZWVkIG1hbnVhbGx5IHNvIHRoYXQgZXZlcnkgcnVuIG9mIHRoZSBjb2RlIHdpbGwgCiMjIHByb2R1Y2UgdGhlIHNhbWUgbnVtYmVycwojIwpzZXQuc2VlZCgxMjM0KQoKbl9zYW1wIDwtIDEwMApjb3Zfc3RyIDwtIHJtdm5vcm0obl9zYW1wLAogICAgICAgICAgICAgICAgICAgbWVhbiA9IHJlcCgwLCBucm93KFJobykpLAogICAgICAgICAgICAgICAgICAgc2lnbWEgPSBjb25zdHJ1Y3RfU2lnbWEoUmhvLCBzaWdtYSwgMC44KSkKCnJlc2lkIDwtIHJlcCgyLjAsIG5fc2FtcCkKCnlfc3RyIDwtIHJub3JtKG5yb3coY292X3N0ciksIG1lYW4gPSBiZXRhMCArIGNvdl9zdHIgJSolIGJldGEsIHNkID0gcmVzaWQpCmRhdF8xIDwtIGRhdGEuZnJhbWUoeV9zdHIsIGNvdl9zdHIsIHJlcCgiU3Ryb25nIiwgbGVuZ3RoKHlfc3RyKSkpCgpjb3Zfc3RyIDwtIHJtdm5vcm0obl9zYW1wLAogICAgICAgICAgICAgICAgICAgbWVhbiA9IHJlcCgwLCBucm93KFJobykpLAogICAgICAgICAgICAgICAgICAgc2lnbWEgPSBjb25zdHJ1Y3RfU2lnbWEoUmhvLCBzaWdtYSwgMC44KSkKeV9zdHIgPC0gcm5vcm0obnJvdyhjb3Zfc3RyKSwgbWVhbiA9IGJldGEwICsgY292X3N0ciAlKiUgYmV0YSwgc2QgPSByZXNpZCkKZGF0XzIgPC0gZGF0YS5mcmFtZSh5X3N0ciwgY292X3N0ciwgcmVwKCJTdHJvbmciLCBsZW5ndGgoeV9zdHIpKSkKCmNvbHVtbl9uYW1lcyA8LSBjKCJ5IiwgcGFzdGUoIngiLCBzZXEoMSwgbGVuZ3RoKGJldGEpKSwgc2VwID0gIiIpLCAiU2NlbmFyaW8iKQpjb2xuYW1lcyhkYXRfMSkgPC0gY29sdW1uX25hbWVzCmNvbG5hbWVzKGRhdF8yKSA8LSBjb2x1bW5fbmFtZXMKYGBgCgojIyMgQ29ycmVsYXRpb24gcGxvdHMgd2l0aCBgcmhvID0gMC44YAoKYGBge3J9CmNvcnJwbG90KGNvcihkYXRfMVssIC1jKDEsMTEpXSkpCmNvcnJwbG90KGNvcihkYXRfMlssIC1jKDEsMTEpXSkpCmNvcnJfMDhfMSA8LSBtYXgoYWJzKGNvcihkYXRfMVssIC1jKDEsMTEpXSkgLSBkaWFnKDEsIDkpKSkKY29ycl8wOF8yIDwtIG1heChhYnMoY29yKGRhdF8yWywgLWMoMSwxMSldKSAtIGRpYWcoMSwgOSkpKQpgYGAKCiMjIyBDb3JyZWxhdGlvbiBwbG90cyB3aXRoIGByaG8gPSAwLjJgCgpgYGB7cn0Kbl9zYW1wIDwtIDEwMApjb3Zfc3RyIDwtIHJtdm5vcm0obl9zYW1wLAogICAgICAgICAgICAgICAgICAgbWVhbiA9IHJlcCgwLCBucm93KFJobykpLAogICAgICAgICAgICAgICAgICAgc2lnbWEgPSBjb25zdHJ1Y3RfU2lnbWEoUmhvLCBzaWdtYSwgMC4yKSkKCnJlc2lkIDwtIHJlcCgyLjAsIG5fc2FtcCkKCnlfc3RyIDwtIHJub3JtKG5yb3coY292X3N0ciksIG1lYW4gPSBiZXRhMCArIGNvdl9zdHIgJSolIGJldGEsIHNkID0gcmVzaWQpCmRhdF8xIDwtIGRhdGEuZnJhbWUoeV9zdHIsIGNvdl9zdHIsIHJlcCgiU3Ryb25nIiwgbGVuZ3RoKHlfc3RyKSkpCgpjb3Zfc3RyIDwtIHJtdm5vcm0obl9zYW1wLAogICAgICAgICAgICAgICAgICAgbWVhbiA9IHJlcCgwLCBucm93KFJobykpLAogICAgICAgICAgICAgICAgICAgc2lnbWEgPSBjb25zdHJ1Y3RfU2lnbWEoUmhvLCBzaWdtYSwgMC4yKSkKeV9zdHIgPC0gcm5vcm0obnJvdyhjb3Zfc3RyKSwgbWVhbiA9IGJldGEwICsgY292X3N0ciAlKiUgYmV0YSwgc2QgPSByZXNpZCkKZGF0XzIgPC0gZGF0YS5mcmFtZSh5X3N0ciwgY292X3N0ciwgcmVwKCJTdHJvbmciLCBsZW5ndGgoeV9zdHIpKSkKCmNvbHVtbl9uYW1lcyA8LSBjKCJ5IiwgcGFzdGUoIngiLCBzZXEoMSwgbGVuZ3RoKGJldGEpKSwgc2VwID0gIiIpLCAiU2NlbmFyaW8iKQpjb2xuYW1lcyhkYXRfMSkgPC0gY29sdW1uX25hbWVzCmNvbG5hbWVzKGRhdF8yKSA8LSBjb2x1bW5fbmFtZXMKCmNvcnJwbG90KGNvcihkYXRfMVssIC1jKDEsMTEpXSkpCmNvcnJwbG90KGNvcihkYXRfMlssIC1jKDEsMTEpXSkpCmNvcnJfMDJfMSA8LSBtYXgoYWJzKGNvcihkYXRfMVssIC1jKDEsMTEpXSkgLSBkaWFnKDEsIDkpKSkKY29ycl8wMl8yIDwtIG1heChhYnMoY29yKGRhdF8yWywgLWMoMSwxMSldKSAtIGRpYWcoMSwgOSkpKQpgYGAKClRoZSBjb3JyZWxhdGlvbiBwbG90cyBtYWtlIGl0IGNsZWFyIHRoYXQgdGhlIG1hZ25pdHVkZSBvZiBhc3NvY2lhdGlvbnMgYW1vbmcgY292YXJpYXRlcyBpcyBtdWNoIGdyZWF0ZXIgd2l0aCBgcmhvID0gMC44YCB0aGFuIGl0IGlzIHdpdGggYHJobyA9IDAuMmAsIGJ1dCBsZXQncyBsb29rIGF0IHRoZSBudW1lcmljYWwgcmVzdWx0cy4KClRoZSBtYXhpbXVtIGNvcnJlbGF0aW9uIGluIGRhdGEgc2V0IDEgKHJobyA9IDAuOCk6IGByIHJvdW5kKGNvcnJfMDhfMSwgMylgCgpUaGUgbWF4aW11bSBjb3JyZWxhdGlvbiBpbiBkYXRhIHNldCAyIChyaG8gPSAwLjgpOiBgciByb3VuZChjb3JyXzA4XzIsIDMpYAoKVGhlIG1heGltdW0gY29ycmVsYXRpb24gaW4gZGF0YSBzZXQgMSAocmhvID0gMC4yKTogYHIgcm91bmQoY29ycl8wMl8xLCAzKWAKClRoZSBtYXhpbXVtIGNvcnJlbGF0aW9uIGluIGRhdGEgc2V0IDIgKHJobyA9IDAuMik6IGByIHJvdW5kKGNvcnJfMDJfMiwgMylgCgpJbiBldmVyeXRoaW5nIHRoYXQgZm9sbG93cywgd2UnbGwgYmUgd29ya2luZyB3aXRoIHRoZSBkYXRhIHNldHMgZ2VuZXJhdGVkIHdoZW4gYHJobyA9IDAuMmAuCgojIyBIb3JzZXNob2UgcHJpb3JzCgpXZSdsbCBiZSB3b3JraW5nIHdpdGggdGhlIGRhdGEgc2V0IGdlbmVyYXRlZCB3aXRoIHJlbGF0aXZlbHkgd2VhayBhc3NvY2lhdGlvbnMgYW1vbmcgY292YXJpYXRlcyBmb3IgdGhlIHJlc3Qgb2YgdGhpcyBleGVyY2lzZSwgYW5kIHdlJ2xsIHN0YXJ0IHdpdGggaG9yc2Vob2UgcHJpb3JzLlteMV0KCmBgYHtyLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRX0KbGlicmFyeShyc3RhbmFybSkKCm9wdGlvbnMobWMuY29yZXMgPSBwYXJhbGxlbDo6ZGV0ZWN0Q29yZXMoKSkKCnNldC5zZWVkKDEyMzQpCgojIyBuIGlzIHRoZSBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiMjIEQgaXMgdGhlIG51bWJlciBvZiBjb3ZhcmlhdGVzCiMjIHAwIGlzIHRoZSBleHBlY3RlZCBudW1iZXIgb2YgaW1wb3J0YW50IGNvdmFyaWF0ZXMKIyMKbiA8LSBucm93KGRhdF8xKQpEIDwtIG5jb2woZGF0XzFbLDI6MTBdKQpwMCA8LSAzCnRhdTAgPC0gcDAvKEQgLSBwMCkgKiAxL3NxcnQobikKcHJpb3JfY29lZmYgPC0gaHMoZ2xvYmFsX3NjYWxlID0gdGF1MCwgc2xhYl9zY2FsZSA9IDEpCgpmaXRfMSA8LSBzdGFuX2dsbSh5IH4geDEgKyB4MiArIHgzICsgeDQgKyB4NSArIHg2ICsgeDcgKyB4OCArIHg5LAogICAgICAgICAgICAgICAgICBkYXRhID0gZGF0XzEsCiAgICAgICAgICAgICAgICAgIHByaW9yID0gcHJpb3JfY29lZmYsCiAgICAgICAgICAgICAgICAgIHJlZnJlc2ggPSAwLAogICAgICAgICAgICAgICAgICBhZGFwdF9kZWx0YSA9IDAuOTk5KQojIyBOb3RlOiBJIHVzZSB0aGUgc2FtZSBwcmlvciBzY2FsZSBoZXJlIGJlY2F1c2UgYm90aCBkYXRhIHNldHMgaGF2ZSB0aGUgc2FtZSBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiMjIGFuZCB0aGUgc2FtZSBleHBlY3RlZCBudW1iZXIgb2YgaW1wb3J0YW50IGNvdmFyaWF0ZXMKIyMKZml0XzIgPC0gc3Rhbl9nbG0oeSB+IHgxICsgeDIgKyB4MyArIHg0ICsgeDUgKyB4NiArIHg3ICsgeDggKyB4OSwKICAgICAgICAgICAgICAgICAgZGF0YSA9IGRhdF8yLAogICAgICAgICAgICAgICAgICBwcmlvciA9IHByaW9yX2NvZWZmLAogICAgICAgICAgICAgICAgICByZWZyZXNoID0gMCwKICAgICAgICAgICAgICAgICAgYWRhcHRfZGVsdGEgPSAwLjk5OSkKCnByZWRpY3RfMSA8LSBwb3N0ZXJpb3JfcHJlZGljdChmaXRfMSkKcHJlZGljdF8yIDwtIHBvc3Rlcmlvcl9wcmVkaWN0KGZpdF8yKQoKZm9yLnBsb3QgPC0gZGF0YS5mcmFtZShPYnNlcnZlZCA9IGRhdF8xJHksIFByZWRpY3RlZCA9IGFwcGx5KHByZWRpY3RfMSwgMiwgbWVhbikpCnAgPC0gZ2dwbG90KGZvci5wbG90LCBhZXMoeCA9IE9ic2VydmVkLCB5ID0gUHJlZGljdGVkKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGUgPSAxLCBpbnRlcmNlcHQgPSAwKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIikgKwogIGdndGl0bGUoIkRhdGEgc2V0IDEiKQpwcmludChwKQoKZm9yLnBsb3QgPC0gZGF0YS5mcmFtZShPYnNlcnZlZCA9IGRhdF8yJHksIFByZWRpY3RlZCA9IGFwcGx5KHByZWRpY3RfMiwgMiwgbWVhbikpCnAgPC0gZ2dwbG90KGZvci5wbG90LCBhZXMoeCA9IE9ic2VydmVkLCB5ID0gUHJlZGljdGVkKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGUgPSAxLCBpbnRlcmNlcHQgPSAwKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIikgKwogIGdndGl0bGUoIkRhdGEgc2V0IDIiKQpwcmludChwKQoKcCA8LSBwbG90KGZpdF8xKSArIGdndGl0bGUoIkVzdGltYXRlcyBmcm9tIGRhdGEgc2V0IDEiKQpwcmludChwKQpwIDwtIHBsb3QoZml0XzIpICsgZ2d0aXRsZSgiRXN0aW1hdGVzIGZyb20gZGF0YSBzZXQgMiIpCnByaW50KHApCgpzdW1tYXJ5KGZpdF8xLCBkaWdpdHMgPSAzKQpzdW1tYXJ5KGZpdF8yLCBkaWdpdHMgPSAzKQoKbmV3X2RhdGEgPC0gZGF0YS5mcmFtZSh4MSA9IDQuMCwgeDIgPSA0LjAsIHgzID0gNC4wLAogICAgICAgICAgICAgICAgICAgICAgIHg0ID0gNC4wLCB4NSA9IDQuMCwgeDYgPSA0LjAsCiAgICAgICAgICAgICAgICAgICAgICAgeDcgPSA0LjAsIHg4ID0gNC4wLCB4OSA9IDQuMCkKCnByZWRpY3RfMSA8LSBwb3N0ZXJpb3JfcHJlZGljdChmaXRfMSwgbmV3X2RhdGEpCnByZWRpY3RfMiA8LSBwb3N0ZXJpb3JfcHJlZGljdChmaXRfMiwgbmV3X2RhdGEpCgpzdW1tYXJpemVfcG9zdGVyaW9yIDwtIGZ1bmN0aW9uKHgsIGNyZWRpYmxlID0gMC45NSwgZGlnaXRzID0gMykgewogIGxvX3AgPC0gKDEuMCAtIGNyZWRpYmxlKS8yLjAKICBoaV9wIDwtIGNyZWRpYmxlICsgbG9fcAogIGNpIDwtIHF1YW50aWxlKHgsIGMobG9fcCwgaGlfcCkpCiAgY2F0KHJvdW5kKG1lYW4oeCksIDMpLCAiICgiLCByb3VuZChjaVsxXSwgMyksICIsIiwgcm91bmQoY2lbMl0sIDMpLCAiKVxuIiwgc2VwID0gIiIpCn0KCnN1bW1hcml6ZV9wb3N0ZXJpb3IocHJlZGljdF8xKQpjYXQoIiAgVHJ1ZSBhbnN3ZXI6ICIsIGJldGEwICsgYXMubWF0cml4KG5ld19kYXRhKSAlKiUgYmV0YSwgIlxuIiwgCiAgICBzZXAgPSAiIikKc3VtbWFyaXplX3Bvc3RlcmlvcihwcmVkaWN0XzIpCmNhdCgiICBUcnVlIGFuc3dlcjogIiwgYmV0YTAgKyBhcy5tYXRyaXgobmV3X2RhdGEpICUqJSBiZXRhLCAiXG4iLCAKICAgIHNlcCA9ICIiKQpgYGAKClRoYXQgbG9va3MgcHJldHR5IGVuY291cmFnaW5nLiBJbiBib3RoIGNhc2VzLCBgeDFgLCBgeDJgLCBhbmQgYHgzYCBhcmUgcGlja2VkIG91dCBhcyBiZWluZyBpbXBvcnRhbnQsIGFuZCB0aGUgZXN0aW1hdGVkIGNvZWZmaWNpZW50cyBmb3IgZWFjaCBhcmUgcmVsYXRpdmVseSBjbG9zZSB0byB0aGUgdmFsdWVzIHVzZWQgdG8gZ2VuZXJhdGUgdGhlIGRhdGEuIFRoZSBldmlkZW5jZSBmb3IgYW4gYXNzb2NpYXRpb24gaW52b2x2aW5nIGFueSBvZiB0aGUgb3RoZXIgY292YXJpYXRlcyBpcyBwcmV0dHkgd2Vhay4KCiMjIFByb2plY3Rpb24gcHJlZGljdGlvbgoKTm93IGxldCdzIHRyeSBwcm9qZWN0aW9uIHByZWRpY3Rpb24gdG8gc2VlIGhvdyBtYW55IHZhcmlhYmxlcyB3ZSdkIGluY2x1ZGUgaW4gYSBtb2RlbCBhbmQgd2hpY2ggb25lcyB0aGV5IGFyZSBpZiB3ZSB3YW50ZWQgZXhwbGljaXRseSB0byBzaW1wbGlmeSB0aGUgbW9kZWwgcmF0aGVyIHRoYW4gZGlyZWN0bHkgaW50ZXJwcmV0aW5nIGNvZWZmaWNpZW50cyBvZiB0aGUgZnVsbCBtb2RlbC4gCgpgYGB7ciwgd2FybmluZyA9IEZBTFNFLCBtZXNzYWdlID0gRkFMU0V9CmxpYnJhcnkocnN0YW5hcm0pCmxpYnJhcnkocHJvanByZWQpCmxpYnJhcnkoYmF5ZXNwbG90KQoKb3B0aW9ucyhtYy5jb3JlcyA9IHBhcmFsbGVsOjpkZXRlY3RDb3JlcygpKQoKcHJvamVjdGlvbl9wcmVkaWN0aW9uIDwtIGZ1bmN0aW9uKGRhdCkgewogIG4gPC0gbnJvdyhkYXQpCiAgRCA8LSBuY29sKGRhdFssMjoxMF0pCiAgcDAgPC0gMwogIHRhdTAgPC0gcDAvKEQgLSBwMCkgKiAxL3NxcnQobikKICBwcmlvcl9jb2VmZiA8LSBocyhnbG9iYWxfc2NhbGUgPSB0YXUwLCBzbGFiX3NjYWxlID0gMSkKICBmaXQgPC0gc3Rhbl9nbG0oeSB+IHgxICsgeDIgKyB4MyArIHg0ICsgeDUgKyB4NiArIHg3ICsgeDggKyB4OSwgZGF0YSA9IGRhdCwgCiAgICAgICAgICAgICAgICAgIHByaW9yID0gcHJpb3JfY29lZmYsIHJlZnJlc2ggPSAwKQoKICB2cyA8LSB2YXJzZWwoZml0LCBtZXRob2QgPSAiZm9yd2FyZCIpCiAgY3ZzIDwtIGN2X3ZhcnNlbChmaXQsIG1ldGhvZCA9ICJmb3J3YXJkIiwgdmVyYm9zZSA9IEZBTFNFKQoKICBwcm9qX3NpemUgPC0gc3VnZ2VzdF9zaXplKGN2cykKICBwcm9qIDwtIHByb2plY3QodnMsIG52ID0gcHJval9zaXplLCBucyA9IDIwMDApCiAgbWNtY19pbnRlcnZhbHMoYXMubWF0cml4KHByb2opLAogICAgICAgICAgICAgICAgIHBhcnMgPSBjKCIoSW50ZXJjZXB0KSIsIG5hbWVzKHZzJHZpbmRbMTpwcm9qX3NpemVdKSwgInNpZ21hIikpCgogIHByZWQgPC0gcHJval9saW5wcmVkKHZzLCB4bmV3ID0gZGF0WywgMjoxMF0sIHluZXcgPSBkYXQkeSwgCiAgICAgICAgICAgICAgICAgICAgICAgbnYgPSBwcm9qX3NpemUsIGludGVncmF0ZWQgPSBUUlVFKQogIGZvci5wbG90IDwtIGRhdGEuZnJhbWUoT2JzZXJ2ZWQgPSBkYXQkeSwKICAgICAgICAgICAgICAgICAgICAgICAgIFByZWRpY3RlZCA9IHByZWQkcHJlZCkKICBwIDwtIGdncGxvdChmb3IucGxvdCwgYWVzKHggPSBPYnNlcnZlZCwgeSA9IFByZWRpY3RlZCkpICsgCiAgICBnZW9tX3BvaW50KCkgKyAKICAgIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpICsKICAgIGdlb21fYWJsaW5lKHNsb3BlID0gMSwgaW50ZXJjZXB0ID0gMCkgKwogICAgbGFicyh4ID0gIk9ic2VydmVkIiwgeSA9ICJQcmVkaWN0ZWQiKQogIHByaW50KHApCgogIHJldHVybihsaXN0KHByb2ogPSBwcm9qLCB2cyA9IHZzLCBjdnMgPSBjdnMsIHByb2pfc2l6ZSA9IHByb2pfc2l6ZSkpCn0KCnN1bW1hcml6ZV9wb3N0ZXJpb3IgPC0gZnVuY3Rpb24oeCwgY3JlZGlibGUgPSAwLjk1LCBkaWdpdHMgPSAzKSB7CiAgbG9fcCA8LSAoMS4wIC0gY3JlZGlibGUpLzIuMAogIGhpX3AgPC0gY3JlZGlibGUgKyBsb19wCiAgY2kgPC0gcXVhbnRpbGUoeCwgYyhsb19wLCBoaV9wKSkKICByZXN1bHRzIDwtIHNwcmludGYoIiUgNS4zZiAoJSA1LjNmLCAlIDUuM2YpXG4iLCBtZWFuKHgpLCBjaVsxXSwgY2lbMl0pCiAgY2F0KHJlc3VsdHMsIHNlcCA9ICIiKQp9CgpzdW1tYXJpemVfcmVzdWx0cyA8LSBmdW5jdGlvbih4LCBjcmVkaWJsZSA9IDAuOTUsIGRpZ2l0cyA9IDMpIHsKICB2YXJzIDwtIGNvbG5hbWVzKGFzLm1hdHJpeCh4KSkKICBmb3IgKGkgaW4gMTpsZW5ndGgodmFycykpIHsKICAgIGxhYmVsIDwtIHNwcmludGYoIiUxMXM6ICIsIHZhcnNbaV0pCiAgICBjYXQobGFiZWwsIHNlcCA9ICIiKQogICAgc3VtbWFyaXplX3Bvc3Rlcmlvcihhcy5tYXRyaXgoeClbLGldKQogIH0KfQoKcHJval8xIDwtIHByb2plY3Rpb25fcHJlZGljdGlvbihkYXRfMSkKcHJval8yIDwtIHByb2plY3Rpb25fcHJlZGljdGlvbihkYXRfMikKCnByaW50KHZhcnNlbF9wbG90KHByb2pfMSRjdnMsIHN0YXRzID0gYygiZWxwZCIsICJybXNlIiksIGRlbHRhcyA9IFRSVUUpKQpwcmludCh2YXJzZWxfcGxvdChwcm9qXzIkY3ZzLCBzdGF0cyA9IGMoImVscGQiLCAicm1zZSIpLCBkZWx0YXMgPSBUUlVFKSkKc3VtbWFyaXplX3Jlc3VsdHMocHJval8xJHByb2opCnN1bW1hcml6ZV9yZXN1bHRzKHByb2pfMiRwcm9qKQpgYGAKCkdvb2QgbmV3cyEgVXNpbmcgcHJvamVjdGlvbiBwcmVkaWN0aW9uIHRvIHNlbGVjdCB0aGUgImltcG9ydGFudCIgY292YXJpYXRlcyBpZGVudGlmaWVkIHRoZSBzYW1lIHRocmVlIGNvdmFyaWF0ZXMgb3VyIGV4YW1pbmF0aW9uIG9mIHRoZSBmdWxsIG1vZGVsIHN1Z2dlc3RlZCB0byB1cyBfKmFuZCpfIHRob3NlIGFyZSB0aGUgdmFyaWFibGVzIHdlIGtub3cgb3JlIGltcG9ydGFudC4KCkJvdHRvbSBsaW5lOiBJZiBhc3NvY2lhdGlvbnMgYW1vbmcgeW91ciBjb3ZhcmlhdGVzIGFyZSBtb2Rlc3QsIGlmIHlvdSBoYXZlIGEgcmVhc29uYWJseSBsYXJnZSBhbW91bnQgb2YgZGF0YSwgXyphbmQgaWYgeW91ciBjb3ZhcmlhdGVzIGluY2x1ZGUgYWxsIG9mIHRob3NlIHRoYXQgaGF2ZSBhIHJlYWwgcmVsYXRpb25zaGlwIHdpdGggdGhlIHJlc3BvbnNlIHZhcmlhYmxlKl8sIHRoZW4geW91IGhhdmUgYSBnb29kIGNoYW5jZSBvZiBkZXRlY3RpbmcgdGhlIHJlbGF0aW9uc2hpcHMgYW5kIG9mIGVzdGltYXRpbmcgdGhlaXIgbWFnbml0dWRlLiBCdXQgd2hhdCBoYXBwZW5zIGlmIHRoZSAidHJ1ZSIgcHJlZGljdG9ycyBhcmUgYW1vbmcgdGhlIGNvdmFyaWF0ZXMgaW5jbHVkZWQgaW4geW91ciBzdHVkeT8KCiMjIFdoYXQgaWYgdGhlICJ0cnVlIiBwcmVkaWN0b3JzIGFyZW4ndCBpbiB0aGUgZGF0YT8KCkluIG91ciBjYXNlIHRoZSB0cnVlIHByZWRpY3RvcnMgYXJlIGB4MWAsIGB4MmAsIGFuZCBgeDNgLiBXZSBrbm93IHRoYXQgdGhleSBpbmZsdWVuY2UgaW5mbHVlbmNlIGB5YCBhbmQgdGhhdCBubyBvdGhlciBjb3ZhcmlhdGVzIGRvLCBiZWNhdXNlIHdlIGdlbmVyYXRlZCB0aGUgZGF0YSB0aGF0IHdheS4gV2hhdCBoYXBwZW5zIGlmIHdlIGZpdCBhIG1vZGVsIHRoYXQgZXhjbHVkZXMgYHgxYCwgYHgyYCwgYW5kIGB4M2AgYW5kIGluY2x1ZGVzIG9ubHkgY292YXJpYXRlcyB0aGF0IHdlIGtub3cgZG9uJ3QgaGF2ZSBhICJyZWFsIiByZWxhdGlvbnNoaXAgdG8gYHlgPwoKYGBge3J9Cm4gPC0gbnJvdyhkYXRfMSkKRCA8LSBuY29sKGRhdF8xWyw1OjEwXSkKcDAgPC0gMwp0YXUwIDwtIHAwLyhEIC0gcDApICogMS9zcXJ0KG4pCnByaW9yX2NvZWZmIDwtIGhzKGdsb2JhbF9zY2FsZSA9IHRhdTAsIHNsYWJfc2NhbGUgPSAxKQoKZml0XzEgPC0gc3Rhbl9nbG0oeSB+IHg0ICsgeDUgKyB4NiArIHg3ICsgeDggKyB4OSwKICAgICAgICAgICAgICAgICAgZGF0YSA9IGRhdF8xLAogICAgICAgICAgICAgICAgICBwcmlvciA9IHByaW9yX2NvZWZmLAogICAgICAgICAgICAgICAgICByZWZyZXNoID0gMCwKICAgICAgICAgICAgICAgICAgYWRhcHRfZGVsdGEgPSAwLjk5OSkKIyMgTm90ZTogSSB1c2UgdGhlIHNhbWUgcHJpb3Igc2NhbGUgaGVyZSBiZWNhdXNlIGJvdGggZGF0YSBzZXRzIGhhdmUgdGhlIHNhbWUgbnVtYmVyIG9mIG9ic2VydmF0aW9ucwojIyBhbmQgdGhlIHNhbWUgZXhwZWN0ZWQgbnVtYmVyIG9mIGltcG9ydGFudCBjb3ZhcmlhdGVzCiMjCmZpdF8yIDwtIHN0YW5fZ2xtKHkgfiB4NCArIHg1ICsgeDYgKyB4NyArIHg4ICsgeDksCiAgICAgICAgICAgICAgICAgIGRhdGEgPSBkYXRfMiwKICAgICAgICAgICAgICAgICAgcHJpb3IgPSBwcmlvcl9jb2VmZiwKICAgICAgICAgICAgICAgICAgcmVmcmVzaCA9IDAsCiAgICAgICAgICAgICAgICAgIGFkYXB0X2RlbHRhID0gMC45OTkpCgpwcmVkaWN0XzEgPC0gcG9zdGVyaW9yX3ByZWRpY3QoZml0XzEpCnByZWRpY3RfMiA8LSBwb3N0ZXJpb3JfcHJlZGljdChmaXRfMikKCmZvci5wbG90IDwtIGRhdGEuZnJhbWUoT2JzZXJ2ZWQgPSBkYXRfMSR5LCBQcmVkaWN0ZWQgPSBhcHBseShwcmVkaWN0XzEsIDIsIG1lYW4pKQpwIDwtIGdncGxvdChmb3IucGxvdCwgYWVzKHggPSBPYnNlcnZlZCwgeSA9IFByZWRpY3RlZCkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fYWJsaW5lKHNsb3BlID0gMSwgaW50ZXJjZXB0ID0gMCkgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpICsKICBnZ3RpdGxlKCJEYXRhIHNldCAxIikKcHJpbnQocCkKCmZvci5wbG90IDwtIGRhdGEuZnJhbWUoT2JzZXJ2ZWQgPSBkYXRfMiR5LCBQcmVkaWN0ZWQgPSBhcHBseShwcmVkaWN0XzIsIDIsIG1lYW4pKQpwIDwtIGdncGxvdChmb3IucGxvdCwgYWVzKHggPSBPYnNlcnZlZCwgeSA9IFByZWRpY3RlZCkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fYWJsaW5lKHNsb3BlID0gMSwgaW50ZXJjZXB0ID0gMCkgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpICsKICBnZ3RpdGxlKCJEYXRhIHNldCAyIikKcHJpbnQocCkKCnAgPC0gcGxvdChmaXRfMSkgKyBnZ3RpdGxlKCJFc3RpbWF0ZXMgZnJvbSBkYXRhIHNldCAxIikKcHJpbnQocCkKcCA8LSBwbG90KGZpdF8yKSArIGdndGl0bGUoIkVzdGltYXRlcyBmcm9tIGRhdGEgc2V0IDIiKQpwcmludChwKQoKc3VtbWFyeShmaXRfMSwgZGlnaXRzID0gMykKc3VtbWFyeShmaXRfMiwgZGlnaXRzID0gMykKCm5ld19kYXRhIDwtIGRhdGEuZnJhbWUoeDQgPSA0LjAsIHg1ID0gNC4wLCB4NiA9IDQuMCwKICAgICAgICAgICAgICAgICAgICAgICB4NyA9IDQuMCwgeDggPSA0LjAsIHg5ID0gNC4wKQoKcHJlZGljdF8xIDwtIHBvc3Rlcmlvcl9wcmVkaWN0KGZpdF8xLCBuZXdfZGF0YSkKcHJlZGljdF8yIDwtIHBvc3Rlcmlvcl9wcmVkaWN0KGZpdF8yLCBuZXdfZGF0YSkKCnN1bW1hcml6ZV9wb3N0ZXJpb3IgPC0gZnVuY3Rpb24oeCwgY3JlZGlibGUgPSAwLjk1LCBkaWdpdHMgPSAzKSB7CiAgbG9fcCA8LSAoMS4wIC0gY3JlZGlibGUpLzIuMAogIGhpX3AgPC0gY3JlZGlibGUgKyBsb19wCiAgY2kgPC0gcXVhbnRpbGUoeCwgYyhsb19wLCBoaV9wKSkKICBjYXQocm91bmQobWVhbih4KSwgMyksICIgKCIsIHJvdW5kKGNpWzFdLCAzKSwgIiwiLCByb3VuZChjaVsyXSwgMyksICIpXG4iLCBzZXAgPSAiIikKfQoKc3VtbWFyaXplX3Bvc3RlcmlvcihwcmVkaWN0XzEpCnN1bW1hcml6ZV9wb3N0ZXJpb3IocHJlZGljdF8yKQpgYGAKClRoZSByZXN1bHRzIGZyb20gYW5hbHlzaXMgb2YgZGF0YSBzZXQgMSBzdWdnZXN0ICh3ZWFrbHkpIHRoYXQgYHg0YCBhbmQgYHg4YCBoYXZlIGEgbmVnYXRpdmUgYXNzb2NpYXRpb24gd2l0aCBgeWAsIGFuZCB0aGF0IGB4OWAgaGFzIGEgcG9zaXRpdmUgYXNzb2NpYXRpb25zLiBUaGF0J3MgcmlnaHQgaW4gdGhlIHNlbnNlIHRoYXQgKGEpIGB4NGAgYW5kIGB4OGAgYXJlIHBvc2l0aXZlbHkgYXNzb2NpYXRlZCB3aXRoIGB4MmAsIHdoaWNoIHdlIGtub3cgaGFzIGEgbmVnYXRpdmUgaW5mbHVlbmNlIG9uIGB5YCBhbmQgKGIpIGB4OWAgaXMgcG9zaXRpdmVseSBhc3NvY2lhdGVkIHdpdGggYm90aCBgeDFgIGFuZCBgeDNgLCB3aGljaCB3ZSBrbm93IGhhdmUgcG9zaXRpdmUgaW5mbHVlbmNlcyBvbiBgeWAuIEluIGRhdGEgc2V0IDIsIGB4NGAgYW5kIGB4NmAgaGF2ZSBhIHdlYWtseSBzdXBwb3J0ZWQgbmVnYXRpdmUgYXNzb2NpYXRpb24gd2l0aCBgeWAsIHdoaWxlIGB4NWAgaGFzIG1vcmUgc3Ryb25nbHkgc3VwcG9ydGVkIHBvc2l0aXZlIGFzc29jaWF0aW9uIHdpdGggYHlgLiBBZ2FpbiwgdGhvc2UgYXNzb2NpYXRpb25zIGFyZSBpbiB0aGUgZGlyZWN0aW9uIHdlJ2QgZXhwZWN0LiBCdXQgdGhlIGFzc29jaWF0aW9ucyBhcmUgd2Vhay4gT3VyIGV2aWRlbmNlIGZvciB0aGVtIGlzIHNoYWt5LgoKV2hhdCBoYXBwZW5zIGlmIHdlIHRyeSB0aGUgcHJvamVjdGlvbiBwcmVkaWN0aW9ucyBtZXRob2QgZm9yIHZhcmlhYmxlIHNlbGVjdGlvbj8KCmBgYHtyfQpwcm9qZWN0aW9uX3ByZWRpY3Rpb24gPC0gZnVuY3Rpb24oZGF0KSB7CiAgbiA8LSBucm93KGRhdCkKICBEIDwtIG5jb2woZGF0KSAtIDEKICBwMCA8LSAzCiAgdGF1MCA8LSBwMC8oRCAtIHAwKSAqIDEvc3FydChuKQogIHByaW9yX2NvZWZmIDwtIGhzKGdsb2JhbF9zY2FsZSA9IHRhdTAsIHNsYWJfc2NhbGUgPSAxKQogIGZpdCA8LSBzdGFuX2dsbSh5IH4geDQgKyB4NSArIHg2ICsgeDcgKyB4OCArIHg5LCBkYXRhID0gZGF0LCAKICAgICAgICAgICAgICAgICAgcHJpb3IgPSBwcmlvcl9jb2VmZiwgcmVmcmVzaCA9IDApCgogIHZzIDwtIHZhcnNlbChmaXQsIG1ldGhvZCA9ICJmb3J3YXJkIikKICBjdnMgPC0gY3ZfdmFyc2VsKGZpdCwgbWV0aG9kID0gImZvcndhcmQiLCB2ZXJib3NlID0gRkFMU0UpCgogIHByb2pfc2l6ZSA8LSBzdWdnZXN0X3NpemUoY3ZzKQogIHByb2ogPC0gcHJvamVjdCh2cywgbnYgPSBwcm9qX3NpemUsIG5zID0gMjAwMCkKICBtY21jX2ludGVydmFscyhhcy5tYXRyaXgocHJvaikpCgogIHByZWQgPC0gcHJval9saW5wcmVkKHZzLCB4bmV3ID0gZGF0WywgLTFdLCB5bmV3ID0gZGF0JHksIAogICAgICAgICAgICAgICAgICAgICAgIG52ID0gcHJval9zaXplLCBpbnRlZ3JhdGVkID0gVFJVRSkKICBmb3IucGxvdCA8LSBkYXRhLmZyYW1lKE9ic2VydmVkID0gZGF0JHksCiAgICAgICAgICAgICAgICAgICAgICAgICBQcmVkaWN0ZWQgPSBwcmVkJHByZWQpCiAgcCA8LSBnZ3Bsb3QoZm9yLnBsb3QsIGFlcyh4ID0gT2JzZXJ2ZWQsIHkgPSBQcmVkaWN0ZWQpKSArIAogICAgZ2VvbV9wb2ludCgpICsgCiAgICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iKSArCiAgICBnZW9tX2FibGluZShzbG9wZSA9IDEsIGludGVyY2VwdCA9IDApICsKICAgIGxhYnMoeCA9ICJPYnNlcnZlZCIsIHkgPSAiUHJlZGljdGVkIikKICBwcmludChwKQoKICByZXR1cm4obGlzdChwcm9qID0gcHJvaiwgdnMgPSB2cywgY3ZzID0gY3ZzLCBwcm9qX3NpemUgPSBwcm9qX3NpemUpKQp9Cgpwcm9qXzEgPC0gcHJvamVjdGlvbl9wcmVkaWN0aW9uKGRhdF8xWywgLWMoMiwzLDQsMTEpXSkKcHJval8yIDwtIHByb2plY3Rpb25fcHJlZGljdGlvbihkYXRfMlssIC1jKDIsMyw0LDExKV0pCgpwcmludCh2YXJzZWxfcGxvdChwcm9qXzEkY3ZzLCBzdGF0cyA9IGMoImVscGQiLCAicm1zZSIpLCBkZWx0YXMgPSBUUlVFKSkKcHJpbnQodmFyc2VsX3Bsb3QocHJval8yJGN2cywgc3RhdHMgPSBjKCJlbHBkIiwgInJtc2UiKSwgZGVsdGFzID0gVFJVRSkpCnN1bW1hcml6ZV9yZXN1bHRzKHByb2pfMSRwcm9qKQpzdW1tYXJpemVfcmVzdWx0cyhwcm9qXzIkcHJvaikKYGBgCgpJIHNhaWQgdGhhdCB0aGUgZGF0YSBzZXQgMSBzdWdnZXN0ZWQgc29tZSByZWxhdGlvbnNoaXBzIHdlYWtseSwgYW5kIHRoZSBwcm9qZWN0aW9uIHByZWRpY3Rpb25zIHNob3cgdGhhdCB0aGUgcmVsYXRpb25zaGlwcyBhcmUgc28gd2VhayB0aGF0IHdlIGRvbid0IGhhdmUgZ29vZCBldmlkZW5jZSBmb3IgaW5jbHVkaW5nIGFueSBvZiB0aGVtLiBUaGVyZSdzIHJlYXNvbmFibGUgZXZpZGVuY2UgZm9yIGluY2x1ZGluZyBvbmUgY292YXJpYXRlLCBgeDVgLCBpbiB0aGUgc2Vjb25kIGRhdGEgc2V0LCB3aGljaCBtZWFucyB3ZSBkb24ndCBoYXZlIGFueSBpbmRpY2F0aW9uIG9mIGFueSBuZWdhdGl2ZSBhc3NvY2lhdGlvbnMgd2l0aCBgeWAuW14yXQoKIyMgQ29uY2x1c2lvbnMKCklmIGFzc29jaWF0aW9ucyBhbW9uZyBvdXIgY292YXJpYXRlcyBhcmUgcmVsYXRpdmVseSB3ZWFrIGFuZCB3ZSBoYXZlIGEgcmVhc29uYWJsZSBhbW91bnQgb2YgZGF0YSwgaXQncyB0ZW1wdGluZyB0byB0aGluayB0aGF0IGFzc29jaWF0aW9ucyB3ZSBkZXRlY3QgaW4gb3VyIGRhdGEgYXJlIGxpa2VseSB0byBiZSByZWFsLiBJdCdzIGNlcnRhaW5seSB0aGUgY2FzZSB0aGF0IHdlJ3JlIGxlc3MgbGlrZWx5IHRvIGJlIG1pc2xlZCB3aGVuIGFzc29jaWF0aW9ucyBhcmUgd2FrZSwgYnV0IHRyeSB0aGlzIHNpbXBsZSB0aG91Z2h0IGV4cGVyaW1lbnQgdG8gcmVhbGl6ZSB3aHkgeW91IHNob3VsZCBiZSBjYXV0aW91cy4KCkltYWdpbmUgdGhhdCB5b3UncmUgaW50ZXJlc3RlZCBpbiBob3cgbGVhZiBtYXNzIHBlciBhcmVhIChMTUEpIHZhcmllcyBhY3Jvc3MgZW52aXJvbm1lbnRhbCBncmFkaWVudHMuIFN1cHBvc2UgeW91J3JlIHdvcmtpbmcgaW4gYSBtb3VudGFpbm91cyByZWdpb24gb2YgdGhlIE5vcnRoIFRlbXBlcmF0ZSBab25lIGFuZCB0aGF0IHlvdSBpbmNsdWRlIGVsZXZhdGlvbiBhcyBhIGNvdmFyaWF0ZSB3aXRob3V0IGluY2x1ZGluZyBhbnkgY292YXJpYXRlcyByZWxhdGVkIHRvIHRlbXBlcmF0dXJlIChlLmcuLCBtZWFuIGFubnVhbCB0ZW1wZXJhdHVyZSwgSmFudWFyeSBtaW5pbXVtIHRlbXBlcmF0dXJlLCBudW1iZXIgb2YgZnJvc3QtZnJlZSBkYXlzKS4gSSBzdXNwZWN0IHlvdTtsbCBmaW5kIGEgcG9zaXRpdmUgYXNzb2NpYXRpb24gYmV0d2VlbiBMTUEgYW5kIGVsZXZhdGlvbiwgYnV0IEkgc2VyaW91c2x5IGRvdWJ0IHRoYXQgaXQncyBlbGV2YXRpb24gX3BlciBzZV8gdGhhdCdzIGltcG9ydGFudC4gUmF0aGVyLCBJIHN1c3BlY3QgdGhhdCB0aGUgcmVsYXRpb25zaGlwLCBhc3N1bWluZyB0aGF0IGl0IGV4aXN0cywgd291bGQgcmVmbGVjdCB0aGUgdGlnaHQgYXNzb2NpYXRpb24gYmV0d2VlbiBkaWZmZXJlbnQgYXNwZWN0cyBvZiB0ZW1wZXJhdHVyZSBhbmQgZWxldmF0aW9uLiBFdmVuIGlmIHlvdSBpbmNsdWRlZCBhIHRlbXBlcmF0dXJlIGNvdmFyaWF0ZSwgZS5nLiwgZnJvc3QtZnJlZSBkYXlzLCBhbmQgeW91IG5vdyBmb3VuZCBhbiBhc3NvY2lhdGlvbiBiZXR3ZWVuIExNQSBhbmQgZnJvc3QtZnJlZSBkYXlzIHdoaWxlIHRoZSBlbGV2YXRpb24gYXNzb2NpYXRpb24gd2FzIG5vdyB3ZWFrIG9yIG5vbi1leGlzdGVudCwgSSdkIHN1Z2dlc3QgdGhhdCB5b3UgdGhpbmsgdmVyeSBjYXJlZnVsbHkgYmVmb3JlIGNvbmNsdWRpbmcgdGhhdCBmcm9zdC1mcmVlIGRheXMgXypyZWFsbHkqXyBtYXR0ZXJzLiBJdCBtaWdodCBiZSBvbmUgb2Ygc2V2ZXJhbCB0ZW1wZXJhdHVyZS1yZWxhdGVkIGNvdmFyaWF0ZXMgdGhhdCBtYXR0ZXIsIG9yIGl0IG1pZ2h0IHNpbXBseSBiZSB0aWdodGx5IGFzc29jaWF0ZWQgd2l0aCBzb21ldGhpbmcgZWxzZSB0aGF0IHJlYWxseSBtYXR0ZXJzLiBJbiB0aGUgYWJzZW5jZSBvZiBhIGNvbnRyb2xsZWQgZXhwZXJpbWVudCBvciBhIHZlcnkgY2FyZWZ1bCBjYXVzYWwgYW5hbHlzaXMse14zfSBhbGwgd2UgY2FuIHJlYWxseSBzYXkgZnJvbSB0aGUgZGF0YSBhbG9uZSBpcyB0aGF0IHRoZXJlJ3MgYW4gYXNzb2NpYXRpb24uICAKSW4gc2hvcnQsIEkgdXJnZSB5b3UgdG8gYmUgdmVyeSBjYXV0aW91cyBpbiBpbnRlcnByZXRpbmcgdGhlIHJlc3VsdHMgb2YgYSBtdWx0aXBsZSByZWdyZXNzaW9uLiBJZiB0aGUgYXNzb2NpYXRpb25zIHlvdSBpZGVudGlmeSBoYXZlIGJlZW4gZm91bmQgcmVwZWF0ZWRseSBpbiBvdGhlciBkYXRhIHNldHMgYW5kIHRoZXJlIGFyZSBnb29kLCBwcmluY2lwbGVkIHJlYXNvbnMgdG8gYmVsaWV2ZSB0aGF0IHRob3NlIGFzc29jaWF0aW9ucyBhcm9zZSB0aHJvdWdoIHRoZSBwcm9jZXNzIHlvdSBhcmUgdHJ5aW5nIHRvIGF0dHJpYnV0ZSB0aGVtIHRvLCBJJ2Qgc2F5IHlvdSdyZSBvbiByZWFzb25hYmx5IHNvbGlkIGdyb3VuZC4gVGhlIGxlc3Mgb2Z0ZW4gdGhlIGFzc29jaWF0aW9ucyB5b3Ugc2VlIGhhdmUgYmVlbiBzZWVuIGJlZm9yZSBhbmQgdGhlIGxlc3MgdGhleSBtaWdodCBoYXZlIGJlZW4gcHJlZGljdGVkIGZyb20gZmlyc3QgcHJpbmNpcGxlcywgdGhlIHNoYWtpZXIgaXMgdGhlIGdyb3VuZCBvbiB3aGljaCB5b3Ugc3RhbmQuIEJ1dCBiZWluZyBjYXV0aW91cyBpbiBpbnRlcnByZXRpbmcgdGhlIG1lYW5pbmcgb2YgYXNzb2NpYXRpb25zIHlvdSBmaW5kLCBkb2Vzbid0IG1lYW4gdGhhdCB0aG9zZSBhc3NvY2lhdGlvbnMgYXJlIHVuaW1wb3J0YW50LiBBc3NvY2lhdGlvbnMgdGhhdCBhcHBlYXIgaW1wb3J0YW50IGluIHlvdXIgYW5hbHlzaXMgYXJlIHdvcnRoIGZ1cnRoZXIgc3R1ZHksIGFuZCB0aGF0IGZ1cnRoZXIgc3R1ZHkgbWF5IGxlYWQgYm90aCB0byBmaW5kaW5nIHRoZW0gcmVwZWF0ZWRseSBhbmQgZGV2ZWxvcGluZyBhbiB1bmRlcnN0YW5kaW5nIG9mIG1lY2hhbmlzbXMgYW5kIHByb2Nlc3NlcyB0aGF0IGNvdWxkIHByb2R1Y2UgdGhlbS4KClteMV06IFlvdSBzaG91bGQgZmluZCBpdCBwcmV0dHkgZWFzeSB0byB0cnkgdGhlIExhc3NvIHVzaW5nIHRoZSBsZWFzdCBzcXVhcmVzIHJlZ3Jlc3Npb24gaW4gYGxtKClgIGlmIHlvdSdyZSBzbyBpbmNsaW5lZC4KCnteMl06IEkgaGF2ZXJuJ3QgdHJpZWQgcnVubmluZyB0aGUgbW9kZWxzIHRoYXQgZXhjbHVkZSBgeDFgLCBgeDJgLCBhbmQgYHgzYCB3aXRoIHRoZSBzdHJvbmcgYXNzb2NpYXRpb25zICgwLjggaW4gdGhlIGNhbGwgdG8gYHJtdm5vcm0oKWApLCBidXQgSSBlbmNvdXJhZ2UgeW91IHRvIHRyeSBpdCBvbiB5b3VyIG93bi4gSSBzdXNwZWN0IHlvdSdsbCBmaW5kIHNvbWUgcmVhc29uYWJseSB3ZWxsIHN1cHBvcnRlZCBhc3NvY2lhdGlvbnMgaWYgeW91IGRvLg==