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

[Ben Bolker]

On 11-04-16 09:29 AM, Douglas Bates 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]))

  I haven't looked carefully, but does this work better for you?

mcmcpvalue <- function(samp,type="fixef")
{
  s <- t(slot(samp,type))
  std <- backsolve(chol(var(s)),
                   cbind(0, t(s)) - colMeans(s),
                   transpose = TRUE)
  sqdist <- colSums(std * std)
  sum(sqdist[-1] > sqdist[1])/nrow(s)
}

mcmcpvalue(markov1)

  I can imagine that the structure of "merMCMC" objects (i.e. the
results of mcmcsamp()) has changed since the mcmcpvalue() function was
written.

  To sort this kind of problem out for yourself, use str(markov1) to see
what is contained in the object.