--- title: "Final Analysis Range limits" author: "Andrew D. Nguyen" email: "anbe642@gmail.com" date: "2019-05-09" output: pdf_document: toc: true toc_depth: 3 editor_options: chunk_output_type: console --- # Load libraries and data ```{r load libraries and data, warning=FALSE, message=FALSE} #maps related library(maps) library(mapdata) library(raster) #library(SDMTools) #analysis related library(gvlma) library(rpart) #library(randomForest) library(MASS) library(MuMIn) options(na.action = "na.fail") library(lme4) #library(MCMCglmm) #figures related library(rpart.plot) library(ggplot2) library(gridExtra) library(dplyr) #Theme for ggplot T<-theme_bw()+theme(text=element_text(size=30),axis.text=element_text(size=30), panel.grid.major=element_blank(), panel.grid.minor.x = element_blank(), panel.grid = element_blank(), legend.key = element_blank(),legend.position="none") library(geosphere) ##loading in physiology data k.dat<-read.csv("../Data/20160726_final_cold_phys_parsed.csv") dim(k.dat) str(k.dat) x<-k.dat%>% filter(State=="Maine") range(x$Lat) range(x$Lon) ``` # Analysis of Distribution ```{r load in data for distribution} full.map<-read.csv("../Data/20151130_maine_ants_distr_megan_modified_ANBE.csv",skip=2) sum(ifelse(full.map$Found_Notfound=="0",1,0)) # absents sum(ifelse(full.map$Found_Notfound=="1",1,0)) #presence dim(full.map) w <- getData('worldclim', var='bio', res=2.5) w2 <- getData('worldclim', var='bio', res=.5,lat=45,lon=-68) # Extract climate values for site lon/lat from bioclim data dbio1 <- extract(w, full.map[,c("Lon","Lat")]) dbio2 <- cbind(full.map, dbio1) names(dbio2)[18:36]<-names(k.dat)[23:41] names(dbio2) dbio2$var<-ifelse(dbio2$Found_Notfound=="0","absent","present") #dbio2$SD/100 #daylength(full.map$Lat,doy=1)#jan 1 #daylength(full.map$Lat,doy=140)#may 19th #daylength(full.map$Lat,doy=170)#june 19th #summer<-daylength(full.map$Lat,doy=200)#july 19 #daylength(full.map$Lat,doy=231) # august 19 #daylength(full.map$Lat,doy=300) # #summary(lm(full.map$Lat~summer)) #plot(full.map$Lat,summer) ``` ## Figure(map) for sampling ```{r figure map for sampling} dbio2$color<-ifelse(full.map$Found_Notfound==0,"black","red")# setting color, black = absent, red = present #making distribution map # black colors are absent #red colors are present plot(w, 1, xlim=c(-76,-67), ylim=c(43,47.5),col="white", axes=F, legend=F,main="",box=FALSE) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont','new hampshire'), add = TRUE) #points(full.map$Lon,full.map$Lat,col=dbio2$color,pch=16,cex=1) #points(full.map$Lon,full.map$Lat,cex=1,col="chartreuse3",pch=16) points(full.map$Lon,full.map$Lat,cex=1.5,col="grey50",pch=16) #points(k.dat$Lon,k.dat$Lat,cex=1,pch=16,col="gray") points(k.dat$Lon,k.dat$Lat,cex=1.6,col="cyan",lwd=3) #Scalebar(-70,43,1,t.cex=.5,scale=50) map.scale(-70,43,ratio=FALSE,lwd=3,relwidth=.1) #rect(-150,25,-72,55,col="white",border="white") #rect(-67,50,-50,25,col="white",border="white") ``` # Decision tree analysis ## rpart, classification ```{r rpart} dbio2$Tmax<-dbio2$Tmax/10 dbio2$SD<-dbio2$SD/1000 dbio2$MAT<-dbio2$MAT/10 dbio2$var<-as.factor(dbio2$var) #cor(dbio2[18:36]) cor(dbio2[18:24]) ##subsetting ones that i have an a priori expectation: mean temp, extreme lows, temp variation, and precipitation summary(lm(Tmin~Tmax+SD,data=dbio2)) ##### sub<-data.frame(cbind(dbio2$MAT,dbio2$Tmin,dbio2$SD,dbio2$TAR,dbio2$ISO,dbio2$MDR,dbio2$AP,dbio2[,31])) names(sub)<-c("MAT","Tmin","SD","TAR","ISO","MDR","AP","PSD") knitr::kable(round(cor(sub),3)) ### sub2<-data.frame(cbind(dbio2$MAT,dbio2$Tmin,dbio2$SD,dbio2$TAR,dbio2$MDR)) names(sub2)<-c("MAT","Tmin","SD","TAR","MDR") knitr::kable(round(cor(sub2),3)) ##### vars<-as.data.frame(cbind(dbio2[,18:36],V1=dbio2[,37])) #all bioclim variables #vars<-as.data.frame(cbind(sub,V1=dbio2[,36])) #vars<-as.data.frame(cbind(pca.clim,V1=dbio2[,36])) #names(vars)[1]<-"V1" form<-as.formula(V1~.) tree.1<-rpart(form,data=vars,control=rpart.control(minsplit=20,cp=0),method="class") printcp(tree.1) #plotcp(tree.1) ###for regression trees par(mfrow=c(1,2)) #rsq.rpart(tree.1) par(mfrow=c(1,1)) rpart.plot(tree.1,type=5,extra=100) plot(tree.1) text(tree.1) ###accuracy #m<-predict(tree.1,vars[-20]) m<-predict(tree.1,vars) m.pre<-ifelse(m[,1]< m[,2],"present","absent") #confusion matrix #following this tutorial #http://eric.univ-lyon2.fr/~ricco/tanagra/fichiers/en_Tanagra_Validation_Croisee_Suite.pdf mc<-table(vars$V1,m.pre);mc sum(ifelse(vars$V1== m.pre,1,0))/nrow(vars) ``` ### try pca ```{r, eval=FALSE, warning=FALSE} nm<-princomp(scale(sub)) nm<-princomp(scale(dbio2[,17:35])) summary(nm) biplot(nm) screeplot(nm) pca.clim<-nm$scores[,1:4] dmo1<-glm(dbio2$Found_Notfound~pca.clim[,1]+pca.clim[,2]+pca.clim[,3],family="binomial") summary(dmo1) plot(pca.clim[,1],dbio2$Found_Notfound) dmo2<-glm(dbio2$Found_Notfound~pca.clim[,2],family="binomial") #abline(dmo2) #dbio2$color<-as.factor(dbio2$color) plot(pca.clim[,1],pca.clim[,2],col=dbio2$color,pch=19,xlab="pc1 (62% variation)",ylab="pc2 (23% variation)") ``` ### random forest package ```{r random forest, eval=FALSE, warning=FALSE} rf_fit<-randomForest(V1~+MAT+ MDR+ISO+SD+Tmax+ Tmin+TAR+TWQ+TDQ+TwarmQ +TminQ+ AP + PWM + PDM +PSD+PWQ+PDQ+PwarmQ +PminQ ,data=vars,ntree=10000,oob_score=TRUE,mtry=5,proximity=TRUE) rf_fit plot(rf_fit,log="y") varImpPlot(rf_fit) MDSplot(rf_fit,vars$V1) #acc<-cbind(rf_fit$predicted,vars$V1)# looking at presence absences for random forest and empirical findings table(rf_fit$predicted,vars$V1)# confusion table sum(ifelse(vars$V1== rf_fit$predicted,1,0))/nrow(vars)# accuracy ``` ## Making maps based on decision tree ```{r color coding by climate variables, eval=FALSE} dbio2$symb<-ifelse(dbio2$Found_Notfound=="1",19,17) #tmax #color gradient #dbio3$Tmax_colgrad <- rbPal(20)[as.numeric(cut(dbio3$Tmax,breaks = 20))] #plot(as.numeric(cut(dbio3$Tmax,breaks = 20)),dbio3$Tmin,col=dbio3$Tmax_colgrad) #cutoff based on classification tree dbio2$Tmax.cut<-ifelse(dbio2$Tmax<24.65,"red","black") dbio2$sd.cut<-ifelse(dbio2$SD>=100.5,"black","red") ##plot plot(w, 5, xlim=c(-72,-67), ylim=c(43,48), axes=F,legend=F,main="",box=FALSE,col=colorRampPalette(c("blue","white"))(10)) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine'), add = TRUE) points(dbio2$Lon,dbio2$Lat,col=dbio2$Tmax.cut,pch=dbio2$symb,cex=1.1) #legend(-69,44,c("Present","Absent"),pch=c(16,17),cex=1,box.col = "white") #points(dbio3$Lon,dbio3$Lat,col=dbio3$sea.cut,pch=dbio3$symb,cex=1.1) #text(dbio3$Lon,dbio3$Lat,labels=dbio3$locality,cex=.5,col=dbio3$color) #plot(dbio3$Tmin,dbio3$Found_Notfound) sea.dat<-subset(dbio2,dbio2$Tmax<24.65) #now with seasonality plot(w, 4, xlim=c(-72,-67), ylim=c(43,48), axes=F,legend=F,main="",box=FALSE,col=colorRampPalette(c("white","orange"))(10)) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine'), add = TRUE) #legend(-69,44,c("Present","Absent"),pch=c(16,17),cex=1,box.col = "white") points(dbio2$Lon,dbio2$Lat,col=dbio2$sd.cut,pch=dbio2$symb,cex=1.1) points(sea.dat$Lon,sea.dat$Lat,col="gray",pch=sea.dat$symb,cex=1.1) #next part of tree, MAT mat.dat<-subset(dbio2,dbio2$seasonality>=100.5) dbio2$MAT.cut<-ifelse(dbio2$MAT<6.55,"black","red") plot(w, 1, xlim=c(-72,-67), ylim=c(43,48), axes=F,legend=F,main="",box=FALSE,col=colorRampPalette(c("white","orange"))(10)) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine'), add = TRUE) #legend(-69,44,c("Present","Absent"),pch=c(16,17),cex=1,box.col = "white") points(dbio2$Lon,dbio2$Lat,col=dbio2$MAT.cut,pch=dbio2$symb,cex=1.1) points(sea.dat$Lon,sea.dat$Lat,col="gray",pch=sea.dat$symb,cex=1.1) points(mat.dat$Lon,mat.dat$Lat,col="gray",pch=mat.dat$symb,cex=1.1) #wettest ``` ### Making maps based on decision tree with cut offs ```{r maps based on prediction, eval=FALSE} #setting the extents extent<-c(-72,-65,42,48) extent<-c(-75,-65,42,48) bew<-crop(w2,extent) #bew<-crop(w,extent) #plot(bew,1) #map("state", c('maine','vermont'), add = TRUE) #id points on map hist(getValues(bew[[1]])) #creat blank mat<-na.omit(bew[[1]]) #blank[na.omit(blank)] <- 250 #fill in cells where mat <40 mat[bew[[1]] < 65] <- 1 mat[bew[[1]] > 65] <- 100 plot(mat,col=c("white","goldenrod"),legend=F,xlim=c(-72,67),ylim=c(43,47.5)) map("state", c('maine','vermont','new hampshire'), add = TRUE) #Designate cut offs #tmax dbio2$Tmax.cut<-ifelse(dbio2$Tmax<24.65,"red","black") #sd dbio2$sd.cut<-ifelse(dbio2$SD>=100.5,"red","black") sea.dat<-subset(dbio2,dbio2$Tmax<24.65) ##mat mat.dat<-subset(dbio2,dbio2$SD>=100.5) dbio2$MAT.cut<-ifelse(dbio2$MAT<6.55,"black","red") ##PwarmQ pwq<-subset(dbio2,dbio2$MAT >= 6.5) dbio2$pwq.cut<-ifelse(dbio2$PwarmQ < 261,"red",ifelse(dbio2$PwarmQ >= 270,"red","black")) #symbosl dbio2$symb<-ifelse(dbio2$Found_Notfound=="1",19,17) lar<-na.omit(bew) ############################ #Change Tmax ############################ lar[bew[[5]] < 246.5] <- 100 lar[bew[[5]] > 246.5] <- 1 #plot plot(lar[[5]],col=c("white","grey75"),legend=T) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont',"new hampshire"), add = TRUE) points(dbio2$Lon,dbio2$Lat,col=dbio2$Tmax.cut,pch=dbio2$symb,cex=1.1) ############################ #-- ############################ #Change SD lar[bew[[4]] < 10050] <- 1 lar[bew[[4]] >= 10050] <- 100 plot(lar[[4]],col=c("white","grey75"),legend=F) #plot(lar[[4]],col=c("red","grey75"),legend=F,add=TRUE) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont',"new hampshire"), add = TRUE) points(dbio2$Lon,dbio2$Lat,col=dbio2$sd.cut,pch=dbio2$symb,cex=1.1,bg=dbio2$sd.cut) points(sea.dat$Lon,sea.dat$Lat,col="grey31",pch=sea.dat$symb,cex=1.1,bg="black") ############################ #-- ############################ #Change MAT lar[bew[[1]] < 65] <- 1 lar[bew[[1]] > 65] <- 100 plot(lar[[1]],col=c("white","grey75"),legend=F) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont',"new hampshire"), add = TRUE) #points(dbio2$Lon,dbio2$Lat,col=dbio2$MAT.cut,pch=dbio2$symb,cex=1.1) points(dbio2$Lon,dbio2$Lat,col=dbio2$color,pch=19,cex=1.1) points(sea.dat$Lon,sea.dat$Lat,col="grey31",pch=sea.dat$symb,cex=1.1) points(mat.dat$Lon,mat.dat$Lat,col="grey31",pch=mat.dat$symb,cex=1.1) ############################ ############################ #Change PWARmQ ############################ lar[bew[[18]] >= 261 & bew[[5]] < 279 ] <- 270 lar[bew[[18]] >=270 ] <- 540 lar[bew[[18]] < 261] <- 135 plot(lar[[18]],col=c("grey80","white","grey66")) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont',"new hampshire"), add = TRUE) points(dbio2$Lon,dbio2$Lat,col=dbio2$pwq.cut,pch=dbio2$symb,cex=1.1) points(sea.dat$Lon,sea.dat$Lat,col="grey31",pch=sea.dat$symb,cex=1.1) points(mat.dat$Lon,mat.dat$Lat,col="grey31",pch=mat.dat$symb,cex=1.1) points(pwq$Lon,pwq$Lat,col="grey31",pch=pwq$symb,cex=1.1) ``` ### Re- Making maps based on decision tree with cut offs ```{r remaking maps, eval=FALSE} ## points related dbio2$coco<-ifelse(dbio2$Found_Notfound=="1","red","black") node2<-subset(dbio2,dbio2$Tmax >= 24.5) node3<-subset(node2,node2$SD <101) node4<-subset(node3,node3$MAT<6.5) #node5<-subset(node4,node4$PwarmQ >= 270) ### Climate Tm<-na.omit(bew[[5]]) Tm[bew[[5]] < 246.5] <- 100 # absent Tm[bew[[5]] > 246.5] <- 1 sd<-na.omit(bew[[4]]) sd[bew[[4]] < 10050] <- 1 sd[bew[[4]] >= 10050] <- 100 # absent #mat mat<-na.omit(mat[[1]]) mat[bew[[1]] < 65] <- 50 mat[bew[[1]] >= 65] <- 500 ##pwarmq #p<-na.omit(mat[[18]]) p<-na.omit(bew[[18]]) p[bew[[18]] < 261] <- -100 p[bew[[18]] >= 261 & bew[[18]] < 270 ] <- 250# absent p[bew[[18]] >=270 ] <- -100 ##Tmax redone plot(Tm,col=c("white","grey75"),legend=F) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont','new hampshire'), add = TRUE) #points(dbio2$Lon,dbio2$Lat,pch=16,col=dbio2$coco) #points(dbio2$Lon,dbio2$Lat,pch=16,col=dbio2$coco) ##SD + Tmax #tmaxsd<-overlay(Tm[[4]],sd[[5]],fun=sum) tmaxsd<-overlay(Tm[[1]],sd[[1]],fun=sum) plot(tmaxsd,col=c("white","grey75"),legend=F) #plot(lar[[4]],col=c("white","black"),legend=F) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont','new hampshire'), add = TRUE) points(dbio2$Lon,dbio2$Lat,pch=16,col=dbio2$coco) points(-73.19686,44.43941,col="blue",pch=16) #points(node2$Lon,node2$Lat,pch=16,col=node2$color) ##Tmax + SD + MAT Tmsdmat<-overlay(tmaxsd,mat,fun=sum) #unique(getValues(Tmsdmat)) plot(Tmsdmat,col=c("white","grey75","pink","grey75"),legend=TRUE) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont',"new hampshire"), add = TRUE) points(dbio2$Lon,dbio2$Lat,pch=16,col=dbio2$coco) #points(node3$Lon,node3$Lat,pch=16,col=node3$color) #points(-73.19686,44.43941,col="blue",pch=16) ##Tmax + SD + MAT PWQ Tmsdmat[Tmsdmat[[1]]==502]<- -1000 # mat presence fin<-overlay(Tmsdmat,p,fun=sum) sort(unique(getValues(fin))) fin[fin[[1]]==51]<-2000 #plot(fin,col=colorRampPalette(c("black", "blue"))(3),legend=TRUE) plot(fin,col=c("pink","pink","grey75","grey75","grey75"),legend=TRUE) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine','vermont',"new hampshire"), add = TRUE) points(dbio2$Lon,dbio2$Lat,pch=16,col=dbio2$coco) #points(node4$Lon,node4$Lat,pch=16,col=node4$color) points(k.dat$Lon,k.dat$Lat,cex=1.6,col="cyan",lwd=3) points(-73.19686,44.43941,col="blue",pch=16) ``` # Analysis of Cold physiology ## Methods figure ```{r methods fig} x<-seq(1:4) y<-c(-5,0,5,25) plot(x,y,ylim=c(-5,25),bty='n',las=1,axes=F,type="n",xlab="Time (hours)",ylab="",xlim=c(0,4)) axis(1,at=c(0,1,2,3,4),las=1) axis(2,at=c(-5,0,5,25),las=1) mtext(side=2,text="Temperature (°C)",line=2) w<-c(25,25,-5,-5,25,25) z.25<-c(0,2,2,3,3,4) lines(z.25,w,lwd=5,lty=6) lines(c(0,1,1),c(5,5,25),lty=3,lwd=5) lines(c(0,1,1),c(0,0,25),lty=2,lwd=5) lines(c(0,1,1),c(-5,-5,25),lty="dotdash",lwd=5) lines(c(1,2,2,3,3,4),w,lwd=5,lty=1) mtext(side=3,"Pre-treatments",padj=1,adj=.05) mtext(side=3,"Chill coma recovery",padj=1,adj=.95) text(1.5,24,label="Recovery") text(2.5,-4,label="Cold shock") ``` ## estimating broad sense (colony level) variance-covariance G matrix ### with manova approach ```{r with manova approach} library(tidyr) o.dat<-read.csv("../Data/20160107_chillcomarecovery_ind_p_colony_dat.csv") #o.dat<-read.csv("../Data/test.csv") str(o.dat) summary(o.dat$Colony);length(summary(o.dat$Colony)) o.dat$treatment_recovery_s<-as.numeric(o.dat$treatment_recovery_s) #19 colonies #converting long to wide wide.o.dat<-spread(o.dat,pretreat_Temp,treatment_recovery_s,) names(wide.o.dat)[3:6]<-c("minusfive","zero","five","twentyfive") head(wide.o.dat) traits<-cbind(wide.o.dat$minusfive,wide.o.dat$zero,wide.o.dat$five,wide.o.dat$twentyfive) dim(wide.o.dat) dim(na.omit(wide.o.dat)) mew<-na.omit(wide.o.dat) dim(mew) #mew<-read.csv("no_nas.csv") gmod2<-manova(cbind(minusfive,zero,five,twentyfive)~Colony,data=mew) summary(gmod2) #calculations for G colony<-as.data.frame(summary(gmod2)$SS[1])/18 error<-as.data.frame((summary(gmod2)$SS[2]))/42 #colony<-as.data.frame(summary(gmod2)$SS[1])/14 #error<-as.data.frame((summary(gmod2)$SS[2]))/32 Gmatrix2<-(colony-error)/3.20401 Gmatrix2 #PCA summary(princomp(Gmatrix2)) summary(princomp(Gmatrix2))$loadings[,1:3] ``` ```{r figures for unpermuted g matrix, include=FALSE, echo=FALSE} par(mfrow=c(1,1)) #gmax 61% plot(c(-5,0,5,25),summary(princomp(Gmatrix2))$loadings[,1],ylim=c(-.7,.7),pch=16,col="blue",cex=2,main="Warmer-cooler trade off (61.70%)",las=1,ylab="Gmax loadings",xlab="Pre-treatment Temperature (°C)",cex.axis=1.1) lines(c(-5,0,5,25),summary(princomp(Gmatrix2))$loadings[,1]) abline(h=0,lty="dotdash",lwd=2) #g2 37% plot(c(-5,0,5,25),summary(princomp(Gmatrix2))$loadings[,2],ylim=c(-1,1),pch=16,col="purple",cex=2,main="Generalist-Specialist trade off (37.91%)",las=1,ylab="g2 loadings",xlab="Pre-treatment Temperature (°C)") lines(c(-5,0,5,25),summary(princomp(Gmatrix2))$loadings[,2]) abline(h=0,lty="dotdash",lwd=2) ``` ### G matrix using lme4 package to fit mixed effects models ```{r mixed effects models to estimate G} o.dat$pretreat_Temp<-as.factor(as.character(o.dat$pretreat_Temp)) m0<-lmer(formula=treatment_recovery_s~1+(pretreat_Temp|Colony),REML=TRUE,data=o.dat) #m1.1<-lmer(formula=treatment_recovery_s~pretreat_Temp+(0+pretreat_Temp|Colony),REML=TRUE,data=o.dat) m1<-lmer(formula=treatment_recovery_s~1+(0+pretreat_Temp|Colony),REML=TRUE,data=o.dat) #o.dat$scale1<-scale(o.dat$treatment_recovery_s) #m1<-lmer(formula=scale1~1+(0+pretreat_Temp|Colony),REML=TRUE,data=o.dat) Gmatrix.pac<-VarCorr(m1)$Colony;Gmatrix.pac Gmatrix.pac[1:4,1:4] #print(Gmatrix.pac,comp=c("Variance","Std.Dev")) print(Gmatrix.pac,comp=c("Std.Dev")) vco<-as.data.frame(VarCorr(m1))[-11,] #CalcICV(Gmatrix.pac[1:4,1:4]) str(Gmatrix.pac) #varcomp(m1) #gmin<-read.csv("Gmatrix_min.csv",header=F) #names(gmin)<-c("-5","0","5","25") #hmin<-princomp(gmin) #summary(hmin) #hmin$loadings[,1:3] h<-princomp(Gmatrix.pac[1:4,1:4]) h$loadings[,1:3] summary(h) plot(c(-5,0,25,5),h$loadings[,1],ylim=c(-.8,.8),pch=16,cex=2) #plot(c(-5,0,25,5),h$loadings[,2],ylim=c(-.8,.8),pch=16,cex=2) ``` ### PC loadings vs temp fig ```{r pc loadings vs temp fig} T<-theme_bw()+theme(text=element_text(size=30),axis.text=element_text(size=30), panel.grid.major=element_blank(), panel.grid.minor.x = element_blank(), panel.grid = element_blank(), legend.key = element_blank(),legend.position="none") n<-as.data.frame(cbind(Gmax=h$loadings[,1],Temp=c(-5,0,25,5))) n$Gmax<-n$Gmax* -1 gplot<-ggplot(n,aes(x=Temp,y=Gmax))+geom_point(size=10)+scale_y_continuous(limit=c(-.6,1),labels=seq(-1,1,.2),breaks=seq(-1,1,.2))+scale_x_continuous(limit=c(-5,25),labels=c(-5,0,5,25),breaks=c(-5,0,5,25))+geom_line(size=2)+xlab("Pre-treatment Temperature (°C)")+ylab(expression(paste("G"[max]," Loadings")))+T+geom_hline(yintercept=0,linetype=6,size=1.5) gplot ``` ## Relationship with cline (stats; ANCOVA) ### basal cold tolerance ```{r cold tolerance analysis, warning=FALSE} k.dat<-read.csv("../Data/20160726_final_cold_phys_parsed.csv") #k.dat<-subset(k.dat,k.dat$Colony != "Canaan") #k.dat<-subset(k.dat,k.dat$Town != "Burlington") #k.dat<-subset(k.dat,k.dat$Colony != "Greenfield") #dim(k.dat) #str(k.dat) k.dat$pretreat_Temp<-as.factor(as.character(k.dat$pretreat_Temp)) #k.dat$pretreat_Temp<-as.numeric(as.character(k.dat$pretreat_Temp)) #stats cold.mod1<-aov(treatment_recovery_s~Tmin*pretreat_Temp+Colony,data=k.dat) cold.mod1<-aov(treatment_recovery_s~MAT*pretreat_Temp+Colony,data=k.dat) #cold.mod1<-lm(treatment_recovery_s~Tmin*pretreat_Temp+Colony#,data=k.dat) summary(cold.mod1) #stepwise aic qc<-stepAIC(cold.mod1,direction="both") summary(qc) ###plotting cold tolerance vs mat ggplot(data=k.dat,aes(y=treatment_recovery_s,x=MAT,colour=factor(pretreat_Temp)))+geom_point()+stat_smooth(method=lm,se=FALSE) #basal dfg<-subset(k.dat,pretreat_Temp=="25") cold.mod101<-lm(treatment_recovery_s~MAT,data=subset(k.dat,pretreat_Temp=="25")) summary(cold.mod101) #induced at 0 assd<-subset(k.dat,pretreat_Temp=="0") assd$hardening<-dfg$treatment_recovery_s-assd$treatment_recovery_s cold.mod102<-lm(hardening~MAT,data=assd) summary(cold.mod102) #mumin subs<-dredge(cold.mod1) head(subs,10) # delta 4 AIC criteria #summary(model.avg(subs, subset = delta < 4)) #top 2 models #summary(model.avg(subs[1:2])) #subs[1] ``` ### hardening ability ```{r hardening ability, warning=FALSE} library(plyr) mew<-ddply(mew,.(Colony),summarize,minusfive=mean(minusfive),zero=mean(zero),five=mean(five),twentyfive=mean(twentyfive)) mew2<-inner_join(mew,k.dat,by="Colony") mew3<-mew2$twentyfive-mew2[,2:5] mew3<-data.frame(mew3[,-4]) names(mew3)<-c("hm5","h0","h5") mew4<-data.frame(cbind(Colony=mew2[,1],MAT=mew2$MAT,Tmin=mew2$Tmin,SD=mew2$SD,mew3)) mew5<-gather(mew4,pretreat_Temp,hardening,hm5:h5) dim(mew5) mew6<-ddply(mew5,.(Colony,pretreat_Temp),summarize,hardening=mean(hardening),Tmin=mean(Tmin),SD=mean(SD),MAT=mean(MAT)) mew6$PT<-rep(c("-5","0","5"),length(mew6$Colony)/3) head(mew6) #stats mew6$PT<-as.factor(as.character(mew6$PT)) #cold.mod8<-aov(hardening~Tmin*PT+Colony,data=mew6) cold.mod8<-aov(hardening~MAT*PT+Colony,data=mew6) #cold.mod8<-lm(hardening~Tmin*PT+Colony,data=mew6) #cold.mod8<-aov(hardening~MAT*PT+Error(Colony/PT),data=mew6) summary(cold.mod8) #stepwise aic qc8<-stepAIC(cold.mod8,direction="both") summary(qc8) mp<-dredge(qc8) summary(model.avg(mp[1:2])) ##Exploring influential points par(mfrow=c(2,2)) plot(qc8) par(mfrow=c(1,1)) ``` ### Figures associated with stats ```{r cold perf v tmin, warning=FALSE} T<-theme_bw()+theme(text=element_text(size=30),axis.text=element_text(size=30), legend.text=element_text(size=28), panel.grid.major=element_blank(), panel.grid.minor.x = element_blank(), panel.grid = element_blank(), legend.key = element_blank(),legend.position=c(.3,.1),legend.title=element_blank()) #legend.direction="horizontal" #ggplot(k.dat,aes(x=Tmin,y=1000/treatment_recovery_s,shape=factor(pretreat_Temp)))+geom_smooth(data=k.dat,method="lm",se=F,aes(colour=factor(pretreat_Temp)),size=1.5)+geom_point(size=4)+T+scale_color_manual(values=c("blue", "purple","red","orange"))+scale_y_continuous("Cold tolerance (1000/CCRT)",limits=c(.7,3.5),breaks=seq(.5,3.5,.5))+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-16.5,-13.5,.5),limits=c(-16.5,-13.5))+T+scale_shape_manual(values=c(15,16,17,18)) k.dat$cold_tol<-max(k.dat$treatment_recovery_s)-k.dat$treatment_recovery_s yy<-ggplot(k.dat,aes(x=Tmin,y=cold_tol,group=factor(pretreat_Temp),colour=factor(pretreat_Temp)))+geom_smooth(data=k.dat,method="lm",se=TRUE,aes(colour=factor(pretreat_Temp)),size=3)+geom_point(size=8)+T+scale_color_manual(values=c("blue", "purple","red","orange"))+scale_y_continuous(expression(paste("Cold Tolerance ( ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,850),breaks=seq(0,850,100))+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-17,-13,1),limits=c(-17,-13))+T+scale_shape_manual(values=c(15,16,17,18)) yy #grid.arrange(yy,beta,ncol=2) #grid.arrange(yy,beta,ncol=2) ``` ### 2 panel fig with basal cold tolerance on top and hardening on bottom ```{r basal cold tolerance and hardening figs, warning=FALSE} #basal cold tolerance ksub<-subset(k.dat,k.dat$pretreat_Temp=="25") ksub$coldplot<-max(ksub$treatment_recovery_s)-ksub$treatment_recovery_s summary(lm(coldplot~Tmin,data=ksub)) #plot(lm(coldplot~Tmin,data=ksub)) mo<-ggplot(ksub,aes(x=Tmin,y=coldplot))+geom_point(size=10)+ geom_smooth(method="lm",se=TRUE,size=4,colour="black")+T+scale_y_continuous(expression(paste("Basal Cold Tolerance","\n"," (25°C; ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,800),breaks=seq(0,800,100))+T+scale_shape_manual(values=c(15,16,17,18))+scale_x_continuous(breaks=seq(-16.5,-13,1),limits=c(-16.5,-13.25))+xlab("") mo #+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-16.5,-13.5,.5),limits=c(-16.5,-13.5))+ zero<-subset(k.dat,k.dat$pretreat_Temp=="0") ksub$hard.zero<-ksub$treatment_recovery_s-zero$treatment_recovery_s #ksub$hard.zero<-zero$treatment_recovery_s-ksub$treatment_recovery_s ######################### nbn3<-plyr::ddply(ksub,.(Colony,pretreat_Temp),summarize,CCRT=mean(hard.zero),Tmin=mean(Tmin),MAT=mean(MAT),SD=mean(SD),MDR=mean(MDR)) #io<-ggplot(nbn3,aes(x=Tmin,y=CCRT))+geom_smooth(method="lm",se=TRUE,size=4,colour="black")+geom_point(size=10,)+scale_y_continuous(expression(paste("Cold Hardening (",CCRT[control], " - ", CCRT[pre-treatment],")")),limits=c(-400,800),breaks=seq(-400,800,200))+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-16.5,-13,1),limits=c(-16.5,-13.25))+scale_shape_manual(values=c(15,16,17,18))+T #io summary(lm(CCRT~Tmin,data=nbn3)) #grid.arrange(io,hbeta,ncol=2) figfor0<-subset(nbn3,nbn3$pretreat_Temp=="0") for0<-ggplot(nbn3,aes(x=Tmin,y=CCRT))+geom_smooth(method="lm",se=TRUE,size=4,colour="black")+geom_point(size=10)+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],")")),limits=c(-400,800),breaks=seq(-400,800,200))+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-16.5,-13.25,1),limits=c(-16.5,-13.25))+scale_shape_manual(values=c(15,16,17,18))+T grid.arrange(mo,for0,nrow=2) #changing y axis mo1<-mo+scale_y_continuous("Cold Tolerance",limits=c(0,800),breaks=seq(0,800,100)) for01<-for0+scale_y_continuous("Cold Hardening",limits=c(-400,800),breaks=seq(-400,800,200)) mo1 for01 grid.arrange(mo1,for01,nrow=2) # cold hardening #io<-ggplot(mew6,aes(x=Tmin,y=hardening,colour=factor(PT)))+geom_smooth(method="lm",se=F,size=4)+geom_point(size=10)+scale_color_manual(values=c("blue", "purple","orange"))+scale_y_continuous(expression(paste("Cold Hardening (",CCRT[control], " - ", CCRT[pre-treatment],")")),limits=c(-400,800),breaks=seq(-400,800,200))+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-16.5,-13.5,.5),limits=c(-16.5,-13.4))+scale_shape_manual(values=c(15,16,17,18))+T mp<-lm(hardening~Tmin + PT, data=mew6) summary(mp) m5<-subset(mew6,mew6$PT=="-5") p0<-subset(mew6,mew6$PT=="0") p5<-subset(mew6,mew6$PT=="5") plot(m5$Tmin,m5$hardening) abline(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="-5"))) abline(a=1823,b=116,col="green") abline(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="0"))) abline(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="5"))) coef(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="-5"))) summary(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="-5")))$coefficients summary(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="5")))$coefficients summary(lm(hardening~Tmin,data=subset(mew6,mew6$PT=="0")))$coefficients io<-ggplot(mew6,aes(x=Tmin,y=hardening,colour=factor(PT)))+geom_abline(slope=72,intercept=1167.35,colour="blue",size=3)+geom_abline(slope=72,intercept=1167-25.16,colour="orange",size=3)+scale_color_manual(values=c("blue", "purple","orange"))+geom_point(size=10)+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],")")),limits=c(-400,800),breaks=seq(-400,800,200))+scale_x_continuous("Local Minimum Temperature (Tmin, °C)",breaks=seq(-16.5,-13.25,1),limits=c(-16.5,-13.25))+T+geom_abline(slope=72,intercept=1167+155.42,colour="purple",size=3) io #grid arrange grid.arrange(mo,io,nrow=2) #changing ``` # lets try figures with MAT ```{r , warning=FALSE} MAT1<-ggplot(ksub,aes(x=MAT,y=coldplot))+geom_point(size=10)+ geom_smooth(method="lm",se=TRUE,size=4,colour="black")+T+scale_y_continuous(expression(paste("Basal Cold Resistance","\n"," (25°C; ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,800),breaks=seq(0,800,100))+T+scale_shape_manual(values=c(15,16,17,18))+scale_x_continuous(breaks=seq(4.5,7,.5),limits=c(4.5,7))+xlab("") MAT1 MAT2<-ggplot(nbn3,aes(x=MAT,y=CCRT))+geom_smooth(method="lm",se=TRUE,size=4,colour="black")+geom_point(size=10)+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],")")),limits=c(-100,800),breaks=seq(0,800,100))+scale_x_continuous("Mean Annual Temperature (MAT, °C)",breaks=seq(4.5,7,.5),limits=c(4.5,7))+scale_shape_manual(values=c(15,16,17,18))+T MAT2 grid.arrange(MAT1,MAT2,nrow=2) #plasticity with seasonality SD1<-ggplot(ksub,aes(x=SD/10,y=coldplot))+geom_point(size=10)+ geom_smooth(method="lm",se=TRUE,size=4,colour="black")+T+scale_y_continuous(expression(paste("Basal Cold Tolerance","\n"," (25°C; ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,800),breaks=seq(0,800,100))+T+scale_shape_manual(values=c(15,16,17,18))+scale_x_continuous("",breaks=seq(9.6,10.15,.15),limits=c(9.6,10.15))#+scale_shape_manual(values=c(15,16,17,18)) SD1 summary(lm(coldplot~SD+MAT,data=ksub)) summary(lm(coldplot~SD+MAT,data=ksub)) SD2<-ggplot(nbn3,aes(x=SD/10,y=CCRT))+geom_point(size=10)+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],")")),limits=c(-100,800),breaks=seq(0,800,100))+stat_smooth(method="lm",formula=y~1,colour="black",size=4)+T+scale_x_continuous("Temperature Seasonality (SD)",breaks=seq(9.6,10.15,.15),limits=c(9.6,10.15))+scale_shape_manual(values=c(15,16,17,18))#+geom_smooth(method="lm",se=TRUE,size=4,colour="black") SD2 grid.arrange(SD1,SD2,nrow=2) ``` ### 3 panel with basal and hardening vs MAT, then basal vs hardening fig ```{r, warning=FALSE} cc2<-ggplot(ksub,aes(x=coldplot,y=hard.zero))+geom_point()+stat_smooth(method="lm")+T+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],")")),limits=c(-100,800),breaks=seq(0,800,100))+scale_x_continuous(expression(paste("Basal Cold Tolerance","\n"," (25°C; ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,800),breaks=seq(0,800,100))+geom_point(size=10)+ geom_smooth(method="lm",se=TRUE,size=4,colour="black") cc2 grobmat<-arrangeGrob(MAT1,MAT2,nrow=2) grid.arrange(grobmat,cc2,ncol=2) ``` ### calculating % improvement ```{r % improvement,eval=FALSE} zero<-zero%>% filter(State!="Vermont") ksubb<-ksub%>% filter(State!="Vermont") nn<-rbind(zero,ksubb[,-46:-48]) nn$delta<-max(nn$treatment_recovery_s)-nn$treatment_recovery_s ggplot(nn,aes(x=Tmin,y=delta,colour=factor(pretreat_Temp)))+geom_point()+stat_smooth(method="lm") #nn$pretreat_Temp<-as.factor(nn$pretreat_Temp) #spread(nn,pretreat_Temp,treatment_recovery_s,) nn2<-melt(nn,id.vars=c("pretreat_Temp","Colony"),measure.vars=c("delta")) wide.nn<-spread(nn2,pretreat_Temp,value)[-19,-2] knitr::kable(wide.nn) round(abs(wide.nn$`0`-wide.nn$`25`)/wide.nn$`25`,3)*100 range(round(abs(wide.nn$`0`-wide.nn$`25`)/wide.nn$`25`,3)*100) ##lets look at the amount of time that has changed itself wide.nn$minutes_change<-abs(wide.nn$`0`-wide.nn$`25`)/60 #range of minutes range(wide.nn$minutes_change) #max(maine$treatment_recovery_s)-maine$treatment_recovery_s #(maine$hard.zero-maine$coldplot2)/maine$coldplot2 ``` ## Figure of cold performance curves (Cold tolerance = 1000/CCRT) **Going to combine with Gmatrix fig to make 2 panel fig** ```{r cold performance curves, warning=FALSE} k.dat$pretreat_Temp<-as.numeric(as.character(k.dat$pretreat_Temp)) #new.dat<-k.dat[k.dat$Colony != "Canaan" & k.dat$Colony != "Saponac",] k.dat$coldplot<-max(k.dat$treatment_recovery_s)-k.dat$treatment_recovery_s new.dat<-k.dat[k.dat$Colony != "Canaan",] new.dat2<-new.dat%>% dplyr::filter(State=="Maine") library(RColorBrewer) bcal<-c(brewer.pal(8,"Dark2"),brewer.pal(12,"Paired")) nocol2<-ggplot(data=new.dat2,aes(x=pretreat_Temp,y=coldplot,group=Colony,col=Colony))+geom_point(size=5)+geom_line(size=1.1)+xlab("Pre-treatment Temperature (°C)")+scale_x_continuous(labels=c(-5,0,5,25),breaks=c(-5,0,5,25))+theme(text=element_text(size=20))+T+scale_y_continuous(expression(paste("Cold Resistance ( ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,850),breaks=seq(0,850,100))+theme(legend.position="none")+scale_colour_manual(values=bcal) nocol2 #boxplots bp<-ggplot(new.dat,aes(y=coldplot,x=factor(pretreat_Temp)))+geom_boxplot(fill="gray")+ylab("")+xlab("Pre-treatment Temperature(°C)")+T+scale_y_continuous(expression(paste("Cold Resistance ( ",CCRT[max], " - ", CCRT[obs],")")),limits=c(0,850),breaks=seq(0,850,100)) #bp ##combining with G matrix fig #grid.arrange(bp,nocol2,gplot,ncol=3) grid.arrange(nocol2,bp,gplot,ncol=3) #grid.arrange(nocol2,gplot,bp,ncol=3) #fixing y axis bpy<-bp+scale_y_continuous("Cold Tolerance",limits=c(0,850),breaks=seq(0,850,100)) bpy ncol3<-nocol2+scale_y_continuous("Cold Tolerance",limits=c(0,850),breaks=seq(0,850,100)) ncol3 grid.arrange(bpy,ncol3,gplot,ncol=3) grid.arrange(ncol3,bpy,gplot,ncol=3) grid.arrange(ncol3,gplot,bpy,ncol=3) ``` # playing iwth k.dat dataset with maps ```{r, eval=FALSE, warning=FALSE} plot(w, 6, xlim=c(-74,-67), ylim=c(43,48), axes=F,legend=F,main="",box=FALSE,col=colorRampPalette(c("blue","white"))(20)) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine'), add = TRUE) #scale basal cold tolerance basal<-subset(k.dat,k.dat$pretreat_Temp=="25") basal$scale<-(scale(1000/basal$treatment_recovery_s,center=TRUE,scale=TRUE)+2.6)/2 ##scaling cold hardening at 0 har<-subset(k.dat,k.dat$pretreat_Temp=="0") har$har<-(scale(basal$treatment_recovery_s-har$treatment_recovery_s)+2.6)/2 ##plotting points basal points(basal$Lon,basal$Lat,pch=19,col="black",cex=basal$scale) points(har$Lon,har$Lat,pch=19,col="black",cex=har$har) #temperature seasonality plot(w, 4, xlim=c(-74,-67), ylim=c(43,48), axes=F,legend=F,main="",box=FALSE,col=colorRampPalette(c("white","black"))(20)) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine'), add = TRUE) points(basal$Lon,basal$Lat,pch=19,col="black",cex=basal$scale) #mdr plot(w, 2, xlim=c(-74,-67), ylim=c(43,48), axes=F,legend=F,main="",box=FALSE,col=colorRampPalette(c("blue","red"))(20)) map("worldHires",c("USA","Canada"),add=TRUE) map("state", c('maine'), add = TRUE) points(basal$Lon,basal$Lat,pch=19,col="black",cex=basal$scale) ``` # Re-analysis with ksub dataset, sets of regressions ```{r , warning=FALSE} new.dat3<-new.dat%>% filter(State!="Vermont") new.dat3$coldplot2<-max(new.dat3$treatment_recovery_s)-new.dat3$treatment_recovery_s #ksub ##basal mm<-lm(treatment_recovery_s~MAT+MDR+ISO+SD+Tmax+Tmin+TAR+TWQ+AP+PWM+PDM+PSD+PWQ+PDQ+PwarmQ+PminQ,data=ksub) summary(mm) summary(stepAIC(mm,direction="both")) ksub2<-data.frame(cbind(ksub$treatment_recovery_s,apply(ksub[,23:41],2,scale))) mm3<-lm(V1~MAT+MDR+ISO+SD+Tmax+Tmin+TAR+TWQ+AP+PWM+PDM+PSD+PWQ+PDQ+PwarmQ+PminQ,data=ksub2) summary(stepAIC(mm3,direction="both")) lapply(ksub2[,-1],function(x){summary(lm(ksub2$V1~x))}) ###hardening ksub3<-data.frame(cbind(ksub$hard.zero,apply(ksub[,23:41],2,scale))) lapply(ksub3[,-1],function(x){summary(lm(ksub3$V1~x))}) ###subsetting only the temp variables #basal mm1<-lm(treatment_recovery_s~MAT+MDR+ISO+SD+Tmax+Tmin+TAR+TWQ,data=ksub) #summary(mm1) summary(stepAIC(mm1,direction="both")) #hardening mm3.1<-lm(V1~MAT+MDR+ISO+SD+Tmax+Tmin+TAR+TWQ,data=ksub2) summary(mm3.1) summary(stepAIC(mm3.1,direction="both")) #correlation between hardening and basal cc2<-ggplot(ksub,aes(x=coldplot,y=hard.zero))+geom_point()+stat_smooth(method="lm")+T cc2 summary(lm(coldplot~hard.zero,data=ksub)) ksub.maine<-ksub%>%filter(State=="Maine") ggplot(ksub.maine,aes(x=coldplot,y=hard.zero))+geom_point()+stat_smooth(method="lm") summary(lm(hard.zero~coldplot,data=ksub.maine)) ksub.low500<-ksub.maine%>%filter(coldplot<500) summary(lm(coldplot~hard.zero,data=ksub.low500)) ksub.above500<-ksub.maine%>%filter(coldplot>500) summary(lm(coldplot~hard.zero,data=ksub.above500)) ksub.maine$coldplot_split<-ifelse(ksub.maine$coldplot<500,"A","B") ggplot(ksub.maine,aes(x=coldplot,y=hard.zero,colour=coldplot_split))+geom_point()+stat_smooth(method="lm") library(segmented) hj<-lm(hard.zero~coldplot,data=ksub.maine) my.seg <- segmented(hj, seg.Z = ~ coldplot, psi = 500) summary(my.seg) my.seg$psi ``` ## redo'ing stats without vt ```{r, warning=FALSE} maine<-ksub%>% dplyr::filter(State!="Vermont") maine$coldplot2<-max(maine$treatment_recovery_s)-maine$treatment_recovery_s summary(lm(hard.zero~MAT,data=maine)) #mean((lm(hard.zero~MAT,data=maine))$residuals) par(mfrow=c(2,2)) plot(lm(hard.zero~MAT,data=maine)) par(mfrow=c(1,1)) #acf(lm(hard.zero~MAT,data=maine)$residuals) #cor.test(maine$MAT,lm(hard.zero~MAT,data=maine)$residuals) #influence.measures(lm(hard.zero~MAT,data=maine)) #gvlma(lm(hard.zero~MAT,data=maine)) #gvlma(lm(hard.zero~MAT+I(MAT^2),data=maine)) summary(lm(hard.zero~MAT+I(MAT^2),data=maine)) #summary(lm(hard.zero~Tmin+I(Tmin^2),data=maine)) #summary(lm(hard.zero~Tmin,data=maine)) hj2<-lm(hard.zero~MAT,data=maine) library(segmented) my.seg <- segmented(hj2, seg.Z = ~ MAT, psi = 5.25) summary(my.seg) my.seg$psi my.seg$coefficients davies.test(my.seg,seg.Z =~MAT) #predicting maine$pred.hard<-predict(my.seg) maine$bp<-ifelse(maine$MAT<5.747,"low","high") ggplot(maine,aes(x=MAT,y=hard.zero,colour=bp))+geom_point()+stat_smooth(method="lm") ##regress predicted values against residuals maine$residualsMAT<-resid(lm(hard.zero~MAT,data=maine)) #p<-predict(lm(hard.zero~MAT,data=maine)) #summary(lm(p~r)) ggplot(maine,aes(x=MAT,y=residualsMAT))+geom_point()+stat_smooth(method="lm") #summary(stepAIC(lm(hard.zero~MAT+I(MAT^2),data=maine),direction="both")) #plot(lm(hard.zero~MAT,data=maine)) #summary(lm(hard.zero~SD,data=maine)) #summary(lm(hard.zero~TAR,data=maine)) #summary(lm(hard.zero~Tmin,data=maine)) #summary(lm(hard.zero~MDR,data=maine)) summary(lm(coldplot~MAT,data=maine)) #summary(lm(coldplot~MAT+I(MAT^2),data=maine)) #summary(lm(coldplot~Tmin+I(Tmin^2),data=maine)) #summary(lm(coldplot~Tmin,data=maine)) ggplot(maine,aes(x=MAT,y=hard.zero))+geom_point()+stat_smooth(method="lm",formula=y~x+I(x^2)) #ggplot(maine,aes(x=MAT,y=hard.zero))+geom_point()+stat_smooth(method="lm",formula=y~x+I(x^2)) #ggplot(maine,aes(x=MAT,y=coldplot))+geom_point()+stat_smooth(method="lm") ##corelation between coldplot and hard.zero summary(lm(hard.zero~coldplot2,data=maine)) maine.low<-maine%>% filter(coldplot2% filter(coldplot2>median(coldplot2)) summary(lm(hard.zero~coldplot2,data=maine.hi)) ###pca temperature variables clim.pca<-princomp(scale(maine[,23:29])) clim.pca<-princomp(scale(maine[,23:33])) summary(clim.pca) clim.pca$loadings[,1:2] summary(lm(maine$hard.zero~clim.pca$scores[,1])) #summary(lm(maine$hard.zero~clim.pca$scores[,1]+I(clim.pca$scores[,1]^2))) #summary(lm(maine$hard.zero~clim.pca$scores[,2])) summary(lm(maine$coldplot~clim.pca$scores[,1])) #summary(lm(maine$coldplot~clim.pca$scores[,2])) ## g matrix estimate main.g<-o.dat%>% filter(Colony!="ted1" & Colony!="ted2" & Colony!="ted8" & Colony!="cent2")%>% droplevels() #levels(main.g$Colony) m1.100<-lmer(formula=treatment_recovery_s~1+(0+pretreat_Temp|Colony),REML=TRUE,data=main.g) summary(m1.100) Gmatrix.pac<-VarCorr(m1.100)$Colony;Gmatrix.pac bba<-Gmatrix.pac[1:4,1:4] maing<-princomp(bba) summary(maing) gload<-data.frame(maing$loadings[,1]) gload$temp<-c(-5,0,25,5) names(gload)[1]<-"loadings" ggplot(gload,aes(x=temp,y=loadings))+geom_point(size=5)+geom_line()+ylim(-.5,.65)+T+geom_hline(yintercept=0) # testing for pretemp treatment effect summary(aov(treatment_recovery_s~factor(pretreat_Temp)+Colony,new.dat3)) TukeyHSD(aov(treatment_recovery_s~factor(pretreat_Temp)+Colony,new.dat3)) summary(lmer(treatment_recovery_s~factor(pretreat_Temp)+(1|Colony),data=new.dat3)) #summary(lmer(treatment_recovery_s~factor(pretreat_Temp)+(1+factor(pretreat_Temp)|Colony),data=new.dat3)) ``` ## Redoing figures without VT ### Figure 4 ```{r, warning=FALSE} #recalculating cold tolerance nocol3<-ggplot(data=new.dat3,aes(x=pretreat_Temp,y=coldplot2,group=Colony,col=Colony))+geom_point(size=5)+geom_line(size=1.1)+xlab("Pre-treatment Temperature (°C)")+scale_x_continuous(labels=c(-5,0,5,25),breaks=c(-5,0,5,25))+theme(text=element_text(size=20))+T+scale_y_continuous(expression(paste("Cold Resistance ( ",CCRT[max], " - ", CCRT[obs],", seconds)")),limits=c(0,600),breaks=seq(0,850,100))+theme(legend.position="none")+scale_colour_manual(values=bcal) nocol3 #boxplots bp3<-ggplot(new.dat3,aes(y=coldplot2,x=factor(pretreat_Temp)))+geom_boxplot(fill="gray")+ylab("")+xlab("Pre-treatment Temperature(°C)")+T+scale_y_continuous(expression(paste("Cold Resistance ( ",CCRT[max], " - ", CCRT[obs],", seconds)")),limits=c(0,600),breaks=seq(0,850,100)) bp3 gplot2<-ggplot(gload,aes(x=temp,y=loadings))+geom_point(size=10)+scale_y_continuous(limit=c(-.6,.63),labels=seq(-1,1,.2),breaks=seq(-1,1,.2))+scale_x_continuous(limit=c(-5,25),labels=c(-5,0,5,25),breaks=c(-5,0,5,25))+geom_line(size=2)+xlab("Pre-treatment Temperature (°C)")+ylab(expression(paste("g"[max]," Loadings")))+T+geom_hline(yintercept=0,linetype=6,size=1.5) gplot2 grid.arrange(nocol3,bp3,gplot2,ncol=3) ``` ### Figure 5 correlation matrix ```{r correlation matrix, warning=FALSE} library(ggcorrplot) # presence absence sites comat<-cor(dbio2[,18:36]) b<-ggcorrplot(comat, hc.order = TRUE,lab=TRUE, type = "upper", outline.col = "white", ggtheme = ggplot2::theme_gray, colors = c("#6D9EC1", "white", "#E46726")) b+theme( axis.title.x = element_blank(), axis.title.y = element_blank(), panel.grid.major = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.ticks = element_blank(),) #common garden sites #comat<-cor(ksub3[,-1]) comat<-cor(maine[,23:33]) comat<-cor(maine[,c("MAT","Tmin","MDR")]) b2<-ggcorrplot(comat, hc.order = TRUE,lab=TRUE, type = "upper", outline.col = "white", ggtheme = ggplot2::theme_gray, colors = c("#6D9EC1", "white", "#E46726")) b2+theme( axis.title.x = element_blank(), axis.title.y = element_blank(), panel.grid.major = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.ticks = element_blank(),) ``` ### Figure 6 ```{r, warning=FALSE} basal<-new.dat3%>% filter(pretreat_Temp=="25") #number of oclonies droplevels(unique(basal$Colony)) MAT1<-ggplot(basal,aes(x=MAT,y=coldplot2))+geom_point(size=10)+ geom_smooth(method="lm",se=TRUE,size=4,colour="black")+T+scale_y_continuous(expression(paste("Basal Cold Resistance","\n"," (25°C; ",CCRT[max], " - ", CCRT[obs],", seconds)")),limits=c(-100,500),breaks=seq(0,800,100))+T+scale_shape_manual(values=c(15,16,17,18))+scale_x_continuous(breaks=seq(4.5,7,.5),limits=c(4.5,7))+xlab("") MAT1 h0<-new.dat3%>% filter(pretreat_Temp=="0") h0$hard.zero<-basal$treatment_recovery_s-h0$treatment_recovery_s MAT2<-ggplot(h0,aes(x=MAT,y=hard.zero))+geom_point(size=10)+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],", seconds)")),limits=c(-200,550),breaks=seq(0,800,100))+scale_x_continuous("Mean Annual Temperature (MAT, °C)",breaks=seq(4.5,7,.5),limits=c(4.5,7))+scale_shape_manual(values=c(15,16,17,18))+T+stat_smooth(method="lm",formula=y~x+I(x^2),size=4,colour="black",span=1.0)#+geom_smooth(method="lm",se=TRUE,size=4,colour="black") MAT2 #summary(lm(hard.zero~MAT+I(MAT^2),data=h0)) grid.arrange(MAT1,MAT2,nrow=2) ## 3rd panel maine$coldplot2<-max(maine$treatment_recovery_s)-maine$treatment_recovery_s cc2<-ggplot(maine,aes(x=coldplot2,y=hard.zero))+geom_point()+stat_smooth(method="lm")+T+scale_y_continuous(expression(paste("Cold Hardening (",CCRT["25°C"], " - ", CCRT["0°C"],", seconds)")),limits=c(-100,550),breaks=seq(0,800,100))+scale_x_continuous(expression(paste("Basal Cold Resistance","\n"," (25°C; ",CCRT[max], " - ", CCRT[obs],", seconds)")),limits=c(0,500),breaks=seq(0,800,100))+geom_point(size=10)+ geom_smooth(method="lm",se=TRUE,size=4,colour="black") cc2 grobmat<-arrangeGrob(MAT1,MAT2,nrow=2) #((MAT1/MAT2)|cc2)+plot_layout(widths = c(.75,1.25)) grid.arrange(grobmat,cc2,ncol=2) ##make cutoff for sara maine$group<-ifelse(maine$coldplot2