Lm

El otro día me pasaron unos datos artificiales para poder probar el ajuste de cierto tipo de modelos. El autor de la simulación construyó tres conjuntos de pares (x,y) y luego los agregó (media de los y agrupando por x) antes de proporcionármelos. ¿Tiene sentido agregar antes de modelar? Incluso sin entrar en el problema del potencial número desigual de observaciones por punto (datos desbalanceados) o las heterogeneidades entre las distintas iteraciones (que nos llevaría al mundo de los modelos mixtos). ...

Ayer se leía en Twitter que "La regresión multinivel debería ser la forma predeterminada de hacer regresión" — Jose Luis Cañadas (@joscani) April 11, 2020 Cabe preguntarse qué pasa si se analizan los mismos datos usando ambas técnicas. Obviamente, hay muchos tipos de datos y supongo que los resultados variarán según qué variante se utilice. Aquí voy a centrarme en unos donde hay medidas repetidas de un factor aleatorio. También voy a situarme en un contexto académico, en el que interesan más las estimaciones de los efectos fijos, que en uno más próximo a mi mundo, la consultoría, donde son más relevantes las estimaciones regularizadas de los efectos aleatorios. ...

Estos días me han preguntado sobre un modelo lineal tal que $y \sim x_1 + \dots$ donde el coeficiente de $x_1$ no se entiende si no es igual a 1. Es como si los datos se creasen de la forma n <- 100 x1 <- rnorm(n) x2 <- rnorm(n) y <- x1 + rnorm(n, .1) + .02 * x2 y se conociese el coeficiente de $x_1$ y no el de $x_2$. Entonces no tiene sentido plantear el modelo lm(y ~ x1 + x2) sino más bien modelo <- lm(y ~ offset(x1) + x2) que hace lo que uno espera. Lo cual se puede comprobar, por ejemplo, comparando ...

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")])

Supongamos que queremos construir un modelo para predecir quién ganará un determinado partido de baloncesto basándonos en datos diversos. Y en un histórico, por supuesto. Podemos utilizar una regresión logística así: set.seed(1234) my.coefs <- -2:2 n <- 200 train.n <- floor(2*n/3) test.error.glm <- function(){ X <- matrix(rnorm(n*5), n, 5) Y <- (0.2 + X %*% my.coefs + rnorm(n)) > 0 train <- sample(1:n, train.n) X <- as.data.frame(X) X$Y <- Y mod.glm <- glm(Y ~ ., data = X[train,], family = binomial) glm.pred <- predict(mod.glm, X[-train,], type = "response") error <- length(glm.pred) - sum(diag(table(glm.pred > 0.5, Y[-train,]))) } errores.glm <- replicate(1000, test.error.glm()) El código anterior hace lo siguiente: Crea las variables aleatorias X (unos predictores) e Y (el resultado de los partidos). Ajusta un modelo logístico a un subconjunto de los datos. Predice sobre el complementario de dichos datos, el conjunto de prueba. Mide el error cometido. Itera el proceso anterior y guarda los errores de clasificación cometidos. Nótese que la variable objetivo es binaria por construcción. ...