2010/11/2 Chitra <cbban...@gmail.com>
>
> yes
>
> >
> >
> > Me too. So you want to do a MC test for Pearson's product-moment
> > correlation, right...?
>
>
So for sample sizes from 3 to about 10 we can use all permutations
[permn(combinat)]- test will be exact! (In our case 7!=5040)
lg<-"lightgreen"
g<-"green"
dg<-"darkgreen"
plot((gamma(1:31)),t="p",main="Suggested tests for r",ylab="Number of
permutations",xlab="n",lwd=2,col=c(rep(lg,10),rep(g,4),rep(dg,17)),log="yx",pch="-",cex=2)
legend(1,range(gamma(4:31))[2],c("exact","MC","cor.test"),col=c(lg,g,dg),pch="-",pt.cex=2)
abline(h=.Machine$integer.max,col=2,lty=3)
We use MC when the number of permutations is very large and we cannot use
them all. Beside, the difference between theoretical distribution for larger
samples >25 will be negligible.
Let's use your data:
> data
Qtot Itot
1 73 684
2 64 451
3 71 378
4 65 284
5 47 179
6 31 117
7 19 69
We get 0.01494540 from cor.test
> cor.test(data[,1],data[,2])
Let's write a function for our test, it might be something like:
cor.test.mc<-function(x,y,n=1e3){
our.data<-cbind(x,y)
if (!is.numeric(our.data[,1]) || !is.numeric(our.data[,2]))
stop("Only numeric variables are allowed.")
l<-length(our.data[,1])
if (l < 3)
stop("At least 3 samples are required.")
DNAME <- paste(deparse(substitute(x)), "and" ,deparse(substitute(y)))
samples<-unique(t(replicate(n,(sample(our.data[,1])))))
loop<-dim(samples)[1]
correlations<-rep(NA,loop)
for(i in 1:loop){
correlations[i]<-cor(our.data[,2],samples[i,])
}
observed<-cor(our.data[,1],our.data[,2])
GE<-sum(correlations>=observed)
LT<-sum(correlations<(-observed))
two.tailed.p<-(GE+LT)/loop
rea<-(loop/gamma(l+1))*100
RVAL <- list(statistic = c(r = observed), p.value = two.tailed.p, method =
"Monte Carlo Pearson's r test" ,
data.name = DNAME,samples=c(" Number of used unique
permutations"=loop),total=c("Percent of all possible
permutations"=round(rea,2)))
class(RVAL) <- "htest"
#But what kind of plot you wish to have - I don't know...
#hist(correlations,col="blue",xlab="r",xlim=c(-1,1),breaks=50)
return(RVAL)
gc()
}
cor.test.mc(data[,1],data[,2])
test<-cor.test.mc(data[,1],data[,2],6e4)
test
test$samples
test$total
#
And that's it. Our p-value is sum of 7/5040 (GE) and 61/5040 (LT).
You may also take a look @ library(MChtest).
Hope this helps!
--
Mi³ego dnia
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org 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.
