[R-sig-ME] binomial fixed-effect p-values by simulation

[Douglas Bates]

On Mon, Aug 25, 2008 at 8:56 AM, Ben Bolker <bolker at ufl.edu> wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
>  My guess, based on Littell et al 2003, would be that
> something like  p %*% V^{-1} %*% p would give you a
> quadratic form that would be chi-squared distributed with
> rank(V) df, or (with an estimated scale parameter) F-distributed
> with (rank(V), n-whatever) df -- or at least nominally
> so, and if you're going to simulate anyway you're going
> to find out how it's _really_ distributed ...
>
> [I have a glmer fit named "zz", and a factor named "status"
> that I want to test all levels == 0 simultaneously]
>
> params <- grep("^status",names(fixef(zz)))
> fixef(zz)[params] %*% solve(vcov(zz)[params,params]) %*% fixef(zz)[params]

(covers face and runs away screaming).

Umm, please don't use grep on names to determine which coefficients
are associated with a given factor.  That's what the terms and assign
attributes are for.

You split the fixed effects according to assign, as in the code for
the anova method.  For example

> gm1 <- glmer(r2 ~ btype + situ + mode + Gender*Anger + (1|id) + (1|item), VerbAgg, binomial)
> str(terms(gm1))
Classes 'terms', 'formula' length 3 r2 ~ btype + situ + mode + Gender * Anger
  ..- attr(*, "variables")= language list(r2, btype, situ, mode, Gender, Anger)
  ..- attr(*, "factors")= int [1:6, 1:6] 0 1 0 0 0 0 0 0 1 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:6] "r2" "btype" "situ" "mode" ...
  .. .. ..$ : chr [1:6] "btype" "situ" "mode" "Gender" ...
  ..- attr(*, "term.labels")= chr [1:6] "btype" "situ" "mode" "Gender" ...
  ..- attr(*, "order")= int [1:6] 1 1 1 1 1 2
  ..- attr(*, "intercept")= int 1
  ..- attr(*, "response")= int 1
  ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
  ..- attr(*, "predvars")= language list(r2, btype, situ, mode, Gender, Anger)
  ..- attr(*, "dataClasses")= Named chr [1:6] "factor" "factor"
"factor" "factor" ...
  .. ..- attr(*, "names")= chr [1:6] "r2" "btype" "situ" "mode" ...
> attr(gm1 at X, "assign")
[1] 0 1 1 2 3 4 5 6

The assign attribute indicates that the first coefficient is the
intercept, the next two are associated with the first-order term
"btype", the next is associated with the first-order term "situ", ...,
the eighth is associated with the second-order term "Gender:Anger".

> Daniel Ezra Johnson wrote:
>> Sorry if this has been covered elsewhere, but if my interest is in
>> testing a single fixed effect _term_ (all coefficients at once) is
>> there an appropriate statistic to simulate for a binomial model?
>>
>> In other words, if I fit a linear model "glmodel" I can simulate one
>> of the F-statistics from anova(glmodel). If there's only one
>> coefficient for the term then F = t^2...
>>
>> If I have a "glmmodel" I can do anova(glmmodel) but I wanted to make
>> sure the F-statistic reported there was a sensible thing to look at
>> since it wasn't quite the square of the z-statistic in the simple
>> case.
>>
>> Maybe it doesn't have an F-distribution but it would still work well
>> as a single-number stand-in for the 'size of a fixed effect' in a
>> simulation...
>>
>> Thanks,
>> D
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.6 (GNU/Linux)
> Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org
>
> iD8DBQFIsroic5UpGjwzenMRAncaAJ9y+Fza6/kU3ya/Bkd699i81/mDPQCfW4DA
> YNHMEw6/cFePwAdvJN6mL7s=
> =x8iG
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>