Como testar e evitar a multicolinearidade no modelo linear misto?

Atualmente, estou executando alguns modelos lineares de efeito misto.

Estou usando o pacote "lme4" em R.

Meus modelos assumem a forma:

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

Antes de executar meus modelos, verifiquei a possibilidade de multicolinearidade entre preditores.

Eu fiz isso por:

Faça um quadro de dados dos preditores

dummy_df <- data.frame(predictor1, predictor2)

Use a função "cor" para calcular a correlação de Pearson entre preditores.

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

Se "correl_dummy_df" fosse maior que 0,80, decidi que o preditor1 e o preditor2 estavam muito correlacionados demais e não foram incluídos nos meus modelos.

Ao fazer algumas leituras, apareceriam maneiras mais objetivas de verificar a multicolinearidade.

Alguém tem algum conselho sobre isso?

O "fator de inflação de variação (VIF)" parece ser um método válido.

O VIF pode ser calculado usando a função "corvif" na embalagem do DEA (não cran). O pacote pode ser encontrado em http://www.highstat.com/book2.htm . O pacote suporta o seguinte livro:

Zuur, AF, Ieno, EN, Walker, N., Saveliev, AA e Smith, GM 2009. Modelos de efeitos mistos e extensões em ecologia com R, 1ª edição. Springer, Nova Iorque.

Parece que uma regra geral é que, se VIF for> 5, a multicolinearidade será alta entre os preditores.

O uso do VIF é mais robusto que a simples correlação de Pearson?

Atualizar

Encontrei um blog interessante em:

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

O blogueiro fornece um código útil para calcular o VIF para modelos do pacote lme4.

Eu testei o código e funciona muito bem. Na minha análise subsequente, descobri que a multicolinearidade não era um problema para meus modelos (todos os valores de VIF <3). Isso foi interessante, já que eu havia encontrado alta correlação de Pearson entre alguns preditores.

Para o cálculo VIF, o usdm também pode ser um pacote (eu preciso instalar o "usdm")

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

Se VIF> 4.0, em geral, presumo que a multicolinearidade remova todas essas variáveis ​​preditoras antes de ajustá-las ao meu modelo

Uma atualização, pois achei esta pergunta útil, mas não consigo adicionar comentários.

O código de Zuur et al. (2009) também está disponível através do material suplementar para uma publicação subsequente (e muito útil) de suas publicações na revista Methods in Ecology and Evolution .

O documento - um protocolo para exploração de dados para evitar problemas estatísticos comuns - fornece conselhos úteis e uma referência muito necessária para justificar os limites de VIF (eles recomendam um limite de 3). O artigo está aqui: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/ completo e o código R está na guia de materiais suplementares (download .zip).

Um guia rápido : a fatores de inflação extrato de variância (VIF) executar sua HighStatLib.r código e usar a função corvif. A função requer um quadro de dados com apenas os preditores (por exemplo, df = data.frame(Dataset[,2:4])se seus dados forem armazenados no conjunto de dados com os preditores nas colunas 2 a 4.

Talvez a qr()função funcione. Se Xé o seu quadro ou matriz de dados, você pode usar qr(X)$pivot. Por exemplo, as qr(X)$pivot= c(1, 2, 4, 5, 7, 8, 3, 6)colunas 3 e 6 são a variável multicolinear.

Para avaliar a multicolinearidade entre preditores ao executar a função draga (pacote MuMIn), inclua a seguinte função max.r como argumento "extra":

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

basta executar a dragagem especificando o número de variáveis ​​preditivas e incluindo a função 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

Isso funciona para os modelos lme4. Para modelos nlme, consulte: https://github.com/rojaff/dredge_mc

O VIF (fator de inflação da variação) pode ser medido simplesmente por:

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