[R] lmer and estimation of p-values: error with mcmcpvalue()

[Christoph Scherber]

Dear Henric,

Thanks, now it works; but how reliable are these estimates? Especially 
with p-values close to 0.05 it is of course important that the range of 
the estimates is not too large. I´ve just run several simulations, each 
of which yielding sometimes quite different p-values.

Best wishes
Christoph


##
Dr. rer. nat. Christoph Scherber
DNPW, Agroecology
University of Goettingen
Waldweg 26
D-37073 Goettingen
http://wwwuser.gwdg.de/~uaoe/Agroecology.html
+49-(0)551-39-8807



Henric Nilsson (Public) schrieb:
> Den Må, 2007-02-12, 13:58 skrev Christoph Scherber:
>> Dear all,
>>
>> I am currently analyzing count data from a hierarchical design, and I?ve
>> tried to follow the suggestions for a correct estimation of p-values as
>> discusssed at R-Wiki
>> (http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests&s=lme%20and%20aov).
>>
>> 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.
>>
>> ##
>> 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) }
>>
>> m1<-lmer(number_pollinators~logpatch+loghab+landscape_diversity+(1|site),quasipoisson)
>>
>> Generalized linear mixed model fit using Laplace
>> Formula: number_pollinators ~ logpatch + loghab + landscape_diversity +
>>      (1 | site)
>>     Data: primula
>>   Family: quasipoisson(log link)
>>     AIC   BIC logLik deviance
>>   84.83 93.75 -37.42    74.83
>> Random effects:
>>   Groups   Name        Variance Std.Dev.
>>   site     (Intercept) 0.036592 0.19129
>>   Residual             1.426886 1.19452
>> number of obs: 44, groups: site, 15
>>
>> Fixed effects:
>>                      Estimate Std. Error t value
>> (Intercept)          -0.4030     0.6857 -0.5877
>> logpatch              0.1091     0.0491  2.2228
>> loghab                0.0875     0.0732  1.1954
>> landscape_diversity   1.0234     0.4850  2.1099
>>
>> Correlation of Fixed Effects:
>>              (Intr) lgptch loghab
>> logpatch     0.091
>> loghab      -0.637 -0.121
>> lndscp_dvrs -0.483 -0.098 -0.348
>>
>>
>> markov1=mcmcsamp(m1,5000)
>> HPDinterval(markov1)
>>
>> mcmcpvalue(as.matrix(markov1)[,1])
> 
> Try `mcmcpvalue(as.matrix(markov1[,1]))'.
> 
> 
> HTH,
> Henric
> 
> 
> 
>> Error in colMeans(samp) : 'x' must be an array of at least two dimensions
>>
>> ______________________________________________
>> R-help at stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
> 
> 
> .
>