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

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.

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

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).

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ą.

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

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

 library(car)
 vif(yourmodel) #this should work for lme4:lmer mixed models.