View source: R/stat_test2.2018ext.R
extended : seems like it's working
1 | test2.2018ext(x, label)
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | ## small test for CRAN submission
dat1 <- matrix(rnorm(60, mean= 1), ncol=2) # group 1 : 30 obs of mean 1
dat2 <- matrix(rnorm(50, mean=-1), ncol=2) # group 2 : 25 obs of mean -1
dmat <- as.matrix(dist(rbind(dat1, dat2))) # Euclidean distance matrix
lab <- c(rep(1,30), rep(2,25)) # corresponding label
test2.2018ext(dmat, lab) # run the code !
## Not run:
## WARNING: computationally heavy.
# Let's compute empirical Type 1 error at alpha=0.05 : should not be on 45 line
niter = 496
pvals1 = rep(0,niter) # 2018 method
pvals2 = rep(0,niter) # extended
for (i in 1:niter){
dat1 = matrix(rnorm(60*3,mean= 1), ncol=3)
dat2 = matrix(rnorm(40*3,mean=-1), ncol=3)
dat = rbind(dat1, dat2)
lab = c(rep(1,60), rep(2,40))
pvals1[i] = test2.2018bbw(as.matrix(dist(dat)), lab)$p.value
pvals2[i] = test2.2018ext(as.matrix(dist(dat)), lab)$p.value
print(paste("iteration ",i," complete..",sep=""))
}
# Visualize the above at multiple significance levels
alphas = seq(from=0.001, to=0.999, length.out=100)
errors1 = rep(0,100)
errors2 = rep(0,100)
for (i in 1:100){
errors1[i] = sum(pvals1<=alphas[i])/niter
errors2[i] = sum(pvals2<=alphas[i])/niter
}
par(mfrow=c(1,2))
plot(alphas, errors1, "b", main="Type 1 : 2018 method",
xlab="alpha", ylab="error", lwd=2)
abline(0,1, lwd=2, col="red")
plot(alphas, errors2, "b", main="Type 1 : extension",
xlab="alpha", ylab="error", lwd=2)
abline(0,1, lwd=2, col="red")
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.