aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorKen Kellner <ken@kenkellner.com>2017-02-17 19:33:12 (GMT)
committerKen Kellner <ken@kenkellner.com>2017-02-17 19:33:12 (GMT)
commit90462ada4e3bd3d77c5195248e21f2fc2fb56449 (patch)
treecfde53a6143e6afb68755e08b8a85fa9ba4fb628
parentc32b75cc40b30777bb7186cd5f7391dc1197dfac (diff)
Add combined figure for talks.
-rw-r--r--figures/fig3_height.R45
-rw-r--r--figures/fig_talk_herbivory.R197
2 files changed, 229 insertions, 13 deletions
diff --git a/figures/fig3_height.R b/figures/fig3_height.R
index b3c35a3..c52c5d6 100644
--- a/figures/fig3_height.R
+++ b/figures/fig3_height.R
@@ -45,9 +45,8 @@ sample.hts2.rab <- (sample.hts.rab-mean(sample.hts.rab))^2
hts.Z.rab <- as.numeric(scale(sample.hts.rab))
hts2.Z.rab <- as.numeric(scale(sample.hts2.rab))
-#########################################################
+########################################################
#Set up figure container
-
tiff(filename="figures/Fig3.tiff",width=84,height=84,units="mm",res=1200, pointsize=8,
compression = "lzw",type='cairo')
@@ -58,6 +57,26 @@ par(mfrow=c(3,1),
require(extrafont)
myfont <- 'Times'
+xshow <- c('n','n',NULL)
+cols <- c('black','black','black')
+
+########################################################
+
+png(filename="figures/herbivory.png",width=6,height=3.8,units="in",res=400, pointsize=10,
+ type='cairo')
+
+par(mfrow=c(1,3),
+ oma = c(4,4,0,0) + 0.1,
+ mar = c(0.5,1,1,1) + 0.1)
+
+#cols <- c('black','firebrick1','dodgerblue')
+
+xshow <- c(NULL,NULL,NULL)
+
+cols<- c(rgb(128,194,0,max=255),
+ rgb(244,124,66,max=255),
+ rgb(241,194,50,max=255))
+
########################################################
#Rabbit sub-figure
@@ -76,14 +95,14 @@ prab.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
mx <- mxrab <- max(prab)
mx.index <-mx.index.rab <- sample.hts.rab[which(prab==mx)]
-plot(sample.hts.rab,prab,type='l',lwd=2,xaxt='n',yaxt='n',ylim=c(0,0.002),
- xlim=c(0,max(sample.hts.rab,na.rm=T)))
+plot(sample.hts.rab,prab,type='l',lwd=2,xaxt=xshow[1],yaxt='n',ylim=c(0,0.002),
+ xlim=c(0,max(sample.hts.rab,na.rm=T)),col=cols[1])
text(15,0.0018,"Rabbit",cex=1.5,adj=c(0,NA))
segments(x0=mx.index[1],y0=mx,x1=mx.index[1],y1=(mx-0.20*mx))
text(mx.index[1],(mx-0.35*mx),round(mx.index[1]))
axis(2,at=c(0,0.001,0.002),labels=c(0,0.001,0.002))
-lines(sample.hts.rab,prab.sum,lwd=2,lty=2)
-legend('topright',legend=c('Oct - April','May - Sept'),lty=c(1,2),lwd=2)
+lines(sample.hts.rab,prab.sum,lwd=2,lty=2,col=cols[1])
+legend('topright',legend=c('Oct - April','May - Sept'),lty=c(1,2),lwd=1)
###################################################
#Deer sub-figure
@@ -103,13 +122,13 @@ pdeer.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
mx <-mxrab<- max(pdeer)
mx.index <-mx.index.rab <- sample.hts.de[which(pdeer==mx)]
-plot(sample.hts.de,pdeer,type='l',lwd=2,xaxt='n',yaxt='n',lty=1,
- xlim=c(0,max(sample.hts.rab,na.rm=T)))
+plot(sample.hts.de,pdeer,type='l',lwd=2,xaxt=xshow[2],yaxt='n',lty=1,
+ xlim=c(0,max(sample.hts.rab,na.rm=T)),col=cols[2])
axis(2,at=c(0,0.10,0.20),labels=c(0,0.10,0.20))
-text(15,0.25,"Deer",cex=1.5,adj=c(0,NA))
+text(7,0.25,"Deer",cex=1.5,adj=c(0,NA))
segments(x0=mx.index[1],y0=mx[1],x1=mx.index[1],y1=(mx-0.05*mx))
text(mx.index[1],(mx-0.15*mx),round(mx.index[1]))
-lines(sample.hts.de,pdeer.sum,lwd=2,lty=2)
+lines(sample.hts.de,pdeer.sum,lwd=2,lty=2,col=cols[2])
##############################################################
#Insect sub-figure
@@ -125,11 +144,11 @@ lin.pred.sum <- ip$tau[1] - (ip$b.species*0 + ip$b.comp*0
pins <- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
pins.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
-plot(sample.hts.rab,pins.sum,type='l',lwd=2,lty=2,ylim=c(0,0.85),yaxt='n',
- xlim=c(0,max(sample.hts.rab,na.rm=T)))
+plot(sample.hts.rab,pins.sum,type='l',lwd=2,lty=2,ylim=c(0,0.85),yaxt='n',xaxt=xshow[3],
+ xlim=c(0,max(sample.hts.rab,na.rm=T)),col=cols[3])
text(15,0.75,"Insect",cex=1.5,adj=c(0,NA))
axis(2,at=c(0,0.40,0.80),labels=c('0','0.40','0.80'))
-lines(sample.hts.rab,pins,lty=1,lwd=2)
+lines(sample.hts.rab,pins,lty=1,lwd=2,col=cols[3])
title(xlab="Seedling Height (cm)",ylab="Probability of Browse Damage",outer=T,line=2.5,
cex.lab=1.5)
diff --git a/figures/fig_talk_herbivory.R b/figures/fig_talk_herbivory.R
new file mode 100644
index 0000000..c36a0b6
--- /dev/null
+++ b/figures/fig_talk_herbivory.R
@@ -0,0 +1,197 @@
+source('function_format_data.R')
+
+#Initial formatting on raw data
+seedling <- format.seedling('data/hee_seedling_master.csv')
+
+#############################################################
+#Seedlings browsed at least once by rabbits, by treatment
+keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid<49)
+browse <- seedling$browseother[keep,]
+browse <- browse + 1
+
+ever.browsed.rabbit <- rep(0,dim(browse)[1])
+for (i in 1:length(ever.browsed.rabbit)){
+ if(!all(is.na(browse[i,]))){
+ test = mean(browse[i,],na.rm=TRUE)
+ if(test>1){ever.browsed.rabbit[i]=1}
+ }
+}
+
+#Get treatment by seedling
+edge <- c(rep(c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0),3))
+harvest <- c(rep(c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),3))
+control <- c(rep(c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1),3))
+
+seedling.covs <- seedling$seedling.data[keep,]
+seed.plotcode <- seedling.covs$plotid
+other.harvest <- other.edge <- other.control <- rep(0,dim(browse)[1])
+for(i in 1:length(ever.browsed.rabbit)){
+ if(harvest[seed.plotcode[i]] == 1){other.harvest[i] <- 1}
+ if(edge[seed.plotcode[i]] == 1){other.edge[i] <- 1}
+ if(control[seed.plotcode[i]] == 1){other.control[i] <- 1}
+}
+
+#Bundle rabbit data
+rab = c(mean(ever.browsed.rabbit[other.control==1]),mean(ever.browsed.rabbit[other.edge==1]),
+ mean(ever.browsed.rabbit[other.harvest==1]))
+
+###############################################
+#Seedlings browsed at least once by insects, by treatment
+
+browse <- seedling$lfdmg[keep,]
+browse <- browse + 1
+
+ever.browsed.insect <- rep(0,dim(browse)[1])
+for (i in 1:length(ever.browsed.insect)){
+ if(!all(is.na(browse[i,]))){
+ test = mean(browse[i,],na.rm=TRUE)
+ if(test>1){ever.browsed.insect[i]=1}
+ }
+}
+
+#Bundle data
+ins = c(mean(ever.browsed.insect[other.control==1]),mean(ever.browsed.insect[other.edge==1]),
+ mean(ever.browsed.insect[other.harvest==1]))
+
+###############################################
+#Seedlings browsed at least once by deer
+
+exclude <- 1 - seedling$plot.data$herbivory
+
+#Plots to keep (deer not excluded, not in shelterwoods)
+plots <- c(1:54)
+keep.plots <- which(plots<49&exclude==0)
+keep.plots <- keep.plots[-2]
+keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid%in%keep.plots)
+browse <- seedling$browsedeer[keep,]
+browse <- browse + 1
+
+ever.browsed.deer <- rep(0,dim(browse)[1])
+for (i in 1:length(ever.browsed.deer)){
+ if(!all(is.na(browse[i,]))){
+ test = mean(browse[i,],na.rm=TRUE)
+ if(test>1){ever.browsed.deer[i]=1}
+ }
+}
+
+seedling.covs <- seedling$seedling.data[keep,]
+seed.sitecode <- seedling.covs$siteid
+seed.plotcode.raw <- seedling.covs$plotid
+seed.plotcode <- rep(1,length(which(seed.plotcode.raw==unique(seed.plotcode.raw)[1])))
+for (i in 2:24){
+ add <- rep(i,length(which(seed.plotcode.raw==unique(seed.plotcode.raw)[i])))
+ seed.plotcode <- c(seed.plotcode,add)
+}
+
+#Get treatment by seedlings
+plot.covs <- seedling$plot.data[keep.plots,]
+
+harvest <- c(1,1,1,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,0,0,0,0)
+edge <- c(0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0)
+control <- c(0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1)
+
+deer.harvest <- deer.edge <- deer.control <- rep(0,dim(browse)[1])
+for(i in 1:length(ever.browsed.deer)){
+ if(harvest[seed.plotcode[i]] == 1){deer.harvest[i] <- 1}
+ if(edge[seed.plotcode[i]] == 1){deer.edge[i] <- 1}
+ if(control[seed.plotcode[i]] == 1){deer.control[i] <- 1}
+}
+
+#Bundle data
+de = c(mean(ever.browsed.deer[deer.control==1]),mean(ever.browsed.deer[deer.edge==1]),
+ mean(ever.browsed.deer[deer.harvest==1]))
+
+############################################################
+#Set up figure container
+
+png(filename="figures/herbivory_transect.png",width=6,height=3.8,units="in",res=400, pointsize=8,
+ type='cairo')
+
+par(mfrow=c(1,2),
+ oma = c(3,0,0,0) + 0.1,
+ mar = c(0,4,1,0.5) + 0.1)
+
+cols <- rev(c(rgb(red=244,green=125,blue=66, maxColorValue=255),
+ rgb(red=241,green=194,blue=50, maxColorValue=255),
+ rgb(red=141,green=213,blue=18, maxColorValue=255)))
+
+#Set up barplot
+barplot(cbind(rab,de,ins),beside=T, ylab="Proportion of Seedlings Browsed",
+ xlab="",names=c('Rabbit','Deer','Insect'),col=cols[1:3],
+ legend.text=c('Forest','Edge','Harvest'),
+ args.legend=list(x=5,y=0.7),ylim=c(0,0.85))
+
+
+#Calculate SE for each
+rab.se <- sqrt(rab*(1-rab))/sqrt(length(ever.browsed.rabbit))
+de.se <- sqrt(de*(1-de))/sqrt(length(ever.browsed.deer))
+ins.se <- sqrt(ins*(1-ins))/sqrt(length(ever.browsed.insect))
+
+structure <- matrix(c(1.5,2.5,3.5,5.5,6.5,7.5,9.5,10.5,11.5),byrow=T,nrow=3)
+
+#Significant differences based on browse model
+sigs <- matrix(c('A','AB','B','A','A','A','A','A','B'),byrow=T,nrow=3)
+
+#Error bars
+wd <- 0.2
+sep <- 0.03
+for (i in 1:3){
+ segments(x0=structure[1,i],y0=rab[i],x1=structure[1,i],y1=rab[i]+1.96*rab.se[i],lwd=1)
+ segments(x0=structure[2,i],y0=de[i],x1=structure[2,i],y1=de[i]+1.96*de.se[i],lwd=1)
+ segments(x0=structure[3,i],y0=ins[i],x1=structure[3,i],y1=ins[i]+1.96*ins.se[i],lwd=1)
+
+ segments(x0=structure[1,i]-wd,y0=rab[i]+1.96*rab.se[i],x1=structure[1,i]+wd,y1=rab[i]+1.96*rab.se[i])
+ segments(x0=structure[2,i]-wd,y0=de[i]+1.96*de.se[i],x1=structure[2,i]+wd,y1=de[i]+1.96*de.se[i])
+ segments(x0=structure[3,i]-wd,y0=ins[i]+1.96*ins.se[i],x1=structure[3,i]+wd,y1=ins[i]+1.96*ins.se[i])
+
+ text(structure[1,i],rab[i]+1.96*rab.se[i]+sep,sigs[1,i])
+ text(structure[2,i],de[i]+1.96*de.se[i]+sep,sigs[2,i])
+ text(structure[3,i],ins[i]+1.96*ins.se[i]+sep,sigs[3,i])
+}
+
+#####################################
+
+par(mar = c(0,4.5,1,0) + 0.1)
+
+#Read in raw data
+pellet.raw <- read.csv('data/hee_pellet.csv',header=T)
+pellet <- pellet.raw[pellet.raw$Pos!='sh',]
+
+#Calculate means/se by treatment and herbivore
+deer.means <- c(mean(pellet$Deer[pellet$Treat=="Matrix"]),
+ mean(pellet$Deer[pellet$Treat=="Edge"]),
+ mean(pellet$Deer[pellet$Treat=="Harvest"],na.rm=TRUE))/5
+
+deer.se <- c(sd(pellet$Deer[pellet$Treat=="Matrix"])/sqrt(length(pellet$Deer[pellet$Treat=="Matrix"])),
+ sd(pellet$Deer[pellet$Treat=="Edge"])/sqrt(length(pellet$Deer[pellet$Treat=="Edge"])),
+ sd(pellet$Deer[pellet$Treat=="Harvest"],na.rm=T)
+ /sqrt(length(pellet$Deer[pellet$Treat=="Harvest"])))
+
+rabbit.means <- c(mean(pellet$Rabbit[pellet$Treat=="Matrix"]),
+ mean(pellet$Rabbit[pellet$Treat=="Edge"]),
+ mean(pellet$Rabbit[pellet$Treat=="Harvest"],na.rm=TRUE))/5
+
+rabbit.se <- c(sd(pellet$Rabbit[pellet$Treat=="Matrix"])/sqrt(length(pellet$Rabbit[pellet$Treat=="Matrix"])),
+ sd(pellet$Rabbit[pellet$Treat=="Edge"])/sqrt(length(pellet$Rabbit[pellet$Treat=="Edge"])),
+ sd(pellet$Rabbit[pellet$Treat=="Harvest"],na.rm=T)
+ /sqrt(length(pellet$Rabbit[pellet$Treat=="Harvest"])))
+
+plot.data <- cbind(rabbit.means,deer.means)
+
+x <- barplot(plot.data,beside=TRUE,names=c('Rabbit','Deer'),
+ ylab=expression("Pellets "~ m^{-2}),ylim=c(0,0.6),col=cols[1:3])
+
+#Plot error bars
+wd=0.2
+for(i in 1:3){
+ segments(x[i,1],rabbit.means[i],x[i,1],rabbit.means[i]+rabbit.se[i])
+ segments(x[i,1]-wd,rabbit.means[i]+rabbit.se[i],x[i,1]+wd,rabbit.means[i]+rabbit.se[i])
+ segments(x[i,2],deer.means[i],x[i,2],deer.means[i]+deer.se[i])
+ segments(x[i,2]-wd,deer.means[i]+deer.se[i],x[i,2]+wd,deer.means[i]+deer.se[i])
+}
+
+#Add significance letters based on randomization test
+text(x=c(1.5,2.5,3.5,5.5,6.5,7.5),y=c(rabbit.means+rabbit.se,deer.means+deer.se)+0.02,
+ c('A','B','B','A','A','A'))
+
+dev.off() \ No newline at end of file