binomial fixed-effect p-values by simulation
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