¿Cómo probar y evitar la multicolinealidad en un modelo lineal mixto?

Actualmente estoy ejecutando algunos modelos lineales de efectos mixtos.

Estoy usando el paquete "lme4" en R.

Mis modelos toman la forma:

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

Antes de ejecutar mis modelos, verifiqué la posible multicolinealidad entre los predictores.

Hice esto por:

Hacer un marco de datos de los predictores

dummy_df <- data.frame(predictor1, predictor2)

Use la función "cor" para calcular la correlación de Pearson entre predictores.

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

Si "correl_dummy_df" era mayor que 0.80, entonces decidí que predictor1 y predictor2 estaban demasiado correlacionados y no estaban incluidos en mis modelos.

Al leer un poco, aparecerían formas más objetivas de verificar la multicolinealidad.

¿Alguien tiene algún consejo sobre esto?

El "Factor de inflación de varianza (VIF)" parece un método válido.

El VIF se puede calcular utilizando la función "corvif" en el paquete AED (no Cran). El paquete se puede encontrar en http://www.highstat.com/book2.htm . El paquete admite el siguiente libro:

Zuur, AF, Ieno, EN, Walker, N., Saveliev, AA y Smith, GM 2009. Modelos y extensiones de efectos mixtos en ecología con R, 1ª edición. Springer, Nueva York.

Parece que una regla general es que si VIF es> 5, entonces la multicolinealidad es alta entre los predictores.

¿Usar VIF es más robusto que la simple correlación de Pearson?

Actualizar

Encontré un blog interesante en:

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

El blogger proporciona un código útil para calcular VIF para modelos del paquete lme4.

He probado el código y funciona muy bien. En mi análisis posterior, descubrí que la multicolinealidad no era un problema para mis modelos (todos los valores VIF <3). Esto fue interesante, dado que anteriormente había encontrado una alta correlación de Pearson entre algunos predictores.

Para el cálculo de VIF, usdm también puede ser un paquete (necesito instalar "usdm")

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

Si VIF> 4.0, generalmente asumo que la multicolinealidad elimina todas esas variables predictoras antes de ajustarlas a mi modelo

Una actualización, ya que esta pregunta me pareció útil pero no puedo agregar comentarios:

El código de Zuur et al. (2009) también está disponible a través del material complementario para su posterior publicación (y muy útil) en la revista Methods in Ecology and Evolution .

El documento, un protocolo para la exploración de datos para evitar problemas estadísticos comunes , proporciona consejos útiles y una referencia muy necesaria para justificar los umbrales de VIF (recomiendan un umbral de 3). El documento está aquí: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/full y el código R está en la pestaña de materiales complementarios (descarga de .zip).

Una guía rápida : para extraer los factores de inflación de varianza (VIF), ejecute su código HighStatLib.r y use la función corvif. La función requiere un marco de datos con solo los predictores (por ejemplo, df = data.frame(Dataset[,2:4])si sus datos se almacenan en el conjunto de datos con los predictores en las columnas 2 a 4.

Quizás la qr()función funcione. Si Xes su marco de datos o matriz, puede usar qr(X)$pivot. Por ejemplo, las qr(X)$pivot= c(1, 2, 4, 5, 7, 8, 3, 6)columnas 3 y 6 son la variable multicolineal.

Para evaluar la multicolinealidad entre predictores al ejecutar la función de dragado (paquete MuMIn), incluya la siguiente función 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))
  }
}

luego simplemente ejecute dragado especificando el número de variables predictoras e incluyendo la función 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

Esto funciona para los modelos lme4. Para modelos nlme ver: https://github.com/rojaff/dredge_mc

El VIF (factor de inflación de varianza) se puede medir simplemente por:

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