VII Jornada Estadística UAM

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

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

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

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\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 , | 2 comentarios

Sobre los ceros en la matriz de precisión

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:

Demostracion.JPG

 

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 selectionBiometrics, 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 lassoThe Annals of Statistics, 1436-1462.

 

 

 

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

Una introducción al paquete CARET

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 | Deja un comentario

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

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:

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

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

Estadística II (Grado en Matemáticas)

Añado a la página de cursos del blog la mayor parte del material que he estado utilizando los dos últimos cursos en la asignatura Estadística II, optativa en el Grado en Matemáticas de la UAM. Dado que es un material de clase que tal vez requeriría más tiempo para ser revisado es posible que aún queden errores e imprecisiones por lo que debe utilizarse con precaución.

Se trata de una segunda asignatura de estadística después de otra obligatoria, en la que se estudia la teoría básica de inferencia. El núcleo central de esta segunda asignatura es el modelo de regresión lineal. Previamente se dan dos temas, uno sobre la distribución normal multivariante que proporciona la base para poder demostrar muchos de los resultados del curso, y otro sobre contrastes de hipótesis no paramétricas (bondad de ajuste, homogeneidad e independencia). El curso termina con una introducción a los métodos de clasificación supervisada en la que se tratan dos clasificadores lineales: el basado en el modelo de regresión logística y la regla de Fisher.

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

Una variable aleatoria que casi siempre toma valores menores que su media

Si la distribución de una variable aleatoria es asimétrica la media puede dar una impresión equivocada de los valores que toma la variable. Dada una probabilidad  p\in (0,1) tan cercana a uno como se quiera, en esta entrada veremos como dar un ejemplo sencillo de variable aleatoria X tal que \mathbb{P}(X < \mathbb{E}(X)) > p, es decir, si elegimos un valor de p cercano a uno, X toma valores que son casi siempre menores que su valor esperado \mathbb{E}(X). Posteriormente veremos que este ejemplo da lugar a un estimador insesgado que se comporta de una manera inesperada. He encontrado el ejemplo en  este blog.

El ejemplo

Sea Z una variable normal estándar y sea \sigma>0. Definimos X=\exp(\sigma Z), es decir, X es una v.a. con distribución log-normal de parámetros 0 y \sigma. La esperanza de esta variable es \mathbb{E}(X)=\exp(\sigma^2/2). Por lo tanto, dado p\in (0,1) arbitrario,

\mathbb{P}(X < \mathbb{E}(X)) = \mathbb{P}(\sigma Z < \sigma^2/2) = \mathbb{P}(Z < \sigma/2) > p,

para un valor de \sigma suficientemente grande.

Estimación de la media de una variable log-normal

Supongamos que X_1,\ldots,X_n es una muestra de v.a.i.i.d. con distribución log-normal de parámetros 0 y \sigma y que queremos estimar la media \mu = \exp(\sigma^2/2) de la población. Consideramos los dos estimadores alternativos siguientes:

  • La media muestral: este es un estimador insesgado que, además, no requiere ningún conocimiento de cuál es la distribución que siguen los datos.
  • Sabiendo que la población es log-normal, se puede usar también el  estimador plug-in que resulta de la expresión \mu = \exp(\sigma^2/2), donde \sigma^2 = \mbox{Var}(\log(X)). Esto lleva a \hat{\mu} =\exp(S^2/2), donde S^2 es la varianza muestral de \log X_1,\ldots,\log X_n. Este estimador tendrá un pequeño sesgo que no debería ser muy importante para tamaños muestrales grandes.

¿Cómo funcionan estos estimadores? Para responder a esta pregunta he simulado muestras muy grandes (de tamaño 10^6) para distintos valores de \sigma. Para cada muestra he calculado el cociente entre cada uno de los dos estimadores anteriores y el valor del parámetro. Los resultados son los que muestra el gráfico:

LogNormal

Los puntos negros son realizaciones de un estimador insesgado, mientras que los blancos corresponden a un estimador que sí tiene sesgo. Así que los resultados parecen ser justo al revés de lo que dice la intuición. El estimador insesgado prácticamente siempre toma valores menores que el parámetro. Sin embargo, su sesgo es cero porque en ocasiones con probabilidad muy pequeña toma un valor muy grande que compensa las numerosas veces en las que el parámetro se infraestima. En la simulación ha ocurrido una vez (el punto negro de la parte superior derecha).

A pesar de no ser insesgado, el estimador plug-in es mucho mejor que la media muestral, para estimar la media de la población. Esto no sorprende porque incorpora el conocimiento de cuál es la distribución de los datos. La media muestral en este caso funciona mal porque hereda el comportamiento de la variable, que está casi siempre por debajo de su media.

Este es el código de R usado en la simulación:

set.seed(100)
n <- 10^6
sigma <- seq(0.1, 5, 0.1)

theta1 = numeric(length(sigma))
theta2 = numeric(length(sigma))
theta = numeric(length(sigma))

for (i in 1:length(sigma)){
z <- rnorm(n)
x <- exp(sigma[i]*z)
theta[i] = exp(sigma[i]^2/2)
theta1[i] = mean(x)
theta2[i] = exp(sd(log(x))^2/2)
}

plot(sigma, theta1/theta, pch=16, xlab='sigma',
ylab='estimador / parámetro', xlim=range(sigma))
abline(h=1)
points(sigma, theta2/theta)
legend('topleft', c('media muestral','plug in'), pch=c(16,1))
Publicado en estadística, probabilidad | Etiquetado | Deja un comentario