Jak przetestować i uniknąć wielokoliniowości w mieszanym modelu liniowym?

25

Obecnie korzystam z modeli liniowych z mieszanym efektem.

Korzystam z pakietu „lme4” w języku R.

Moje modele mają postać:

model <- lmer(response ~ predictor1 + predictor2 + (1 | random effect))

Przed uruchomieniem moich modeli sprawdziłem możliwą wielokoliniowość między predyktorami.

Zrobiłem to przez:

Utwórz ramkę danych predyktorów

dummy_df <- data.frame(predictor1, predictor2)

Użyj funkcji „cor”, aby obliczyć korelację Pearsona między predyktorami.

correl_dummy_df <- round(cor(dummy_df, use = "pair"), 2) 

Jeśli „correl_dummy_df” był większy niż 0,80, to zdecydowałem, że predyktor1 i predyktor2 były zbyt silnie skorelowane i nie zostały uwzględnione w moich modelach.

Czytając, pojawiłyby się bardziej obiektywne sposoby sprawdzenia wielokoliniowości.

Czy ktoś ma jakieś porady na ten temat?

„Współczynnik inflacji wariancji (VIF)” wydaje się być jedną z prawidłowych metod.

VIF można obliczyć za pomocą funkcji „corvif” w pakiecie AED (non-cran). Pakiet można znaleźć na stronie http://www.highstat.com/book2.htm . Pakiet obsługuje następującą książkę:

Zuur, AF, Ieno, EN, Walker, N., Saveliev, AA & Smith, GM 2009. Modele z efektami mieszanymi i rozszerzenia w ekologii z R, 1. edycja. Springer, Nowy Jork.

Wygląda na to, że ogólna zasada jest taka, że ​​jeśli VIF jest> 5, to wielokoliniowość jest wysoka między predyktorami.

Czy używanie VIF jest bardziej niezawodne niż prosta korelacja Pearsona?

Aktualizacja

Znalazłem interesującego bloga pod adresem:

http://hlplab.wordpress.com/2011/02/24/diagnosing-collinearity-in-lme4/

Blogger udostępnia użyteczny kod do obliczania VIF dla modeli z pakietu lme4.

Przetestowałem kod i działa świetnie. W mojej późniejszej analizie odkryłem, że wielokoliniowość nie stanowiła problemu dla moich modeli (wszystkie wartości VIF <3). Było to interesujące, biorąc pod uwagę, że wcześniej znalazłem wysoką korelację Pearsona między niektórymi predyktorami.

mjburns
źródło
6
(1) AEDPakiet został wycofany ; zamiast tego tylko source("http://www.highstat.com/Book2/HighstatLibV6.R")dla corviffunkcji. (2) Mam nadzieję, że udzielę prawdziwej odpowiedzi, ale (a) Uważam, że VIF bierze pod uwagę wielokoliniowość (np. Możesz mieć trzy predyktory, z których żaden nie ma silnych korelacji par, ale kombinacja liniowa A i B jest silnie skorelowana z C ) i (b) mam poważne zastrzeżenia co do mądrości porzucania terminów kolinearnych; patrz Graham Ecology 2003, doi: 10.1890 / 02-3114
Ben Bolker
Dzięki Ben. Zaktualizowałem powyższy post, aby uwzględnić Twoje sugestie.
mjburns
@BenBolker, czy możesz krótko wyjaśnić, dlaczego jesteś przeciwny porzucaniu terminów kolinearnych? Doceniam odniesienie, ale może też polubić wersję Cliff Notes. Dzięki!
Bajcz
korekta w odpowiedzi Bena .. URL tohttp://highstat.com/Books/BGS/GAMM/RCodeP2/HighstatLibV6.R
Manoj Kumar

Odpowiedzi:

10

Do obliczeń VIF usdm może być również pakietem (muszę zainstalować „usdm”)

library(usdm)
df = # Data Frame
vif(df)

Jeśli VIF> 4.0, ogólnie zakładam, że wielokoliniowość usuwa wszystkie te zmienne predykcyjne przed ich dopasowaniem do mojego modelu

Sohsum
źródło
Trochę addendum, którego można użyć do filtrowania zmiennych, takich jak wykluczenie wszystkich, które pokazują korelację powyżej .4jako vifcor(vardata,th=0.4). Podobnie możesz użyć, vifstep(vardata,th=10)aby odrzucić wszystkie większe niż 10.
SIslam
Nie działa dla HLM
Mox
7

Aktualizacja, ponieważ uznałem to pytanie za przydatne, ale nie mogę dodawać komentarzy -

Kod Zuur i in. (2009) jest również dostępny poprzez materiał uzupełniający do kolejnej (i bardzo przydatnej) publikacji w czasopiśmie Methods in Ecology and Evolution .

Artykuł - Protokół eksploracji danych w celu uniknięcia typowych problemów statystycznych - zawiera przydatne porady i bardzo potrzebne odniesienie do uzasadnienia progów VIF (zalecają próg 3). Artykuł znajduje się tutaj: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/full, a kod R znajduje się w zakładce materiałów dodatkowych (pobieranie .zip).

Krótki przewodnik : aby wyodrębnić współczynniki inflacji wariancji (VIF), uruchom ich kod HighStatLib.r i użyj funkcji corvif. Ta funkcja wymaga ramki danych zawierającej tylko predyktory (więc na przykład, df = data.frame(Dataset[,2:4])jeśli dane są przechowywane w zestawie danych z predyktorami w kolumnach od 2 do 4).

lityczny
źródło
1

Może qr()funkcja zadziała. Jeśli Xtwoja ramka danych lub macierz, możesz użyć qr(X)$pivot. Na przykład qr(X)$pivot= c(1, 2, 4, 5, 7, 8, 3, 6)wtedy kolumny 3 i 6 są zmienną wielokoliniową.

JunhuiLi
źródło
1

Aby ocenić wielokoliniowość między predyktorami podczas uruchamiania funkcji pogłębiarki (pakiet MuMIn), dołącz następującą funkcję maks. R jako argument „dodatkowy”:

max.r <- function(x){
  corm <- cov2cor(vcov(x))
  corm <- as.matrix(corm)
  if (length(corm)==1){
    corm <- 0
    max(abs(corm))
  } else if (length(corm)==4){
  cormf <- corm[2:nrow(corm),2:ncol(corm)]
  cormf <- 0
  max(abs(cormf))
  } else {
    cormf <- corm[2:nrow(corm),2:ncol(corm)]
    diag(cormf) <- 0
    max(abs(cormf))
  }
}

następnie po prostu uruchom pogłębianie, określając liczbę zmiennych predykcyjnych i włączając funkcję max.r:

options(na.action = na.fail)
Allmodels <- dredge(Fullmodel, rank = "AIC", m.lim=c(0, 3), extra= max.r) 
Allmodels[Allmodels$max.r<=0.6, ] ##Subset models with max.r <=0.6 (not collinear)
NCM <- get.models(Allmodels, subset = max.r<=0.6) ##Retrieve models with max.r <=0.6 (not collinear)
model.sel(NCM) ##Final model selection table

Działa to dla modeli lme4. Modele nlme patrz: https://github.com/rojaff/dredge_mc

r.jaffe
źródło
1

VIF (współczynnik inflacji wariancji) można zmierzyć po prostu przez:

 library(car)
 vif(yourmodel) #this should work for lme4:lmer mixed models.
Hassan.JFRY
źródło