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

[Douglas Bates]

The short answer (and it has to be short as I am typing this in a
hotel lobby in Dortmund) is that I have changed the format of the
value of the mcmcsamp function and not correspondingly changed the
function to calculate p-values.


On Fri, Aug 8, 2008 at 8:32 PM, David Afshartous
<dafshartous at med.miami.edu> wrote:
>
> 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)
> }
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>