Pregunta sobre regresión logística

Quiero ejecutar una regresión logística binaria para modelar la presencia o ausencia de conflicto (variable dependiente) a partir de un conjunto de variables independientes durante un período de 10 años (1997-2006), y cada año tiene 107 observaciones. Mis independientes son:

• degradación de la tierra (categórica para 2 tipos de degradación);
• aumento de la población (0- no; 1-sí);
• tipo de sustento (0 - tipo uno; 1 - tipo dos);
• densidad de población (tres niveles de densidad);
• NDVI continuo (productividad máxima de verduras);
• NDVI t - 1 (disminución de verduras respecto al año anterior - 0 - no; 1 - sí) y$_{t-1}$
• y NDVI t - 2 (disminución de las verduras de dos años anteriores - 0- no; 1- sí).$_{t-2}$

Soy bastante nuevo en todo esto, este es un proyecto que mi profesor me ha dado, por lo que agradecería algún consejo u orientación. Ya he probado la multicoliniaridad.

Esencialmente, mis datos se dividen en 107 unidades de observación (regiones espaciales) que cubren 10 años (1070 en total) y para cada unidad de observación, da un valor de 'instantánea' de las condiciones de las variables independientes en ese momento dentro de esa unidad ( región). ¿Quiero saber cómo configurar mi regresión logística (o tabla) para reconocer los 107 valores de cada año por separado para que se puedan evaluar los cambios temporales del NDVI entre diferentes unidades de años?

¡Este es realmente un problema extremadamente sofisticado y una pregunta difícil de su profesor!

En términos de cómo organizar sus datos, un rectángulo de 1070 x 10 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 ajustar un modelo, la función glm () como sugiere @ gui11aume hará lo básico ...

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

... pero esto tiene el problema de que trata el "país" (supongo que tiene el país como sus 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 se permite la autocorrelación.

Puede abordar el primer problema con un modelo de efectos mixtos lineales generalizados como, por ejemplo, en el paquete lme4 de Bates et al 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 último problema restante es la autocorrelación en sus 10 períodos. Básicamente, sus 10 puntos de datos en cada país no valen tanto como si fueran 10 puntos distribuidos independientes e identificados elegidos al azar. No conozco una solución de software ampliamente disponible para la autocorrelación en los residuos de un modelo multinivel con una respuesta no normal. Ciertamente no está implementado en lme4. Otros pueden saber más que yo.

Este tutorial es completo.

En R, debe preparar sus datos, digamos variable dataen a data.frame, cuya primera columna es su variable 0-1 (conflicto) y las otras columnas son los predictores. Para las variables categóricas, debe asegurarse de que sean de tipofactor . Para asegurarse de que la columna 3, por ejemplo, tenga esta propiedad, puede aplicarla mediante data[,3] <- as.factor(data[,3]).

Entonces es solo cuestión de

glm(data, family="binomial")

Esto supone implícitamente que tiene un modelo aditivo y le da los valores estimados. Para obtener una salida más completa, con prueba para parámetros individuales, puede hacer

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