[R-sig-ME] p values with lmer (mcmcpvalue error)

[Douglas Bates]

On Sat, Apr 16, 2011 at 4:58 AM, Kaitlyn Gaynor <kaitlyng at gmail.com> wrote:
> Hello,

> I am currently using lmer to make a fairly large model with 9 fixed
> factors (8 of them with 2-4 levels, and 1 of them continuous) and 2
> random factors, >1000 observations.  I realize there is a lot of
> debate about p-values for parameter estimates but they’ve been
> requested by a journal so I’m trying to calculate them.  I have been
> trying to obtain p-values using the mcmcpvalue function:

> fit<-lmer(vig~Act+Height+FD+IGE+Pos.For+PR+Kin+PDDO+Rank+(1|Grp/Subj),
> data=vigilance)
> mcmcpvalue <- function(samp)
> { std <- backsolve(chol(var(samp)),
> cbind(0, t(samp)) - colMeans(samp),
> transpose = TRUE)
> sqdist <- colSums(std * std)
> sum(sqdist[-1] > sqdist[1])/nrow(samp) }
> markov1 <- mcmcsamp(fit, 5000)
> HPDinterval(markov1)

> mcmcpvalue(as.matrix(markov1[,1]))
> mcmcpvalue(as.matrix(markov1[,2]))
> ...
>
> HPDinterval generates CIs for each treatment, but mcmcpvalue function
> produces the following error messages:
>
> Error in markov1[, 1] : object of type 'S4' is not subsettable
> Error in as.matrix(markov1[, 1]) :
>  error in evaluating the argument 'x' in selecting a method for
> function 'as.matrix'
>
> Any help would be greatly appreciated!

With over 1000 observations I would feel comfortable quoting p-values
from likelihood ratio tests.  First switch to the maximum likelihood
(ML) estimation method rather than REML estimation then fit the model
without a particular term and evaluate the likelihood ratio (which is
somewhat confusingly produced by the anova() function).  That is

fitML <- update(fit, REML=FALSE)
fit1 <- update(fit, . ~ . - Act)
anova(fit1, fitML)
...