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)