[R-sig-ME] variance inflation factor
Martin Maechler
maechler at stat.math.ethz.ch
Fri Feb 23 22:47:30 CET 2018
>>>>> 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