Estadística

Sobre la peculiarísima implementación del modelo lineal en (pseudo-)scikit-learn

Si ejecutas

import numpy as np
from sklearn.linear_model import LinearRegression

n = 1000
X = np.random.rand(n, 2)

Y = np.dot(X, np.array([1, 2])) + 1 + np.random.randn(n) / 2
reg = LinearRegression().fit(X, Y)

reg.intercept_
reg.coef_

se obtiene más o menos lo esperado. Pero si añades una columna linealmente dependiente,

X = np.column_stack((X, 1 * X[:,1]))

ocurren cosas de la más calamitosa especie:

Y = np.dot(X, np.array([1, 2, 1])) + 1 + np.random.randn(n) / 2
reg = LinearRegression().fit(X, Y)
reg.coef_
# array([ 9.89633991e-01, -1.63740303e+14,  1.63740303e+14])

Comentarios:

  • Diríase que la implementación del modelo lineal en scikit-learn no es la que se estudia por doquier (la prima, la inversa, etc.); sospecho que convierte el error cuadrático en una función que depende de los coeficientes y se la pasan a un optimizador (más o menos) genérico.
  • Supongo que la implementación actual pasa todos las pruebas unitarias.
  • Y sospecho, además, que las pruebas unitarias no las ha planteado un estadístico.
  • Así que tal vez los de scikit-learn no saben que tienen problemas de colinealidad; y si alguien se lo ha comentado, igual no han comprendido el issue.
  • Para la predicción, el problema arriba apuntado no es tal. Aun con coeficientes desaforados y siempre que no haya problemas de precisión numérica, tanto da hacer las cosas como todo el mundo o implementando ocurrencias como la anterior.
  • Pero para todo lo demás (p.e., impacto de variables, etc.), la implementación es de traca y no vale para maldita de Dios la cosa.
  • Aunque a estas alturas de siglo, ¿quién en su sano juicio usa el modelo lineal básico?
  • Además, en la práctica, no hay problemas de colinealidad (ni aproximada o, mucho menos, exacta). Hummm…
  • ¿O sí? Mi colega Cañadas ha escrito una entrada en su blog sobre la codificación de variables categóricas donde denuncia cómo en Python las funciones más habituales crean por defecto columnas linealmente dependientes por defecto (por no omitir el primer nivel). O sea, en Python, si no andas con cuidado, acabas con la suela llena de kk de perro.

Estacionalidad semanal de la mortalidad

Continúo con esto.

(El pico de primeros de agosto corresponde a la ola de calor).

Resumen:

  • Acusada variación intrasemanal (con un intevalo de variación del 5%, que es mucho).
  • Los intervalos de confianza de los nuevos modelos contienen mucho mejor las observaciones reales; muchos picos observados dejan de parecer anómalos para quedar dentro de las franjas de variación esperadas.

Mortalidad y domingos

Es sabido que nacen menos niños en domingo, efecto, parece de la planificación de partos. Tengo cierto indicio, que voy a ver cuándo (y si) puedo corroborar, de que también hay menos fallecimientos. Mírese

que es una gráfica de mortalidad diaria en España y en la que, contumazmente, los días 2 de junio (domingo), 9 de junio (domingo) y 16 de junio (domingo), las observaciones quieren salirse de las bandas que para que no trotaran libérrimamente por el plano cartesiano construyó el bueno de Monsieur Poisson.

Modelización de retrasos: una aplicación del análisis de supervivencia

En vigilancia epidemiológica contamos eventos (p.e., muertes o casos de determinadas enfermedades). Lo que pasa es que el caso ocurrido en el día 0 puede notificarse con un retraso de 1, 2, 3… o incluso más días. En algunas aplicaciones, incluso semanas.

¿Cómo estimar el número de casos ocurridos el día 0 el día, p.e., 5?

Se puede aplicar el análisis de la supervivencia donde el evento muerte se reinterpreta como notificación. El el día 0 todos los sujetos están vivos y, poco a poco, van cayendo. Como en los consabidos modelos/gráficos de Kaplan-Meier,

Bayes no había previsto esto

Muestreo. Se trata de seleccionar unas unidades experimentales (proceso caro) y tratar de estimar una proporción (p.e.) en la población total.

Existen técnicas para estimar el valor N mínimo para garantizar cierto margen de error. Pero dichas técnicas requieren conocer (algo d-) el resultado del experimento para estimar N (p.e. una estimación de la proporción que cabe esperar).

Circulus in demonstrando.

Bayes. Ve examinando unidades y actualiza tus intervalos de credibilidad hasta que tengan la anchura solicitada.

Causalidad. Atribución. Madrid Central.

Si hay algo inaprensible, es la causalidad. No la que entiende Maripili, claro, sino esta. Pero vivimos en tiempos de tremendamente polémicas y presuntamente potentísimas y causas eficientes. Verbigracia, la desigualdad… y Madrid Central:

Argumentas en términos causales cuando esperas que te lea Maripili. Entre gente seria solemos hablar más bien de atribución. Lo de la atribución consiste en tratar de repartir un efecto entre posibles causas potenciales. Como típicamente no hay criterio definitivo, en la práctica funciona así:

Modelos GARCH (o: no me cuentes tu vida, dame el pxxx modelo generativo y ya)

Los modelos GARCH son otra de esas cosas de las que oyes hablar y como nunca se convierten en problemas de los de carne en el asador, preocupan poco y ocupan menos (más allá de que sabes que se trata de modelos similares a los de series temporales de toda la vida donde la varianza varía de cierta forma a lo largo del tiempo). Pero comienzas a leer cosas como esta y no te enteras de nada: solo hay letras y llamadas a funciones oscuras y oscurantistas.

¿Y si quitamos el puntico de arriba a la izquierda?

Esta entrada es una illustración de otra de no hace mucho, Análisis de la discontinuidad + polinomios de grado alto = … Mirad:

Se ha hecho un análisis de la discontinuidad usando parábolas a ambos lados del punto de corte. Y la discontinuidad no es pequeña. Pero me juego un buen cacho de lo que quede de mi reputación a que mucho de ella la explica el puntico de arriba a la izquierda.