power.R

analysis <- function(x,y) {
	#cor.test(x,y,method="kendall")$p.value
	cor.test(x,y,method="spearman")$p.value
	#unlist(summary(aov(lm(y ~ x))))["Pr(>F)1"]
}

pvals <- c()
pvals.binned <- c()
for(reps in 1:200) {
	year <- 1:65
	x <- 11.2 + (0.8 / 65) * year + rnorm(year, 0, 0.5)
	y <- 15.5 - (0.3 / 0.8) * x + rnorm(x, 0, 0.3)
	xb <- apply(matrix(x, ncol=13), 2, mean)
	yb <- apply(matrix(y, ncol=13), 2, max)
	pvals <- c(pvals, analysis(x,y))
	pvals.binned <- c(pvals.binned, analysis(xb,yb))
}

par(mfrow=c(2,2))
hist(pvals, xlim=c(0,1), breaks=20, ylim=c(0,200))
plot(y ~ x, main="Typical data raw")
hist(pvals.binned, xlim=c(0,1), breaks=20, ylim=c(0,200))
plot(yb ~ xb, main="Typical data binned")