Comment tester et éviter la multicolinéarité dans un modèle linéaire mixte?

26

J'utilise actuellement des modèles linéaires à effets mixtes.

J'utilise le package "lme4" dans R.

Mes modèles prennent la forme:

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

Avant d'exécuter mes modèles, j'ai vérifié la possible multicolinéarité entre les prédicteurs.

Je l'ai fait par:

Faire une trame de données des prédicteurs

dummy_df <- data.frame(predictor1, predictor2)

Utilisez la fonction «cor» pour calculer la corrélation de Pearson entre les prédicteurs.

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

Si "correl_dummy_df" était supérieur à 0,80, alors j'ai décidé que le prédicteur1 et le prédicteur2 étaient trop fortement corrélés et ils n'étaient pas inclus dans mes modèles.

En faisant quelques lectures, il semblerait que des moyens plus objectifs de vérifier la multicolinéarité.

Quelqu'un a-t-il des conseils à ce sujet?

Le «facteur d'inflation de la variance (VIF)» semble être une méthode valable.

Le VIF peut être calculé en utilisant la fonction "corvif" dans le package AED (non-cran). Le package peut être trouvé à http://www.highstat.com/book2.htm . Le package prend en charge le livre suivant:

Zuur, AF, Ieno, EN, Walker, N., Saveliev, AA & Smith, GM 2009. Modèles d'effets mixtes et extensions en écologie avec R, 1ère édition. Springer, New York.

Il semble qu'une règle générale soit que si VIF est> 5, la multicolinéarité est élevée entre les prédicteurs.

L'utilisation de VIF est-elle plus robuste que la simple corrélation de Pearson?

Mise à jour

J'ai trouvé un blog intéressant sur:

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

Le blogueur fournit du code utile pour calculer VIF pour les modèles à partir du package lme4.

J'ai testé le code et cela fonctionne très bien. Dans mon analyse ultérieure, j'ai trouvé que la multicolinéarité n'était pas un problème pour mes modèles (toutes les valeurs VIF <3). C'était intéressant, étant donné que j'avais précédemment trouvé une forte corrélation de Pearson entre certains prédicteurs.

mjburns
la source
6
(1) Le AEDpackage a été supprimé ; à la place, juste source("http://www.highstat.com/Book2/HighstatLibV6.R")pour la corviffonction. (2) J'espère apporter une vraie réponse, mais (a) je crois que VIF prend en compte la multicolinéarité (par exemple, vous pouvez avoir trois prédicteurs, dont aucun n'a de fortes corrélations par paires, mais la combinaison linéaire de A et B est fortement corrélée avec C ) et b) j'ai de fortes réserves quant à l'opportunité de supprimer les termes colinéaires; voir Graham Ecology 2003, doi: 10.1890 / 02-3114
Ben Bolker
Merci Ben. J'ai mis à jour mon message ci-dessus pour inclure vos suggestions.
mjburns
@BenBolker, pouvez-vous expliquer très brièvement pourquoi vous êtes contre l'abandon des termes colinéaires? J'apprécie la référence mais je pourrais aussi aimer une version de Cliff Notes. Merci!
Bajcz
correction dans la réponse de Ben .. l'URL esthttp://highstat.com/Books/BGS/GAMM/RCodeP2/HighstatLibV6.R
Manoj Kumar

Réponses:

10

Pour le calcul VIF, usdm peut également être un package (j'ai besoin d'installer "usdm")

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

Si VIF> 4.0, je suppose généralement que la multicolinéarité supprime toutes ces variables prédictives avant de les insérer dans mon modèle

Sohsum
la source
Un petit addenda que vous pouvez utiliser pour filtrer des variables comme exclure tout ce qui montre la corrélation ci-dessus .4comme vifcor(vardata,th=0.4). De même, vous pouvez utiliser vifstep(vardata,th=10)pour éliminer tous les plus de 10.
SIslam
Ne fonctionne pas pour HLM
Mox
7

Une mise à jour, car j'ai trouvé cette question utile mais je ne peux pas ajouter de commentaires -

Le code de Zuur et al. (2009) est également disponible via le matériel supplémentaire à une publication ultérieure (et très utile) de la leur dans la revue Methods in Ecology and Evolution .

Le document - Un protocole d'exploration des données pour éviter les problèmes statistiques courants - fournit des conseils utiles et une référence bien nécessaire pour justifier les seuils VIF (ils recommandent un seuil de 3). Le document est ici: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/full et le code R est dans l'onglet des documents supplémentaires (téléchargement .zip).

Un guide rapide : pour extraire les facteurs d'inflation de la variance (VIF), exécutez leur code HighStatLib.r et utilisez la fonction corvif. La fonction nécessite une trame de données avec uniquement les prédicteurs (donc, par exemple, df = data.frame(Dataset[,2:4])si vos données sont stockées dans Dataset avec les prédicteurs des colonnes 2 à 4.

lithique
la source
1

Peut-être que la qr()fonction fonctionnera. Si Xc'est votre trame de données ou matrice, vous pouvez utiliser qr(X)$pivot. Par exemple, qr(X)$pivot= c(1, 2, 4, 5, 7, 8, 3, 6)alors les colonnes 3 et 6 sont la variable multicollinéaire.

JunhuiLi
la source
1

Pour évaluer la multicolinéarité entre les prédicteurs lors de l'exécution de la fonction de dragage (package MuMIn), incluez la fonction max.r suivante comme argument "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))
  }
}

exécutez simplement dredge en spécifiant le nombre de variables prédictives et en incluant la fonction 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

Cela fonctionne pour les modèles lme4. Pour les modèles nlme, voir: https://github.com/rojaff/dredge_mc

r.jaffe
la source
1

Le VIF (variance inflation factor) peut être mesuré simplement par:

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