aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorKen Kellner <ken.kellner@gmail.com>2016-04-20 21:39:34 (GMT)
committerKen Kellner <ken.kellner@gmail.com>2016-04-20 21:39:34 (GMT)
commitb2206bef81fefa2c2ee7a34a664063cd6394ce82 (patch)
tree0278acfd1e9e9545a9ed8f1350aae0f46e91c03e
parent0f8bbccec26c1eb46febd9db511f470c924511ea (diff)
Fix error in calculation of sample seedling heights for figure. Fix error in assigning plots to treatments in deer analysis. Update manuscript to reflect corrections. Delete old manuscript.
-rw-r--r--analysis_deerbrowse.R2
-rw-r--r--analysis_rabbitbrowse.R2
-rw-r--r--figures/figure_height.R85
-rw-r--r--figures/figure_pellets.R2
-rw-r--r--figures/figure_transect.R2
-rw-r--r--tables/table_modeloutput_cat.R65
6 files changed, 96 insertions, 62 deletions
diff --git a/analysis_deerbrowse.R b/analysis_deerbrowse.R
index 7885977..8de62ef 100644
--- a/analysis_deerbrowse.R
+++ b/analysis_deerbrowse.R
@@ -53,6 +53,8 @@ for (i in 1:dim(surv)[1]){
ht.raw[i,(j-1)] <- ht.raw[i,(j-2)]
}}}}
ht.raw[,8] <- ht.raw[,7]
+ht.raw.de <- ht.raw
+
ht <- (ht.raw - mean(ht.raw,na.rm=TRUE)) / sd(ht.raw,na.rm=TRUE)
ht2.raw <- (ht.raw - mean(ht.raw,na.rm=TRUE))^2
diff --git a/analysis_rabbitbrowse.R b/analysis_rabbitbrowse.R
index adffa4c..ed6f1e5 100644
--- a/analysis_rabbitbrowse.R
+++ b/analysis_rabbitbrowse.R
@@ -41,6 +41,8 @@ for (i in 1:dim(surv)[1]){
ht.raw[i,(j-1)] <- ht.raw[i,(j-2)]
}}}}
ht.raw[,8] <- ht.raw[,7]
+ht.raw.rab <- ht.raw
+
ht <- (ht.raw - mean(ht.raw,na.rm=TRUE)) / sd(ht.raw,na.rm=TRUE)
ht2.raw <- (ht.raw - mean(ht.raw,na.rm=TRUE))^2
diff --git a/figures/figure_height.R b/figures/figure_height.R
index 9177250..ee76c4b 100644
--- a/figures/figure_height.R
+++ b/figures/figure_height.R
@@ -2,19 +2,27 @@
## Figure: Height effects on Browse Intensity ##
################################################
-#Run analysis code up cleanup of height data
-sample.hts <- seq(10,max(ht.raw,na.rm=T),0.2)
+sample.hts.de <- sort(as.numeric(na.omit(as.vector(ht.raw.de))))
-sample.hts2 <- (sample.hts-mean(sample.hts))^2
+sample.hts2.de <- (sample.hts.de-mean(sample.hts.de))^2
-hts.Z <- as.numeric(scale(sample.hts))
+hts.Z.de <- as.numeric(scale(sample.hts.de))
-hts2.Z <- as.numeric(scale(sample.hts2))
+hts2.Z.de <- as.numeric(scale(sample.hts2.de))
-op <- par()
+########################################################
+
+sample.hts.rab <- sort(as.numeric(na.omit(as.vector(ht.raw.rab))))
+
+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))
+
+#########################################################
-########################################
#Separate figures
require(extrafont)
@@ -29,10 +37,10 @@ par(mfrow=c(3,1),
###
rp <- browserabbit.output.cat$mean
-lin.pred <- rp$tau[1] - (rp$b.distance*0 + rp$b.species*0 +rp$b.comp*0 + rp$b.ht*hts.Z
- + rp$b.ht2*hts2.Z+rp$b.species)
-lin.pred.sum <- rp$tau[1] - (rp$b.distance*0 + rp$b.species*0 + rp$b.comp*0 + rp$b.ht*hts.Z
- + rp$b.ht2*hts2.Z + rp$b.season +rp$b.species)
+lin.pred <- rp$tau[1] - (rp$b.species*0 +rp$b.comp*0 + rp$b.ht*hts.Z.rab
+ + rp$b.ht2*hts2.Z.rab+rp$b.species)
+lin.pred.sum <- rp$tau[1] - (rp$b.species*0 + rp$b.comp*0 + rp$b.ht*hts.Z.rab
+ + rp$b.ht2*hts2.Z.rab + rp$b.season +rp$b.species)
#Probability of being browsed at all by rabbits
prab<- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
@@ -40,22 +48,23 @@ prab.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
#Max value
mx <-mxrab<- max(prab)
-mx.index <-mx.index.rab <- sample.hts[which(prab==mx)]
+mx.index <-mx.index.rab <- sample.hts.rab[which(prab==mx)]
-plot(sample.hts,prab,type='l',lwd=2,xaxt='n',yaxt='n')
-text(15,0.003,"Rabbit",cex=1.5,adj=c(0,NA))
-segments(x0=mx.index,y0=mx,x1=mx.index,y1=(mx-0.25*mx))
-text(mx.index,(mx-0.35*mx),mx.index)
-axis(2,at=c(0,0.004,0.008),labels=c(0,0.004,0.008))
-lines(sample.hts,prab.sum,lwd=2,lty=2)
-legend('topleft',legend=c('Oct - April','May - Sept'),lty=c(1,2),lwd=2)
+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)))
+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),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)
###
dp <- browsedeer.output.cat$mean
-lin.pred <- dp$tau[1] - (dp$b.distance*0 + dp$b.species*0 + dp$b.comp*0 + dp$b.ht*hts.Z
- + dp$b.ht2*hts2.Z+dp$b.species)
-lin.pred.sum <- dp$tau[1] - (dp$b.distance*0 + dp$b.species*0 + dp$b.comp*0 + dp$b.ht*hts.Z
- + dp$b.ht2*hts2.Z + dp$b.season +dp$b.species)
+lin.pred <- dp$tau[1] - (dp$b.species*0 + dp$b.comp*0 + dp$b.ht*hts.Z.de
+ + dp$b.ht2*hts2.Z.de + dp$b.species)
+lin.pred.sum <- dp$tau[1] - (dp$b.species*0 + dp$b.comp*0 + dp$b.ht*hts.Z.de
+ + dp$b.ht2*hts2.Z.de + dp$b.season +dp$b.species)
#Probability of being browsed at all by deer
pdeer<- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
@@ -63,30 +72,32 @@ pdeer.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
#Max Value
mx <-mxrab<- max(pdeer)
-mx.index <-mx.index.rab <- sample.hts[which(pdeer==mx)]
+mx.index <-mx.index.rab <- sample.hts.de[which(pdeer==mx)]
-plot(sample.hts,pdeer,type='l',lwd=2,xaxt='n',yaxt='n',lty=1)
+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)))
axis(2,at=c(0,0.10,0.20),labels=c(0,0.10,0.20))
-text(15,mx-0.25*mx,"Deer",cex=1.5,adj=c(0,NA))
-segments(x0=mx.index,y0=mx,x1=mx.index,y1=(mx-0.05*mx))
-text(mx.index,(mx-0.12*mx),mx.index)
-lines(sample.hts,pdeer.sum,lwd=2,lty=2)
+text(15,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),mx.index[1])
+lines(sample.hts.de,pdeer.sum,lwd=2,lty=2)
###
ip <- browseinsect.output.cat$mean
-lin.pred <- ip$tau[1] - (ip$b.distance*0 + ip$b.species*0 + ip$b.comp*0
- + ip$b.ht*hts.Z +ip$b.species)
-lin.pred.sum <- ip$tau[1] - (ip$b.distance*0 + ip$b.species*0 + ip$b.comp*0
- + ip$b.ht*hts.Z + ip$b.season*1 + ip$b.species)
+lin.pred <- ip$tau[1] - (ip$b.species*0 + ip$b.comp*0
+ + ip$b.ht*hts.Z.rab +ip$b.species)
+lin.pred.sum <- ip$tau[1] - (ip$b.species*0 + ip$b.comp*0
+ + ip$b.ht*hts.Z.rab + ip$b.season*1 + ip$b.species)
#Probability of being browsed at all by insects
pins <- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
pins.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
-plot(sample.hts,pins.sum,type='l',lwd=2,ylim=c(0,0.35),lty=2,yaxt='n')
-text(15,0.15,"Insect",cex=1.5,adj=c(0,NA))
-axis(2,at=c(0,0.15,0.30),labels=c('0','0.15','0.30'))
-lines(sample.hts,pins,lty=1,lwd=2)
+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)))
+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)
title(xlab="Seedling Height (cm)",ylab="Probability of Browse Damage",outer=T,line=2.5,
cex.lab=1.5)
diff --git a/figures/figure_pellets.R b/figures/figure_pellets.R
index af3dd57..84d2eb3 100644
--- a/figures/figure_pellets.R
+++ b/figures/figure_pellets.R
@@ -74,7 +74,7 @@ par(mgp=c(2.5,1,0),mar = c(2.5,3.7,1.3,0.5) + 0.1)
x <- barplot(plot.data,beside=TRUE,names=c('Rabbit','Deer'),
ylab=expression("Pellets "~ m^{-2}),ylim=c(0,0.6),
legend.text=c('Forest','Edge','Harvest'),
- args.legend=list('topright'))
+ args.legend=list(x='topright'))
wd=0.2
for(i in 1:3){
diff --git a/figures/figure_transect.R b/figures/figure_transect.R
index 9b94cd5..f300fcf 100644
--- a/figures/figure_transect.R
+++ b/figures/figure_transect.R
@@ -111,7 +111,7 @@ par(mgp=c(2.5,1,0),mar = c(2.5,3.5,1.3,0.5) + 0.1)
barplot(cbind(rab,de,ins),beside=T, ylab="Proportion of Seedlings Browsed",
xlab="",names=c('Rabbit','Deer','Insect'),
legend.text=c('Forest','Edge','Harvest'),
- args.legend=list(x=5,y=0.7),ylim=c(0,0.85))
+ args.legend=list(x='topleft'),ylim=c(0,0.85))
box()
rab.se <- sqrt(rab*(1-rab))/sqrt(length(ever.browsed.rabbit))
diff --git a/tables/table_modeloutput_cat.R b/tables/table_modeloutput_cat.R
index 1b7a546..ce95ec8 100644
--- a/tables/table_modeloutput_cat.R
+++ b/tables/table_modeloutput_cat.R
@@ -2,44 +2,62 @@
## Cumulative Link Model Output Table ##
########################################
+########################################
+## Cumulative Link Model Output Table ##
+########################################
+
#Decimal places to include
rd <- 2
-#Format deer model output
-td <- browsedeer.output.cat
-include <- c(8,9,10)
-mn1 <-unlist(td$mean)[include]
-q2.5 <- round(unlist(td$q2.5),rd)[include]
-q97.5 <- round(unlist(td$q97.5),rd)[include]
-ci1 <- paste("(",sprintf("%.2f",q2.5),",",sprintf("%.2f",q97.5),')',sep="")
-f1 <- sprintf("%.2f",round(unlist(td$f),rd))[include]
-ov1 <- unlist(td$overlap0)[include]
-
#Format rabbit model output (1 missing value)
tr <- browserabbit.output.cat
-include <- c(8,9,10)
-mn2 <-unlist(tr$mean)[include]
+include <- c(4:17)
+mn1 <-unlist(tr$mean)[include]
q2.5 <- round(unlist(tr$q2.5),rd)[include]
q97.5 <- round(unlist(tr$q97.5),rd)[include]
+ci1 <- paste("(",sprintf("%.2f",q2.5),",",sprintf("%.2f",q97.5),')',sep="")
+f1 <- sprintf("%.2f",round(unlist(tr$f),rd))[include]
+ov1 <- unlist(tr$overlap0)[include]
+
+#Format deer model output
+td <- browsedeer.output.cat
+include <- c(4:17)
+mn2 <-unlist(td$mean)[include]
+q2.5 <- round(unlist(td$q2.5),rd)[include]
+q97.5 <- round(unlist(td$q97.5),rd)[include]
ci2 <- paste("(",sprintf("%.2f",q2.5),",",sprintf("%.2f",q97.5),')',sep="")
-f2 <- sprintf("%.2f",round(unlist(tr$f),rd))[include]
-ov2 <- unlist(tr$overlap0)[include]
+f2 <- sprintf("%.2f",round(unlist(td$f),rd))[include]
+ov2 <- unlist(td$overlap0)[include]
-#Format insect model output (2 missing values)
+#Format insect model output (1 missing value)
ti <- browseinsect.output.cat
-include <- c(8,9,10)
+include <- c(4:17)
mn3 <-unlist(ti$mean)[include]
+mn3[8] <- NA
q2.5 <- round(unlist(ti$q2.5),rd)[include]
q97.5 <- round(unlist(ti$q97.5),rd)[include]
ci3 <- paste("(",sprintf("%.2f",q2.5),",",sprintf("%.2f",q97.5),')',sep="")
+ci3 <- c(ci3[1:7],'-',ci3[9:14])
f3 <- sprintf("%.2f",round(unlist(ti$f),rd))[include]
+f3 <- c(f3[1:7],'-',f3[9:14])
ov3 <- unlist(ti$overlap0)[include]
+ov3 <- c(ov3[1:7],TRUE,ov3[9:14])
#Combine into raw data frame
ov <- cbind(ov1,ov1,ov1,ov2,ov2,ov2,ov3,ov3,ov3)
+
+ov[which(!ov,arr.ind=T)] <- TRUE
+
+#ov <- cbind(ov1,ov1,ov2,ov2,ov3,ov3)
tab <- data.frame(mn1,ci1,f1,mn2,ci2,f2,mn3,ci3,f3)
-row.names(tab) <- c('Edge-Forest','Harvest-Forest','Harvest-Edge')
+#tab <- data.frame(mn1,ci1,mn2,ci2,mn3,ci3)
+row.names(tab) <- c('Plot SD','Subplot SD','Seedling SD',
+ 'White Oak','Competition','Growing Season','Height','Height\\textsuperscript{2}',
+ 'Edge','Harvest','Harvest-Edge',
+ 'Prior Rabbit','Prior Deer','Prior Insect')
colnames(tab) <- rep(c('Mean','95% CI','\\textit{f}'),3)
+#colnames(tab) <- rep(c('Mean','95% CI'),3)
+
#Correct spacing of mean values so decimal point locations match in LaTeX
hold <- tab
@@ -69,14 +87,15 @@ for (i in 1:length(ov[,1])){
#Manually comment out cline
require(Hmisc)
capture.output(
- latex(tab,title='Contrast',file="",
- cgroup = c('Deer','Rabbit','Insect'), n.cgroup = c(3,3,3),
+ latex(tab,title='Parameter',file="",
+ rgroup = c("Random", "Fixed"), n.rgroup=c(3, nrow(tab)-3),
+ cgroup = c('Rabbit','Deer','Insect'), n.cgroup = c(3,3,3),
booktabs=TRUE,
col.just = strsplit("lcclcclcc", "")[[1]],
- caption="Differences between transect positions based on cumulative link models where
- position on transect was categorical instead of continous.",
- insert.bottom = '{\\scriptsize *95\\% CI does not include 0}',
- label = "modeloutputcat",
+ caption="Parameter estimates from the cumulative link models of browse damage (index) by white-tailed deer,
+ eastern cottontail rabbits, and insects on oak seedlings.",
+ #insert.bottom = '{\\scriptsize *95\\% CI does not include 0}',
+ label = "modeloutput",
size = 'scriptsize'),
file="paper/table_modeloutput_cat.tex")