Pregunta acerca de la regresión logística

Question

Quiero correr una regresión logística binaria para el modelo de la presencia o ausencia de conflicto (variable dependiente) a partir de un conjunto de variables independientes sobre un período de 10 años (1997-2006), con cada año que tiene 107 observaciones. Mi independientes son:

• la degradación de la tierra (categóricos para 2 tipos de degradación);
• el incremento de la población (0 - no, 1-sí);
• los medios de subsistencia de tipo (0 - tipo uno; 1 - tipo dos);
• la densidad de la población (tres niveles de densidad);
• NDVI continua (máx. verduras de la productividad);
• NDVI$_{t-1}$ (descenso en la veg que el año anterior - 0 - no, 1 -sí) y
• y NDVI$_{t-2}$ (descenso en la veg a partir de los dos años anteriores - 0 - no, 1 - sí).

Soy bastante nuevo en todo esto - este es un proyecto de mi profesor me ha dado - y entonces yo estaría agradecido de algún consejo o guía. He probado de multicolliniarity ya.

Básicamente mis datos se divide en 107 unidades de observación espacial (regiones) que abarca 10 años (1070 en total) y para cada unidad de observación que le da ser una 'instantánea' del valor de las condiciones de las variables independientes en ese momento dentro de esa unidad (región). Quiero saber como configurar mi regresión logística (o tabla) para reconocer el 107 valores de cada año por separado, de modo que el temporal de NDVI cambios entre los distintos unidad de años puede ser evaluado?

Answer 1

Esto es realmente un muy sofisticado problema y una difícil pedir a tu profesor!

En términos de cómo organizar los datos, una 1070 x 10 rectángulo está bien. Por ejemplo, en R:

> conflict.data <- data.frame(
+ confl = sample(0:1, 1070, replace=T),
+ country = factor(rep(1:107,10)),
+ period = factor(rep(1:10, rep(107,10))),
+ landdeg = sample(c("Type1", "Type2"), 1070, replace=T),
+ popincrease = sample(0:1, 1070, replace=T),
+ liveli =sample(0:1, 1070, replace=T),
+ popden = sample(c("Low", "Med", "High"), 1070, replace=T),
+ NDVI = rnorm(1070,100,10),
+ NDVIdecl1 = sample(0:1, 1070, replace=T),
+ NDVIdecl2 = sample(0:1, 1070, replace=T))
> head(conflict.data)
  confl country period landdeg popincrease liveli popden     NDVI NDVIdecl1 NDVIdecl2
1     1       1      1   Type1           1      0    Low 113.4744         0         1
2     1       2      1   Type2           1      1   High 103.2979         0         0
3     0       3      1   Type2           1      1    Med 109.1200         1         1
4     1       4      1   Type2           0      1    Low 112.1574         1         0
5     0       5      1   Type1           0      0   High 109.9875         0         1
6     1       6      1   Type1           1      0    Low 109.2785         0         0
> summary(conflict.data)
     confl           country         period     landdeg     popincrease         liveli        popden         NDVI          NDVIdecl1        NDVIdecl2     
 Min.   :0.0000   1      :  10   1      :107   Type1:535   Min.   :0.0000   Min.   :0.0000   High:361   Min.   : 68.71   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.0000   2      :  10   2      :107   Type2:535   1st Qu.:0.0000   1st Qu.:0.0000   Low :340   1st Qu.: 93.25   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :1.0000   3      :  10   3      :107               Median :1.0000   Median :1.0000   Med :369   Median : 99.65   Median :1.0000   Median :0.0000  
 Mean   :0.5009   4      :  10   4      :107               Mean   :0.5028   Mean   :0.5056              Mean   : 99.84   Mean   :0.5121   Mean   :0.4888  
 3rd Qu.:1.0000   5      :  10   5      :107               3rd Qu.:1.0000   3rd Qu.:1.0000              3rd Qu.:106.99   3rd Qu.:1.0000   3rd Qu.:1.0000  
 Max.   :1.0000   6      :  10   6      :107               Max.   :1.0000   Max.   :1.0000              Max.   :130.13   Max.   :1.0000   Max.   :1.0000  
                  (Other):1010   (Other):428                                                                                                              
> dim(conflict.data)
[1] 1070   10

Para el ajuste de un modelo, el glm() la función como @gui11aume se sugiere hacer las cosas básicas...

mod <- glm(confl~., family="binomial", data=conflict.data)
anova(mod)

... pero esto tiene el problema de que se trata de "el país" (estoy asumiendo que usted tiene el país como su 107 unidades) como un efecto fijo, mientras que un efecto aleatorio es más apropiado. También trata el período como un factor simple, no hay autocorrelación permitido.

Usted puede abordar el primer problema con una generalizada modelo lineal de efectos mixtos, como por ejemplo en Bates et al lme4 paquete en R. Hay una buena introducción a algunos aspectos de esto aquí. Algo como

library(lme4)
mod2 <- lmer(confl ~ landdeg + popincrease + liveli + popden + 
    NDVI + NDVIdecl1 + NDVIdecl2 + (1|country) +(1|period), family=binomial,
    data=conflict.data)
summary(mod2)

sería un paso adelante.

Ahora su última problema restante es de autocorrelación a través de sus 10 períodos. Básicamente, sus 10 puntos de datos en cada país, no valen tanto como si fueran 10 elegidos al azar independiente y identicall puntos distribuidos. Yo no soy consciente de que una solución de software ampliamente disponible para la autocorrelación en los residuales de un modelo multinivel con un no-respuesta Normal. Ciertamente no es implementado en lme4. Otros pueden saber más que yo.

Answer 2

Este tutorial es muy completo.

En R, es necesario preparar los datos, decir que la variable data en data.frame, la primera columna de la que es su 0-1 variable (conflicto) y las demás columnas son los predictores. Para las variables categóricas, usted debe asegurarse de que son de tipo factor. Para asegurarse de que la columna 3, por ejemplo, tiene esta propiedad, usted puede aplicar por data[,3] <- as.factor(data[,3]).

Entonces es sólo una cuestión de

glm(data, family="binomial")

Implícitamente se supone que tiene un modelo aditivo y le da los valores estimados. Para obtener una visión más integral de la salida, con la prueba de los parámetros individuales, usted puede hacer

summary(glm(data, family="binomial"))