Caminos aleatorios

¿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 incluso 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.

Anuncios

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

Una introducción al paquete CARET

Publicado el 11/08/2017 por JRB

El paquete CARET (classification and regression training, Kuhn (2016)) incluye una serie de funciones que actúan como interfaz para decenas de métodos complejos de clasificación y regresión. Utilizar este paquete en lugar de las funciones originales de los métodos presenta dos ventajas:

Permite utilizar un código unificado para aplicar reglas de clasificación muy distintas, implementadas en diferentes paquetes.
Es más fácil poner en práctica algunos procedimientos usuales en problemas de clasificación. Por ejemplo, hay funciones específicas para dividir la muestra en datos de entrenamiento y datos de test o para ajustar parámetros mediante validación cruzada.

He tenido que usar este paquete para elaborar algunos ejemplos de un artículo y, de paso, he escrito una introducción muy básica con un resumen de lo que he aprendido. La podéis encontrar en este enlace.

Una “chuleta” que resulta útil para recordar lo básico del uso de Caret se puede encontrar en este enlace.

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

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

Publicado el 23/07/2017 por JRB

La Facultad de Ciencias de la Universidad Autónoma de Madrid (UAM) y la Unidad de Bioestadística del Instituto IMDEA Alimentación han organizado conjuntamente el curso de formación continua en Estadística Aplicada con R (tercera edición) que se va a impartir entre septiembre y diciembre de 2017 en Madrid de forma presencial

El programa está compuesto por seis módulos que se pueden cursar de forma independiente:

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

Todos los módulos se impartirán en las aulas de informática de la Facultad de Ciencias de la UAM. Más información en el siguiente enlace: http://tinyurl.com/ydgo79bf .