[R-sig-ME] binomial fixed-effect p-values by simulation
Ben Bolker
bolker at ufl.edu
Mon Aug 25 19:28:35 CEST 2008
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
Douglas Bates wrote:
> On Mon, Aug 25, 2008 at 8:56 AM, Ben Bolker <bolker at ufl.edu> wrote:
> 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".
>
Sorry about the running-away-screaming part.
I'm just happy that you don't disagree that vehemently with the
statistical (as opposed to the R-programming) part of my advice ...
thanks
Ben
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.6 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org
iD8DBQFIsuvCc5UpGjwzenMRAi0sAKCDy5ktn0LdgUNYSUCstEbpEqMa5gCfUzlt
gOv6ckTkE8V7dSMzPildRBI=
=vTxi
-----END PGP SIGNATURE-----
More information about the R-sig-mixed-models
mailing list