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

[Ben Bolker]

-----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-----