cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,1,0,1,0,0,0,0,1,1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1,1,1,1,0,1,1,1)
dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0, 0,1,1,1,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0)
florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,0,0,1,0)
mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,0,1,1,1,1,0,0,1,1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1)
data<-cbind(cerana,dorsata,florea,mellifera)
raw_score<-sum(rowSums(data[,1:4])^2)
raw_score
n<-nrow(data)
randomisations<-NULL
for(i in 1:100000){
	d1<-sample(data[,1],n,replace=F)
	d2<-sample(data[,2],n,replace=F)
	d3<-sample(data[,3],n,replace=F)
	d4<-sample(data[,4],n,replace=F)
	shuffled_score<-sum((d1+d2+d3+d4)^2) 
	randomisations<-c(randomisations,shuffled_score)
	} 
hist(randomisations)
abline(v=raw_score,lwd=2,col=2)
text(528,12000,"observed",srt=270,col="red")
equals<-length(which(randomisations==raw_score))
below<-length(which(randomisations<raw_score))+0.5*equals
above<-length(which(randomisations>raw_score))+0.5*equals
text(480,16000,paste("less than = ",below,sep=""))
text(560,16000,paste("greater than = ",above,sep=""),xpd=NA)
text(560,12000,paste("P = ",round(above/(above+below),4)),col="blue")