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

[Kaitlyn Gaynor]

That works, thanks!  Though it only provides p-values for each factor,
not for each level of each factor.  This may be sufficient, but in the
past I've seen p-values presented along with the estimates and
t-values for each level of each factor (ex. for Act=Feed, Act=Rest ...
, not just for Act)  Any way to do this?

On 4/16/11, Douglas Bates <bates at stat.wisc.edu> wrote:
> 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)
> ...
>