Вопрос о логистической регрессии

Я хочу запустить бинарную логистическую регрессию, чтобы смоделировать наличие или отсутствие конфликта (зависимой переменной) из набора независимых переменных в течение 10-летнего периода (1997-2006 гг.), Причем каждый год имеет 107 наблюдений. Мои независимые:

• деградация земель (категорически для 2 типов деградации);
• увеличение численности населения (0 - нет; 1 - да);
• тип средств к существованию (0 - тип один; 1 - тип два);
• плотность населения (три уровня плотности);
• NDVI непрерывный (максимальная производительность овощей);
• NDVI t - 1 (снижение вегетационного периода по сравнению с предыдущим годом - 0 - нет; 1-да) и$_{t-1}$
• и NDVI t - 2 (снижение вегетации за последние два года - 0 - нет; 1 - да).$_{t-2}$

Я довольно новичок во всем этом - это проект, который мне дал лектор, - и поэтому я был бы благодарен за некоторые советы или рекомендации. Я уже проверил на мультиколлинеарность уже.

По сути, мои данные разделены на 107 единиц наблюдения (пространственных областей), охватывающих 10 лет (всего 1070), и для каждой единицы наблюдения это дает «снимок» значения условий независимых переменных в это время в пределах этой единицы ( область). Я хочу знать, как настроить мою логистическую регрессию (или таблицу) для распознавания 107 значений каждого года отдельно, чтобы можно было оценить временные изменения NDVI между различными единичными годами?

На самом деле это чрезвычайно сложная проблема и сложная задача вашего лектора!

С точки зрения того, как вы организовываете свои данные, прямоугольник 1070 x 10 - это хорошо. Например, в 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

Для подгонки модели, функция glm (), как предлагает @ gui11aume, сделает основы ...

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

... но проблема в том, что он рассматривает "страну" (я предполагаю, что у вас есть страна в качестве 107 единиц) как фиксированный эффект, тогда как случайный эффект более уместен. Он также рассматривает период как простой фактор, автокорреляция не допускается.

Вы можете обратиться к первой проблеме , связанной с обобщенной линейной модели смешанных эффектов , как, например , в Бейтс и др - х lme4 пакет в R. Там хороший введение в некоторые аспекты этого здесь . Что-то вроде

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

было бы шагом вперед.

Теперь вашей последней оставшейся проблемой является автокорреляция между вашими 10 периодами. По сути, ваши 10 точек данных по каждой стране не стоят столько, сколько если бы они были 10 случайно выбранными независимыми и одинаково распределенными точками. Мне не известно о широко доступном программном решении для автокорреляции в остатках многоуровневой модели с ненормальным откликом. Конечно, это не реализовано в lme4. Другие могут знать больше, чем я.

Этот урок является всеобъемлющим.

В R вам нужно подготовить свои данные, скажем, переменную dataв a data.frame, первый столбец которой - ваша переменная 0-1 (конфликт), а другие столбцы - предикторы. Для категориальных переменных вы должны убедиться, что они имеют тип factor. Чтобы убедиться, что столбец 3, скажем, имеет это свойство, вы можете принудительноdata[,3] <- as.factor(data[,3]) .

Тогда это только вопрос

glm(data, family="binomial")

Это подразумевает, что у вас есть аддитивная модель, и дает вам оценочные значения. Чтобы получить более полный вывод, с тестом на отдельные параметры, вы можете сделать

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