aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorKen Kellner <ken.kellner@gmail.com>2016-04-20 20:19:13 (GMT)
committerKen Kellner <ken.kellner@gmail.com>2016-04-20 20:19:13 (GMT)
commit0f8bbccec26c1eb46febd9db511f470c924511ea (patch)
treea7d306d0c3abf9233c3d33250ef091998fbb535c
parent826dd49b1b27fd516e786f2be7d57b02fa6efb9c (diff)
Correct edge/harvest classification error. Add supplements.
-rw-r--r--analysis_deerbrowse.R4
-rw-r--r--analysis_rabbitbrowse.R2
-rw-r--r--figures/figure_pellets.R2
-rw-r--r--figures/figure_transect.R22
-rw-r--r--supplements/supplement_s1.Rmd64
-rw-r--r--supplements/supplement_s2.Rmd114
6 files changed, 200 insertions, 8 deletions
diff --git a/analysis_deerbrowse.R b/analysis_deerbrowse.R
index 3c6fb56..7885977 100644
--- a/analysis_deerbrowse.R
+++ b/analysis_deerbrowse.R
@@ -83,8 +83,8 @@ nplots <- length(keep.plots)
plot.covs <- seedling$plot.data[keep.plots,]
-edge <- c(rep(c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0),3),rep(0,6))
-harvest <- c(rep(c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),3),rep(0,6))
+edge <- c(rep(c(0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0),3),rep(0,6))[keep.plots]
+harvest <- c(rep(c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),3),rep(0,6))[keep.plots]
plot.sitecode <- plot.covs$siteid
diff --git a/analysis_rabbitbrowse.R b/analysis_rabbitbrowse.R
index 4025c80..adffa4c 100644
--- a/analysis_rabbitbrowse.R
+++ b/analysis_rabbitbrowse.R
@@ -153,4 +153,6 @@ require(jagsUI)
browserabbit.output.cat <- jags(data=jags.data,inits=inits,parameters.to.save=params,model.file=modFile,
n.chains=3,n.iter=12000,n.burnin=8000,n.thin=40,parallel=TRUE)
+browserabbit.output.cat <- update(browserabbit.output.cat,n.iter=4000,n.thin=40)
+
save(browserabbit.output.cat,file="output/browserabbit_output_cat.Rdata")
diff --git a/figures/figure_pellets.R b/figures/figure_pellets.R
index 1b029d2..af3dd57 100644
--- a/figures/figure_pellets.R
+++ b/figures/figure_pellets.R
@@ -69,7 +69,7 @@ require(extrafont)
myfont <- 'Times'
pdf(file="paper/figure_pellets.pdf",width=5,height=3.9,family=myfont,pointsize=12)
-par(mgp=c(3,1,0),mar = c(2.3,4.2,1.3,0.5) + 0.1)
+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),
diff --git a/figures/figure_transect.R b/figures/figure_transect.R
index d7a8a5d..9b94cd5 100644
--- a/figures/figure_transect.R
+++ b/figures/figure_transect.R
@@ -106,12 +106,12 @@ require(extrafont)
myfont <- "Times"
pdf(file="paper/figure_transect.pdf",width=5,height=3.9,family=myfont,pointsize=12)
-par(mgp=c(3,1,0),mar = c(2.5,4.2,1.3,0.5) + 0.1)
+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.8))
+ args.legend=list(x=5,y=0.7),ylim=c(0,0.85))
box()
rab.se <- sqrt(rab*(1-rab))/sqrt(length(ever.browsed.rabbit))
@@ -124,10 +124,22 @@ 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)
+sigs <- matrix(c('A','AB','B','A','A','A','A','A','B'),byrow=T,nrow=3)
+
+wd <- 0.2
+sep <- 0.03
for (i in 1:3){
- segments(x0=structure[1,i],y0=rab[i]-1.96*rab.se[i],x1=structure[1,i],y1=rab[i]+1.96*rab.se[i],lwd=2)
- segments(x0=structure[2,i],y0=de[i]-1.96*de.se[i],x1=structure[2,i],y1=de[i]+1.96*de.se[i],lwd=2)
- segments(x0=structure[3,i],y0=ins[i]-1.96*ins.se[i],x1=structure[3,i],y1=ins[i]+1.96*ins.se[i],lwd=2)
+ 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])
}
dev.off()
diff --git a/supplements/supplement_s1.Rmd b/supplements/supplement_s1.Rmd
new file mode 100644
index 0000000..38c5541
--- /dev/null
+++ b/supplements/supplement_s1.Rmd
@@ -0,0 +1,64 @@
+---
+title: "Supplement S1: Generation of Competition Index"
+author: "Kenneth F. Kellner"
+date: "April 20, 2016"
+output: pdf_document
+---
+
+For each visit to each subplot during the study, we measured five interspecific competition-related variables:
+
+1. Percent cover of herbaceous vegetation (ordinal scale 0-6)
+2. Percent cover of woody vegetation < 2 m in height (ordinal scale 0-6)
+3. Density of woody stems 0-50 cm in height
+4. Density of woody stems 50-100 cm in height
+5. Density of woody stems >100 cm in height
+
+These variables were stored in data frame *comp*:
+
+```{r, eval=FALSE}
+head(comp)
+```
+
+```{r, echo=FALSE}
+comp <- read.csv('../data/competition.csv',header=TRUE)[,c(1,5:10)]
+compdisp <- comp[c(1,2,3,5,6,7),]
+knitr::kable(head(compdisp))
+```
+
+The ordinal percent cover variables were converted back to a percent scale by assigning them to the midpoint of the ordinal category percent range:
+
+```{r}
+conv <- c(0.5,3,15,37.5,62.5,85,97.5)
+comp$herb <- conv[comp$herb+1]
+comp$woody <- conv[comp$woody+1]
+```
+
+We then ran a principal components analysis on the five variables:
+
+```{r}
+pc <- prcomp(formula=~herb+woody+stems10.50+stems50.100+stems100,data=comp[,3:7],
+ center=TRUE,scale=TRUE,na.action=na.exclude)
+summary(pc)
+```
+
+Based on the broken-stick criterion, we retained on the first two principal components, which together explained 77% of the variability. We rotated then the principal components using varimax:
+
+```{r}
+v <- varimax(pc$rotation[,1:2])
+v
+```
+
+We interpreted the first principal component to be a general metric of woody competition, and the second to be a metric of herbaceous competition. Since we were primarily interested in woody plant competition with oak in this study, we selected the first principal component as our index of competition.
+
+Finally, we calculated the competition index for each subplot using the loadings we obtained above:
+
+```{r}
+compindex <- (apply(comp[,3:7],2,scale) %*% (-1*v$loadings))[,1]
+```
+
+```{r,echo=FALSE}
+comp <- read.csv('../data/competition.csv',header=TRUE)[,c(1,5:10)]
+knitr::kable(head(cbind(comp,compindex)[c(1,2,3,5,6,7),]))
+```
+
+Larger values of *compindex* correspond to plots with a greater amount of woody competition. \ No newline at end of file
diff --git a/supplements/supplement_s2.Rmd b/supplements/supplement_s2.Rmd
new file mode 100644
index 0000000..1e01381
--- /dev/null
+++ b/supplements/supplement_s2.Rmd
@@ -0,0 +1,114 @@
+---
+title: "Supplement S2: BUGS Model Code"
+author: "Kenneth F. Kellner"
+date: "April 20, 2016"
+output: pdf_document
+---
+
+```{r, eval=FALSE}
+model {
+ model {
+
+ #Likelihood
+
+ #Plot random effect
+ for (i in 1:nplots){
+ plot.mean[i] ~ dnorm(0,plot.tau)
+ }
+
+ #Subplot random effect
+ for (i in 1:nsubplots){
+ subplot.mean[i] ~ dnorm(0,subplot.tau)
+ }
+
+ for (i in 1:nseedlings){
+
+ #Seedling random effect
+ seed.mean[i] ~ dnorm(0,seed.tau)
+
+ #Iterate through samples
+ for (j in 1:nsamples[i]){
+
+ #Linear predictor
+ mu[i,j] <- seed.mean[i]
+ + plot.mean[seed.plotcode[i]]
+ + subplot.mean[seed.sitecode[i]]
+ + b.species*species[i]
+ + b.comp*comp[seed.plotcode[i],j]
+ + b.ht*ht[i,j] + b.ht2*ht2[i,j]
+ + b.season*season[j]
+ + b.deer*browse.deer[i,j]
+ + b.insect*browse.insect[i,j]
+ + b.rabbit*browse.rabbit[i,j]
+ + b.edge*edge[seed.plotcode[i]]
+ + b.harvest*harvest[seed.plotcode[i]]
+
+ #Cumulative logistic probabilities
+ logit(Q[i,j,1]) <- tau[1] - mu[i,j]
+ p[i,j,1] <- Q[i,j,1]
+
+ for (k in 2:3){
+ logit(Q[i,j,k]) <- tau[k] - mu[i,j]
+ p[i,j,k] <- Q[i,j,k] - Q[i,j,k-1]
+ }
+
+ #Corresponding probabilities
+ p[i,j,4] <- 1 - Q[i,j,3]
+
+ #Model observed browse severity data as categorical
+ browse[i,j] ~ dcat(p[i,j,1:4])
+
+ #Calculate expected value
+ ev[i,j] <- 1*p[i,j,1] + 2*p[i,j,2] + 3*p[i,j,3] + 4*p[i,j,4]
+
+ #Absolute residual for real datapoint
+ res[cucount[i,j]] <- abs(browse[i,j] - ev[i,j])
+
+ #Simulate new datapoint
+ browse.new[i,j] ~ dcat(p[i,j,1:4])
+
+ #Absolute residual for simulated datapoint
+ res.new[cucount[i,j]] <- abs(browse.new[i,j] - ev[i,j])
+
+ }}
+
+ #Derived quantities
+
+ #Difference between harvest and edge
+ diff.treat <- b.harvest - b.edge
+
+ #Posterior predictive check stuff
+ fit <- sum(res[])
+ fit.new <- sum(res.new[])
+
+ #Threshold priors
+ for(k in 1:3){
+ tau0[k] ~ dnorm(0,.01)
+ }
+ tau[1:3] <- sort(tau0[1:3])
+
+ #Priors
+
+ plot.tau <- pow(site.sd,-2)
+ plot.sd ~ dunif(0,100)
+
+ subplot.tau <- pow(plot.sd,-2)
+ subplot.sd ~ dunif(0,100)
+
+ seed.tau <- pow(seed.sd,-2)
+ seed.sd ~ dunif(0,100)
+
+ b.edge ~ dnorm(0,0.01)
+ b.harvest ~ dnorm(0,0.01)
+ b.species ~ dnorm(0,0.01)
+ b.ht ~ dnorm(0,0.01)
+ b.ht2 ~ dnorm(0,0.01)
+ b.comp ~ dnorm(0,0.01)
+ b.season ~ dnorm(0,0.01)
+ b.deer ~ dnorm(0,0.01)
+ b.insect ~ dnorm(0,0.01)
+ b.rabbit ~ dnorm(0,0.01)
+
+}
+
+ ``` \ No newline at end of file