# How to test and avoid multicollinearity in mixed linear model?

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.

• (1) The AED package has been discontinued; instead, just source("http://www.highstat.com/Book2/HighstatLibV6.R") for the corvif function. (2) Hope to provide a real answer, but (a) I believe VIF takes multicollinearity into account (e.g. you may have three predictors, none of which has strong pairwise correlations, but the linear combination of A and B is strongly correlated with C) and (b) I have strong reservations about the wisdom of dropping collinear terms; see Graham Ecology 2003, doi:10.1890/02-3114 Jan 22, 2014 at 21:27
• Thanks Ben. I've updated my above post to include your suggestions. Jan 24, 2014 at 0:09
• @BenBolker, can you elaborate very briefly why you are against dropping collinear terms? I appreciate the reference but might also like a Cliff Notes version. Thanks! Mar 24, 2017 at 16:39
• correction in the Ben's response.. the URL is http://highstat.com/Books/BGS/GAMM/RCodeP2/HighstatLibV6.R Jul 24, 2018 at 21:26

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

• A bit addendum you can use thresold to filter variables like exclude all that show correlation above .4 as vifcor(vardata,th=0.4). Likewise you can use vifstep(vardata,th=10) to discard all greater than 10. Jan 22, 2017 at 8:22
• Doesn't work for HLM
– Mox
Jul 24, 2018 at 22:02

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.

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

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


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