Skip to content

large discrepancy between p-values from mcmc-sampling and t-test

2 messages · Martin Eklund, David Duffy

Original
Martin Eklund
# 
Dear everyone,

My question concerns a large difference between p-values based on mcmc-sampling from a lmer model and those based on t-tests. I am fully aware of the fact that p-values based on t-tests are anti-conservative and that the mcmc-based ones therefore make more sense. However, the following example makes me slightly confused:

Data:

strata,method,y
1,A,1.023
2,A,1.051
3,A,1.025
4,A,1.03
5,A,1.055
6,A,1.116
7,A,1.108
8,A,1.245
9,A,0.983
10,A,1.174
1,B,1.086
2,B,1.074
3,B,1.095
4,B,1.13
5,B,1.089
6,B,1.186
7,B,1.083
8,B,1.293
9,B,1.015
10,B,1.225


Code:

library(lme4)
library(languageR)
test <- read.csv(file="testData.csv")	# Assuming the data above has been saved in a file called testData.csv
t.test(test$y[1:10], test$y[11:20], paired=T)
wilcox.test(test$y[1:10], test$y[11:20], paired=T)
M <- lmer(y ~ method+(1|strata), data=test)
pvalues <- pvals.fnc(M)
pvalues$fixed
barplot(test$y[1:10]-test$y[11:20])
binom.test(9, 10)	# 9 "successes" out of 10 trials, as seen in the barplot from the previous line

Note the significant results in all tests, except for the mcmc-based p-values in the lmer model. Even the sign test in the end result in a significant result, despite have relatively bad power (compared to the t-test and Wilcoxon test). My questions are: What may cause this large discrepancy in p-values based on the mcmc-sampling and the other tests? Can it simply be down to the small sample size? Would you in this case still trust the mcmc-sampling despite the fact that all other tests disagree with the p-value based on the mcmc?

Thank you very much!

Martin.

P.S. I'm not entirely sure if this is the right forum for this question. I would appreciate any pointers about where to ask this question if it is off-topic here. Thanks!



R version 2.15.0 (2012-03-30)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] sv_SE.ISO8859-1/sv_SE.ISO8859-1/sv_SE.ISO8859-1/C/sv_SE.ISO8859-1/sv_SE.ISO8859-1

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

other attached packages:
[1] nlme_3.1-104     languageR_1.4    lme4_0.999375-42 Matrix_1.0-6    
[5] lattice_0.20-6  

loaded via a namespace (and not attached):
[1] grid_2.15.0   stats4_2.15.0 tcltk_2.15.0  tools_2.15.0
1 day later
David Duffy
#
On Mon, 28 May 2012, Martin Eklund wrote:

            

      
      
      
    

        
You are relying on mcmcsamp here, which like everything else has been 
affected by various revisions of the lme4 code.  I think some versions 
were giving wrong answers, which I'm pretty sure is the case here (I 
double-checked this running in JAGS, where the 95%HPD looks more 
sensible: est=0.047, 0.016-0.075).

Bugs model:

  model
    {
       for( k in 1 : P ) {
          for( i in 1 : N ) {
             Y[i , k] ~ dnorm(m[i , k], tau1)
             m[i , k] <- mu + T[i , k] * phi / 2 + delta[i]
             T[i , k] <- 2*k - 3
          }
       }
       for( i in 1 : N ) {
          delta[i] ~ dnorm(0.0, tau2)
       }
       tau1 ~ dgamma(0.001, 0.001) sigma1 <- 1 / sqrt(tau1)
       tau2 ~ dgamma(0.001, 0.001) sigma2 <- 1 / sqrt(tau2)
       mu ~ dnorm(0.0, 1.0E-6)
       phi ~ dnorm(0.0, 1.0E-6)
    }