[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