En el primer número de la novísima revista Spanish Journal of Statistics aparece un artículo con un título tentador: Recovering income distributions from aggregated data via micro-simulations.
Es decir, un artículo que nos puede permitir, por ejemplo, muestrear lo que la AEAT llama rendimientos a partir de lo que publica (aquí):
Uno de los métodos de los que sostienen el ignominioso a mí me funciona está basado en el modelo
[El título de esta entrada tiene un + delante porque ya escribí sobre el asunto tiempo atrás.]
Con la excusa de la reciente publicación del paquete ivreg (para el ajuste de modelos con variables instrumentales, por si el contexto no lo hace evidente), he mirado a ver quién estaba construyendo y ajustando modelos generativos menos triviales que los míos (véase el enlace anterior) para que quede más claro de qué va la cosa. Porque la explicación típica, que adopta formas no muy distintas de
Imputación (que es algo en lo que muy a regañadientes estoy trabajando estos días).
Si de verdad tienes que imputar datos en una tabla (y solo en ese caso), solo hay un criterio: construye un modelo para predecir los valores faltantes en función del resto y reemplaza el NA por la su predicción.
El modelo puede ser tan tonto como
lm(my_col ~ 1, na.rm = T)
que resulta en la popular estrategia de reemplazar los NAs por la media del resto de las observaciones. Cambiando lm por otras cosas funciones más molonas y la fórmula por otras más complejas en que intervengan otras columnas se obtienen métodos más potentes. Se pueden usar GAMs (como en mtsdi) o random forests (como en missForest), pero la idea está clara. Es solo la naturaleza del problema la que nos invita a decantarnos por una u otra opción.
Tres situaciones. La primera:
n <- 20
y <- 15
test <- prop.test(y, n, p = .5)
test$p.value
# [1] 0.04417134
test$conf.int
# 0.5058845 0.9040674La segunda:
n <- 200
y <- 115
test <- prop.test(y, n, p = 0.5)
test$p.value
#[1] 0.04030497
test$conf.int
# 0.5032062 0.6438648Y la tercera:
n <- 2000
y <- 1046
test <- prop.test(y, n, p = 0.5)
test$p.value
#[1] 0.0418688
test$conf.int
# 0.5008370 0.5450738En resumen:
La pregunta: ¿qué circunstancia es más favorable? Una respuesta, aquí.
He intentado replicar los resultados de la entrada de ayer con GAM (vía mgcv) así (véase el enlace anterior para la definición de los datos):
library(mgcv)
modelo_gam <- gam(
y ~ x + s(id, bs = "re"),
data = datos,
method = "REML",
family = "poisson")Y nada:
Error in gam(y ~ x + s(id, bs = "re"), data = datos, method = "REML", : Model has more coefficients than data
Sí, ya sé que tengo más variables que observaciones. Pero, ¿no es para eso que estoy usando efectos aleatorios?
Construyo unos datos (artificiales, para conocer la verdad):
n <- 10000
x1 <- rnorm(n)
x2 <- rnorm(n)
probs <- -2 + x1 + x2
probs <- 1 / (1 + exp(-probs))
y <- sapply(probs, function(p) rbinom(1, 1, p))
dat <- data.frame(y = y, x1 = x1, x2 = x2)Construyo un modelo de clasificación (logístico, que hoy no hace falta inventar, aunque podría ser cualquier otro):
summary(glm(y ~ x1 + x2, data = dat, family = binomial))
#Call:
#glm(formula = y ~ x1 + x2, family = binomial, data = dat)
#
#Deviance Residuals:
# Min 1Q Median 3Q Max
#-2.2547 -0.5967 -0.3632 -0.1753 3.3528
#
#Coefficients:
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -2.05753 0.03812 -53.97 <2e-16 ***
#x1 1.01918 0.03386 30.10 <2e-16 ***
#x2 1.00629 0.03405 29.55 <2e-16 ***
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#(Dispersion parameter for binomial family taken to be 1)
#
# Null deviance: 9485.2 on 9999 degrees of freedom
#Residual deviance: 7373.4 on 9997 degrees of freedom
#AIC: 7379.4
#
#Number of Fisher Scoring iterations: 5Correcto.
Abundo sobre mi entrada del otro día. Usando números aleatorios hirsutos,
n <- 200
x <- runif(n)
plot(cumsum(x - .5), type = "l")produce
mientras que
library(randtoolbox)
s <- sobol(n, 1, scrambling = 3)
plot(cumsum(s - .5), type = "l")genera
que tiene un cariz totalmente distinto.
Contemplando y comparando
y
se me han venido a la mente los adjetivos hirsuto y pocholo para calificar las respectivas formas de aleatoriedad que representan. La primera es el resultado del habitual
n <- 200
x <- runif(n)
y <- runif(n)
plot(x, y, pch = 16)mientras que la segunda exige el más sofisticado
library(randtoolbox)
s <- sobol(n, 2, scrambling = 3)
x <- s[,1]
y <- s[,2]
plot(x, y, pch = 16)Se ve que Sobol quería rellenar más armoniosamente el espacio. Me temo que, al hablar de aleatoriedad, muchos de nosotros también (p.e., esto).
Tengo una entrada perpetuamente pendiente que se pospone, entre otras cosas, porque aún no he encontrado una manera satisfactoria para muestrear histogramas. Una de las vías sería dar con (y ajustar) una distribución subyacente que generase unos histogramas similares.
Hoy voy a contar un ejemplo de cómo puede fallar tal estrategia.
Por un lado he bajado datos de la distribución de renta en España del INE:
Por otro, me he dejado convencer temporalmente de que la distribución de Burr podría ser conveniente para modelar la distribución de ingresos de los hogares (Wikipedia dixit!).
Uno de los proyectos en los que estoy trabajando últimamente está relacionado con un problema de optimización no lineal: tengo un modelo (o una familia de modelos) no lineales con una serie de parámetros, unos datos y se trata de lo que no mercería más explicación: encontrar los que minimizan cierta función de error.
Tengo implementadas dos vías:
nls porque esa función me queda muy corta).Ambas tienen sus ventajas y desventajas. La una es rápida y la otra no; la una me da poca información sobre los parámetros y la otra, mucha; una me permite introducir mucha información a priori y la otra casi nada, etc.
