assoc <- function (poly.frame, polylist,response,permutations) 
{
tstat.max <- vector(mode="numeric",length=permutations) #to hold the max t from each permutation
tstat.perm <- vector(mode="numeric",length=length(polylist)) #to hold the t from each lm at each permutation
attach(poly.frame)
results <- numeric(length(polylist))

names(results) <- polylist

#first find tstat for each polymorphism
                for(j in polylist) {
                if (nlevels(get(j))>1) {
			results[j] <- cor.test(get(response),get(j))$statistic[[1]]
			#print(c("t statistic of correlation for ",j,cor.test(get(response),get(j))$statistic[[1]]),quote=false)
                        #tstat <- cor.test(get(response),get(j))$statistic#old
                        #print(c(j,tstat) , quote=F) #old 
			}
                else
                        print(c(j,"singularity"),quote=F)}


# now do permutations
for(i in 1:permutations) {#do the number of permutations listed in "permutations"
        perm <- sample(eval(parse(text=c("poly.frame$",response)))) #permute
                for(j in polylist) { #for each permutation regress on each marker
                if (nlevels(get(j))>1) tstat.perm[j] <-
		abs(cor.test(perm,get(j))$statistic)}
tstat.max[i] <- max(tstat.perm,na.rm=T) } #for each permutation take largest tstat
#tstat.max <- sort(tstat.max)
detach(poly.frame)

c(sort(abs(round(results,2))),"the experiment-wise P 0.1 threshold for t is: ",quantile(tstat.max,0.9),".05 ", quantile(tstat.max,.95)," .01 ",quantile(tstat.max,0.99)," .005",quantile(tstat.max,.995))

}

perm.sem <- function (poly.frame, polylist,response,permutations) { #function to calculate SEM by permutation
	attach(poly.frame)
	means <- array(dim=c(2,length(polylist),permutations),dimnames=list(gt=c("Cvi","Ler"),poly=polylist,perm=1:permutations))
	colnames(means) <- polylist
	for(i in 1:permutations) {#do the number of permutations listed in "permutations"
        perm <- sample(eval(parse(text=c("poly.frame$",response)))) #permute
        	for (p in polylist) {
        		means[,p,i] <- tapply(perm,get(p),mean,na.rm=T)
        		}
        }
        detach(poly.frame)
        apply(means,c(1,2),sd,na.rm=T)
      }


make.graph <- function(poly.frame,polylist,response,col=NULL,permutations=1000) {
	p.means <- matrix(nrow=2,ncol=length(polylist))
	colnames(p.means) <- polylist
	rownames(p.means) <- c("Cvi","Ler")
	if (col == "white") col <- c("white","black")
	p.sems <- perm.sem(poly.frame,polylist,response,permutations)
	for (p in polylist) {
		p.means[,p] <- tapply(poly.frame[,response],poly.frame[,p],mean,na.rm=T)
		}
		x <- barplot(p.means,beside=T,ylim=c(0,14),#max(p.means+p.sems,na.rm=T)),
			#col=col,main=response,density=c(-1,15),
			col=c(col,"white"),main=response,
			names.arg = unlist(strsplit(polylist,split=".",fixed=T))[seq(2,length(polylist)*2,2)],
			cex.names=1.5,cex.axis=1.5,cex.main=1.5)
		arrows(x,p.means+p.sems,x,p.means-p.sems,angle=90,code=3,length=.025)
		}
		
	
	
polys <- read.csv("PHYB_ALL_P1-P7 041607 for paper.csv")

polys$White.resid <-  NA

tmp.resid <- resid(lm(polys$White~polys$Blue+polys$FarRed))

for (i in names(tmp.resid)) { #this is necessary because there is a missing value in FarRed
	polys$White.resid[as.numeric(i)] <- tmp.resid[as.character(i)]
	}
		
summary(polys)

names(polys)

responses <- names(polys)[c(4:7,10)]

names(polys)[11:17] <- c("site.1","site.3", "site.4","site.5","site.7","site.8","site.12") # Change to Daniele's numbering system

sites <- names(polys)[c(11:17)]

for (i in responses){
  print(i)
  print(assoc(polys,sites,i,10000))
      }
      
sites <- sites[c(1,2,3,5,7)]
col.list <- unlist(strsplit("white blue red darkred grey30",split=" "))
names(col.list) <- responses

pdf(file="newbarcharts.pdf",width=10,height=7.5)
par(mfrow=c(2,3))
for(i in responses) {
	make.graph(polys,sites,i,col.list[i],10000)
	}
dev.off()	
	
#for paper, use grayscale:
pdf(file="newbarcharts.pdf",width=7.09,height=10.5)
par(mfrow=c(3,2))
for(i in responses) {
	make.graph(polys,sites,i,"grey50",10000)
	}
dev.off()

#check to see if entire effect can be explained by site 3
attach(polys)

newData<- data.frame(Red,site.1,site.3,site.12)
newData <- newData[complete.cases(newData),]
summary(newData)
detach(polys)
attach(newData)

lm1 <- lm(Red ~ site.1)
summary(lm1)
lm3 <- lm(Red ~ site.3)
summary(lm3)
lm12 <- lm(Red ~ site.12)
summary(lm12)

lm13 <- lm(Red ~ site.1 + site.3)
summary(lm13)
anova(lm3,lm13)

lm3.12 <- lm(Red ~ site.3 + site.12)
summary(lm3.12)
anova(lm3.12,lm3)

anova (lm1,lm13)
anova(lm12,lm3.12)

library(boot)

#try bootstrap estimation of R2
boot.lm <- function(data,indices,formula,dep) {
	data <- data[indices,]
	if(nlevels(data[,dep])<2) return(NA)
	else summary(lm(formula=formula,data=data))$adj.r.squared
	}
	
boot.lm1 <- boot(newData,boot.lm,1000,formula=Red~site.1,dep="site.1")

#large confidence intervals!



resid1 <- resid(lm1)
resid3 <- resid(lm3)
resid12 <- resid(lm12)

lm3r1 <- lm(resid1~site.3)
summary(lm3r1)

lm1r3 <- lm(resid3~site.1)
summary(lm1r3)

lm3r12 <- lm(resid12~site.3)
summary(lm3r12)

lm12r3 <- lm(resid3~site.12)
summary(lm12r3)