aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorKen Kellner <ken.kellner@gmail.com>2016-09-29 20:24:18 (GMT)
committerKen Kellner <ken.kellner@gmail.com>2016-09-29 20:24:18 (GMT)
commit4fd101c69808fccb9920ceec29e19ea7eed6ab7b (patch)
tree18986a930c2bc17529903440a58ba913dce7249f
parent988b170b1899752d2c15e969fcce37c01a033d63 (diff)
Clean up repo, get rid of old code, add example data.
-rw-r--r--.gitignore2
-rw-r--r--LICENSE.md21
-rw-r--r--analysis_browse.R (renamed from analysis_deerbrowse.R)69
-rw-r--r--analysis_insectbrowse.R158
-rw-r--r--analysis_pellet_counts.R40
-rw-r--r--analysis_rabbitbrowse.R160
-rw-r--r--data/example_competition.csv11
-rw-r--r--data/example_pellet.csv11
-rw-r--r--data/example_seedling_master.csv11
-rw-r--r--figures/fig1_transect.R (renamed from figures/figure_transect.R)31
-rw-r--r--figures/fig2_pellets.R (renamed from figures/figure_pellets.R)48
-rw-r--r--figures/fig3_height.R (renamed from figures/figure_height.R)127
-rw-r--r--figures/table1_model_output.R (renamed from tables/table_modeloutput_cat.R)10
-rw-r--r--function_format_data.R (renamed from script_format_data.R)77
-rw-r--r--model_browse.R26
-rw-r--r--supplements/supplement_s2.Rmd112
-rw-r--r--tables/table_modeloutput.R89
17 files changed, 295 insertions, 708 deletions
diff --git a/.gitignore b/.gitignore
index d2198e5..b3a1013 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,7 +2,7 @@
.Rhistory
.RData
oak-herbivory.Rproj
-data*
+hee_*
output*
*tiff
*.bbl
diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000..a862114
--- /dev/null
+++ b/LICENSE.md
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2016 Kenneth F. Kellner
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/analysis_deerbrowse.R b/analysis_browse.R
index 8de62ef..aaeaeef 100644
--- a/analysis_deerbrowse.R
+++ b/analysis_browse.R
@@ -1,22 +1,44 @@
-########################################################
-## Analysis of deer browse intensity on oak seedlings ##
-## Ordinal logit model ##
-########################################################
+###################################################
+## Analysis of browse intensity on oak seedlings ##
+## Ordinal logit model ##
+###################################################
-source('script_format_data.R')
+source('function_format_data.R')
#Initial formatting on raw data
-seedling <- format.seedling('data/seedlingmaster.csv')
-exclude <- 1 - seedling$plot.data$herbivory
-
-#Plots to keep (deer not excluded, not in shelterwoods)
+seedling <- format.seedling('data/hee_seedling_master.csv')
plots <- c(1:54)
-keep.plots <- which(plots<49&exclude==0)
-keep.plots <- keep.plots[-2]
+###########################################################
+#Deer browse analysis
+exclude <- 1 - seedling$plot.data$herbivory
+keep.plots <- which(plots<49&exclude==0) #exclude shelterwoods & plots deer were excluded from
+keep.plots <- keep.plots[-2]
#Only keep seedlings that "established"
keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid%in%keep.plots)
+#Browse value (scale to 1-4 instead of 0-3)
+browse <- seedling$browsedeer[keep,]
+browse <- browse + 1
+
+############################################################
+#Rabbit browse analysis
+#Only keep seedlings that "established" & were not in shelterwoods
+keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid<49)
+
+browse <- seedling$browseother[keep,]
+
+############################################################
+#Insect browse analysis
+#Only keep seedlings that "established" & were not in shelterwoods
+keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid<49)
+
+browse <- seedling$lfdmg[keep,]
+
+
+############################################################
+#Common to all analyses
+
#Get survival data as baseline for when browsing could occur
surv <- seedling$surv.sprout[keep,]
nsamples <- numeric(dim(surv)[1])
@@ -60,10 +82,6 @@ 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
ht2 <- (ht2.raw-mean(ht2.raw,na.rm=TRUE))/sd(ht2.raw,na.rm=TRUE)
-#Browse value (scale to 1-4 instead of 0-3)
-browse <- seedling$browsedeer[keep,]
-browse <- browse + 1
-
#Previous browse
browse.deer.raw <- seedling$browsedeer[keep,]
browse.deer <- cbind(rep(0,dim(browse.deer.raw)[1]),browse.deer.raw[,1:7])
@@ -126,7 +144,6 @@ jags.data <- c('browse','nseedlings','nplots','nsites','nsamples'
,'ht','ht2'
,'species'
,'edge','harvest'
- #,'distance','distance2'
,'comp'
,'cucount'
,'season'
@@ -136,12 +153,14 @@ jags.data <- c('browse','nseedlings','nplots','nsites','nsamples'
################################
#Model file
+#Remove ht2 parameter for insects
modFile <- 'model_browse.R'
################################
#Parameters to save
+#Comment out ht2 parameter for insects
params <- c('tau','site.sd','plot.sd','seed.sd'
,'b.species'
@@ -149,7 +168,6 @@ params <- c('tau','site.sd','plot.sd','seed.sd'
,'b.season'
,'b.ht','b.ht2'
,'b.edge','b.harvest','diff.treat'
- #,'b.distance','b.distance2'
,'b.rabbit','b.deer','b.insect'
,'fit','fit.new'
)
@@ -166,10 +184,17 @@ inits <- function(){
require(jagsUI)
-#Switch to only deer
-
+#Deer analysis
browsedeer.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)
-browsedeer.output.cat <- update(browsedeer.output.cat,n.iter=4000,n.thin=40)
-
+ n.chains=3,n.iter=16000,n.burnin=12000,n.thin=40,parallel=TRUE)
save(browsedeer.output.cat,file="output/browsedeer_output_cat.Rdata")
+
+#Rabbit analysis
+browserabbit.output.cat <- jags(data=jags.data,inits=inits,parameters.to.save=params,model.file=modFile,
+ n.chains=3,n.iter=16000,n.burnin=12000,n.thin=40,parallel=TRUE)
+save(browserabbit.output.cat,file="output/browserabbit_output_cat.Rdata")
+
+#Insect analysis
+browseinsect.output.cat <- jags(data=jags.data,inits=inits,parameters.to.save=params,model.file=modFile,
+ n.chains=3,n.iter=16000,n.burnin=12000,n.thin=40,parallel=TRUE)
+save(browseinsect.output.cat,file="output/browseinsect_output_cat.Rdata")
diff --git a/analysis_insectbrowse.R b/analysis_insectbrowse.R
deleted file mode 100644
index 84945ec..0000000
--- a/analysis_insectbrowse.R
+++ /dev/null
@@ -1,158 +0,0 @@
-##########################################################
-## Analysis of insect browse intensity on oak seedlings ##
-## Ordinal logit model ##
-##########################################################
-
-source('script_format_data.R')
-
-#Initial formatting on raw data
-seedling <- format.seedling('data/seedlingmaster.csv')
-
-#Only keep seedlings that "established" & were not in shelterwoods
-keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid<49)
-
-#Get survival data as baseline for when browsing could occur
-surv <- seedling$surv.sprout[keep,]
-nsamples <- numeric(dim(surv)[1])
-for (i in 1:dim(surv)[1]){
- nsamples[i] <- length(na.omit(surv[i,]))
- if(0%in%surv[i,]){nsamples[i] <- nsamples[i]-1}
-}
-
-#Seedling-level covariates
-nseedlings <- dim(surv)[1]
-seedling.covs <- seedling$seedling.data[keep,]
-seed.sitecode <- seedling.covs$siteid
-seed.plotcode <- seedling.covs$plotid
-species <- seedling.covs$species
-age <- seedling.covs$age
-
-#Format height data
-ht.raw <- as.matrix(cbind(seedling$height[,1],seedling$height[,1],seedling$height[,2],seedling$height[,2],seedling$height[,3],
- seedling$height[,3],seedling$height[,4],seedling$height[,4]))
-ht.raw <- ht.raw[keep,]
-ht.raw[is.na(ht.raw[,1])&age==1,1] <- mean(ht.raw[age==1,1],na.rm=TRUE)
-ht.raw[is.na(ht.raw[,1])&age==0,1] <- mean(ht.raw[age==0,1],na.rm=TRUE)
-
-for (i in 1:dim(surv)[1]){
- for (j in 2:8){
- if(!is.na(surv[i,j])){
- if(is.na(ht.raw[i,(j-1)])){
- ht.raw[i,(j-1)] <- ht.raw[i,(j-2)]
- }}}}
-ht.raw[,8] <- ht.raw[,7]
-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
-
-ht2 <- (ht2.raw-mean(ht2.raw,na.rm=TRUE))/sd(ht2.raw,na.rm=TRUE)
-
-#Browse - simplify to presence/absence for now
-browse <- seedling$lfdmg[keep,]
-
-browse.deer.raw <- seedling$browsedeer[keep,]
-browse.deer <- cbind(rep(0,dim(browse.deer.raw)[1]),browse.deer.raw[,1:7])
-browse.deer[which(is.na(browse.deer),arr.ind=T)] <- 0
-browse.deer[which(browse.deer>1,arr.ind=T)] <- 1
-
-browse.rabbit.raw <- seedling$browseother[keep,]
-browse.rabbit <- cbind(rep(0,dim(browse.rabbit.raw)[1]),browse.rabbit.raw[,1:7])
-browse.rabbit[which(is.na(browse.rabbit),arr.ind=T)] <- 0
-browse.rabbit[which(browse.rabbit>1,arr.ind=T)] <- 1
-
-browse.insect.raw <- seedling$lfdmg[keep,]
-browse.insect <- cbind(rep(0,dim(browse.insect.raw)[1]),browse.insect.raw[,1:7])
-browse.insect[which(is.na(browse.insect),arr.ind=T)] <- 0
-browse.insect[which(browse.insect>1,arr.ind=T)] <- 1
-
-#Final response variable (all browse events)
-browse <- browse + 1
-
-#Format plot-level variables
-nplots <- 48
-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))
-#shelter <- c(rep(0,48),rep(1,6))
-plot.sitecode <- seedling$plot.data$siteid
-exclude <- 1 - seedling$plot.data$herbivory
-distance <- seedling$plot.data$distanceZ
-distance[4] <- 0
-distance2 <- seedling$plot.data$distance2Z
-distance2[4] <- 0
-
-#Competition
-comp <- seedling$comp.data[,1,]
-comp[which(is.na(comp),arr.ind=TRUE)] <- 0
-
-#Site level variables
-nsites <- 12
-
-#Cumulative sample count index
-cucount = matrix(NA, nseedlings, 8)
-index=1
-for(i in 1:nseedlings){
- for(j in 1:nsamples[i]){
- cucount[i,j] = index
- index = index+1
- }}
-
-#Season vector
-#summer = 1
-#initial,oct11,may12,oct12,may13,oct13,may14,oct14
-season <- c(1,1,0,1,0,1,0,1)
-
-###############################
-
-#Bundle data for JAGS
-
-jags.data <- c('browse','nseedlings','nplots','nsites','nsamples'
- ,'seed.plotcode','seed.sitecode','plot.sitecode'
- ,'ht'#,'ht2'
- ,'species'
- ,'edge','harvest'
- #,'distance','distance2'
- ,'comp'
- ,'cucount'
- ,'season'
- ,'browse.deer','browse.insect','browse.rabbit'
-)
-
-################################
-
-#Model file
-
-modFile <- 'model_browse.R'
-
-################################
-
-#Parameters to save
-
-params <- c('tau','site.sd','plot.sd','seed.sd'
- ,'b.species'
- ,'b.comp'
- ,'b.season'
- ,'b.ht','b.ht2'
- ,'b.edge','b.harvest','diff.treat'
- #,'b.distance','b.distance2'
- ,'b.rabbit','b.deer','b.insect'
- ,'fit','fit.new'
-)
-
-################################
-
-inits <- function(){
- list(
- "tau0" =c(-0.5,0,0.5)
- )
-}
-
-#Run analysis
-
-require(jagsUI)
-
-#Switch to only insects
-
-browseinsect.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)
-
-save(browseinsect.output.cat,file="output/browseinsect_output_cat.Rdata")
diff --git a/analysis_pellet_counts.R b/analysis_pellet_counts.R
new file mode 100644
index 0000000..4c21fc7
--- /dev/null
+++ b/analysis_pellet_counts.R
@@ -0,0 +1,40 @@
+###########################################################
+## Pellet Counts ~ Harvest Treatment: Randomization Test ##
+###########################################################
+
+#Read in raw data
+pellet.raw <- read.csv('data/hee_pellet.csv',header=T)
+pellet <- pellet.raw[pellet.raw$Pos!='sh',]
+
+#Actual test statistic
+actual <- aov(pellet$Deer~pellet$Treat)
+obtainedF <- summary(actual)[[1]]$'F value'[1]
+
+#Actual pairwise test stats
+pair.actual <- TukeyHSD(actual)
+obtained.HvE <- pair.actual$"pellet$Treat"[1,1]
+obtained.MvE <- pair.actual$"pellet$Treat"[2,1]
+obtained.MvH <- pair.actual$"pellet$Treat"[3,1]
+
+#Conduct randomization test
+nsims <- 1000
+
+dep <- pellet$Deer
+testF <- testHvE <- testMvE <- testMvH <- vector(length=nsims)
+
+for (i in 1:nsims){
+
+ indep1 <- sample(pellet$Treat[pellet$Year==2012])
+ indep2 <- sample(pellet$Treat[pellet$Year==2013])
+ indep3 <- sample(pellet$Treat[pellet$Year==2014])
+ indep <- unlist(list(indep1,indep2,indep3))
+
+ test <- aov(dep~indep)
+ testF[i] <- summary(test)[[1]]$'F value'[1]
+
+ pair.test <- TukeyHSD(test)
+ testHvE[i] <- pair.test$"indep"[1,1]
+ testMvE[i] <- pair.test$"indep"[2,1]
+ testMvH[i] <- pair.test$"indep"[3,1]
+
+} \ No newline at end of file
diff --git a/analysis_rabbitbrowse.R b/analysis_rabbitbrowse.R
deleted file mode 100644
index ed6f1e5..0000000
--- a/analysis_rabbitbrowse.R
+++ /dev/null
@@ -1,160 +0,0 @@
-##########################################################
-## Analysis of rabbit browse intensity on oak seedlings ##
-## Ordinal logit model ##
-##########################################################
-
-source('script_format_data.R')
-
-#Initial formatting on raw data
-seedling <- format.seedling('data/seedlingmaster.csv')
-
-#Only keep seedlings that "established" & were not in shelterwoods
-keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid<49)
-
-#Get survival data as baseline for when browsing could occur
-surv <- seedling$surv.sprout[keep,]
-nsamples <- numeric(dim(surv)[1])
-for (i in 1:dim(surv)[1]){
- nsamples[i] <- length(na.omit(surv[i,]))
- if(0%in%surv[i,]){nsamples[i] <- nsamples[i]-1}
-}
-
-#Seedling-level covariates
-nseedlings <- dim(surv)[1]
-seedling.covs <- seedling$seedling.data[keep,]
-seed.sitecode <- seedling.covs$siteid
-seed.plotcode <- seedling.covs$plotid
-species <- seedling.covs$species
-age <- seedling.covs$age
-
-#Format height data
-ht.raw <- as.matrix(cbind(seedling$height[,1],seedling$height[,1],seedling$height[,2],seedling$height[,2],seedling$height[,3],
- seedling$height[,3],seedling$height[,4],seedling$height[,4]))
-ht.raw <- ht.raw[keep,]
-ht.raw[is.na(ht.raw[,1])&age==1,1] <- mean(ht.raw[age==1,1],na.rm=TRUE)
-ht.raw[is.na(ht.raw[,1])&age==0,1] <- mean(ht.raw[age==0,1],na.rm=TRUE)
-
-for (i in 1:dim(surv)[1]){
- for (j in 2:8){
- if(!is.na(surv[i,j])){
- if(is.na(ht.raw[i,(j-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
-
-ht2 <- (ht2.raw-mean(ht2.raw,na.rm=TRUE))/sd(ht2.raw,na.rm=TRUE)
-
-#Browse - simplify to presence/absence for now
-browse <- seedling$browseother[keep,]
-
-
-browse.deer.raw <- seedling$browsedeer[keep,]
-browse.deer <- cbind(rep(0,dim(browse.deer.raw)[1]),browse.deer.raw[,1:7])
-browse.deer[which(is.na(browse.deer),arr.ind=T)] <- 0
-browse.deer[which(browse.deer>1,arr.ind=T)] <- 1
-
-browse.insect.raw <- seedling$lfdmg[keep,]
-browse.insect <- cbind(rep(0,dim(browse.insect.raw)[1]),browse.insect.raw[,1:7])
-browse.insect[which(is.na(browse.insect),arr.ind=T)] <- 0
-browse.insect[which(browse.insect>1,arr.ind=T)] <- 1
-
-browse.rabbit.raw <- seedling$browseother[keep,]
-browse.rabbit <- cbind(rep(0,dim(browse.rabbit.raw)[1]),browse.rabbit.raw[,1:7])
-browse.rabbit[which(is.na(browse.rabbit),arr.ind=T)] <- 0
-browse.rabbit[which(browse.rabbit>1,arr.ind=T)] <- 1
-
-#Final response variable (all browse events)
-browse <- browse + 1
-
-#Format plot-level variables
-nplots <- 48
-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))
-#shelter <- c(rep(0,48),rep(1,6))
-plot.sitecode <- seedling$plot.data$siteid
-exclude <- 1 - seedling$plot.data$herbivory
-distance <- seedling$plot.data$distanceZ
-distance[4] <- 0
-distance2 <- seedling$plot.data$distance2Z
-distance2[4] <- 0
-
-#Competition
-comp <- seedling$comp.data[,1,]
-comp[which(is.na(comp),arr.ind=TRUE)] <- 0
-
-#Site level variables
-nsites <- 12
-
-#Cumulative sample count index
-cucount = matrix(NA, nseedlings, 8)
-index=1
-for(i in 1:nseedlings){
- for(j in 1:nsamples[i]){
- cucount[i,j] = index
- index = index+1
- }}
-
-#Season vector
-#summer = 1
-#initial,oct11,may12,oct12,may13,oct13,may14,oct14
-season <- c(1,1,0,1,0,1,0,1)
-
-###############################
-
-#Bundle data for JAGS
-
-jags.data <- c('browse','nseedlings','nplots','nsites','nsamples'
- ,'seed.plotcode','seed.sitecode','plot.sitecode'
- ,'ht','ht2'
- ,'species'
- ,'edge','harvest'
- #,'distance','distance2'
- ,'comp'
- ,'cucount'
- ,'season'
- ,'browse.deer','browse.insect','browse.rabbit'
-)
-
-################################
-
-#Model file
-
-modFile <- 'model_browse.R'
-
-################################
-
-#Parameters to save
-
-params <- c('tau','site.sd','plot.sd','seed.sd'
- ,'b.species'
- ,'b.comp'
- ,'b.season'
- ,'b.ht','b.ht2'
- ,'b.edge','b.harvest','diff.treat'
- #,'b.distance','b.distance2'
- ,'b.rabbit','b.deer','b.insect'
- ,'fit','fit.new'
-)
-
-inits <- function(){
- list(
- "tau0" =c(-0.5,0,0.5)
- )
-}
-
-################################
-
-#Run analysis
-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/data/example_competition.csv b/data/example_competition.csv
new file mode 100644
index 0000000..f6c1661
--- /dev/null
+++ b/data/example_competition.csv
@@ -0,0 +1,11 @@
+"Date","Unit","Treatment","Site","Plotid","herb","woody","stems10.50","stems50.100","stems100","Dens1","Dens2","Dens3","Dens4","DensMean","Notes"
+"7/14/2011",3,"Clearcut","1HC",1,1,3,6,1,0,4,1,12,3,5,""
+"7/18/2011",6,"Clearcut","4HC",29,1,1,0,0,0,95,94,92,93,93.5,""
+"7/19/2011",6,"Clearcut","2HC",21,3,5,11,2,4,9,10,0,0,4.75,""
+"7/11/2011",3,"Clearcut","3HN",10,2,2,1,0,0,1,0,0,0,0.25,""
+"7/19/2011",6,"Clearcut","2NC",23,3,1,1,0,0,0,0,0,0,0,""
+"7/18/2011",6,"Clearcut","3NN",28,0,0,0,0,0,94,91,93,93,92.75,""
+"7/13/2011",3,"Clearcut","3NC",11,1,1,0,0,0,29,9,13,66,29.25,""
+"7/11/2011",3,"Clearcut","4HC",13,3,1,0,0,0,93,91,91,90,91.25,""
+"7/19/2011",6,"Clearcut","2HN",22,3,3,2,1,1,94,81,55,73,75.75,""
+"7/11/2011",3,"Clearcut","4NC",15,3,2,3,0,0,2,0,0,3,1.25,""
diff --git a/data/example_pellet.csv b/data/example_pellet.csv
new file mode 100644
index 0000000..081f371
--- /dev/null
+++ b/data/example_pellet.csv
@@ -0,0 +1,11 @@
+"Year","Date","Unit","Pos","Treat","Direction","T1","T2","T3","T4","T5","Deer","Rabbit","Notes"
+2012,"3/14/2012",3,"sh","Shelter","center","0","0","0","0","0",0,0,""
+2012,"3/13/2012",9,"sh","Shelter","up","0","0","0","0","1",0,1,""
+2013,"3/11/2013",3,"E","Edge","up","0","0","0","0","0",0,0,""
+2012,"3/14/2012",3,"C","Harvest","down","0","0","0","0","0",0,0,""
+2012,"3/13/2012",6,"E","Edge","up","1r","0","0","0","0",0,1,""
+2012,"3/13/2012",6,"M","Matrix","down","0","0","0","0","0",0,0,""
+2012,"3/13/2012",9,"E","Edge","up","0","0","0","0","0",0,0,""
+2012,"3/13/2012",9,"C","Harvest","up","0","0","0","0","0",0,0,""
+2013,"3/11/2013",3,"M","Matrix","up","0","0","0","0","0",0,0,""
+2012,"3/13/2012",9,"E","Edge","center","2D","0","0","0","2D",4,0,"out of forest - thicket"
diff --git a/data/example_seedling_master.csv b/data/example_seedling_master.csv
new file mode 100644
index 0000000..40fc89f
--- /dev/null
+++ b/data/example_seedling_master.csv
@@ -0,0 +1,11 @@
+"Unit","Treatment","Plot","Pos","Species","Age","Plant.date","Date","Surv","Ht","Rt","Lfdmg","Browse"
+3,"Clearcut","3HC","A3","W",0,"5/9/2011","7/11/2011",0,NA,NA,NA,""
+3,"Clearcut","2HN","A6","W",0,"5/9/2011","7/13/2011",0,NA,NA,NA,""
+3,"Clearcut","4HC","A1","W",0,"5/9/2011","7/13/2011",1,9.3,2.94,0,"0"
+3,"Clearcut","3HN","A8","B",0,"5/9/2011","7/11/2011",1,9,3.01,0,"0"
+3,"Clearcut","3NN","A2","B",0,"5/9/2011","7/13/2011",1,7,2.2,0,"0"
+3,"Clearcut","1NC","C3","B",0,"5/9/2011","7/14/2011",0,NA,NA,NA,""
+3,"Clearcut","3NN","F2","B",1,"5/23/2011","7/13/2011",1,50.5,4.52,0,"0"
+3,"Clearcut","3HC","D3","W",1,"5/23/2011","7/11/2011",1,42.1,6.75,0,"0"
+3,"Clearcut","1NC","A7","W",0,"5/9/2011","7/14/2011",1,13,3.33,0,"0"
+3,"Clearcut","3HN","C8","B",0,"5/9/2011","7/11/2011",1,5.9,2.56,0,"0"
diff --git a/figures/figure_transect.R b/figures/fig1_transect.R
index 22459ac..ba12aea 100644
--- a/figures/figure_transect.R
+++ b/figures/fig1_transect.R
@@ -1,12 +1,14 @@
-#####################################################
-## Figure: Proportion Browsed by Transect Position ##
-#####################################################
+#######################################################
+## Figure 1: Proportion Browsed by Transect Position ##
+#######################################################
-source('script_format_data.R')
+source('function_format_data.R')
#Initial formatting on raw data
-seedling <- format.seedling('data/seedlingmaster.csv')
+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
@@ -19,6 +21,7 @@ for (i in 1:length(ever.browsed.rabbit)){
}
}
+#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))
@@ -32,15 +35,16 @@ for(i in 1:length(ever.browsed.rabbit)){
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,]))){
@@ -49,10 +53,13 @@ for (i in 1:length(ever.browsed.insect)){
}
}
+#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)
@@ -80,6 +87,7 @@ for (i in 2:24){
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)
@@ -93,9 +101,12 @@ for(i in 1:length(ever.browsed.deer)){
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
tiff(filename="figures/Fig1.tiff",width=84,height=84,units="mm",res=1200, pointsize=8,
compression = "lzw",type='cairo')
@@ -106,6 +117,7 @@ require(extrafont)
myfont <- "Times"
#pdf(file="paper/figure_transect.pdf",width=5,height=3.9,family=myfont,pointsize=12)
+#Set up barplot
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",
@@ -114,18 +126,17 @@ barplot(cbind(rab,de,ins),beside=T, ylab="Proportion of Seedlings Browsed",
args.legend=list(x='topleft'),ylim=c(0,0.85))
box()
-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))
-
+#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){
diff --git a/figures/figure_pellets.R b/figures/fig2_pellets.R
index 486f6ae..82b4efd 100644
--- a/figures/figure_pellets.R
+++ b/figures/fig2_pellets.R
@@ -1,44 +1,12 @@
-################################################
-## Figure: Pellet Counts by Transect Position ##
-################################################
-
-pellet.raw <- read.csv('data/pellet.csv',header=T)
+##################################################
+## Figure 2: Pellet Counts by Transect Position ##
+##################################################
+#Read in raw data
+pellet.raw <- read.csv('data/hee_pellet.csv',header=T)
pellet <- pellet.raw[pellet.raw$Pos!='sh',]
-actual <- aov(pellet$Deer~pellet$Treat)
-obtainedF <- summary(actual)[[1]]$'F value'[1]
-
-pair.actual <- TukeyHSD(actual)
-obtained.HvE <- pair.actual$"pellet$Treat"[1,1]
-obtained.MvE <- pair.actual$"pellet$Treat"[2,1]
-obtained.MvH <- pair.actual$"pellet$Treat"[3,1]
-
-nsims <- 1000
-
-dep <- pellet$Deer
-
-testF <- testHvE <- testMvE <- testMvH <- vector(length=nsims)
-
-for (i in 1:nsims){
-
- indep1 <- sample(pellet$Treat[pellet$Year==2012])
- indep2 <- sample(pellet$Treat[pellet$Year==2013])
- indep3 <- sample(pellet$Treat[pellet$Year==2014])
- indep <- unlist(list(indep1,indep2,indep3))
-
- test <- aov(dep~indep)
- testF[i] <- summary(test)[[1]]$'F value'[1]
-
- pair.test <- TukeyHSD(test)
- testHvE[i] <- pair.test$"indep"[1,1]
- testMvE[i] <- pair.test$"indep"[2,1]
- testMvH[i] <- pair.test$"indep"[3,1]
-
-}
-
-##########################################################################
-
+#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
@@ -59,6 +27,7 @@ rabbit.se <- c(sd(pellet$Rabbit[pellet$Treat=="Matrix"])/sqrt(length(pellet$Rabb
plot.data <- cbind(rabbit.means,deer.means)
+#Set up plot container
tiff(filename="figures/Fig2.tiff",width=84,height=84,units="mm",res=1200, pointsize=8,
compression = "lzw",type='cairo')
@@ -69,6 +38,7 @@ require(extrafont)
myfont <- 'Times'
#pdf(file="paper/figure_pellets.pdf",width=5,height=3.9,family=myfont,pointsize=12)
+#Set up barplot
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'),
@@ -76,6 +46,7 @@ x <- barplot(plot.data,beside=TRUE,names=c('Rabbit','Deer'),
legend.text=c('Forest','Edge','Harvest'),
args.legend=list(x='topright'))
+#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])
@@ -86,6 +57,7 @@ for(i in 1:3){
box()
+#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'))
diff --git a/figures/figure_height.R b/figures/fig3_height.R
index fafa51c..b3c35a3 100644
--- a/figures/figure_height.R
+++ b/figures/fig3_height.R
@@ -1,56 +1,79 @@
-################################################
-## Figure: Height effects on Browse Intensity ##
-################################################
-
-
+##################################################
+## Figure 3: Height effects on Browse Intensity ##
+##################################################
+
+#Read in model ouptut
+load('output/browsedeer_output_cat.Rdata')
+load('output/browseinsect_output_cat.Rdata')
+load('output/browserabbit_output_cat.Rdata')
+
+#Raw seedling height data
+source('function_format_data.R')
+seedling <- format.seedling('data/hee_seedling_master.csv')
+
+seedling.covs <- seedling$seedling.data
+age <- seedling.covs$age
+ht.raw <- as.matrix(cbind(seedling$height[,1],seedling$height[,1],seedling$height[,2],
+ seedling$height[,2],seedling$height[,3],
+ seedling$height[,3],seedling$height[,4],seedling$height[,4]))
+ht.raw[is.na(ht.raw[,1])&age==1,1] <- mean(ht.raw[age==1,1],na.rm=TRUE)
+ht.raw[is.na(ht.raw[,1])&age==0,1] <- mean(ht.raw[age==0,1],na.rm=TRUE)
+
+#########################################
+#Seedling heights - deer exposed only
+
+exclude <- 1 - seedling$plot.data$herbivory
+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)
+ht.raw.de <- ht.raw[keep,]
+
+#Sort and normalize heights
sample.hts.de <- sort(as.numeric(na.omit(as.vector(ht.raw.de))))
-
sample.hts2.de <- (sample.hts.de-mean(sample.hts.de))^2
-
hts.Z.de <- as.numeric(scale(sample.hts.de))
-
hts2.Z.de <- as.numeric(scale(sample.hts2.de))
########################################################
+#Seedling heights - rabbit and insect exposed
+keep <- which(seedling$surv.sprout[,1]==1&seedling$seedling.data$plotid<49)
+ht.raw.rab <- ht.raw[keep,]
+#Sort and normalize heights
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
+#Set up figure container
tiff(filename="figures/Fig3.tiff",width=84,height=84,units="mm",res=1200, pointsize=8,
compression = "lzw",type='cairo')
-require(extrafont)
-#myfont <- "CM Roman"
-myfont <- 'Times'
-
-#pdf(file="paper/figure_height.pdf",width=3.9,height=3.9,family=myfont,pointsize=12)
-
par(mfrow=c(3,1),
oma = c(4,4,0,0) + 0.1,
mar = c(0.5,1,1,1) + 0.1)
-###
+require(extrafont)
+myfont <- 'Times'
+
+########################################################
+#Rabbit sub-figure
+
+#Get predicted browse probability by height and season
rp <- browserabbit.output.cat$mean
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)
+ + rp$b.ht2*hts2.Z.rab+rp$b.species) #winter
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)
+ + rp$b.ht2*hts2.Z.rab + rp$b.season +rp$b.species) #summer
-#Probability of being browsed at all by rabbits
-prab<- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
-prab.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
+#Probability of being browsed (any severity) by rabbit
+prab<- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
+prab.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
#Max value
-mx <-mxrab<- max(prab)
+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),
@@ -62,14 +85,17 @@ 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)
-###
+###################################################
+#Deer sub-figure
+
+#Get predicted probabilities by height and season
dp <- browsedeer.output.cat$mean
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
+#Probability of being browsed (any severity) by deer
pdeer<- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
pdeer.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
@@ -85,14 +111,17 @@ 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)
-###
+##############################################################
+#Insect sub-figure
+
+#Get predicted probabilities by height and season
ip <- browseinsect.output.cat$mean
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
+#Probability of being browsed (any severity) by insects
pins <- 1 - (exp(lin.pred) / (exp(lin.pred) + 1))
pins.sum <- 1 - (exp(lin.pred.sum) / (exp(lin.pred.sum) + 1))
@@ -105,40 +134,4 @@ 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)
-dev.off()
-Sys.setenv(R_GSCMD = "C:/Program Files/gs/gs9.19/bin/gswin64c.exe")
-embed_fonts("paper/figure_height.pdf")
-
-######################################
-
-#Combined figure
-tiff(filename="figures/Fig3.tiff",width=84,height=84,units="mm",res=1200, pointsize=8,
- compression = "lzw",type='cairo')
-require(extrafont)
-#font_install('fontcm')
-#loadfonts()
-myfont <- "CM Roman"
-#pdf(file="../dissertation/figures/fig3-2.pdf",width=3.9,height=3.9,family=myfont,pointsize=12)
-
-par(mgp=c(3,1,0),mar = c(5,4.5,2,2) + 0.1)
-plot(sample.hts,pdeer/max(pdeer),type='l',lwd=2,xlab="Seedling Height (cm)"
- ,ylab="Scaled Browse Probability",ylim=c(0,1))
-#text(20,mx-0.15*mx,"Insect",cex=2)
-segments(x0=mx.index.deer,y0=1,x1=mx.index.deer,y1=0.93)
-text(mx.index.deer,.9,round(mx.index.deer))
-
-lines(sample.hts,prab/max(prab),type='l',lwd=2,lty=2)
-#text(20,mx-0.15*mx,"Insect",cex=2)
-segments(x0=mx.index.rab,y0=1,x1=mx.index.rab,y1=0.87)
-text(mx.index.rab,.84,round(mx.index.rab))
-
-lines(sample.hts,pins/max(pins),type='l',lwd=2,lty=3)
-#text(20,mx-0.15*mx,"Insect",cex=2)
-segments(x0=mx.index.ins,y0=1,x1=mx.index.ins,y1=0.93)
-#text(mx.index.ins,.9,mx.index.ins)
-
-legend('topleft',lty=c(1,2,3),legend=c('Deer','Rabbit','Insect'),bty='n')
-
-dev.off()
-Sys.setenv(R_GSCMD = "C:/Program Files/gs/gs9.18/bin/gswin64c.exe")
-embed_fonts("../dissertation/figures/fig3-2.pdf")
+dev.off() \ No newline at end of file
diff --git a/tables/table_modeloutput_cat.R b/figures/table1_model_output.R
index ce95ec8..5e4947e 100644
--- a/tables/table_modeloutput_cat.R
+++ b/figures/table1_model_output.R
@@ -1,10 +1,6 @@
-########################################
-## Cumulative Link Model Output Table ##
-########################################
-
-########################################
-## Cumulative Link Model Output Table ##
-########################################
+###########################################
+## Table 1: Cumulative Link Model Output ##
+###########################################
#Decimal places to include
rd <- 2
diff --git a/script_format_data.R b/function_format_data.R
index 8a50cd8..fb88af0 100644
--- a/script_format_data.R
+++ b/function_format_data.R
@@ -1,22 +1,23 @@
-
+##########################################################
+## Function to format raw seedling growth/survival data ##
+##########################################################
format.seedling <- function(eh.file){
-#Read in file
+#Read in data files
raw <- read.csv(eh.file,header=TRUE,na.strings="")
-
-coords <- read.csv('data/plotcoords.csv',header=TRUE)
-
-comp <- read.csv('data/competition.csv',header=TRUE)
+coords <- read.csv('data/hee_plotcoords.csv',header=TRUE)
+comp <- read.csv('data/hee_competition.csv',header=TRUE)
#Tier 1: Site (not unit), 15 total (3*4+3)
#Tier 2: Plot (4 per clearcut site, 2 per shelterwood, 12*4+3*2 = 54 total)
#Tier 3: Individual seedlings x time
+########################################################################
#Site info
-
siteid <- 1:15
unit <- c(rep(3,4),rep(6,4),rep(9,4),3,6,9)
+
#treatments
exclosure <- c(rep(1:4,3),5,5,5)
opening <- c(rep(c(1,1,0,0),3),0,0,0)
@@ -24,10 +25,11 @@ edge <- c(rep(c(0,0,1,0),3),0,0,0)
matrix <- c(rep(c(0,0,0,1),3),0,0,0)
shelter <- c(rep(0,12),1,1,1)
+#Combine site data
site.data <- data.frame(siteid,unit,exclosure,opening,edge,matrix,shelter)
+########################################################################
#Plot info
-
unqplots <- as.character(sort(unique(raw$Plot)))
plotunit <- c(rep(3,16),rep(6,16),rep(9,16),3,3,6,6,9,9)
plot.siteid <- c(c(gl(4,4)),c(gl(4,4))+4,c(gl(4,4))+8,13,13,14,14,15,15)
@@ -41,6 +43,7 @@ distance2Z <- as.numeric(scale((coords[,4] - mean(coords[,4],na.rm=TRUE))^2))
#SW = 1 NE = 0
aspect <- as.numeric(coords[,5]=="SW")
canopy <- as.numeric(scale(comp$DensMean[1:54]))
+canopy2 <- comp$DensMean[1:54]/100
#Generate competition PCA
@@ -71,31 +74,37 @@ for (i in 1:4){
index2 = index2 + 2
}
-plot.data <- data.frame(unit=plotunit,siteid=plot.siteid,plotid,code,herbivory,competition,distanceZ,distance2Z,aspect,canopy)
+#Combine plot data
+plot.data <- data.frame(unit=plotunit,siteid=plot.siteid,plotid,code,herbivory,competition,distanceZ,distance2Z,aspect,canopy,canopy2)
+####################################################################################
#Seedling info
+#Get plot ID for each seedling
seed.plotid <- vector(length=dim(raw)[1])
compare <- paste(as.character(plot.data$unit),as.character(plot.data$code),sep="")
-
for (i in 1:dim(raw)[1]){
hold <- paste(as.character(raw$Unit[i]),as.character(raw$Plot[i]),sep="")
seed.plotid[i] <- which(hold==compare)
}
+#Get site ID for each seedling
seed.siteid <- vector(length=dim(raw)[1])
for (i in 1:length(seed.siteid)){
seed.siteid[i] <- plot.data$siteid[seed.plotid[i]]
}
+#Get seedling species and convert to binary
species <- rep(0,dim(raw)[1])
species[raw$Species=='W'] <- 1
+#Get mean/sd height values for 0-0 and 1-0 seedlings
mean0 <- suppressWarnings(mean(as.double(as.vector(raw$Ht[raw$Age==0])),na.rm=TRUE))
sd0 <- suppressWarnings(sd(as.double(as.vector(raw$Ht[raw$Age==0])),na.rm=TRUE))
mean1 <- suppressWarnings(mean(as.double(as.vector(raw$Ht[raw$Age==1])),na.rm=TRUE))
sd1 <- suppressWarnings(sd(as.double(as.vector(raw$Ht[raw$Age==1])),na.rm=TRUE))
+#Calculate initial height Z-scores
initialhtZ <- numeric(length(species))
for(i in 1:length(species)){
if(raw$Age[i]==0){
@@ -106,18 +115,21 @@ for(i in 1:length(species)){
}
initialhtZ[which(is.na(initialhtZ),arr.ind=TRUE)] <- 0
-
+#Combine seedling data
seedling.data <- data.frame(siteid=seed.siteid,plotid=seed.plotid,pos=raw$Pos,species,age=raw$Age,initialhtZ=initialhtZ,
initialht=raw$Ht)
-#Seedling survival matrix
+#############################################################
+#Observation-level data
+#Seedling survival matrix
+#Set up matrix
ind <- as.vector(sapply(names(raw),
function(i){
strsplit(i,split='[.]')[[1]][1]
}))
-
surv <- raw[,which(ind=='Surv')]
+
#Fix values in last survival column
surv[which(surv[,ncol(surv)]=='?'),ncol(surv)] <- 0
surv <- suppressWarnings(matrix(as.numeric(as.matrix(surv)),nrow=nrow(surv),ncol=ncol(surv)))
@@ -131,14 +143,11 @@ for (i in 1:nrow(surv)){
}
}
-#Identify seedlings that re-appeared
-
+#Identify seedlings that re-appeared after "dying" and correct
reappear <- rep(0,nrow(surv))
for (i in 1:nrow(surv)){
-
if(any(is.na(surv[i,]))){
-
-
+
st <- min(which(surv[i,]==0))
if(any(!is.na(surv[i,(st+1):ncol(surv)]))){
@@ -159,23 +168,21 @@ for (i in 1:nrow(surv)){
}
}
+#Extract notes for each observation
notes <- raw[,which(ind=='Notes')]
notes <- cbind(rep(NA,nrow(notes)),notes)
names(notes) <- c('Notes',paste('Notes.',1:(ncol(surv)-1),sep=""))
notes <- as.matrix(notes)
#Identify resprouts
-
sprout.keywords <- c('sprout','resprout','sprout?','sprout/alive','weakly sprouting',
'sprout; black?','sprouting','sprout?; dug')
sprout <- matrix(0,ncol=ncol(notes),nrow=nrow(notes))
sprout[which(notes%in%sprout.keywords,arr.ind=TRUE)] <- 1
-#which(rowSums(sprout)>0)
-
-#Survival including sprouting
+#Survival matrix including sprouting
surv.sprout <- surv
for (i in 1:nrow(surv)){
@@ -186,23 +193,19 @@ for (i in 1:nrow(surv)){
surv.sprout[i,1:ncol(surv)] <- 1
}
}
-
}
-#Identify problematic seedlings
+#Identify problematic seedlings with code below
#t3 <- which(reappear==1)
#t3 <- t3[!which(reappear==1)%in%which(rowSums(sprout)>0)]
#surv[t3,] # should be empty
-#Leaf damage
-
+#Get leaf damage values
lfdmg <- raw[,which(ind=='Lfdmg')]
lfdmg <- suppressWarnings(matrix(as.numeric(as.matrix(lfdmg)),nrow=nrow(lfdmg),ncol=ncol(lfdmg)))
-
-#Sample dates
+#Get sample dates
dateraw <- data.frame(raw[,which(ind=='Date')])
-
sample.dates <- data.frame()
plant.date <- vector(length=length(siteid))
@@ -231,32 +234,27 @@ for (i in 1:ncol(sample.dates)){
}
-#Height
-
+#Get seedling height
heightraw <- data.frame(raw[,which(ind=='Ht')])
height <- suppressWarnings(matrix(as.numeric(as.matrix(heightraw)),nrow=nrow(heightraw),ncol=ncol(heightraw)))
-#Height growth
-
+#Get seedling height growth
htgrowth <- matrix(NA,nrow=nrow(height),ncol=(ncol(height)-1))
for (i in 1:ncol(htgrowth)){
htgrowth[,i] <- height[,i+1] - height[,i]
}
-#Root collar Diameter
-
+#Get seedling root collar Diameter (rcd)
rcdraw <- data.frame(raw[,which(ind=='Rt')])
rcd <- suppressWarnings(matrix(as.numeric(as.matrix(rcdraw)),nrow=nrow(rcdraw),ncol=ncol(rcdraw)))
-#Root collar diameter growth
-
+#Get seedling root collar diameter growth
rcdgrowth <- matrix(NA,nrow=nrow(rcd),ncol=(ncol(rcd)-1))
for (i in 1:ncol(rcdgrowth)){
rcdgrowth[,i] <- rcd[,i+1] - rcd[,i]
}
-#Browse damage
-
+#Browse damage from mammals
browseraw <- data.frame(raw[,which(ind=='Browse')])
browsedeer <- matrix(NA,nrow=nrow(browseraw),ncol=ncol(browseraw))
@@ -276,6 +274,7 @@ browse[which((browseraw=='1r'|browseraw=='1R'|browseraw=='1d'|browseraw=='1D'),a
browse[which((browseraw=='2r'|browseraw=='2R'|browseraw=='2d'|browseraw=='2D'),arr.ind=TRUE)] <- 2
browse[which((browseraw=='3r'|browseraw=='3R'|browseraw=='3d'|browseraw=='3D'),arr.ind=TRUE)] <- 3
+#Combine all values in output
output <- list(site.data=site.data,plot.data=plot.data,seedling.data=seedling.data,
comp.data=comp.data,
surv=surv,surv.sprout=surv.sprout,
diff --git a/model_browse.R b/model_browse.R
index 326ff00..53e2ad0 100644
--- a/model_browse.R
+++ b/model_browse.R
@@ -1,33 +1,44 @@
+#####################################
+## Browse damage on oak seedlings ##
+## Hierarchical Ordinal Regression ##
+#####################################
+
+## @knitr model.browse
model {
#Likelihood
+ #Site random effect
for (i in 1:nsites){
site.mean[i] ~ dnorm(0,site.tau)
}
+ #Plot random effect
for (i in 1:nplots){
plot.mean[i] ~ dnorm(0, plot.tau)
}
+ #Iterate over seedlings
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]]
+ site.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.ht*ht[i,j]
+ + b.ht2*ht2[i,j] #Comment out for insect analysis
+ b.season*season[j]
+ b.deer*browse.deer[i,j]
+ b.insect*browse.insect[i,j]
+ b.rabbit*browse.rabbit[i,j]
- #+ b.distance*distance[seed.plotcode[i]]
- #+ b.distance2*distance2[seed.plotcode[i]]
+ b.edge*edge[seed.plotcode[i]]
+ b.harvest*harvest[seed.plotcode[i]]
@@ -64,29 +75,22 @@ model {
diff.treat <- b.harvest - b.edge
-
#Derived quantities
fit <- sum(res[])
fit.new <- sum(res.new[])
#Priors
-
grand.mean ~ dunif(-100,100)
-
site.tau <- pow(site.sd,-2)
site.sd ~ dunif(0,100)
-
plot.tau <- pow(plot.sd,-2)
plot.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.distance ~ dnorm(0,0.01)
- #b.distance2 ~ dnorm(0,0.01)
b.species ~ dnorm(0,0.01)
b.ht ~ dnorm(0,0.01)
b.ht2 ~ dnorm(0,0.01)
diff --git a/supplements/supplement_s2.Rmd b/supplements/supplement_s2.Rmd
index 9275607..2c01204 100644
--- a/supplements/supplement_s2.Rmd
+++ b/supplements/supplement_s2.Rmd
@@ -10,111 +10,11 @@ Kenneth F. Kellner and Robert K. Swihart
Department of Forestry and Natural Resources, Purdue University
email: kkellner@purdue.edu
+```{r, echo=FALSE}
+library(knitr)
+read_chunk('../model_browse.R')
+```
```{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
+<<model.browse>>
+``` \ No newline at end of file
diff --git a/tables/table_modeloutput.R b/tables/table_modeloutput.R
deleted file mode 100644
index f8e8bfb..0000000
--- a/tables/table_modeloutput.R
+++ /dev/null
@@ -1,89 +0,0 @@
-########################################
-## Cumulative Link Model Output Table ##
-########################################
-
-#Decimal places to include
-rd <- 2
-
-#Format deer model output
-td <- browsedeer.output
-include <- c(1,2,3,7,10,15,16,18,17,11,12,8,9)
-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
-include <- c(1,2,3,7,9,14,15,17,16,10,11,8)
-mn2 <-c(unlist(tr$mean)[include],NA)
-q2.5 <- round(unlist(tr$q2.5),rd)[include]
-q97.5 <- round(unlist(tr$q97.5),rd)[include]
-ci2 <- c(paste("(",sprintf("%.2f",q2.5),",",sprintf("%.2f",q97.5),')',sep=""),'-')
-f2 <- c(sprintf("%.2f",round(unlist(tr$f),rd))[include],'-')
-ov2 <- c(unlist(tr$overlap0)[include],TRUE)
-
-#Format insect model output (2 missing values)
-ti <- browseinsect.output
-include <- c(1,2,3,7,9,13,14,16,15,10,8)
-mn3 <-c(unlist(ti$mean)[include],NA)
-mn3 <- c(mn3[1:10],NA,mn3[11:12])
-q2.5 <- round(unlist(ti$q2.5),rd)[include]
-q97.5 <- round(unlist(ti$q97.5),rd)[include]
-ci3 <- c(paste("(",sprintf("%.2f",q2.5),",",sprintf("%.2f",q97.5),')',sep=""),'-')
-ci3 <- c(ci3[1:10],'-',ci3[11:12])
-f3 <- c(sprintf("%.2f",round(unlist(ti$f),rd))[include],'-')
-f3 <- c(f3[1:10],'-',f3[11:12])
-ov3 <- c(unlist(ti$overlap0)[include],TRUE)
-ov3 <- c(ov3[1:10],TRUE,ov3[11:12])
-
-#Combine into raw data frame
-ov <- cbind(ov1,ov1,ov1,ov2,ov2,ov2,ov3,ov3,ov3)
-tab <- data.frame(mn1,ci1,f1,mn2,ci2,f2,mn3,ci3,f3)
-row.names(tab) <- c('Plot SD','Subplot SD','Seedling SD',
- 'Competition','White Oak','Growing Season','Prior Deer','Prior Rabbit','Prior Insect',
- 'Height','Height\\textsuperscript{2}','Distance','Distance\\textsuperscript{2}')
-colnames(tab) <- rep(c('Mean','95% CI','\\textit{f}'),3)
-
-#Correct spacing of mean values so decimal point locations match in LaTeX
-hold <- tab
-for (i in 1:length(ov[,1])){
- for (j in c(1,4,7)){
-
- if(is.na(hold[i,j])){tab[i,j] = '~~~-'; next}
-
- val <- sprintf("%.2f",hold[i,j])
- nval <- sprintf("%.2f",-hold[i,j])
-
- if((i<4)){tab[i,j] <- paste('~',val,sep=""); next}
-
- if(!ov[i,j]){
- if(hold[i,j]<0){tab[i,j] <- paste('-$',nval,'^{*}$',sep="")
- } else {tab[i,j] <- paste('~$',val,'^{*}$',sep="")}
-
- } else {
- if(hold[i,j]<0){tab[i,j] <- paste(val,sep="")
- } else {tab[i,j] <- paste('~',val,sep="")}
- }
- }
-
-}
-
-#Write LaTeX table code to output file
-#Manually comment out cline
-require(Hmisc)
-capture.output(
- latex(tab,title='Parameter',file="",
- rgroup = c("Random", "Fixed"), n.rgroup=c(3, nrow(tab)-3),
- cgroup = c('Deer','Rabbit','Insect'), n.cgroup = c(3,3,3),
- booktabs=TRUE,
- col.just = strsplit("lcclcclcc", "")[[1]],
- 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.tex")
-