Skip to content

Collinearity tests (e.g. VIF) for glmmTMB package

9 messages · John Fox, Ben Bolker, Williamson, Michael +3 more

Original
John Fox
# 
Dear Mike,

I'm not sufficiently familiar with the objects produced by glmmTMB() to answer definitively, and I'm also not entirely sure why you want to check for collinearity, but maybe the following would help:

You can used vcov() to return the variances and covariances of coefficients in the various parts of the "glmmTMB" model. For example:

---------------- snip ------------
Conditional model:
            (Intercept)        sppPR        sppDM       sppEC-A      sppEC-L     sppDES-L       sppDF       minedno
(Intercept)  0.04245503 -0.012754751 -0.013349646 -0.0125136751 -0.013436038 -0.013225977 -0.01391389 -0.0305911919
sppPR       -0.01275475  0.077687602  0.011642383  0.0119168647  0.011903658  0.011843477  0.01185466  0.0013084323
sppDM       -0.01334965  0.011642383  0.020980164  0.0117137251  0.011868129  0.011728938  0.01171587  0.0015986374
sppEC-A     -0.01251368  0.011916865  0.011713725  0.0404883426  0.011904829  0.011680709  0.01185958  0.0009042868
sppEC-L     -0.01343604  0.011903658  0.011868129  0.0119048294  0.017500122  0.011744878  0.01192195  0.0016761527
sppDES-L    -0.01322598  0.011843477  0.011728938  0.0116807092  0.011744878  0.016968986  0.01186668  0.0015556516
sppDF       -0.01391389  0.011854661  0.011715873  0.0118595830  0.011921947  0.011866683  0.02370581  0.0021442905
minedno     -0.03059119  0.001308432  0.001598637  0.0009042868  0.001676153  0.001555652  0.00214429  0.0350573728

Zero-inflation model:
               zi~(Intercept)     zi~sppPR     zi~sppDM   zi~sppEC-A   zi~sppEC-L  zi~sppDES-L     zi~sppDF   zi~minedno
zi~(Intercept)     0.08027669 -0.055011989 -0.064230942 -0.056164325 -0.064230942 -0.066122481 -0.064230942 -0.028881293
zi~sppPR          -0.05501199  0.157151941  0.060172003  0.062766076  0.060172003  0.059563719  0.060172003 -0.009287683
zi~sppDM          -0.06423094  0.060172003  0.122669211  0.060357133  0.061653087  0.061956976  0.061653087  0.004639967
zi~sppEC-A        -0.05616432  0.062766076  0.060357133  0.135723657  0.060357133  0.059862868  0.060357133 -0.007546778
zi~sppEC-L        -0.06423094  0.060172003  0.061653087  0.060357133  0.122669211  0.061956976  0.061653087  0.004639967
zi~sppDES-L       -0.06612248  0.059563719  0.061956976  0.059862868  0.061956976  0.123808814  0.061956976  0.007497634
zi~sppDF          -0.06423094  0.060172003  0.061653087  0.060357133  0.061653087  0.061956976  0.122669211  0.004639967
zi~minedno        -0.02888129 -0.009287683  0.004639967 -0.007546778  0.004639967  0.007497634  0.004639967  0.043632782

---------------- snip ------------

In this case, there are two components to the model -- the conditional model and the zero-inflation model -- and I believe that they are independent, so you should be able to eliminate the intercept from each and compute VIFs for the other coefficients:

---------------- snip ------------
sppPR    sppDM  sppEC-A  sppEC-L sppDES-L    sppDF  minedno 
1.154418 1.918674 1.340247 2.317812 2.344363 1.767413 1.006961
zi~sppPR    zi~sppDM  zi~sppEC-A  zi~sppEC-L zi~sppDES-L    zi~sppDF  zi~minedno 
   1.503986    1.699895    1.614313    1.699895    1.707801    1.699895    1.079338

---------------- snip ------------

Of course, it would be nice to automate this and to compute generalized VIFs for terms with more than one coefficient, but I don't see where to recover the necessary information about the structures of the component models from the "glmmTMB" object.

I'm cc'ing Ben Bolker in case he has something to add (or correct).

I hope this helps,
 John

--------------------------------------
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
Web: socialsciences.mcmaster.ca/jfox/
Ben Bolker
# 
I'm not sure what you else need to know about the component structures?
structures of the component models from the "glmmTMB" object.

   I've been meaning to write a vif.glmmTMB method, but I was planning
to just add a "component" argument to make the user choose: the names of
the vcov() components are "cond" and "zi" (and there could be a "disp"
component if there's a non-trivial dispersion model ...)
On 2019-05-09 3:28 p.m., Fox, John wrote:
John Fox
# 
Hi Ben,

Take a look at the vif.merMod() method from the car package:
function(mod, ...) {
  if (any(is.na(fixef(mod)))) 
    stop ("there are aliased coefficients in the model")
  v <- as.matrix(vcov(mod))
  assign <- attr(model.matrix(mod), "assign")
  if (names(fixef(mod)[1]) == "(Intercept)") {
    v <- v[-1, -1]
    assign <- assign[-1]
  }
  else warning("No intercept: vifs may not be sensible.")
  terms <- labels(terms(mod))
  n.terms <- length(terms)
  if (n.terms < 2) stop("model contains fewer than 2 terms")
  R <- cov2cor(v)
  detR <- det(R)
  result <- matrix(0, n.terms, 3)
  rownames(result) <- terms
  colnames(result) <- c("GVIF", "Df", "GVIF^(1/(2*Df))")
  for (term in 1:n.terms) {
    subs <- which(assign == term)
    result[term, 1] <- det(as.matrix(R[subs, subs])) *
      det(as.matrix(R[-subs, -subs])) / detR
    result[term, 2] <- length(subs)
  }
  if (all(result[, 2] == 1)) result <- result[, 1]
  else result[, 3] <- result[, 1]^(1/(2 * result[, 2]))
  result
}

So for each model component, I'd need the "assign" attribute from the model matrix. My relatively cursory examination of objects produced by glmmTMB() didn't turn up whether (or where) this information is stored.

I agree that for glmmTMB() VIFs could be computed componentwise, with the user specifying which component(s) and perhaps a default of "all". 

Best,
 John
4 days later
Williamson, Michael
# 
Thanks for your input, and apologies for the delay in reply.

I just wanted to test to make sure there was no correlation between my explanatory variables.

That's really helpful much appreciated.

Mike

-----Original Message-----
From: Ben Bolker <bbolker at gmail.com> 
Sent: 09 May 2019 20:31
To: Fox, John <jfox at mcmaster.ca>; Williamson, Michael <michael.williamson at kcl.ac.uk>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: Collinearity tests (e.g. VIF) for glmmTMB package



  I'm not sure what you else need to know about the component structures?
structures of the component models from the "glmmTMB" object.

   I've been meaning to write a vif.glmmTMB method, but I was planning to just add a "component" argument to make the user choose: the names of the vcov() components are "cond" and "zi" (and there could be a "disp"
component if there's a non-trivial dispersion model ...)
On 2019-05-09 3:28 p.m., Fox, John wrote:
Daniel Lüdecke
# 
Hi John and Ben,

iirc, model.matrix() for glmmTMB objects returns an assign-attribute,
however, I think you can only return the model matrix for the
count-component, not for the zero-inflated component. I don't know if there
might be situations where you have models that have different variables in
the ZI-formula as compared to the count-formula, but in such cases, VIF
can't be calculated.

Best
Daniel

-----Urspr?ngliche Nachricht-----
Von: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] Im
Auftrag von Fox, John
Gesendet: Freitag, 10. Mai 2019 00:31
An: Ben Bolker <bbolker at gmail.com>; Williamson, Michael
<michael.williamson at kcl.ac.uk>
Cc: r-sig-mixed-models at r-project.org
Betreff: Re: [R-sig-ME] Collinearity tests (e.g. VIF) for glmmTMB package

Hi Ben,

Take a look at the vif.merMod() method from the car package:
function(mod, ...) {
  if (any(is.na(fixef(mod)))) 
    stop ("there are aliased coefficients in the model")
  v <- as.matrix(vcov(mod))
  assign <- attr(model.matrix(mod), "assign")
  if (names(fixef(mod)[1]) == "(Intercept)") {
    v <- v[-1, -1]
    assign <- assign[-1]
  }
  else warning("No intercept: vifs may not be sensible.")
  terms <- labels(terms(mod))
  n.terms <- length(terms)
  if (n.terms < 2) stop("model contains fewer than 2 terms")
  R <- cov2cor(v)
  detR <- det(R)
  result <- matrix(0, n.terms, 3)
  rownames(result) <- terms
  colnames(result) <- c("GVIF", "Df", "GVIF^(1/(2*Df))")
  for (term in 1:n.terms) {
    subs <- which(assign == term)
    result[term, 1] <- det(as.matrix(R[subs, subs])) *
      det(as.matrix(R[-subs, -subs])) / detR
    result[term, 2] <- length(subs)
  }
  if (all(result[, 2] == 1)) result <- result[, 1]
  else result[, 3] <- result[, 1]^(1/(2 * result[, 2]))
  result
}

So for each model component, I'd need the "assign" attribute from the model
matrix. My relatively cursory examination of objects produced by glmmTMB()
didn't turn up whether (or where) this information is stored.

I agree that for glmmTMB() VIFs could be computed componentwise, with the
user specifying which component(s) and perhaps a default of "all". 

Best,
 John
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Marya Verdel
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING
Dimitris Rizopoulos
# 
Hi Daniel,

In the development version of GLMMadaptive 
(https://github.com/drizopoulos/GLMMadaptive) that can fit similar 
models as glmmTMB, I have put a VIF() method (based on car::vif) that 
can calculate VIFs also for zero-inflated models. E.g.,

library("GLMMadaptive")
fm <- mixed_model(y ~ time + sex, random = ~ 1 | id, data = <your_data>, 
family = zi.negative.binomial(), zi_fixed = ~ sex, zi_random = ~ 1 | id)

VIF(fm) # for the fixed effects of the non-zero part
VIF(fm, type = "zi_fixed") # fixed effects zero part


Best,
Dimitris
On 5/14/2019 12:25 PM, Daniel L?decke wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      John Fox
    

    

      
      #
    

  


  

      
Dear Mike,

            

      
      
      
    

        
If that's all you want to do the you could just check directly, but "collinearity" in a mixed model is more complicated and vif() defines it using the correlations of the coefficients. In a linear model fit by OLS the two are the same but not more generally.

Best,
 John

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      John Fox
    

    

      
      #
    

  


  

      
Dear Daniel,

            

      
      
      
    

        
There is no intrinsic reason why the predictors in the two parts of the model have to be the same.

As I mentioned in a related response, it's not the correlations among the columns of the model matrix but among the coefficients that figures in VIFs computed for models other than linear models fit by OLS. In glmmTMB, the correlations among the coefficients will differ in different parts of the model even if the model matrices are the same.

Best,
 John

            

      
      
      
    

        

      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      d@iuedecke m@iii@g oii uke@de
    

    

      
      #
    

  


  

      
Dear John,

sorry for confusion. I wasn't thinking of calculating correlations from the
model-matrix, but of the "assign"-attribute. That's what you use
model.matrix() for, and glmmTMB is missing the option to get the
model-matrix for the zero-inflated part, hence you need to workaround this
in order to get the assignments.

Best
Daniel

-----Urspr?ngliche Nachricht-----
Von: Fox, John <jfox at mcmaster.ca> 
Gesendet: Dienstag, 14. Mai 2019 14:58
An: Daniel L?decke <d.luedecke at uke.de>; 'Ben Bolker' <bbolker at gmail.com>;
'Williamson, Michael' <michael.williamson at kcl.ac.uk>
Cc: r-sig-mixed-models at r-project.org
Betreff: RE: [R-sig-ME] Collinearity tests (e.g. VIF) for glmmTMB package

Dear Daniel,

            

      
      
      
    

        
might

            

      
      
      
    

        
ZI-

            

      
      
      
    

        
There is no intrinsic reason why the predictors in the two parts of the
model have to be the same.

As I mentioned in a related response, it's not the correlations among the
columns of the model matrix but among the coefficients that figures in VIFs
computed for models other than linear models fit by OLS. In glmmTMB, the
correlations among the coefficients will differ in different parts of the
model even if the model matrices are the same.

Best,
 John

            

      
      
      
    

        
model

            

      
      
      
    

        
"disp"

            

      
      
      
    

        

      
      
    

        

      
      
    

        
Uwe

            

      
      
      
    

        

      
      
    

        
--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Marya Verdel
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING