probit - modelos lineales generalizados en r
Modelos logísticos multinomiales multinivel en R (6)
Problema: necesito estimar un conjunto de modelos multiniveles logísticos multinivel y no puedo encontrar un paquete R adecuado. ¿Cuál es el mejor paquete de R para estimar tales modelos? STATA 13 recientemente agregó esta función a sus modelos de efectos combinados de múltiples niveles, por lo que la tecnología para estimar dichos modelos parece estar disponible.
Detalles: varias preguntas de investigación requieren la estimación de modelos de regresión logística multinomial en los que la variable de resultado es categórica. Por ejemplo, los biólogos podrían estar interesados en investigar qué tipo de árboles (p. Ej., Pinos, arces, robles) se ven más afectados por la lluvia ácida. Los investigadores de mercado podrían estar interesados si existe una relación entre la edad de los clientes y la frecuencia de las compras en Target, Safeway o Walmart. Estos casos tienen en común que la variable de resultado es categórica (no ordenada) y las regresiones logísticas multinomiales son el método de estimación preferido. En mi caso, estoy investigando las diferencias en los tipos de migración humana, con la variable de resultado (mig) codificada 0 = no migrada, 1 = migración interna, 2 = migración internacional. Aquí hay una versión simplificada de mi conjunto de datos:
migDat=data.frame(hhID=1:21,mig=rep(0:2,times=7),age=ceiling(runif(21,15,90)),stateID=rep(letters[1:3],each=7),pollution=rep(c("high","low","moderate"),each=7),stringsAsFactors=F)
hhID mig age stateID pollution
1 1 0 47 a high
2 2 1 53 a high
3 3 2 17 a high
4 4 0 73 a high
5 5 1 24 a high
6 6 2 80 a high
7 7 0 18 a high
8 8 1 33 b low
9 9 2 90 b low
10 10 0 49 b low
11 11 1 42 b low
12 12 2 44 b low
13 13 0 82 b low
14 14 1 70 b low
15 15 2 71 c moderate
16 16 0 18 c moderate
17 17 1 18 c moderate
18 18 2 39 c moderate
19 19 0 35 c moderate
20 20 1 74 c moderate
21 21 2 86 c moderate
Mi objetivo es estimar el impacto de la edad (variable independiente) en las probabilidades de (1) migrar internamente frente a no migrar, (2) migrar internacionalmente versus no migrar, (3) migrar internamente versus migrar internacionalmente. Una complicación adicional es que mis datos operan a diferentes niveles de agregación (por ejemplo, la contaminación opera a nivel estatal) y también me interesa predecir el impacto de la contaminación del aire (contaminación) en las probabilidades de embarcarse en un tipo particular de movimiento.
Soluciones torpes: se podría estimar un conjunto de modelos de regresión logística separados al reducir el conjunto de datos para cada modelo a solo dos tipos de migración (por ejemplo, Modelo 1: solo casos codificados mig = 0 y mig = 1; Modelo 2: solo casos codificados mig = 0 y mig = 2; Modelo 3: solo casos codificados mig = 1 y mig = 2). Este modelo simple de regresión logística multinivel podría estimarse con lme4, pero este enfoque es menos ideal porque no tiene en cuenta adecuadamente el impacto de los casos omitidos. Una segunda solución sería ejecutar modelos de multinivel logísticos multinivel en MLWiN hasta R usando el paquete R2MLwiN. Pero como MLWiN no es de código abierto y el objeto generado es difícil de usar, preferiría evitar esta opción. Sobre la base de una búsqueda exhaustiva en Internet, parece que hay cierta demanda de tales modelos, pero no estoy al tanto de un buen paquete R. Por lo tanto, sería genial si algunos expertos que han ejecutado tales modelos pudieran proporcionar una recomendación y si hubiera más de un paquete tal vez indiquen algunas ventajas / desventajas. Estoy seguro de que dicha información sería un recurso muy útil para múltiples usuarios de R. ¡¡Gracias!!
Mejor raphael
Aquí hay una implementación (no la mía). Acabo de trabajar en este código. Además, de esta manera realmente sabrás lo que está pasando bajo el capó.
http://www.nhsilbert.net/docs/rcode/multilevel_multinomial_logistic_regression.R
En general, existen dos formas de ajustar un modelo multinomial de una variable categórica con grupos J: (1) Estimar simultáneamente los contrastes J-1; (2) Estimación de un modelo logit separado para cada contraste.
¿Producen estos dos métodos los mismos resultados? No, pero los resultados son a menudo similares.
¿Qué método es mejor? Simultáneamente, la adaptación es más precisa (vea a continuación una explicación de por qué)
¿Por qué alguien usaría modelos logit separados entonces? (1) el paquete lme4 no tiene rutina para ajustar simultáneamente modelos multinomiales y no hay otro paquete R de varios niveles que pueda hacer esto. Por lo tanto, los modelos logit separados son actualmente la única solución práctica si alguien quiere estimar modelos multinomiales multinivel en R. (2) Como han argumentado algunos estadísticos poderosos (Begg y Gray, 1984; Allison, 1984, p. 46-47), logit separado Los modelos son mucho más flexibles ya que permiten la especificación independiente de la ecuación del modelo para cada contraste.
¿Es legítimo usar modelos logit separados? Sí, con algunos descargos de responsabilidad. Este método se denomina "Aproximación de Begg y Gray". Begg y Gray (1984, p. 16) demostraron que este "método individualizado es altamente eficiente". Sin embargo, existe cierta pérdida de eficiencia y la aproximación de Begg y Gray produce errores estándar más grandes (Agresti 2002, p. 274). Como tal, es más difícil obtener resultados significativos con este método y los resultados pueden considerarse conservadores. Esta pérdida de eficiencia es mínima cuando la categoría de referencia es grande (Begg y Gray, 1984; Agresti 2002). Los paquetes R que emplean la aproximación de Begg y Gray (no de varios niveles) incluyen mlogitBMA (Sevcikova y Raftery, 2012).
¿Por qué una serie de modelos logit individuales es imprecisa? En mi ejemplo inicial tenemos una variable ( migration ) que puede tener tres valores A (sin migración), B (migración interna), C (migración internacional). Con una sola variable predictiva x (edad), los modelos multinomiales se parametrizan como una serie de contrastes binomiales como sigue (Long y Cheng, 2004 p. 277):
Eq. 1: Ln(Pr(B|x)/Pr(A|x)) = b0,B|A + b1,B|A (x)
Eq. 2: Ln(Pr(C|x)/Pr(A|x)) = b0,C|A + b1,C|A (x)
Eq. 3: Ln(Pr(B|x)/Pr(C|x)) = b0,B|C + b1,B|C (x)
Para estos contrastes deben mantenerse las siguientes ecuaciones:
Eq. 4: Ln(Pr(B|x)/Pr(A|x)) + Ln(Pr(C|x)/Pr(A|x)) = Ln(Pr(B|x)/Pr(C|x))
Eq. 5: b0,B|A + b0,C|A = b0,B|C
Eq. 6: b1,B|A + b1,C|A = b1,B|C
El problema es que estas ecuaciones (ecuación 4-6) en la praxis no se mantienen exactamente porque los coeficientes se estiman en base a muestras ligeramente diferentes ya que solo se utilizan los casos de los dos grupos contrastantes y se omiten los casos del tercer grupo. Los programas que estiman simultáneamente los contrastes multinomiales aseguran que la ec. 4-6 espera (Long y Cheng, 2004 p. 277). No sé exactamente cómo funciona este modelo de resolución "simultánea". ¿Tal vez alguien pueda proporcionar una explicación? El software que realiza el ajuste simultáneo de modelos multinomiales multinivel incluye MLwiN (Steele 2013, p. 4) y STATA (comando xlmlogit, Pope, 2014).
Referencias:
Agresti, A. (2002). Análisis de datos categóricos (2ª ed.). Hoboken, NJ: John Wiley & Sons.
Allison, PD (1984). Análisis de la historia del evento. Thousand Oaks, CA: Sage Publications.
Begg, CB, y Gray, R. (1984). Cálculo de parámetros de regresión logística policotómica utilizando regresiones individualizadas. Biometrika, 71 (1), 11-18.
Long, SJ, y Cheng, S. (2004). Modelos de regresión para resultados categóricos. En M. Hardy y A. Bryman (Eds.), Manual de análisis de datos (pp. 258-285). Londres: SAGE Publications, Ltd.
Papa, R. (2014). En el centro de atención: Conozca el nuevo comando xlmlogit de Stata. Stata News, 29 (2), 2-3.
Sevcikova, H., y Raftery, A. (2012). Estimación del modelo logit multinomial utilizando la aproximación de Begg & Gray.
Steele, F. (2013). Módulo 10: Modelos de un solo nivel y multinivel para conceptos de respuestas nominales. Bristol, Reino Unido, Centro para el modelado multinivel.
Estoy sorprendido de que esta técnica se considere "estándar" y "equivalente", aunque podría ser una buena solución práctica. (Supongo que sería mejor revisar las referencias de Allison y Dobson & Barnett). Para el caso multinomial simple (sin agrupaciones, medidas repetidas, etc.) Begg y Gray (1984) proponen utilizar logits binomiales K-1 en una categoría de referencia como una aproximación (aunque buena) en muchos casos para logit multinomial en toda regla. Demuestran cierta pérdida de eficiencia cuando usan una categoría de referencia única, aunque es pequeña para los casos en que se usa una categoría de línea de base de alta frecuencia como referencia. Agresti (2002: p. 274) proporciona un ejemplo donde hay un pequeño aumento en los errores estándar, incluso cuando la categoría de referencia constituye más del 70% de 219 casos en un ejemplo de cinco categorías.
Tal vez no sea un gran problema, pero no veo cómo mejoraría la aproximación al agregar una segunda capa de aleatoriedad.
Referencias
Agresti, A. (2002). Análisis de datos categóricos. Hoboken NJ: Wiley.
Begg, CB, y Gray, R. (1984). Cálculo de parámetros de regresión logística policotómica utilizando regresiones individualizadas. Biometrika, 71 (1), 11-18.
Estoy tratando el mismo problema y una posible solución que encontré parece recurrir al equivalente de poisson (loglinear / count) del modelo logístico multinomial como se describe en esta lista de mailinglist , estas buenas diapositivas o en Agresti (2013: 353-356). Por lo tanto, debería ser posible utilizar la función glmer(... family=poisson) del paquete lme4 con cierta agregación de los datos.
Referencia:
Agresti, A. (2013) Análisis de datos categóricos. Hoboken, Nueva Jersey: Wiley.
Te recomendaré usar el paquete "mlogit"
Una pregunta anterior, pero creo que recientemente surgió una opción viable es brms , que utiliza el programa Bayesian Stan para ejecutar el modelo. Por ejemplo, si desea ejecutar una regresión logística multinomial en los datos del iris :
b1 <- brm (Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width,
data=iris, family="categorical",
prior=c(set_prior ("normal (0, 8)")))
Y para obtener una regresión ordinal (no es apropiado para el iris , por supuesto), cambiaría family="categorical" a family="acat" (o cratio o sratio , dependiendo del tipo de regresión ordinal que desee) y asegúrese de que la variable dependiente está ordered .
Aclaración según el comentario de Raphael: esta llamada de brm compila su fórmula y argumentos en el código de Stan . Stan lo compila en C ++ y usa el compilador de C ++ de su sistema, que es necesario. En una Mac, por ejemplo, puede que necesite instalar las Herramientas de desarrollo gratuitas para obtener C ++. No estoy seguro acerca de Windows. Linux debería tener instalado C ++ por defecto.)
Aclaración según el comentario de brms : brms maneja fácilmente modelos multinivel y también usa la fórmula R (1 | groupvar) para agregar una intersección de grupo (aleatoria) para un grupo, (1 + foo | groupvar) para agregar una intersección aleatoria y pendiente, etc.