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.

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:

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 at gmail.com]
Sent: Thursday, May 9, 2019 3:31 PM
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?

  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 at 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 at 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 at kcl.ac.uk<mailto:michael.williamson at kcl.ac.uk>




	[[alternative HTML version deleted]]