[R-sig-ME] Anova()'s Wald chisquare test for predictors with df > 1 in a logistic GLMM

[Fox, John]

Dear Juho,

Actually, for a main-effect-only model, there is no distinction between the "type-II" and "type-III" Wald tests computed by Anova(). It's only when there are terms (beyond the constant) related by marginality that the distinction arises.

Best,
 John 

> -----Original Message-----
> From: Juho Kristian Ruohonen [mailto:juho.kristian.ruohonen using gmail.com]
> Sent: Thursday, April 18, 2019 3:31 PM
> To: Fox, John <jfox using mcmaster.ca>
> Cc: Paul Johnson <paul.johnson using glasgow.ac.uk>; Shaul Oreg via R-sig-mixed-
> models <r-sig-mixed-models using r-project.org>
> Subject: Re: [R-sig-ME] Anova()'s Wald chisquare test for predictors with df > 1
> in a logistic GLMM
> 
> I was talking type-II tests of a main-effects model all along, which I apologize
> for not specifying at the outset. Anyway, this is now resolved. Many thanks to
> John for weighing in and to Paul for the pedagogical code snippets! (And here I
> I was assuming that the multi-df Wald test would be as simple as the single-
> parameter one. Joke's on me.)
> 
> Best,
> 
> J
> 
> 
> to 18. huhtik. 2019 klo 19.58 Fox, John (jfox using mcmaster.ca
> <mailto:jfox using mcmaster.ca> ) kirjoitti:
> 
> 
> 	Dear Paul,
> 
> 	The Wald() function suffices for the simple case you use, where there
> are no interactions. The trick in the more complex case of "type-II" tests for
> models with terms related by marginality is to generate the correct restriction
> matrix, which doesn't simply pick out certain coefficients.
> 
> 	Best,
> 	 John
> 
> 	  -------------------------------------------------
> 	  John Fox, Professor Emeritus
> 	  McMaster University
> 	  Hamilton, Ontario, Canada
> 	  Web: http::/socserv.mcmaster.ca/jfox
> <http://socserv.mcmaster.ca/jfox>
> 
> 	> On Apr 18, 2019, at 12:18 PM, Paul Johnson
> <paul.johnson using glasgow.ac.uk <mailto:paul.johnson using glasgow.ac.uk> >
> wrote:
> 	>
> 	>> I'm trying to figure out how the Wald chi-square test is performed
> for a
> 	>> nominal predictor with more than two levels in the context of a
> logistic
> 	>> GLMM.
> 	>
> 	> Keep in mind that Wald tests and LRTs for GLMMs come with health
> warnings, e.g.
> 	> https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#tests-
> of-effects-i.e.testing-that-several-parameters-are-simultaneously-zero
> 	>
> 	>> But how are they calculated for a multi-level nominal predictor?
> 	>>
> 	>> My (un)educated guess is that for each level of the nominal
> variable, the
> 	>> z-score is squared and the sum of these squares compared to the
> right-tail
> 	>> probability of the chi-squared distribution with DF equal to the
> number of
> 	>> levels of the predictor minus one. And indeed this square-the-z
> approach
> 	>> seems to correctly reproduce the results of Anova() for predictors
> with a
> 	>> single degree of freedom. But I can't make it reproduce the results
> of
> 	>> Anova() for predictors with more than one level. Hence my
> question: how is
> 	>> the test statistic calculated?
> 	>
> 	> When the null hypothesis is that more than one fixed effect is zero
> (as with a factor with >2 levels), the covariances of the fixed effects have to be
> be taken into account. The Wald statistic is
> 	>
> 	>    t(R %*% b) %*% solve(R %*% vc %*% t(R)) %*% (R %*% b)
> 	>
> 	> where R is a restriction matrix that picks out the contrasts you're
> interested in, b is the fixed effect estimates vector, and vc is the estimated
> covariance matrix of the fixed effects (and we’re assuming that as is usual the
> null hypothesis values of b are all zero). In the case where you’re testing a
> multilevel categorical variable, R %*% b is simply the vector of estimates for
> the k-1 non-reference levels and R %*% vc %*% t(R) is the k-1 X k-1 block of
> the covariance matrix for those k-1 fixed effects. If the covariances were all
> zero (e.g. by design — but I don’t think they can be for multilevel factors), the
> Wald statistic formula above would reduce to your sum of squares of z
> method.
> 	>
> 	> Below is an binomial glmer example with a Wald test function I
> found on the web years ago. It gives the same result as car::Anova but it’s
> easier to see the inner workings.
> 	>
> 	> Best wishes,
> 	> Paul
> 	>
> 	>
> 	> # wald z-test function
> 	> # (adapted from Stijn Ruiter's function
> http://stijnruiter.nl/blog/?p=309 - broken link)
> 	>
> 	> Wald<-
> 	>  function(object, # lme4 or glmmTMB fit
> 	>           R, # restrictions matrix. null hypothesis: R %*% fixef(object) = q
> 	>           q = NULL, # null values. by default q is a vector of zeroes
> 	>           gmmTMB.model = "cond")
> 	>  {
> 	>    require(lme4)
> 	>    if (!is.matrix(R)) stop("Restrictions must be a matrix")
> 	>    if(is.null(q)) q <- rep(0, nrow(R))
> 	>    b <- fixef(object)
> 	>    if(class(b) == "fixef.glmmTMB") b <- b[[gmmTMB.model]]
> 	>    vc <- vcov(object)
> 	>    if("vcov.glmmTMB" %in% class(vc)) vc <- vc[[gmmTMB.model]]
> 	>    w <- t(R %*% b - q) %*% solve(R %*% vc %*% t(R)) %*% (R %*% b -
> q)
> 	>    pw <- pchisq(w[1], length(q), lower.tail = FALSE)
> 	>    cat("*************\n* Wald Test *\n*************\n")
> 	>    cat("lme4 fixed effects:\n")
> 	>    print(fixef(object))
> 	>    cat("\nRestrictions:\n")
> 	>    print(R)
> 	>    cat("\nq = ",q)
> 	>    cat("\nChi-square:", round(w[1],3), " df = ", length(q))
> 	>    cat("\nProb x>chisq:", round(pw, 5), "\n")
> 	>    return(pw)
> 	>  }
> 	>
> 	>
> 	> # first example from ?glmer
> 	> library(lme4)
> 	> library(car)
> 	> (gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
> 	>              data = cbpp, family = binomial))
> 	>
> 	>
> 	> R <- cbind(0, diag(3))
> 	> # ...but ideally you want to generate the restriction matrix
> automatically
> 	>
> 	> Wald(gm1, R = R)
> 	> anova(gm1, update(gm1, ~ . - period))
> 	> Anova(gm1)
> 	>
> 	>
> 	> _______________________________________________
> 	> R-sig-mixed-models using r-project.org <mailto:R-sig-mixed-models using r-
> project.org>  mailing list
> 	> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
>