How to test and avoid multicollinearity in mixed linear model?

Question

I am currently running some mixed effect linear models.

I am using the package "lme4" in R.

My models take the form:

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

Before running my models, I checked for possible multicollinearity between predictors.

I did this by:

Make a dataframe of the predictors

dummy_df <- data.frame(predictor1, predictor2)

Use the "cor" function to calculate Pearson correlation between predictors.

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

If "correl_dummy_df" was greater than 0.80, then I decided that predictor1 and predictor2 were too highly correlated and they were not included in my models.

In doing some reading, there would appear more objective ways to check for multicollinearity.

Does anyone have any advice on this?

The "Variance Inflation Factor (VIF)" seems like one valid method.

VIF can be calculated using the function "corvif" in the AED package (non-cran). The package can be found at http://www.highstat.com/book2.htm. The package supports the following book:

Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. 2009. Mixed effects models and extensions in ecology with R, 1st edition. Springer, New York.

Looks like a general rule of thumb is that if VIF is > 5, then multicollinearity is high between predictors.

Is using VIF more robust than simple Pearson correlation?

Update

I found an interesting blog at:

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

The blogger provides some useful code to calculate VIF for models from the lme4 package.

I've tested the code and it works great. In my subsequent analysis, I've found that multicollinearity was not an issue for my models (all VIF values < 3). This was interesting, given that I had previously found high Pearson correlation between some predictors.

Answer 1

For VIF calculation usdm can also be package ( I need to install "usdm")

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

If VIF > 4.0 then I generally assume multicollinearity remove all those Predictor Variables before fitting them into my model

Answer 2

An update, since I found this question useful but can't add comments -

The code from Zuur et al. (2009) is also available via the supplementary material to a subsequent (and very useful) publication of their's in the journal Methods in Ecology and Evolution.

The paper - A protocol for data exploration to avoid common statistical problems - provides useful advice and a much needed reference for justifying VIF thresholds (they recommend a threshold of 3). The paper is here: http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2009.00001.x/full and the R code is in the supplementary materials tab (.zip download).

A quick guide: to extract variance inflation factors (VIF) run their HighStatLib.r code and use the function corvif. The function requires a data frame with just the predictors (so, for example, df = data.frame(Dataset[,2:4]) if your data are stored in Dataset with the predictors in columns 2 to 4.

Answer 3

VIF (variance inflation factor) can be measured simply by:

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

Answer 4

Maybe qr() function will work. If X is your data frame or matrix, you can use qr(X)$pivot. For example, qr(X)$pivot= c(1, 2, 4, 5, 7, 8, 3, 6), then column 3 and 6 is the multicollinear variable.

Answer 5

To assess multicollinearity between predictors when running the dredge function (MuMIn package), include the following max.r function as the "extra" argument:

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

then simply run dredge specifying the number of predictor variables and including the max.r function:

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

This works for lme4 models. For nlme models see: https://github.com/rojaff/dredge_mc