Estoy tratando de hacer un modelo de regresión lineal múltiple para correlacionar una dependiente $Y$ variable (normalmente distribuida) frente a un conjunto de 642 variables. Estas 642 variables codifican la presencia o ausencia de un determinado grupo químico en una molécula, por lo que sólo pueden suponer 1 o 0 como valor. El tamaño del conjunto de datos es de 1300 experimentos (uno $Y$ ) y 642 $X$ variables. El problema que tengo es que las 642 variables están puestas en 1 sólo el 0,3% de las veces, es decir, mi matriz X de 1300 x 642 está formada por 2600 1 y 832000 0 . Nada más (excepto el $Y$ por supuesto). ¿La alta frecuencia de 0 ¿crear un problema para la implementación de MLR en R? ¿Existe un método estadístico más adecuado para tratar este caso?
Respuestas
¿Demasiados anuncios?
@JoeDanger ha proporcionado una buena respuesta (+1); debería tener cierta preocupación por la posibilidad de sobreajuste . Una regla general en la regresión múltiple es que se necesitan al menos 10 observaciones por variable (véase aquí: reglas generales sobre el tamaño mínimo de la muestra para la regresión múltiple para más información sobre este tema). Con 642 variables, eso significa que quiere 6.420 observaciones, mientras que sólo tiene 1.300 (~20%). Esto es sólo una regla general y no hay nada mágico en ese número, pero debe tenerlo en cuenta.
Además de sus puntos, no hay nada en el algoritmo que R utiliza para ajustar el modelo que será lanzado por tener en su mayoría 0 s, si esa es su principal preocupación. Sin embargo, tendrá un bajo poder para detectar los efectos de esas variables en la respuesta. Se puede pensar en un análisis de regresión con una variable ficticia como un Prueba t y, por lo tanto, se produce una situación en la que los grupos son de tamaño desigual. He discutido esta cuestión aquí: cómo-se-debe-interpretar-la-comparación-de-medias-de-diferentes-tamaños-de-muestra que puede ayudarle a comprender mejor el asunto. Una regresión múltiple como su situación será análoga.
Siempre que cada grupo químico aparezca en al menos dos experimentos puedes tener éxito.
Una buena manera de explorar los problemas es con una simulación. Lo que importa es la cantidad de cualquier variación aleatoria en los valores de $Y$ se compara con el tamaño de los efectos: cuando esa variación es relativamente pequeña, se puede estimar con bastante fiabilidad los efectos y saber que son reales; cuando esa variación es relativamente grande, hay que olvidarse de ella.
Exploremos utilizando datos simulados.
set.seed(17)
x.0 <- replicate(642, sample(c(rep(0, 1298), rep(1, 2))))
x.1 <- t(replicate(1300, sample(c(rep(0, 641), 1))))
x <- x.0 + x.1
sum(x) # Need to be around 2600
par(mfrow=c(1,2))
plot(apply(x, 1, sum), xlab="Case", ylab="Groups per case", cex=0.6)
plot(apply(x, 2, sum), xlab="Group", ylab="Cases per group", cex=0.6)
Los gráficos confirman que cada caso incluye al menos un grupo químico y que cada grupo químico aparece en al menos dos casos (por lo que es posible evaluar su importancia). Hay $1300$ casos, $642$ variables, y $2584$ en el conjunto de datos, lo suficientemente cerca de la descripción en el problema.
La función sim genera los verdaderos valores de $Y$ y añade una variación normal iid de una desviación estándar especificada sd . Ejecuta la regresión lineal múltiple lm entonces traza (a) un histograma de los valores "observados" de $Y$ y (b) un histograma de los valores p para valores p entre 0 y p.limit . (Sus parámetros opcionales se pasan a los dos hist llamadas). Reducción de p.limit nos permite ampliar el rango de los valores p más pequeños para verlos mejor.
sim <- function(beta, sd, p.limit=1, ...) {
d <- data.frame(x, y=x %*% beta + rnorm(dim(x)[1], sd=sd))
fit <- lm(y ~ ., data=d)
hist(d$y, xlab="Y", main="Histogram of Y", ...)
p <- summary(fit)$coefficients[, 4]
hist(p[p <= p.limit], xlab="p", main="Histogram of p-values", ...)
}Sin efectos
Ahora es fácil explorar los escenarios. Comience con el caso en el que hay no efectos:
beta <- rep(0, dim(x)[2])
sim(beta, 1, freq=FALSE)
Como era de esperar, los valores p se distribuyen de forma aproximadamente uniforme en todo su rango posible. Esto significa que se debe esperar alrededor de $.05 \times 642 = 32$ de ellos para ser menos de $5$ y alrededores $.01 \times 642 = 6$ para ser menos que $1$ % a pesar de que ninguno de las variables tiene alguna relación con $Y$ en absoluto: tiene un problema con los falsos positivos. Una forma de hacer frente a esto es la corrección de Bonferroni, que estipula que sólo los valores p inferiores a $0.05/642$ debe ser considerada significativa en el en general $5$ % de nivel. (Hay métodos mejores, pero ésta es una forma fácil de ilustrar el papel que desempeña cualquier corrección de comparaciones múltiples).
Los efectos en las siguientes simulaciones se distribuyen normalmente con una desviación estándar de $1$ , centrado en $0$ .
Efectos fuertes
Una desviación estándar en los errores de sólo $0.1$ es bastante pequeño en comparación con estos efectos. Veamos qué tal funcionan los mínimos cuadrados ordinarios con la corrección de Bonferroni:
set.seed(17)
beta <- rnorm(dim(x)[2])
sim(beta, .1, p.limit=0.05/642)
Funciona maravillosamente bien. La mayoría de los efectos son tan grandes (en comparación con la variación aleatoria) que prácticamente todos ellos se consideran significativos.
Efectos débiles
Pero no des por sentado que tus datos serán así de bonitos. Veamos quizás una situación más realista en la que los tamaños del efecto son comparables a la variación aleatoria:
sim(beta, 1, p.limit=0.05/642)
Lamentablemente, sólo un par de docenas de efectos pueden considerarse significativos aquí. Simplemente no se tiene el poder de detectar la mayoría de ellos.
Conclusiones
-
Debe protegerse contra los falsos positivos adoptando algún tipo de procedimiento de comparaciones múltiples.
-
Si tiene un gran conjunto de datos con efectos fuertes, la regresión múltiple por mínimos cuadrados ordinarios puede funcionar estupendamente. De lo contrario, puede descubrir que los datos simplemente no son suficientes para detectar los efectos que puedan existir.
No creará problemas, aunque sí lo hará la inclusión de las 642 variables en su regresión. Sólo tienes 1300 observaciones para cada variable y si incluyes todas las variables estarás pidiendo a los datos que respondan a 642 preguntas (preguntando si esta variable es significativa o no), lo cual es demasiado en mi opinión. Algo del orden de 10 variables sería más apropiado para datos de este tamaño.
También es importante señalar que todas sus variables explicativas son variables ficticias que indican la inclusión o no de un grupo químico. Esto no causa ningún problema, pero debe tener en cuenta que los coeficientes estimados en su regresión le darán la media de Y para esa variable en particular. Por ejemplo, si realizo una regresión de la altura sobre una constante y una variable ficticia que indica si alguien es hombre o no (1=hombre, 0=mujer), el coeficiente del término constante sólo me dirá la altura media de todas las mujeres y la suma del coeficiente del término constante y la variable ficticia de género me dirá la altura media de todos los hombres de la muestra. Así, si su Y es, por ejemplo, el coste del experimento y una de sus variables ficticias indica la presencia de bencina, su constante + coeficiente le dirá el coste medio de un experimento con bencina.
A mí me preocuparía sobre todo que el rendimiento de su modelo se erosionara muy rápidamente a medida que añadiera más y más variables explicativas. Sin embargo, la relación entre ceros y unos no es un problema.
