Extracting p-values from an glmmML fitted model

spunky

spunky

Can't make spagetti
#1
#1
So, I downloaded the glmmML R package to get p-values for both fixed and random effects but I don't know how to extract said p-values.

Example:

Code: 
id <- factor(rep(1:20, rep(5, 20)))
y <- rbinom(100, prob = rep(runif(20), rep(5, 20)), size = 1)
x <- rnorm(100)
dat <- data.frame(y = y, x = x, id = id)
mod<-glmmML(y ~ x, data = dat, cluster = id)

>mod
Call:  glmmML(formula = y ~ x, data = dat, cluster = id) 


               coef se(coef)      z Pr(>|z|)
(Intercept) -1.3840   0.9609 -1.440    0.150
x            0.6809   0.3459  1.968    0.049

Scale parameter in mixing distribution:  3.436 gaussian 
Std. Error:                              0.6324 

        LR p-value for H_0: sigma = 0:  9.567e-10 

Residual deviance: 96.2 on 97 degrees of freedom        AIC: 102.2
I've been using my usual bag of tricks with str(mod), unclass(mod), coef(mod), etc. but, although I can get the fixed effects coefficients. I'd really like to get what's on the Pr(>|z|) column and that p value of 9.567e-10 for the variance component....

....

Help? :)
 
Dason

Dason

Ambassador to the humans
#2
#2
Code: 
> mod

Call:  glmmML(formula = y ~ x, data = dat, cluster = id)


                coef se(coef)       z Pr(>|z|)
(Intercept) -0.37000   0.4103 -0.9017    0.367
x            0.04571   0.2414  0.1893    0.850

Scale parameter in mixing distribution:  1.432 gaussian
Std. Error:                              0.4267

        LR p-value for H_0: sigma = 0:  0.0004959

Residual deviance: 126.3 on 97 degrees of freedom       AIC: 132.3
> class(mod)
[1] "glmmML"
> print.glmmML
function (x, digits = max(3, getOption("digits") - 3), na.print = "",
    ...)
{
    cat("\\nCall: ", deparse(x$call), "\\n\\n")
    savedig <- options(digits = digits)
    on.exit(options(savedig))
    coef <- x$coefficients
    se <- x$coef.sd
    tmp <- cbind(coef, se, coef/se, signif(1 - pchisq((coef/se)^2,
        1), digits - 1))
    dimnames(tmp) <- list(names(coef), c("coef", "se(coef)",
        "z", "Pr(>|z|)"))
    cat("\\n")
    prmatrix(tmp)
    cat("\\nScale parameter in mixing distribution: ", x$sigma,
        x$prior, "\\n")
    cat("Std. Error:                             ", x$sigma.sd,
        "\\n")
    pv <- 0.5 * pchisq(x$cluster.null.deviance - x$deviance,
        df = 1, lower.tail = FALSE)
    cat("\\n        LR p-value for H_0: sigma = 0: ", pv, "\\n")
    if (x$boot) {
        cat("\\n Bootstrap p-value for H_0: sigma = 0: ", x$bootP,
            "(", x$boot, ")\\n")
    }
    cat("\\nResidual deviance:", format(signif(x$deviance, digits)),
        "on", x$df.residual, "degrees of freedom", "\\tAIC:",
        format(signif(x$aic, digits)), "\\n")
}
<bytecode: 0x0000000036be2a48>
<environment: namespace:glmmML>
> # Note the lines starting with pv <-
> pv <- 0.5 * pchisq(mod$cluster.null.deviance - mod$deviance, df = 1, lower.tail = FALSE)
> pv
[1] 0.0004959339
Now that you know where to look I leave the exercise of extracting the individual p values to the reader.
 
Last edited:
spunky

spunky

Can't make spagetti
#5
#5
You must log in or register to reply here.