[R-sig-ME] variance inflation factor

Fox, John jfox at mcmaster.ca
Sat Feb 24 01:17:18 CET 2018


Dear Martin,

I've already handled this in the development version of car, where I introduced a vif.merMod() method. The general point may still be valid.

Best,
 John

> -----Original Message-----
> From: Martin Maechler [mailto:maechler at stat.math.ethz.ch]
> Sent: Friday, February 23, 2018 4:48 PM
> To: Fox, John <jfox at mcmaster.ca>
> Cc: Ben Bolker <bbolker at gmail.com>; r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] variance inflation factor
> 
> >>>>> Fox, John <jfox at mcmaster.ca>
> >>>>>     on Thu, 22 Feb 2018 17:34:59 +0000 writes:
> 
>     > Hi Ben,
>     >> -----Original Message----- From: Ben Bolker
>     >> [mailto:bbolker at gmail.com] Sent: Thursday, February 22,
>     >> 2018 12:31 PM To: Fox, John <jfox at mcmaster.ca> Cc: Toni
>     >> Hernandez-Matias <ahmatias at gmail.com>;
>     >> r-sig-mixed-models at r- project.org Subject: Re: [R-sig-ME]
>     >> variance inflation factor
>     >>
>     >> OK. I guess we could also change lme4 (hard to see why it
>     >> would break anything, and I can check downstream packages
>     >> ...)
> 
>     > That's might open a can of worms. I suppose,
>     > alternatively, you could introduce a cov2cor() method for
>     > "dpoMatrix" objects, which probably would be safer.
> 
>     > Best, John
> 
> I strongly agree with John here.
> 
> Initial note:
> Look at the examples on the  ?lmer   help page, or just run
>   example(lmer)
> 
> to see how things work already, basically
> 
>    V <- vcov( lmer(....)  )
>    R <- as(V, "corMatrix")
> 
> Note that with a "dpoMatrix" it is implicitly known that you have a positive
> (semi)definite symmetric matrix and that *is* useful information for subsequent
> methods / algorithms... even though that information may rarely be made used
> of in current R code outside the Matrix package.
> 
> I did not look into 'car',  but just using lme4, cov2cor()  already works :
> 
> > require(lme4)
> 
> > (fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy))
> Linear mixed model fit by REML ['lmerMod']
> Formula: Reaction ~ Days + (Days | Subject)
>    Data: sleepstudy
> REML criterion at convergence: 1743.628
> Random effects:
>  Groups   Name        Std.Dev. Corr
>  Subject  (Intercept) 24.740
>           Days         5.922   0.07
>  Residual             25.592
> Number of obs: 180, groups:  Subject, 18 Fixed Effects:
> (Intercept)         Days
>      251.41        10.47
> > (V <- vcov(fm1))
> 2 x 2 Matrix of class "dpoMatrix"
>             (Intercept)      Days
> (Intercept)   46.574562 -1.451096
> Days          -1.451096  2.389463
> > cov2cor(V)
> 2 x 2 Matrix of class "dpoMatrix"
>             (Intercept)       Days
> (Intercept)   1.0000000 -0.1375534
> Days         -0.1375534  1.0000000
> > as(V, "corMatrix")
> 2 x 2 Matrix of class "corMatrix"
>             (Intercept)       Days
> (Intercept)   1.0000000 -0.1375534
> Days         -0.1375534  1.0000000
> > sessionInfo()
> R version 3.4.3 Patched (2018-02-19 r74280)
> Platform: x86_64-pc-linux-gnu (64-bit)
> Running under: Fedora 26 (Twenty Six)
> 
> [.........................]
> 
> attached base packages:
> [1] graphics  grDevices datasets  stats     utils     methods   base
> 
> other attached packages:
> [1] lme4_1.1-16    Matrix_1.2-12  fortunes_1.5-4 sfsmisc_1.1-2
> 
> loaded via a namespace (and not attached):
>  [1] minqa_1.2.4     MASS_7.3-48     compiler_3.4.3  tools_3.4.3
>  [5] Rcpp_0.12.15    splines_3.4.3   nlme_3.1-131.1  grid_3.4.3
>  [9] nloptr_1.0.4    lattice_0.20-35
> >
> ----------------------
> 
> So, I remain a bit puzzled on why you observe the problems you observe.
> May 'car' needs to import some stuff from Matrix or lme4  so these methods do
> work there, too ?
> 
> Martin Maechler
> ETH Zurich
> 
>     >> On Thu, Feb 22, 2018 at 12:24 PM, Fox, John <jfox at mcmaster.ca> wrote:
>     >> > Hi Ben,
>     >> >
>     >> >> -----Original Message-----
>     >> >> From: Ben Bolker [mailto:bbolker at gmail.com]
>     >> >> Sent: Thursday, February 22, 2018 11:53 AM
>     >> >> To: Fox, John <jfox at mcmaster.ca>
>     >> >> Cc: Toni Hernandez-Matias <ahmatias at gmail.com>; r-sig-mixed-
> models at r-
>     >> >> project.org
>     >> >> Subject: Re: [R-sig-ME] variance inflation factor
>     >> >>
>     >> >>   I don't know if there's a good reason that vcov() is a "dpoMatrix"
>     >> >> (from the Matrix package) rather than a plain old base-R matrix, but
>     >> >> as.matrix() converts it harmlessly (and shouldn't have any ill
>     >> >> effects on something that's already a matrix ...)
>     >> >
>     >> > Yes, I know that, but it prevents the default method for car::vif() from
> working
>     >> because cov2cor() checks for a matrix:
>     >> >
>     >> >> vif(mod)
>     >> > Error in cov2cor(v) : 'V' is not a square numeric matrix In addition:
>     >> > Warning message:
>     >> > In vif.default(mod) :
>     >> > Error in cov2cor(v) : 'V' is not a square numeric matrix
>     >> >
>     >> > It should be easy to get around this -- I'll probably just create a
>     >> > local, unexported version of cov2cor() in the car package, something
>     >> > like
>     >> >
>     >> > cov2cor <- function(V) stats::cov2cor(as.matrix(V))
>     >> >
>     >> > Best,
>     >> >  John
>     >> >
>     >> >>
>     >> >> On Thu, Feb 22, 2018 at 11:01 AM, Fox, John <jfox at mcmaster.ca>
> wrote:
>     >> >> > Dear Toni,
>     >> >> >
>     >> >> >> -----Original Message-----
>     >> >> >> From: Toni Hernandez-Matias [mailto:ahmatias at gmail.com]
>     >> >> >> Sent: Thursday, February 22, 2018 10:01 AM
>     >> >> >> To: Fox, John <jfox at mcmaster.ca>
>     >> >> >> Cc: r-sig-mixed-models at r-project.org
>     >> >> >> Subject: Re: [R-sig-ME] variance inflation factor
>     >> >> >>
>     >> >> >> Dear John,
>     >> >> >>
>     >> >> >> Thank you very much for your message.
>     >> >> >>
>     >> >> >> Did you know a bibliographic reference for this method,
>     >> >> >> specifically in the case of generalized linear mixed models?
>     >> >> >
>     >> >> > Sorry, I don't, but don't see why the general idea wouldn't apply
>     >> >> > here (as Ben
>     >> >> mentioned in his response). By the way, the more general vif()
>     >> >> function in the car package works for lme() in nlme but not for
>     >> >> lmer() or glmer() in lme4 (because the result returned by vcov() for
>     >> >> a "merMod" object isn't of class "matrix"). We should probably make it
>     >> work.
>     >> >> >
>     >> >> > Best,
>     >> >> >  John
>     >> >> >
>     >> >> >>
>     >> >> >> Thank you again,
>     >> >> >>
>     >> >> >> Antonio
>     >> >> >>
>     >> >> >> On Thu, Feb 22, 2018 at 3:51 PM, Fox, John <jfox at mcmaster.ca
>     >> >> >> <mailto:jfox at mcmaster.ca> > wrote:
>     >> >> >>
>     >> >> >>
>     >> >> >>       Dear Toni,
>     >> >> >>
>     >> >> >>       Yes, that looks reasonable, although I'm not sure why it's
>     >> >> >> necessary to test for more than one fixed-effect intercept at the
>     >> >> >> beginning of the fixed-effect coefficient vector and how an
>     >> >> >> intercept could be named "Intercept" rather than "(Intercept)".
>     >> >> >> Assuming that there's a reason that escapes me, here's a
>     >> >> >> simplified
>     >> >> >> version:
>     >> >> >>
>     >> >> >>       vif.lme <- function (fit) {
>     >> >> >>         ## adapted from rms::vif
>     >> >> >>         v <- vcov(fit)
>     >> >> >>         nam <- names(fixef(fit))
>     >> >> >>         ## exclude intercepts
>     >> >> >>         ns <- sum(nam == "Intercept" | nam == "(Intercept)")
>     >> >> >>         if (ns > 0) v <- v[-(1:ns), -(1:ns), drop = FALSE]
>     >> >> >>         diag(solve(cov2cor(v)))
>     >> >> >>       }
>     >> >> >>
>     >> >> >>       I hope this helps,
>     >> >> >>        John
>     >> >> >>
>     >> >> >>
>     >> >> >>       -----------------------------
>     >> >> >>       John Fox, Professor Emeritus
>     >> >> >>       McMaster University
>     >> >> >>       Hamilton, Ontario, Canada
>     >> >> >>       Web: socialsciences.mcmaster.ca/jfox/
>     >> >> >> <http://socialsciences.mcmaster.ca/jfox/>
>     >> >> >>
>     >> >> >>
>     >> >> >>
>     >> >> >>
>     >> >> >>       > -----Original Message-----
>     >> >> >>       > From: R-sig-mixed-models
>     >> >> >> [mailto:r-sig-mixed-models-bounces at r-
>     >> >> >> <mailto:r-sig-mixed-models-bounces at r->
>     >> >> >>       > project.org <http://project.org> ] On Behalf Of Toni
>     >> >> >> Hernandez- Matias
>     >> >> >>       > Sent: Thursday, February 22, 2018 7:28 AM
>     >> >> >>       > To: r-sig-mixed-models at r-project.org
>     >> >> >> <mailto:r-sig-mixed-models at r- project.org>
>     >> >> >>       > Subject: [R-sig-ME] variance inflation factor
>     >> >> >>       >
>     >> >> >>       > Dear all,
>     >> >> >>       >
>     >> >> >>       > Could some help me to calculate Variance Inflation Factors
>     >> >> >> of models fitted
>     >> >> >>       > with glmer and lmer ?
>     >> >> >>       >
>     >> >> >>       > I found the method I copy below but I am not sure it is correct:
>     >> >> >>       >
>     >> >> >>       > vif.lme <- function (fit) {
>     >> >> >>       >      ## adapted from rms::vif
>     >> >> >>       >      v <- vcov(fit)
>     >> >> >>       >      nam <- names(fixef(fit))
>     >> >> >>       >      ## exclude intercepts
>     >> >> >>       >      ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
>     >> >> >>       >      if (ns > 0) {
>     >> >> >>       >          v <- v[-(1:ns), -(1:ns), drop = FALSE]
>     >> >> >>       >          nam <- nam[-(1:ns)] }
>     >> >> >>       >      d <- diag(v)^0.5
>     >> >> >>       >      v <- diag(solve(v/(d %o% d)))
>     >> >> >>       >      names(v) <- nam
>     >> >> >>       >      v }
>     >> >> >>       >
>     >> >> >>       >
>     >> >> >>       >
>     >> >> >>       > ##la estimamos para el modelo 1
>     >> >> >>       >
>     >> >> >>       > vif.lme(mod1)
>     >> >> >>       >
>     >> >> >>       >
>     >> >> >>       > Thank you very much in advance,
>     >> >> >>       >
>     >> >> >>       > Antonio
>     >> >> >>       >
>     >> >> >>       >
>     >> >> >>       > --
>     >> >> >>       >
>     >> >> >>
>     >> *********************************************************
>     >> >> >>       >
>     >> >> >>       > Antonio Hernandez Matias
>     >> >> >>       >
>     >> >> >>       > Equip de Biologia de la Conservació
>     >> >> >>       > Departament de Biologia Evolutiva, Ecología i Ciències
>     >> >> >> Ambientals Facultat de
>     >> >> >>       > Biologia  i Institut de Recerca de la Biodiversitat (IRBio)
> Universitat de
>     >> >> >>       > Barcelona (UB) Av. Diagonal, 643
>     >> >> >>       > Barcelona      08028
>     >> >> >>       > Spain
>     >> >> >>
>     >> >> >>       > Telephone: +34-934035857 <tel:%2B34-934035857>
>     >> >> >> <+34%20934%2003%2058%2057>
>     >> >> >>       > FAX: +34-934035740 <tel:%2B34-934035740>
>     >> >> >> <+34%20934%2003%2057%2040>
>     >> >> >>       > e-mail: ahernandezmatias at ub.edu
>     >> >> >> <mailto:ahernandezmatias at ub.edu>
> 
>     > _______________________________________________
>     > R-sig-mixed-models at r-project.org mailing list
>     > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



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