[R-sig-ME] Collinearity tests (e.g. VIF) for glmmTMB package
Fox, John
j|ox @end|ng |rom mcm@@ter@c@
Tue May 14 14:58:15 CEST 2019
Dear Daniel,
> -----Original Message-----
> From: Daniel Lüdecke [mailto:d.luedecke using uke.de]
> Sent: Tuesday, May 14, 2019 6:26 AM
> To: Fox, John <jfox using mcmaster.ca>; 'Ben Bolker' <bbolker using gmail.com>;
> 'Williamson, Michael' <michael.williamson using kcl.ac.uk>
> Cc: r-sig-mixed-models using r-project.org
> Subject: AW: [R-sig-ME] Collinearity tests (e.g. VIF) for glmmTMB package
>
> 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.
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
>
> 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://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> _______________________________________________
> R-sig-mixed-models using 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
More information about the R-sig-mixed-models
mailing list