[R-sig-ME] mcmcpvalue and error message

[Eva Frei]

Dear all,

I am currently analyzing data from a mixed model with lmer.
I`ve tried to follow the suggestions for a correct estimation of p-values as discussed at R-Wiki (http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests&s=lme%20and%20aov).

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)
}

Example that works from D. Bates on RWiki:
fm1Adg <- lmer(adg ~ InitWt*Treatment - 1 + (1|Block), AvgDailyGain)
AdgS1 <- mcmcsamp(fm1Adg, 50000)
library(coda)
HPDinterval(AdgS1)
mcmcpvalue(as.matrix(AdgS1)[, 6:8])

In my example, I must use the HPDinterval of the library(lme4) because of S4 classes I think, and that my Output of HPDinterval is splitted into subsets $fixef, $ST and $sigma.
And when I try to extract information exspecially from one subset with
@fixef, that gives me not an mcmcpvalue, but the following error message:

My example:
library(lme4)
fm1Marie<-lmer(a ~ covar + b*c + (1|d) + (1|d:e),na.action=na.omit)
MarieS1<-mcmcsamp(fm1Marie, 50000)
HPDinterval(MarieS1) #gives me an Output with subsets
HPDinterval(MarieS1 at fixef)
str(MarieS1 at fixef)
mcmcpvalue(as.matrix(MarieS1 at fixef)[3:5,])

Error in chol(val(samp)): error in evaluating the argument 'x' in selecting a method for function 'chol'. 


My questions:
Is this because I use a newer version of HPDinterval than D. Bates?

Is this because the matrix is not of positive infinites in my case? 

Can I write an other mcmcpvalue that allows non-positive values in the matrix?

How I can even so became mcmcpvalues for my fixed effects?

Or can I suppress the slots $ST and $sigma? I`m only interessted in the p-values of my covariable covar and the fixed effects b, c, b:c!

Thanks for any advice!


Marieantoine

PhD. student, Switzerland



>>
>> However, I have the problem that my model only consists of parameters
>> with just 1 d.f. (intercepts, slopes), so that the "mcmcpvalue" function
>> defined below obviously produces error messages.
>>
>> How can I proceed in estimating the p-values, then?
>>
>> I very much acknowledge any suggestions.
>>
>> Best regards
>> Christoph.

-- 
GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT!