[R-sig-ME] Collinearity diagnostics for (mixed) multinomial models
John Fox
j|ox @end|ng |rom mcm@@ter@c@
Mon Feb 28 01:23:48 CET 2022
Dear Juho,
I've now had a chance to think about this problem some more, and I
believe that the approach I suggested is correct. I also had an
opportunity to talk the problem over a bit with Georges Monette, who
coauthored the paper that introduced generalized variance inflation
factors (GVIFs). On the other hand, the results produced by
performance::check_collinearity() for multinomial logit models don't
seem to be correct (see below).
Here's an example, using the nnet::multinom() function to fit a
multinomial logit model, with alternative parametrizations of the LHS of
the model:
--------- snip -----------
> library(nnet) # for multinom()
> library(carData) # for BEPS data set
> # alternative ordering of the response levels:
> BEPS$vote1 <- factor(BEPS$vote, levels=c("Labour", "Liberal
Democrat", "Conservative"))
> levels(BEPS$vote)
[1] "Conservative" "Labour" "Liberal Democrat"
> levels(BEPS$vote1)
[1] "Labour" "Liberal Democrat" "Conservative"
> m <- multinom(vote ~ . - vote1, data=BEPS)
# weights: 33 (20 variable)
initial value 1675.383740
iter 10 value 1345.935273
iter 20 value 1150.956807
iter 30 value 1141.921662
iter 30 value 1141.921661
iter 30 value 1141.921661
final value 1141.921661
converged
> m1 <- multinom(vote1 ~ . - vote, data=BEPS)
# weights: 33 (20 variable)
initial value 1675.383740
iter 10 value 1280.439304
iter 20 value 1165.513772
final value 1141.921662
converged
> rbind(coef(m), coef(m1)) # compare coefficients
(Intercept) age economic.cond.national
economic.cond.household
Labour 0.9515214 -0.021913989 0.5575707
0.15839096
Liberal Democrat 1.4119306 -0.016810735 0.1810761
-0.01196664
Liberal Democrat 0.4604567 0.005102666 -0.3764928
-0.17036682
Conservative -0.9514466 0.021912305 -0.5575644
-0.15838744
Blair Hague Kennedy Europe
political.knowledge
Labour 0.8371764 -0.90775585 0.2513436 -0.22781308
-0.5370612
Liberal Democrat 0.2937331 -0.82217625 0.6710567 -0.20004624
-0.2034605
Liberal Democrat -0.5434408 0.08559455 0.4197027 0.02776465
0.3336068
Conservative -0.8371670 0.90778068 -0.2513735 0.22781092
0.5370545
gendermale
Labour 0.13765774
Liberal Democrat 0.12640823
Liberal Democrat -0.01125898
Conservative -0.13764849
> c(logLik(m), logLik(m1)) # same fit to the data
[1] -1141.922 -1141.922
> # covariance matrices for coefficients:
> V <- vcov(m)
> V1 <- vcov(m1)
> cbind(colnames(V), colnames(V1)) # compare
[,1] [,2]
[1,] "Labour:(Intercept)" "Liberal
Democrat:(Intercept)"
[2,] "Labour:age" "Liberal Democrat:age"
[3,] "Labour:economic.cond.national" "Liberal
Democrat:economic.cond.national"
[4,] "Labour:economic.cond.household" "Liberal
Democrat:economic.cond.household"
[5,] "Labour:Blair" "Liberal
Democrat:Blair"
[6,] "Labour:Hague" "Liberal
Democrat:Hague"
[7,] "Labour:Kennedy" "Liberal
Democrat:Kennedy"
[8,] "Labour:Europe" "Liberal
Democrat:Europe"
[9,] "Labour:political.knowledge" "Liberal
Democrat:political.knowledge"
[10,] "Labour:gendermale" "Liberal
Democrat:gendermale"
[11,] "Liberal Democrat:(Intercept)"
"Conservative:(Intercept)"
[12,] "Liberal Democrat:age" "Conservative:age"
[13,] "Liberal Democrat:economic.cond.national"
"Conservative:economic.cond.national"
[14,] "Liberal Democrat:economic.cond.household"
"Conservative:economic.cond.household"
[15,] "Liberal Democrat:Blair" "Conservative:Blair"
[16,] "Liberal Democrat:Hague" "Conservative:Hague"
[17,] "Liberal Democrat:Kennedy" "Conservative:Kennedy"
[18,] "Liberal Democrat:Europe" "Conservative:Europe"
[19,] "Liberal Democrat:political.knowledge"
"Conservative:political.knowledge"
[20,] "Liberal Democrat:gendermale"
"Conservative:gendermale"
> int <- c(1, 11) # remove intercepts
> colnames(V)[int]
[1] "Labour:(Intercept)" "Liberal Democrat:(Intercept)"
> colnames(V1)[int]
[1] "Liberal Democrat:(Intercept)" "Conservative:(Intercept)"
> V <- V[-int, -int]
> V1 <- V1[-int, -int]
> age <- c(1, 10) # locate age coefficients
> colnames(V)[age]
[1] "Labour:age" "Liberal Democrat:age"
> colnames(V1)[age]
[1] "Liberal Democrat:age" "Conservative:age"
> V <- cov2cor(V) # compute coefficient correlations
> V1 <- cov2cor(V1)
> # compare GVIFs:
> c(det(V[age, age])*det(V[-age, -age])/det(V),
+ det(V1[age, age])*det(V1[-age, -age])/det(V1))
[1] 1.046232 1.046229
--------- snip -----------
For curiosity, I applied car::vif() and
performance::check_collinearity() to these models to see what they would
do. Both returned the wrong answer. vif() produced a warning, but
check_collinearity() didn't:
--------- snip -----------
> car::vif(m1)
age economic.cond.national economic.cond.household
15.461045 22.137772 16.693877
Blair Hague Kennedy
14.681562 7.483039 15.812067
Europe political.knowledge gender
6.502119 4.219507 2.313885
Warning message:
In vif.default(m1) : No intercept: vifs may not be sensible.
> performance::check_collinearity(m)
# Check for Multicollinearity
Low Correlation
Term VIF Increased SE Tolerance
age 1.72 1.31 0.58
economic.cond.national 1.85 1.36 0.54
economic.cond.household 1.86 1.37 0.54
Blair 1.63 1.28 0.61
Hague 1.94 1.39 0.52
Kennedy 1.70 1.30 0.59
Europe 2.01 1.42 0.50
political.knowledge 1.94 1.39 0.52
gender 1.78 1.33 0.56
> performance::check_collinearity(m1)
# Check for Multicollinearity
Low Correlation
Term VIF Increased SE Tolerance
age 1.19 1.09 0.84
economic.cond.national 1.42 1.19 0.70
economic.cond.household 1.32 1.15 0.76
Blair 1.50 1.22 0.67
Hague 1.30 1.14 0.77
Kennedy 1.19 1.09 0.84
Europe 1.34 1.16 0.75
political.knowledge 1.30 1.14 0.77
gender 1.23 1.11 0.81
--------- snip -----------
I looked at the code for vif() and check_collinearity() to see where
they went wrong. Both failed to handle the two intercepts in the model
correctly -- vif() thought there was no intercept and
check_collinearity() just removed the first intercept but not the second.
In examining the code for check_collinearity(), I discovered a couple of
additional disconcerting facts. First, part of the code seems to be
copied from vif.default(). Second, as a consequence,
check_collinearity() actually computes GVIFs rather than VIFs (and
doesn't reference either the Fox and Monette paper introducing GVIFs or
the car package) but doesn't seem to understand that, and, for example,
takes the squareroot of the GVIF (reported in the column marked
"Increased SE") rather than the 2p root (when there are p > 1
coefficients in a term).
Here's the relevant code from the two functions (where . . . denotes
elided lines) -- the default method for vif() and .check_collinearity(),
which is called by check_collinearity.default():
--------- snip -----------
> car:::vif.default
function (mod, ...)
{
. . .
v <- vcov(mod)
assign <- attr(model.matrix(mod), "assign")
if (names(coefficients(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)
. . .
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)
}
. . .
}
> performance:::.check_collinearity
function (x, component, verbose = TRUE)
{
v <- insight::get_varcov(x, component = component, verbose = FALSE)
assign <- .term_assignments(x, component, verbose = verbose)
. . .
if (insight::has_intercept(x)) {
v <- v[-1, -1]
assign <- assign[-1]
}
else {
if (isTRUE(verbose)) {
warning("Model has no intercept. VIFs may not be sensible.",
call. = FALSE)
}
}
. . .
terms <- labels(stats::terms(f[[component]]))
. . .
n.terms <- length(terms)
if (n.terms < 2) {
if (isTRUE(verbose)) {
warning(insight::format_message(sprintf("Not enough model
terms in the %s part of the model to check for multicollinearity.",
component)), call. = FALSE)
}
return(NULL)
}
R <- stats::cov2cor(v)
detR <- det(R)
. . .
for (term in 1:n.terms) {
subs <- which(assign == term)
. . .
result <- c(result, det(as.matrix(R[subs, subs])) *
det(as.matrix(R[-subs, -subs]))/detR)
. . .
}
. . .
}
--------- snip -----------
So, the upshot of all this is that you should be able to do what you
want, but not with either car::vif() or
performance::check_collinearity(). Instead, either write your own
function or do the computations in a script.
There's also a lesson here about S3 default methods: The fact that a
default method returns a result rather than throwing an error or a
warning doesn't mean that the result is the right answer.
I hope this helps,
John
On 2022-02-26 3:45 p.m., Juho Kristian Ruohonen wrote:
> Dear John W,
>
> Thank you very much for the tip-off! Apologies for not responding earlier
> (gmail apparently decided to direct your email right into the junk folder).
> I am very pleased to note that the package you mention does indeed work
> with *brms* multinomial models! Thanks again!
>
> Best,
>
> Juho
>
> pe 25. helmik. 2022 klo 19.23 John Willoughby (johnwillec using gmail.com)
> kirjoitti:
>
>> Have you tried the check_collinearity() function in the performance
>> package? It's supposed to work on brms models, but whether it will work on
>> a multinomial model I don't know. It works well on mixed models generated
>> by glmmTMB().
>>
>> John Willoughby
>>
>>
>> On Fri, Feb 25, 2022 at 3:01 AM <r-sig-mixed-models-request using r-project.org>
>> wrote:
>>
>>> Send R-sig-mixed-models mailing list submissions to
>>> r-sig-mixed-models using r-project.org
>>>
>>> To subscribe or unsubscribe via the World Wide Web, visit
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>> or, via email, send a message with subject or body 'help' to
>>> r-sig-mixed-models-request using r-project.org
>>>
>>> You can reach the person managing the list at
>>> r-sig-mixed-models-owner using r-project.org
>>>
>>> When replying, please edit your Subject line so it is more specific
>>> than "Re: Contents of R-sig-mixed-models digest..."
>>>
>>>
>>> Today's Topics:
>>>
>>> 1. Collinearity diagnostics for (mixed) multinomial models
>>> (Juho Kristian Ruohonen)
>>>
>>> ----------------------------------------------------------------------
>>>
>>> Message: 1
>>> Date: Fri, 25 Feb 2022 10:23:25 +0200
>>> From: Juho Kristian Ruohonen <juho.kristian.ruohonen using gmail.com>
>>> To: John Fox <jfox using mcmaster.ca>
>>> Cc: "r-sig-mixed-models using r-project.org"
>>> <r-sig-mixed-models using r-project.org>
>>> Subject: [R-sig-ME] Collinearity diagnostics for (mixed) multinomial
>>> models
>>> Message-ID:
>>> <
>>> CAG_dBVfZr1-P7Q3kbE8TGPm-_2sJixdGCHCtWM9Q9PEnd8ftZw using mail.gmail.com>
>>> Content-Type: text/plain; charset="utf-8"
>>>
>>> Dear John (and anyone else qualified to comment),
>>>
>>> I fit lots of mixed-effects multinomial models in my research, and I
>> would
>>> like to see some (multi)collinearity diagnostics on the fixed effects, of
>>> which there are over 30. My models are fit using the Bayesian *brms*
>>> package because I know of no frequentist packages with multinomial GLMM
>>> compatibility.
>>>
>>> With continuous or dichotomous outcomes, my go-to function for
>> calculating
>>> multicollinearity diagnostics is of course *vif()* from the *car*
>> package.
>>> As expected, however, this function does not report sensible diagnostics
>>> for multinomial models -- not even for standard ones fit by the *nnet*
>>> package's *multinom()* function. The reason, I presume, is because a
>>> multinomial model is not really one but C-1 regression models (where C
>> is
>>> the number of response categories) and the *vif()* function is not
>> designed
>>> to deal with this scenario.
>>>
>>> Therefore, in order to obtain meaningful collinearity metrics, my present
>>> plan is to write a simple helper function that uses *vif() *to calculate
>>> and present (generalized) variance inflation metrics for the C-1
>>> sub-datasets to which the C-1 component binomial models of the overall
>>> multinomial model are fit. In other words, it will partition the data
>> into
>>> those C-1 subsets, and then apply *vif()* to as many linear regressions
>>> using a made-up continuous response and the fixed effects of interest.
>>>
>>> Does this seem like a sensible approach?
>>>
>>> Best,
>>>
>>> Juho
>>>
>>>
>>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
web: https://socialsciences.mcmaster.ca/jfox/
More information about the R-sig-mixed-models
mailing list