[R-sig-ME] binomial fixed-effect p-values by simulation
Douglas Bates
bates at stat.wisc.edu
Mon Aug 25 17:23:05 CEST 2008
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
>
More information about the R-sig-mixed-models
mailing list