Colinealidad y posterioris
En esta entrada voy a crear un conjunto de datos donde dos variables tienen una correlación muy alta, ajustar un modelo de regresión y obtener la siguiente representación de la distribución a posteriori de los coeficientes,
donde se aprecia el efecto de la correlación entre x1 y x2.
El código,
library(mvtnorm)
library(rstan)
library(psych)
n <- 100
corr_coef <- .9
x <- rmvnorm(n, c(0, 0),
sigma = matrix(c(1, corr_coef, corr_coef, 1), 2, 2))
plot(x)
x1 <- x[,1]
x2 <- x[,2]
x3 <- runif(n) - 0.5
y <- 1 + .4 * x1 - .2 * x2 + .1 * x3 + rnorm(n, 0, .1)
summary(lm(y ~ x1 + x2 + x3))
stan_code <- "
data {
int N;
vector[N] y;
vector[N] x1;
vector[N] x2;
vector[N] x3;
}
parameters {
real a;
real a1;
real a2;
real a3;
real sigma;
}
model {
a ~ cauchy(0,10);
a1 ~ cauchy(0,2.5);
a2 ~ cauchy(0,2.5);
a3 ~ cauchy(0,2.5);
y ~ normal(a + a1 * x1 + a2 * x2 + a3 * x3, sigma);
}"
datos_stan <- list(
N = n,
y = y,
x1 = x1,
x2 = x2,
x3 = x3
)
fit2 <- stan(model_code = stan_code,
data = datos_stan,
iter = 10000, warmup = 2000,
chains = 2, thin = 4)
res <- as.data.frame(fit2)
pairs.panels(res[, c("a", "a1", "a2", "a3", "sigma")])
