[R-sig-ME] HPDinterval, and mcmcpvalue from R-Wiki discussion

[David Afshartous]

All,

I am trying to apply HPDinterval and mcmcpvalue as shown in the discussion
at R-Wiki  
((http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests&s=lme%20and%
20aov).

Three questions:
1) why am I getting the error message for HPDinterval given that I'm
following the exact steps in the discussion?
2) is it normal that my p-value (p=.72) differs that much from the one in
the discussion (p=.46)?
3) why isn't mcmcpvalue written to work on a single column? (see below)

Thanks,
David



> sessionInfo()
R version 2.7.1 (2008-06-23)
i386-apple-darwin8.10.1

locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] SASmixed_0.4-4     coda_0.13-2        Hmisc_3.4-3        nlme_3.1-89
lme4_0.999375-24   Matrix_0.999375-11
[7] lattice_0.17-8 

loaded via a namespace (and not attached):
[1] cluster_1.11.11 grid_2.7.1      tools_2.7.1
 data(AvgDailyGain, package = "SASmixed")
 (fm1Adg <- lmer(adg ~ InitWt*Treatment - 1 + (1|Block), AvgDailyGain))

> AdgS1 <- mcmcsamp(fm1Adg, 50000)
> library(coda)
> HPDinterval(AdgS1)
Error in UseMethod("HPDinterval") :
  no applicable method for "HPDinterval"

> mcmcpvalue(as.matrix(AdgS1)[, 6:8])
[1] 0.7231
> mcmcpvalue(as.matrix(AdgS1)[, 1])
Error in base::colMeans(x, na.rm = na.rm, dims = dims) :
  'x' must be an array of at least two dimensions
> 



mcmcpvalue <- function(samp)
{
   ## elementary version that creates an empirical p-value for the
   ## hypothesis that the columns of samp have mean zero versus a
   ## general multivariate distribution with elliptical contours.

   ## differences from the mean standardized by the observed
   ## variance-covariance factor
   std <- backsolve(chol(var(samp)),
                    cbind(0, t(samp)) - colMeans(samp),
                    transpose = TRUE)
   sqdist <- colSums(std * std)
   sum(sqdist[-1] > sqdist[1])/nrow(samp)
}