Regresión

Análisis de la discontinuidad + polinomios de grado alto = ...

Una técnica que, al parecer, es muy del gusto de los economistas es lo del análisis de la discontinuidad. Es como todo lo que tiene que ver con causalImpact pero usando técnicas setenteras (regresiones independientes a ambos lados del punto de corte).

Si a eso le sumas que las regresiones pueden ser polinómicas con polinomios de alto grado… pasan dos cosas:

  • Tienes una probabilidad alta de obtener un resultado significativo, i.e., publicable.
  • Pero que se deba solo al ruido producido por el método (corte discreto, inestabilidad polinómica, etc.).

Es decir, la habitual chocolatada que algunos llaman ciencia (cierto, algunos dirán que mala ciencia, pero que, ¡ah!, nos cobran al mismo precio que la buena).

¿Quitar variables no significativas?

Contexto: modelos de regresión con de varias a muchas variables. Muy particularmente cuando interesa la predicción.

Pseudoproblema: ¿quitamos las variables no significativas?

Los manualitos (muy queridos de enseñantes, porque les dan reglas sencillitas; muy queridos también de los aprendientes, por el mismo motivo) rezan que sí. Se quitan y a otra cosa.

La regla adulta es:

  • Si el coeficiente es grande y tiene el signo correcto, ¡enhorabuena!
  • Si el coeficiente es pequeño, la variable no hace ni bien ni mal. Y hay más motivos para dejarla que para quitarla.
  • Pero si el coeficiente es grande y el signo es contrario a lo que cabría esperar (p.e., a más gripe menos fallecidos, a más capacidad económica menos compra media, etc.), ¡ah!, toca volver a replantear el modelo seriamente.

Nota: en lo anterior no he usado la palabra significativo. Si alguien quiere traducir grande y pequeño en términos de la ocurrencia de hace ochenta años de un inglés que sostenía que el tabaco era sano, allá él.

Tres de seis consejos para mejorar las regresiones

Por si alguien se lo perdió, están aquí. De los seis, mencionaré tres que me están resultando muy útiles en un proyecto actual.

De todos ellos, el que más a rajatabla sigo es el primero: ajustar muchos modelos. Pudiera parecer trampa: buscar y rebuscar por si sale algo. Sin embargo, es una técnica que plantearse como una manera de familiarizarse y aprender la estructura de los datos. Los modelos (explicativos, como los que justifican esta entrada) no dejan de ser resúmenes de conjuntos de datos y no es sino ajustando diversos modelos que uno aprende si, por ejemplo, un coeficiente varía por año o provincia.

Modelos directos, inversos y en los que tanto da

Continúo con esto que concluí con una discusión que me negué a resolver sobre la geometría de los errores.

Que es la manera de entender que los problemas directos e inversos no son exactamente el mismo. Digamos que no es una medida invariante frente a reflexiones del plano (que es lo que hacemos realmente al considerar el modelo inverso).

¿Pero y si medimos la distancia (ortogonal) entre los puntos $latex (x,y)$ y la curva $latex y = f(x)$ (o, equivalentemente, $latex x = f^{-1}(x)$)? Entonces daría (o debería dar) lo mismo.

Polinomios monótonos

R

Recibí un mensaje el otro día sobre polinomios monótonos. Mejor dicho, sobre el ajuste de datos usando polinomios monótonos. Frente a un modelo del tipo y ~ x (x e y reales) donde la relación entre las dos variables es

  • manifiestamente no lineal y
  • necesariamente monótina, p.e., creciente (por consideraciones previas),

cabe considerar ajustar un polinomio monótono, i.e., realizar una regresión polinómica con la restricción adicional de que el polinomio de ajuste resultante sea monótono.

Un problema inverso de regresión

He estado pensando qué tipo de ejercicios de estadística (y modelos estadísticos) plantear a mis alumnos del máster de data science de la UTAD.

Así que les he dado unos datos, los X, relativamente grandes (y sin problemas de colinealidad y similares) y les voy a pedir que me construyan la y de manera que los coeficientes obtenidos sean, aproximadamente, iguales a unos dados. A ver qué tal se les da.

El problema de la estimación inversa

Supongamos que tenemos unos niños de los que sabemos las edades $latex x_i$ y las alturas $latex y_i$. Supongamos además que podemos estimar las segundas en función de las primeras con un modelo lineal clásico

$$ y_i \sim N(a_0 + a_1 x_1, \sigma).$$

Este modelo nos permite, dada una edad, estimar la altura y los correspondientes intervalos de confianza. Pero, dada una altura, ¿qué nos dice de la edad? Este es el problema conocido como de la estimación inversa.

Bajo el capó del particionamiento recursivo basado en modelos

Una de las mayores contrariedades de estar sentado cerca de alguien que es más matemático que un servidor (de Vds., no de silicio) es que oye siempre preguntar por qué. Una letanía de preguntas me condujo a leer papelotes que ahora resumo.

Primero, unos datos:

set.seed(1234)

n <- 100

x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)

y <- 0.3 + 0.2 * x1 + 0.5 * (x2 > 0) + 0.2 * rnorm(n)

Luego, un modelo:

modelo <- lm(y ~ x1)
summary(modelo)

# Call:
#   lm(formula = y ~ x1)
#
# Residuals:
#   Min      1Q  Median      3Q     Max
# -0.9403 -0.2621  0.0420  0.2299  0.6877
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.55632    0.03364  16.538  < 2e-16 ***
#   x1           0.21876    0.03325   6.579 2.34e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3323 on 98 degrees of freedom
# Multiple R-squared:  0.3063,  Adjusted R-squared:  0.2992
# F-statistic: 43.28 on 1 and 98 DF,  p-value: 2.341e-09

Pocos que no entiendan cómo se han generado los datos advertirían lo malo de su especificación: hemos omitido una variable explicativa cuyo efecto ha ido a incrementar el error de manera que los tests habituales de bondad de ajuste no advierten.

Medianas ponderadas en R

La mediana de 1:3 es 2. Pero puede ser que queramos dar a 1:3 los pesos 2, 1, 2. En ese caso, el cálculo de la mediana sigue siendo sencillo (y sigue siendo 2). Pero la situación puede complicarse más.

Mientras los pesos sean enteros, todavía pueden usarse trucos:

x <- 1:3
pesos <- c(2,1,2)
median(rep(x, times = pesos ))

¿Pero qué hacemos cuando hay pesos fraccionarios? Bueno, en realidad, podemos ordenar:

n <- 1000

x <- runif(n)
pesos <- runif(n)
o <- order(x)
x.o <- x[o]
pesos.o <- pesos[o]
x.o[min(which(cumsum(pesos.o) > .5 * sum(pesos.o)))]

Pero me parece más limpio usar el paquete quantreg:

Regresión por cuantiles en R y SAS

Hace un tiempo, con la aburridora perspectiva de un largo viaje en metro hasta mi casa ensombreciendo mi futuro más inminente, decidí regalarme algún tipo de amena lectura. A tal fin, imprimí un articulillo que, bajo la perspectiva de SAS, me introducía a una técnica que se vino a mí como por azar. O, bajo otro punto de vista, una técnica que, también por azar, había esquivado hasta tal fecha un encontronazo con mi husmeadora curiosidad.