[R-sig-ME] HPDinterval, and mcmcpvalue from R-Wiki discussion
Douglas Bates
bates at stat.wisc.edu
Sun Aug 10 21:16:24 CEST 2008
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
>
More information about the R-sig-mixed-models
mailing list