The pari from angeella

pari's People

Contributors

Stargazers

Watchers

pari's Issues

Bug in permTest

Hi Angela,

There seems to be an important bug in your calculation of two-sample Welch statistics: the test statistics you report are inconsistent (and much larger) than the ones obtained by t.test. See the example below with m=2 tests using n=100 observations, with no signal.

> n <- 100
> m <- 2
> groups <- rep(c(0, 1), each = n/2)
> 
> mat <- matrix(rnorm(m*n), m, n) 
> t.test(mat[1, ]~groups)

	Welch Two Sample t-test

data:  mat[1, ] by groups
t = 1.4507, df = 92.122, p-value = 0.1503
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1055058  0.6772781
sample estimates:
mean in group 0 mean in group 1 
      0.1640547      -0.1218315 

> out <- pARI::permTest(X = mat, label = groups, B = 10, rand = FALSE)
> out$pv
[1] 1.239105e-10 3.401634e-01
> out$Test
[1] 7.2533998 0.9592307

In order to spot this type of errors, I highly recommend testing your code, see https://r-pkgs.org/tests.html.

In this case a minimal test is to ensure that your function reproduces the results of stats::t.test. This is what we do in the sansSouci package, see e.g. https://github.com/pneuvial/sanssouci/blob/develop/tests/testthat/test_rowWelchTests.R

Recommend Projects