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

[Reinhold Kliegl]

The test applies to two-level factors. You could convert factors with
n  levels to  a set of n-1 predictors and enter the predictors instead
of the factor in the model. Then you can apply an anova() to test each
degree of freedom separately.
Reinhold Kliegl

On Sat, Apr 16, 2011 at 6:36 PM, Kaitlyn Gaynor <kaitlyng at gmail.com> wrote:
> 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)
>> ...
>>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>