Pytanie o regresję logistyczną

Chcę uruchomić binarną regresję logistyczną, aby modelować obecność lub brak konfliktu (zmienna zależna) z zestawu zmiennych niezależnych w okresie 10 lat (1997-2006), przy czym każdego roku ma 107 obserwacji. Moi niezależni to:

• degradacja gruntów (kategoryczna dla 2 rodzajów degradacji);
• wzrost populacji (0 - nie; 1 - tak);
• typ źródła utrzymania (0 - typ jeden; 1 - typ drugi);
• gęstość zaludnienia (trzy poziomy gęstości);
• Ciągłe NDVI (maksymalna wydajność warzyw);
• NDVI (spadek warzyw w stosunku do poprzedniego roku - 0 - nie; 1-tak) i$_{t-1}$
• i NDVI (spadek wegetacji od dwóch lat - 0 - nie; 1 - tak).$_{t-2}$

Jestem całkiem nowy - to projekt, który dał mi mój wykładowca - i dlatego byłbym wdzięczny za porady lub wskazówki. Testowałem już pod kątem wielokoliczności.

Zasadniczo moje dane są podzielone na 107 jednostek obserwacji (regiony przestrzenne) obejmujących 10 lat (łącznie 1070) i ​​dla każdej jednostki obserwacji daje wartość „migawki” warunków zmiennych niezależnych w tym czasie w obrębie tej jednostki ( region). Chcę wiedzieć, jak skonfigurować moją regresję logistyczną (lub tabelę) w celu osobnego rozpoznawania 107 wartości każdego roku, aby można było ocenić czasowe zmiany NDVI między różnymi latami jednostkowymi?

To jest naprawdę bardzo skomplikowany problem i trudne pytanie od wykładowcy!

Jeśli chodzi o sposób organizacji danych, prostokąt 1070 x 10 jest w porządku. Na przykład w 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

Aby dopasować model, funkcja glm (), jak sugeruje @ gui11aume, zrobi podstawy ...

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

... ale ma to problem polegający na tym, że traktuje „kraj” (zakładam, że masz kraj jako swoje 107 jednostek) jako stały efekt, podczas gdy efekt losowy jest bardziej odpowiedni. Traktuje także okres jako prosty czynnik, nie dopuszcza się autokorelacji.

Pierwszy problem można rozwiązać za pomocą uogólnionego liniowego modelu efektów mieszanych, jak np. Pakiet lme4 Batesa i in. W R. Ładne wprowadzenie do niektórych aspektów tego tutaj . Coś jak

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

byłby krok naprzód.

Teraz twoim ostatnim pozostałym problemem jest autokorelacja przez 10 okresów. Zasadniczo twoje 10 punktów danych w każdym kraju nie jest tak warte, jakby były 10 losowo wybranymi niezależnymi i identycznie rozmieszczonymi punktami. Nie znam powszechnie dostępnego oprogramowania do autokorelacji w resztkach modelu wielopoziomowego z nietypową odpowiedzią. Na pewno nie jest zaimplementowany w lme4. Inni mogą wiedzieć więcej niż ja.

Ten samouczek jest wyczerpujący.

W R musisz przygotować swoje dane, powiedzmy zmienną dataw data.frame, której pierwsza kolumna to zmienna 0-1 (konflikt), a pozostałe kolumny to predyktory. W przypadku zmiennych jakościowych należy upewnić się, że są one typu factor. Aby upewnić się, że kolumna 3, powiedzmy, ma tę właściwość, możesz ją egzekwować data[,3] <- as.factor(data[,3]).

To tylko kwestia

glm(data, family="binomial")

Zakłada to domyślnie, że masz model addytywny i podaje wartości szacunkowe. Aby uzyskać bardziej kompleksowy wynik, możesz wykonać test poszczególnych parametrów

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