[R-sig-ME] Collinearity tests (e.g. VIF) for glmmTMB package

D. Rizopoulos d@r|zopou|o@ @end|ng |rom er@@mu@mc@n|
Tue May 14 12:32:57 CEST 2019


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:
> 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 using r-project.org] Im
> Auftrag von Fox, John
> Gesendet: Freitag, 10. Mai 2019 00:31
> An: Ben Bolker <bbolker using gmail.com>; Williamson, Michael
> <michael.williamson using kcl.ac.uk>
> Cc: r-sig-mixed-models using 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:
> 
>> car:::vif.merMod
> 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
> 
>> -----Original Message-----
>> From: Ben Bolker [mailto:bbolker using gmail.com]
>> Sent: Thursday, May 9, 2019 3:31 PM
>> To: Fox, John <jfox using mcmaster.ca>; Williamson, Michael
>> <michael.williamson using kcl.ac.uk>
>> Cc: r-sig-mixed-models using 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?
>>
>>>   but I don't see where to recover the necessary information about
>>> the
>> 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:
>>> 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 ------------
>>>
>>>> library(glmmTMB)
>>>> example("glmmTMB")
>>>> v <- vcov(m3)
>>>> v
>>> 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 ------------
>>>
>>>> diag(solve(cov2cor(v[[1]][-1, -1])))
>>>     sppPR    sppDM  sppEC-A  sppEC-L sppDES-L    sppDF  minedno
>>> 1.154418 1.918674 1.340247 2.317812 2.344363 1.767413 1.006961
>>>
>>>> diag(solve(cov2cor(v[[2]][-1, -1])))
>>>     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/
>>>
>>>
>>>
>>>> -----Original Message-----
>>>> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-
>>>> project.org] On Behalf Of Williamson, Michael via
>>>> R-sig-mixed-models
>>>> Sent: Thursday, May 9, 2019 11:13 AM
>>>> To: r-sig-mixed-models using r-project.org
>>>> Subject: [R-sig-ME] Collinearity tests (e.g. VIF) for glmmTMB
>>>> package
>>>>
>>>> Good Afternoon,
>>>>
>>>> I've been running a few generalised linear mixed models on my data.
>>>> Due to convergence issues, down to the size of the data set, I was
>>>> recommended to switch to the glmmTMB package from the glmer
>>>> function in lme4.. The models are running much better now with no
>>>> more convergence issues.
>>>>
>>>> I'm looking to test the collinearity of my models, but the VIF
>>>> function in the car package does not work with the glmmTMB package.
>>>> Does anyone know of any packages or functions that can be used to
>>>> calculate collinearity from model outputs generated by glmmTMB?
>>>>
>>>> Many thanks,
>>>>
>>>> Mike Williamson
>>>>
>>>> Email:
>>>> michael.williamson using kcl.ac.uk<mailto:michael.williamson using kcl.ac.uk>
>>>>
>>>>
>>>>
>>>>
>>>> 	[[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models using r-project.org mailing list
>>>> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C949964ad1c114eb5f8fd08d6d8568a0e%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C636934263533396201&sdata=uSPgGHD4Kd%2FBi%2F%2BmCRYDeSzdKiHGaeqLxVuqzMlKbPU%3D&reserved=0
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C949964ad1c114eb5f8fd08d6d8568a0e%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C636934263533396201&sdata=uSPgGHD4Kd%2FBi%2F%2BmCRYDeSzdKiHGaeqLxVuqzMlKbPU%3D&reserved=0
> 
> --
> 
> _____________________________________________________________________
> 
> Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg | https://eur01.safelinks.protection.outlook.com/?url=www.uke.de&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C949964ad1c114eb5f8fd08d6d8568a0e%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C636934263533396201&sdata=7TTw4wZTddjd%2BvLuu3Lw1MpPzFI5CurHZiliTy10ke8%3D&reserved=0
> Vorstandsmitglieder: Prof. Dr. Burkhard Göke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Prölß, Marya Verdel
> _____________________________________________________________________
> 
> SAVE PAPER - THINK BEFORE PRINTING
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C949964ad1c114eb5f8fd08d6d8568a0e%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C636934263533396201&sdata=uSPgGHD4Kd%2FBi%2F%2BmCRYDeSzdKiHGaeqLxVuqzMlKbPU%3D&reserved=0
> .
> 

-- 
Dimitris Rizopoulos
Professor of Biostatistics
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014
Web (personal): http://www.drizopoulos.com/
Web (work): http://www.erasmusmc.nl/biostatistiek/
Blog: http://iprogn.blogspot.nl/


More information about the R-sig-mixed-models mailing list