Caminos aleatorios

Las componentes lineales de una regla de clasificación cuadrática

Publicado el 04/09/2018 por JRB

¿Cuándo una regla lineal es equivalente a una cuadrática?

Las reglas de clasificación cuadráticas son muy flexibles y permiten distinguir poblaciones que difieren no solo en las medias sino también en la estructura de covarianzas. De hecho, si queremos clasificar entre dos poblaciones normales con matrices de covarianzas distintas, la regla óptima (o regla Bayes) es cuadrática.

Pero incluso en situaciones en las que la regla de clasificación óptima es cuadrática, es posible que no se pierda nada (en términos de error de clasificación) al sustituirla por una regla lineal más simple. Consideremos el ejemplo siguiente:

En este ejemplo tenemos dos muestras de entrenamiento procedentes de poblaciones normales con matrices de covarianzas diferentes. La regla de clasificación óptima es cuadrática y corresponde a la línea discontinua del gráfico. Está claro que si la cambiamos por una regla lineal no perdemos nada en error de clasificación y que la complejidad que introducen los términos cuadráticos es superflua en este caso. Podemos usar una regla lineal, mucho más simple y fácil de interpretar.

Aproximación de reglas cuadráticas mediante productos de hiperplanos

En un artículo que se publicará proximamente hemos considerado el problema de cómo aproximar una regla cuadrática mediante productos de hiperplanos, con la idea de identificar situaciones como la del ejemplo anterior. La ventaja es que las aproximaciones solo dependen de combinaciones lineales de las variables, lo que permite definir reglas de clasificación más fáciles de interpretar. Además, es posible saber cuándo el producto es parecido a la función cuadrática original. En el ejemplo anterior el resultado de aplicar nuestro método a los datos disponibles es el siguiente:

optima-prod

En este caso solo uno de los dos hiperplanos del producto aproximante (las líneas continuas) es necesario para obtener una regla lineal con propiedades equivalentes a la regla cuadrática original. En otros casos pueden ser necesarios los dos hiperplanos o incluso otros productos de hiperplanos complementarios. Consideremos el ejemplo siguiente:

El gráfico de arriba a la derecha corresponde a la aproximación a la regla optima basada en un único producto de hiperplanos. Su error de clasificación va a ser bastante mejor que el obtenido mediante una regla lineal (por ejemplo Fisher) al precio de considerar dos funciones lineales en lugar de una. Incluso nos podemos acercar más a la regla cuadrática óptima si consideramos un segundo producto de hiperplanos (gráfico de abajo a la derecha). Para ello tendríamos que tener en cuenta cuatro combinaciones lineales de las variables originales.

La idea es que combinando productos de hiperplanos se obtiene una regla de clasificación lineal a trozos que combina las ventajas de las reglas cuadráticas (error de clasificación) y de las lineales (facilidad de interpretación). Los hiperplanos que intervienen en los productos constituyen la estructura lineal de la regla cuadrática.

El método se aplica a cualquier regla cuadrática, no necesariamente a las óptimas. Usando técnicas de regularización se han propuesto muchas reglas cuadráticas incluso en casos de alta dimensión en los que no es posible estimar las matrices de covarianzas. Un artículo reciente en el que se mencionan muchas de ellas es Qin (2018).

¿Es muy difícil determinar la estructura lineal de una regla cuadrática?

Si la regla cuadrática consiste en clasificar $x$ en una de las poblaciones si y solo si $f(x)>0$ , donde $f(x)=x'Ax+2a'x+c$ , entonces para calcular los productos de hiperplanos basta diagonalizar la matriz $B=aa'-cA$ , que depende de manera muy sencilla de la matriz $A$ , el vector $a$ y el escalar $c$ . Y esto para cualquier dimensión. De manera similar a como ocurre en el análisis de componentes principales, si muchos de los autovalores de $B$ son pequeños es posible aproximar bien la regla cuadrática con un número pequeño de productos de hiperplanos. De esta manera también se puede considerar nuestro método como una forma de reducir la dimensión del problema.

Si el lector tiene curiosidad por saber por qué la matriz $B$ es justo la que hace falta diagonalizar para encontrar los productos de hiperplanos aproximantes tendrá que leer nuestro artículo 🙂

Referencia

Berrendero, J.R. and Cárcamo, J. (To appear). Linear components of quadratic classifiers. Advances in Data Analysis and Classification. Preprint | doi

Anuncios

Publicado en estadística | Etiquetado clasificación supervisada, investigación | Deja un comentario

El Grand Slam de las revistas estadísticas

Publicado el 16/08/2018 por JRB

Según wikipedia el término Grand Slam procede de juegos de cartas de la familia del whist y en este contexto ya se usaba en el siglo XVII. Desde 1930 se usa en golf para referirse al hecho de ganar los cuatro principales torneos en la misma temporada. Al parecer el periodista Alan Gould fue el primero en usar este término referido al tenis en un artículo de 1933 en el periódico The Reading Eagle:

GrandSlam

Pero, ¿cuál sería el Grand Slam de las revistas de estadística? Si a los investigadores en estadística les preguntaran cuáles son las revistas más importantes, ¿cuál sería su respuesta?

Una encuesta

Recientemente se han publicado los resultados de 2018 del ranking de Shanghái. Junto con la lista de mejores universidades (ya suficientemente comentada en diversos foros) se publican también los resultados de una encuesta en la que se pregunta a científicos de distintas áreas, entre otras cosas, por las revistas de su especialidad que consideran más prestigiosas. Me ha parecido interesante mirar los resultados porque reflejan la importancia de las revistas de una forma más cualitativa que los controvertidos índices de impacto. En el área de estadística son cuatro las revistas que aparecen, lo que podríamos llamar el Grand Slam de las revistas estadísticas. Estos son los resultados:

shanghai

Si comparamos con el ranking de revistas en función del índice de impacto hay ciertas discrepancias. Annals, JASA y Royal B son (todavía) Q1, aunque no son ciertamente las tres revistas con mayor índice de impacto. Por otra parte, hay muchos años en los que Biometrika no aparece como Q1, con lo que ello implica en procesos de acreditación, solicitudes de sexenios, etc.

Biometrika

Y sin embargo, a quién no le gustaría publicar en una revista fundada por Pearson, Weldon y Galton a principios del siglo XX y en la que han aparecido contribuciones fundamentales en la historia de la metodología estadística. Así cuenta David Cox los comienzos de Biometrika en un artículo publicado para celebrar su centenario:

biometrika

En 1908 Gosset publicó en Biometrika su artículo sobre la distribución t de Student. Más de un siglo después miles de científicos siguen usando a diario esta distribución. Me cuesta imaginar un artículo que contribuya más al impacto de una revista científica.

Publicado en estadística | Etiquetado historia, profesión | Deja un comentario

¿Con qué probabilidad existe el estimador de máxima verosimilitud en regresión logística?

Publicado el 21/06/2018 por JRB

Un nuevo resultado sobre el estimador de máxima verosimilitud en regresión logística

En una entrada anterior de este blog se comentaba un conocido resultado sobre la existencia del estimador de máxima verosimilitud (EMV) en el modelo de regresión logística. El EMV existe si y solo si las nubes de puntos de las variables regresoras correspondientes a los dos grupos de la variable respuesta ( $Y=0$ e $Y=1$ ) no se pueden separar linealmente mediante un hiperplano (la separación no tiene por qué ser estricta).

En un artículo publicado recientemente en arxiv, Emmanuel Candes y Pragya Sur han profundizado en esta situación cuando las variables regresoras tienen distribución normal. Si $n$ es el tamaño muestral y $p$ es el número de variables regresoras, en el artículo se investiga con qué probabilidad podemos esperar que exista un hiperplano separador (y por lo tanto el EMV no exista) cuando $n\to\infty$ y $p/n\to\kappa$ , con $0<\kappa < 1$ . Esta probabilidad es relevante en la práctica porque, dada la gran cantidad de datos a nuestra disposición, muchas veces se puede presentar la posibilidad de ajustar el modelo de regresión logística con un número elevado de variables regresoras.

El resultado es sorprendente ya que existe un brusco cambio de régimen en función de si $\kappa$ está por encima o por debajo de cierto umbral. Si está por debajo, el EMV existe en el límite con probabilidad uno; mientras que si está por encima, el EMV existe con probabilidad cero. El umbral solo depende del valor del término independiente del modelo $\beta_0$ y del límite de las varianzas $\gamma=\lim_{n\to\infty} Var(\beta'x)$ , donde $\beta$ es el vector de coeficientes de la regresión y $x$ es el vector de variables regresoras.

Una pequeña simulación

Tenía curiosidad por comprobar la relevancia empírica de este resultado y, en particular, si se notaba claramente el cambio de régimen para tamaños muestrales moderados. Para ello, he hecho una pequeña simulación y comparto aquí los resultados.

He considerado el caso más sencillo en el que $\beta_0=\beta=0$ , es decir, la variable respuesta es independiente de las variables regresoras. En este caso, puede comprobarse que el umbral que determina el cambio de régimen es $\kappa=1/2$ . Por lo tanto, si el número de variables regresoras está por debajo de la mitad del tamaño muestral debemos esperar probabilidades altas de existencia del EMV mientras que si está por encima se espera una probabilidad de existencia baja. Hay que decir que el umbral crítico es una función decreciente de $\gamma$ y por lo tanto si fuese $\beta\neq 0$ el umbral sería aún menor.

Para $n=200$ he generado muestras con $p$ variables regresoras independientes con distribución normal estándar, para $p=2,\ldots,200$ y, en cada caso, he comprobado si son linealmente separables (doy más detalles al final sobre esta comprobación). Este experimento lo he repetido 100 veces. En el siguiente gráfico se muestra la proporción de veces en las que el EMV existe en función del número de variables regresoras:

Existencia200

Los resultados coinciden totalmente con lo esperado. Hasta una dimensión cercana a 100 el EMV existe siempre. En un entorno de 100 la proporción de veces en las que existe desciende bruscamente a cero y se mantiene en cero a partir de valores de la dimensión ligeramente por encima de 100.

A continuación he repetido el experimento para tamaños muestrales aún menores ( $n=25, 50, 75$ y $100$ ) con los resultados siguientes:

Existencia25a100

Vemos que ya para $n=25$ el efecto del cambio de régimen se nota claramente. Incluso para tamaños muestrales pequeños no es recomendable ajustar un modelo de regresión logística con un número de variables regresoras por encima del umbral que describen Candes y Pragya en su artículo. Desgraciadamente este umbral no es conocido porque depende en general de los valores de los parámetros que queremos estimar.

¿Qué ocurre para variables regresoras con otras distribuciones?

El teorema se cumple bajo condiciones más generales, al menos en el caso en que $y$ es independiente de $x$ . De hecho, bajo la condición $\beta_0=\beta=0$ hay resultados generales ya en los años sesenta debidos a Tom Cover [Cover, 1965]. A continuación aparece el resultado que he obtenido simulando variables regresoras independientes con distribución exponencial de media uno:

Existencia200exp

Comprobación de la separabilidad lineal

Para comprobar si existe un hiperplano separador he seguido el enfoque basado en optimización lineal propuesto por Konis (2007) Para ello, he usado la función separator incluida en el paquete de R safeBinaryRegression. Como la función está un poco escondida, yo la he copiado directamente en mi código, que se puede encontrar aquí por si alguien quiere repetir o modificar lo que yo he hecho.

Referencias

Cover, T. M. (1965). Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition. IEEE transactions on electronic computers, 3, 326-334.

Candes, E. y Pragya, S. (2018). The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression. https://arxiv.org/abs/1804.09753

Konis, K., B. (2007). Linear programming algorithms for detecting separated data in binary logistic regression models. PhD thesis, University of Oxford.

Publicado en estadística | Etiquetado máxima verosimilitud, regresión logística | Deja un comentario

Curso de estadística aplicada con R (cuarta edición)

Publicado el 23/05/2018 por JRB

Los Departamentos de Matemáticas y Ecología de la Facultad de Ciencias de la Universidad Autónoma de Madrid organizan, en colaboración con el departamento de Bioestadística de GEICAM, la cuarta edición del curso Estadística aplicada con el software R. El programa está compuesto por diferentes módulos que pueden ser cursados en su totalidad o separadamente, dependiendo de los diferentes intereses y conocimientos de los alumnos. Los módulos se impartirán entre septiembre y noviembre de 2018 en la Facultad de Ciencias de la UAM:

Introducción a R (este módulo se impartirá en dos turnos, el 6-7 y el 20-21 de septiembre)
Métodos de regresión y análisis multivariante con R (3 a 5 de octubre)
Métodos de Regresión Avanzados para la Investigación en Ciencias Naturales con R (16 a 18 de octubre).
Estadística Aplicada a la Investigación Biomédica con R (22 a 26 de octubre).
Modelos Mixtos / Jerárquicos / Multinivel con R (20 a 22 de noviembre).
Técnicas Estadísticas de Data Mining con R (26 a 30 de noviembre).

Más información en este enlace.

Publicado en estadística, R | Etiquetado cursos | Deja un comentario

VII Jornada Estadística UAM

Publicado el 12/02/2018 por JRB

jornada

El próximo día 23 de febrero de 2018, viernes, tendremos en la sala 520 del Departamento de Matemáticas (módulo 17 de la Facultad de Ciencias) la VII jornada estadística UAM, con diferentes charlas de miembros del departamento y de algunos profesores invitados de otras universidades. Este es el programa completo de la jornada:

jornada-programa

Publicado en estadística | Etiquetado seminarios | Deja un comentario

Regresión logística y datos con grupos linealmente separables

Publicado el 29/11/2017 por JRB

Recientemente, los participantes en una encuesta realizada por la página web kaggle afirmaban que el modelo de regresión logística era la técnica estadística más utilizada en ciencia de datos. Este modelo proporciona una regla lineal de clasificación. Sin embargo, si tratamos de ajustar un modelo de regresión logística a datos con dos grupos linealmente separables (en principio, los datos que mejor se adaptan a reglas de clasificación lineales) nos encontraremos con problemas. En esta entrada se explica por qué.

El modelo

Se observan datos $(y_1,x_1),\ldots, (y_n,x_n)$ , donde $y_i =0$ o $y_i=1$ es la variable respuesta y $x_i\in \mathbb{R}^p$ es un vector de $p$ variables regresoras. Consideramos que estos datos siguen el modelo:

$\mathbb{P}(Y_i=1 | x_i) = F(\beta_0 + \beta'x_i)$

donde $\beta_0\in \mathbb{R}$ , $\beta\in \mathbb{R}^p$ y $F$ es una función de distribución estrictamente creciente, con $0<F(x)<1$ para todo $x\in\mathbb{R}$ .

El caso más importante viene dado por el modelo de regresión logística, para el que $F(x)=(1+e^{-x})^{-1}$ . Si $F(x)$ es la función de distribución normal estándar, el modelo se llama probit. Los parámetros $\beta_0\in \mathbb{R}$ y $\beta\in \mathbb{R}^p$ se suelen estimar usando el método de máxima verosimilitud.

Muestras separadas linealmente

Supongamos que existen $\bar\beta_0\in \mathbb{R}$ y $\bar\beta\in \mathbb{R}^p$ tales que $\bar\beta_0 + \bar\beta' x_i >0$ siempre que $y_i=1$ y $\bar\beta_0 + \bar\beta' x_i <0$ siempre que $y_i=0$ . Esto equivale a decir que los datos correspondientes a las dos clases pueden ser separados linealmente mediante un hiperplano. Un ejemplo de esta situación es la muestra siguiente (los puntos negros corresponden a los ceros y los rojos a los unos):

ejemplo

Si tratamos de ajustar con R un modelo de regresión logística a estos datos obtenemos el siguiente resultado:

output

Vemos un aviso según el cual parece que el algoritmo de optimización utlizado para maximizar la verosimilitud no llega a converger. ¿Por qué surge este problema?

No existe el estimador de máxima verosimilitud

Vamos a ver que si los grupos son separables linealmente el estimador de máxima verosimilitud no existe: dado un vector de parámetros cualquiera al que corresponde un cierto valor de la función de verosimilitud, siempre es posible encontrar otro vector para el que la función de verosimilitud es aún mayor.

La función de verosimilitud que hay que maximizar para estimar los parámetros del modelo es:

$L(\beta_0,\beta) = \prod_{i:\, y_i=1} F(\beta_0 +\beta'x_i) \prod_{i:\, y_i=0} (1-F(\beta_0 +\beta'x_i))$ .

Siempre se verifica $L(\beta_0,\beta)<1$ .

Sean $\bar\beta_0\in \mathbb{R}$ y $\bar\beta\in \mathbb{R}^p$ los valores que definen el hiperplano separador y definamos la sucesión $\beta_n = n\bar\beta$ .

Si $y_i=1$ se tiene que

$\lim_{n\to\infty} F(\bar\beta_0 +\beta_n' x_i) = 1,$

mientras que si $y_i=0$ se verifica

$\lim_{n\to\infty} F(\bar\beta_0 +\beta_n' x_i) = 0.$

Como consecuencia $\lim_{n\to\infty} L(\bar\beta_0,\beta_n)=1$ , es decir, haciendo $n$ suficientemente grande el valor de la verosimilitud se puede aproximar todo lo que queramos a su supremo.

Más información

En Silvapulle (1981) se caracterizan las situaciones para las que existen los estimadores de máxima verosimilitud. En particular se demuestra que si las dos clases están solapadas (no son linealmente separables), entonces el estimador de máxima verosimilitud sí que existe. La situación en la que hay más de dos clases se analiza en Albert y Anderson (1984). Ambos artículos han sido las referencias básicas para escribir esta entrada.

Código

Este es el código de R que he usado para generar los datos y ajustar el modelo:

library(MASS)
n0 <- 100
n1 <- 100
mu0 <- c(0,0)
mu1 <- c(5,5)
sigma <- diag(c(1,1))
y <- factor(c(rep(0, n0), rep(1, n1)))
set.seed(100)
x0 <- mvrnorm(n0, mu0, sigma)
x1 <- mvrnorm(n1, mu1, sigma)
datos <- data.frame(rbind(x0, x1), y)
plot(datos[,-3], col=datos$y, bty='l', pch=16)
glm(y ~ ., data = datos, family = binomial)

Referencias

Albert, A., Anderson, J. A. (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71, 1–10.

Silvapulle, M. (1981). On the Existence of Maximum Likelihood Estimators for the Binomial Response Models. Journal of the Royal Statistical Society, Series B, 43, 310-313.

Publicado en estadística | Etiquetado máxima verosimilitud, regresión logística | 2 comentarios

Sobre los ceros en la matriz de precisión

Publicado el 25/10/2017 por JRB

En esta entrada se explica cómo interpretar la aparición de ceros en algunas posiciones de la matriz de precisión de un vector aleatorio, que es la inversa de su matriz de covarianzas.

Un ejemplo

Supongamos que $\epsilon_1$ , $\epsilon_2$ y $\epsilon_3$ son v.a.i.i.d. con distribución normal estándar y definamos las variables

$X_1=\epsilon_1$ ,

$X_2=0.8X_1 + \epsilon_2$ ,

$X_3 = 0.8X_1 +\epsilon_3$ .

Es muy fácil calcular la matriz de covarianzas $\Sigma$ del vector $X:=(X_1,X_2,X_3)$ :

$\Sigma = \left(\begin{array}{ccc} 1 & 0.8 & 0.8 \\ 0.8 & 1.64 & 0.64 \\ 0.8 & 0.64 & 1.64 \end{array}\right).$

Ninguna de las covarianzas es igual a cero, lo que indica que existen relaciones lineales entre cada par de variables. Esto es obvio si consideramos la relación entre $X_1$ y $X_2$ o entre $X_1$ y $X_3$ a partir de las correspondientes definiciones. Sin embargo, aunque también existe una relación lineal entre $X_2$ y $X_3$ , es una relación lineal inducida por la influencia que $X_1$ ejerce sobre ambas. Una vez eliminado el efecto de $X_1$ no habría relación alguna entre $X_2$ y $X_3$ . De hecho, dado que el vector $X$ es normal, las variables $X_2$ y $X_3$ son independientes si condicionamos al valor de $X_1$ , lo que implica que una vez conocido el valor de $X_1$ , el conocimiento de $X_2$ no aporta ninguna información para predecir $X_3$ , y recíprocamente.

En la matriz $\Sigma$ están mezcladas las relaciones lineales directas e indirectas entre las variables por lo que a partir de $\Sigma$ no podemos distinguir entre ellas. Sin embargo, existe una manera fácil de detectar cuándo la única relación lineal entre dos variables es la inducida por la influencia del resto. El método tiene que ver con la inversa de $\Sigma$ , que a veces se llama matriz de precisión $\Omega = \Sigma^{-1}$ . En nuestro ejemplo tenemos:

$\Omega = \Sigma^{-1} = \left(\begin{array}{ccc} 2.28 & -0.8 & -0.8 \\ -0.8 & 1 & 0 \\ -0.8 & 0 & 1 \end{array}\right).$

Observamos que en las posiciones correspondientes a las variables $X_2$ y $X_3$ aparecen ceros. Esto no es casualidad ya que siempre que en una posición de la matriz de precisión aparece un cero, la única relación lineal entre las variables correspondientes a la fila y la columna de esa posición es la inducida por el resto de variables. Por lo tanto, la detección de variables condicionalmente incorreladas requiere únicamente identificar los ceros de la matriz de precisión.

Demostración

A continuación se da una demostración de las afirmaciones anteriores. Se asume que existen todas las matrices inversas que aparecen:

Más información

Los problemas de estimación cuando se impone que algunas variables son condicionalmente incorreladas (lo que equivale como hemos visto a imponer que algunas entradas de la matriz de precisión valen cero) fueron originalmente estudiados por Dempster (1972) en un artículo clásico (más de mil citas).

La estimación de matrices de precisión, especialmente en el caso de vectores de alta dimensión para los que la matriz de precisión tiene muchos ceros (es sparse), ha recibido cierta atención en la literatura estadística reciente y es la base de los llamados modelos gráficos gausianos (GGM).

Algunas referencias relacionadas:

Dempster, A. P. (1972). Covariance selection. Biometrics, 157-175.

Lauritzen, S. L. (1996). Graphical models. Clarendon Press.

Meinshausen, N., y Bühlmann, P. (2006). High-dimensional graphs and variable selection with the lasso. The Annals of Statistics, 1436-1462.

Publicado en estadística | Etiquetado análisis multivariante, correlación | Deja un comentario