Ants and Ecosystem Functions in HF-HeRE
JAK and AME
This code handles raw data for ants, nitrogen, decomposition and soil respiration from ant subplots in the Harvard Forest Hemlock Removal Experiment. The data are cleaned, transformed and analyzed, and final figures for publication are produced.
The following packages must be installed for the code to work: plyr ggplot2 reshape vegan gdata plspm segmented
The code takes the following .csv data files as inputs: 2014resinmaster.csv compare-all-saplings.csv Hemlock Ant Decomp Master_with formulas.csv hf106-03-shrub.csv hf107-01-light.csv hf160-03-ants-updated-for-2014.csv Simes Ants Exclosures Carbon flux data 160-05 through 2008.csv simeslightuse-annotated.csv understory-veg2014.csv veg2014.csv veg2014-predict.csv These files should all be in the same directory, which should be set to the R working directory prior to running the code.
Chunk 1: Load packages, read ant, nitrogen, decomp and carbon data, clean it and put it in generally consistent format.
require(plyr)## Loading required package: plyrrequire(ggplot2)## Loading required package: ggplot2require(grid)## Loading required package: gridrequire(reshape)## Loading required package: reshape
##
## Attaching package: 'reshape'
##
## The following objects are masked from 'package:plyr':
##
## rename, round_anyrequire(vegan)## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.0-10require(gdata)## Loading required package: gdata
## gdata: Unable to locate valid perl interpreter
## gdata:
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata:
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.
##
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.
##
## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.
##
## Attaching package: 'gdata'
##
## The following object is masked from 'package:stats':
##
## nobs
##
## The following object is masked from 'package:utils':
##
## object.sizerequire(plspm)## Loading required package: plspm
## Loading required package: amap
## Loading required package: diagram
## Loading required package: shape
## Loading required package: tester
## Loading required package: turnerrequire(segmented)## Loading required package: segmented## Warning: package 'segmented' was built under R version 3.1.2###############################################################################
#Ants
###############################################################################
ant.master<-read.table("hf160-03-ants-updated-for-2014.csv",sep=",",header=T)
#Fill in blank values with 1, pool cases where there are multiple lines
#for the same species in the same trap.
ant.master[is.na(ant.master)]<-1
ant.master<-ddply(ant.master,
.(Date,Treatment,Block,Plot,Subplot,SubplotT,Subfamily,Genus,Species,SpecCode),
summarise,Count=sum(Count))
#Classify date and add year column
ant.master$Date<-as.character(ant.master$Date)
ant.master$Date<-as.Date(ant.master$Date,"%Y%m%d")
ant.master$Year<-as.numeric(format(ant.master$Date,"%Y"))
ant.master$SubplotT <- revalue(ant.master$SubplotT,
c("control" = "Control", "disturbance" = "Disturbance Control",
"exclosure" = "Ants Excluded"))
###############################################################################
#Nitrogen
###############################################################################
resin.master<-read.table("2014resinmaster.csv",header=T,sep=",")
resin.master<-resin.master[c("Month.extracted..end.of.quarter.","year","quarter",
"block","canopy","plot","group","treatment",
"ug.NH4.gresin.day","ug.NO3.gresin.day")]
colnames(resin.master)<-c("month","year","quarter","block","canopy","plot",
"group","treatment","NH4","NO3")
#Convert month extracted to date of measurement
resin.master$month<-paste(resin.master$month,"-01",sep="")
resin.master$month<-as.Date(resin.master$month,"%Y-%m-%d")
#Average A/B subplots and split by ion
NO3master<-ddply(resin.master,.(month,year,quarter,block,canopy,plot,treatment),
summarise,NO3=mean(NO3,na.rm=TRUE))
NH4master<-ddply(resin.master,.(month,year,quarter,block,canopy,plot,treatment),
summarise,NH4=mean(NH4,na.rm=TRUE))
###############################################################################
#Carbon
###############################################################################
carbon.master<-read.table("Simes Ants Exclosures Carbon flux data 160-05 through 2008.csv",
sep=",",header=T)
#generate "run" numbers by clumping all observations are taken in <5 days.
#Incomplete runs are removed.
carbon.master$Run<-carbon.master$Date
levels(carbon.master$Run)<-c(0,1,1,0,2,2,3,3,0,0,0,0,0,6,6,6,8,8,7,9,9,7,10,
11,11,10,13,13,12,12,14,14,16,16,15,15)
carbon.master<-carbon.master[-which(carbon.master$Run==0),]
###############################################################################
#Decomp
###############################################################################
decomp.master<-read.table("Hemlock Ant Decomp Master_with formulas.csv",
sep=",",header=TRUE,comment.char="")
#Find and delete all observations from CWT,
#all observations from addition subplots,
#and the handful of observations with no plot ID column.
todelete=which(decomp.master$Site=="CWT")
decomp.master<-decomp.master[-todelete,]
todelete=which(decomp.master$Ants=="addition")
decomp.master<-decomp.master[-todelete,]
todelete=which(is.na(decomp.master$Plot..))
decomp.master<-decomp.master[-todelete,]
#get rid of excess columns, clean up names, change some ID values for ease of use
row.names(decomp.master)<-NULL
decomp.master<-decomp.master[-c(1,7,10:23)]
colnames(decomp.master)<-c("Date","Months","Plot","Subplot",
"SubplotT","Litter","Treatment","Remaining")
decomp.master$Block=ifelse(decomp.master$Plot>3&decomp.master$Plot<8,"Ridge","Valley")
decomp.master$Months<-as.numeric(as.character(decomp.master$Months))
decomp.master$Treatment<-factor(decomp.master$Treatment)
decomp.master$Treatment<-factor(decomp.master$Treatment,
levels(decomp.master$Treatment)[c(2,3,4,1)])
decomp.master$Date<-as.Date(decomp.master$Date,format="%m.%d.%Y")
#Split file into lignin and cellulose files
lig<-which(decomp.master$Litter=="Lignin")
lignin.master<-decomp.master[lig,]
cell.master<-decomp.master[-lig,]
row.names(lignin.master)<-NULLChunk 2: Transforation and derived values
###############################################################################
#Ants
###############################################################################
#Generate dataframe with richness and abundance for each observation
ant.observations<-ddply(ant.master,
.(Year,Block,Plot,Treatment,SubplotT,Subfamily,Genus,Species,SpecCode),
summarise,Count=sum(Count))
ant.richness.count<-ddply(ant.observations,.(Year,Plot,SubplotT),summarise,
total=sum(Count),richness=length(SpecCode), .drop=FALSE)
ant.richness.count$Block<-ifelse(ant.richness.count$Plot>3&ant.richness.count$Plot<8,"Ridge","Valley")
ant.richness.count$Treatment<-ifelse(ant.richness.count$Plot==1|ant.richness.count$Plot==5,"G",
ifelse(ant.richness.count$Plot==2|ant.richness.count$Plot==4,
"L",ifelse(ant.richness.count$Plot>6,"HD","HC")))
ant.richness.count$Treatment<-factor(ant.richness.count$Treatment,levels=c("HC","G","L","HD"))
#add pre/post/overlap
ant.richness.count$invasion <- ifelse(ant.richness.count$Year < 2010, "pre",
ifelse(ant.richness.count$Year == 2010, "during", "post"))
#Calculate abundance of each genus in each observation
ant.genus<-ddply(ant.observations,.(Year,Treatment,SubplotT,Genus),
summarise,abundance=sum(Count), .drop=F)
ant.genus.totals<-ddply(ant.observations,.(Year,Treatment,SubplotT),summarise,total=sum(Count))
ant.genus$total<-rep(ant.genus.totals$total,each=9)
#Put data in wide form. Block and Treatment excluded to allow add.missing to work.
ant.master.wide<-cast(ant.master,Year+Plot+SubplotT~SpecCode,value="Count",
fun.aggregate=sum,add.missing=T)
ant.master.wide$Year<-as.numeric(as.character(ant.master.wide$Year))
ant.master.wide$Plot<-as.numeric(as.character(ant.master.wide$Plot))
Block<-ifelse(ant.master.wide$Plot>3&ant.master.wide$Plot<8,"Ridge","Valley")
Treatment<-ifelse(ant.master.wide$Plot==1|ant.master.wide$Plot==5,"G",
ifelse(ant.master.wide$Plot==2|ant.master.wide$Plot==4,"L",
ifelse(ant.master.wide$Plot>6,"HD","HC")))
ant.master.wide<-cbind(Block,Treatment,ant.master.wide)
#Using this, generate hill numbers
raw.div<-ant.master.wide[,!(names(ant.master.wide)%in%c("Block","Treatment","Year","Plot","SubplotT"))]
H<-diversity(raw.div,index="shannon")
ant.div<-cbind(ant.master.wide[c("Block","Treatment","Year","Plot","SubplotT")],H)
ant.div$eH<-exp(ant.div$H)
ant.div$Treatment<-factor(ant.div$Treatment,levels=c("HC","G","L","HD"))
ant.div$invasion <- ifelse(ant.div$Year < 2010, "pre", ifelse(ant.div$Year == 2010, "during", "post"))
###############################################################################
#Nitrogen
###############################################################################
#Put data in wide form with the three subplot treatments as individual columns
NO3master<-cast(NO3master,month+year+quarter+block+canopy+plot~treatment,value="NO3")
NH4master<-cast(NH4master,month+year+quarter+block+canopy+plot~treatment,value="NH4")
#generate new columns for exclosure-control and disturbance-control
NO3master$D.C<-NO3master$D-NO3master$C
NO3master$X.C<-NO3master$X-NO3master$C
NH4master$D.C<-NH4master$D-NH4master$C
NH4master$X.C<-NH4master$X-NH4master$C
#Put data back in long form with a line for each value
NO3master<-as.data.frame(NO3master)
NH4master<-as.data.frame(NH4master)
NO3master<-melt(NO3master,id.vars=c("month","year","quarter","block","canopy","plot"),measure.vars=c("D.C","X.C"))
NH4master<-melt(NH4master,id.vars=c("month","year","quarter","block","canopy","plot"),measure.vars=c("D.C","X.C"))
#split data mid-2010 for pre/post invasion
NH4master$invasion<-as.factor(ifelse(NH4master$month<as.Date("2010-04-02"),"Pre-adelgid","Post-adelgid"))
NO3master$invasion<-as.factor(ifelse(NO3master$month<as.Date("2010-04-02"),"Pre-adelgid","Post-adelgid"))
#refactor for plotting- only need to do one because only need one legend
NH4master$variable <- revalue(NH4master$variable, c("D.C" = "Disturbance Control", "X.C" = "Ants Excluded"))
###############################################################################
#Carbon
###############################################################################
carbon.rel<-ddply(carbon.master,.(Run,Block,Plot,Treatment,SubplotT),summarise,co2=mean(co2.flux),DOY=DOY[1],Date=Date[1])
carbon.rel$Block<-ifelse(carbon.rel$Block==1,"Ridge","Valley")
#Transform to difference from control
carbon.rel<-cast(carbon.rel,Run+Block+Plot+Treatment+DOY+Date~SubplotT,value="co2")
carbon.rel$D.C<-carbon.rel$disturbed-carbon.rel$control
carbon.rel$X.C<-carbon.rel$exclosure-carbon.rel$control
carbon.rel<-as.data.frame(carbon.rel)
carbon.rel<-melt(carbon.rel, id.vars=c("Run","Block","Plot","Treatment","DOY","Date"),measure.vars=c("D.C","X.C"))
carbon.rel$Date<-(as.Date(carbon.rel$Date))
###############################################################################
#Decomp
###############################################################################
#further cleanup
lignin.master<-ddply(lignin.master,.(Date,Months,Plot,Subplot,SubplotT,Treatment,Block),summarise,Remaining=mean(Remaining,na.rm=TRUE))
cell.master<-cell.master[-6]
#Add a series of rows indicating 100% remaining at 9-22-2008
chunk<-cell.master[which(cell.master$Months==3),]
chunk$Months=0
chunk$Date<-as.Date("09.22.2008",format="%m.%d.%Y")
chunk$Remaining=100
ligmaster.extended<-rbind(chunk,lignin.master)
cellmaster.extended<-rbind(chunk,cell.master)
#Add grouping variable, run regression
#input includes 100% initial values for cellulose but not lignin
cellmaster.extended$ptreat<-paste(cellmaster.extended$Plot,cellmaster.extended$SubplotT)
lignin.master$ptreat<-paste(lignin.master$Plot,lignin.master$SubplotT)
cell.fitlist<-dlply(cellmaster.extended,.(ptreat),lm,formula=Remaining~Months)
lig.fitlist<-dlply(lignin.master,.(ptreat),lm,formula=Remaining~Months)
#Create a function to extract slope and SE of slope from these lists.
extractfun <- function(m) {
cf <- coef(m)
tinfo <- summary(m)$coefficients[2,2]
data.frame(m = cf[2], slope.se = tinfo)
}
#Generate dataframes with m and b (slope and intercetp) for each subplot type
cell.coefs<-ldply(cell.fitlist,extractfun)
lig.coefs<-ldply(lig.fitlist,extractfun)
Treatment<-rep(c("Girdled","Logged","Hemlock Control","Logged","Girdled",
"Hemlock Control","Hardwood Control","Hardwood Control"),each=3)
SubplotT<-rep(c("C","D","X"),8)
Block<-rep(c("v","v","v","r","r","r","r","v"),each=3)
cell.coefs<-cbind(Treatment,SubplotT,Block,cell.coefs)
lig.coefs<-cbind(Treatment,SubplotT,Block,lig.coefs)Chunk 3: ANOVAs for ants, nitrogen, decomp and carbon
###############################################################################
#Ants
###############################################################################
summary(aov(data=ant.richness.count,total~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 60.17 60.17
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 2408.3 802.8 3.509 0.165
## Residuals 3 686.3 228.8
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 1484.1 742.0 6.26 0.0231 *
## Treatment:SubplotT 6 604.8 100.8 0.85 0.5662
## Residuals 8 948.3 118.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 3575 3575 44.53 2.64e-10 ***
## Residuals 191 15336 80
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=ant.richness.count,richness~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 1.5 1.5
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 161.81 53.94 11.25 0.0386 *
## Residuals 3 14.39 4.80
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 33.95 16.977 4.461 0.050 *
## Treatment:SubplotT 6 18.05 3.008 0.790 0.602
## Residuals 8 30.44 3.806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 71.11 71.11 42.98 5.02e-10 ***
## Residuals 191 316.00 1.65
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=ant.div,eH~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 1.404 1.404
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 75.64 25.212 15.78 0.0242 *
## Residuals 3 4.79 1.598
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 14.25 7.124 2.920 0.112
## Treatment:SubplotT 6 11.52 1.920 0.787 0.604
## Residuals 8 19.52 2.440
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 27.25 27.247 33.02 3.56e-08 ***
## Residuals 191 157.62 0.825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#pre/post adelgid, using != to include "during" in both categories
summary(aov(data=ant.richness.count[which(ant.richness.count$invasion != "post"),],
total~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 14.01 14.01
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 791.4 263.81 4.462 0.125
## Residuals 3 177.4 59.12
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 574.5 287.26 7.103 0.0168 *
## Treatment:SubplotT 6 357.1 59.52 1.472 0.2987
## Residuals 8 323.5 40.44
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 1586 1586.2 47.85 5.33e-10 ***
## Residuals 95 3149 33.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=ant.richness.count[which(ant.richness.count$invasion != "post"),],
richness~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 0.2083 0.2083
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 55.69 18.564 3.74 0.154
## Residuals 3 14.89 4.964
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 16.62 8.308 6.647 0.0199 *
## Treatment:SubplotT 6 12.18 2.031 1.624 0.2566
## Residuals 8 10.00 1.250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 51.34 51.34 40.28 7.45e-09 ***
## Residuals 95 121.06 1.27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=ant.div[which(ant.div$invasion != "post"),],
eH~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 0.194 0.194
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 25.989 8.663 4.279 0.132
## Residuals 3 6.074 2.025
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 8.480 4.240 3.773 0.0701 .
## Treatment:SubplotT 6 7.674 1.279 1.138 0.4204
## Residuals 8 8.991 1.124
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 24.47 24.47 34.96 5.27e-08 ***
## Residuals 95 66.50 0.70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=ant.richness.count[which(ant.richness.count$invasion != "pre"),],
total~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 11.41 11.41
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 2628 876.1 2.036 0.287
## Residuals 3 1291 430.4
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 939.1 469.6 2.214 0.172
## Treatment:SubplotT 6 776.9 129.5 0.611 0.718
## Residuals 8 1696.4 212.0
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 96 96.27 0.907 0.343
## Residuals 95 10088 106.19summary(aov(data=ant.richness.count[which(ant.richness.count$invasion != "pre"),],
richness~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 3.008 3.008
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 153.96 51.32 21.71 0.0155 *
## Residuals 3 7.09 2.36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 22.05 11.025 1.986 0.199
## Treatment:SubplotT 6 17.82 2.969 0.535 0.769
## Residuals 8 44.40 5.550
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 1.84 1.837 1.249 0.267
## Residuals 95 139.76 1.471summary(aov(data=ant.div[which(ant.div$invasion != "pre"),],
eH~Treatment*SubplotT+Error(Block/Treatment/SubplotT)+Year))##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 1.973 1.973
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 73.90 24.63 41.78 0.00603 **
## Residuals 3 1.77 0.59
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 8.874 4.437 1.737 0.236
## Treatment:SubplotT 6 11.847 1.974 0.773 0.613
## Residuals 8 20.437 2.555
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 1 0.79 0.7876 1.14 0.288
## Residuals 95 65.64 0.6909###############################################################################
#Nitrogen
###############################################################################
#Models on Full data
summary(aov(data=NO3master,value~canopy*variable+Error(block/canopy/variable)+month))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 0.6232 0.6232
##
## Error: block:canopy
## Df Sum Sq Mean Sq F value Pr(>F)
## canopy 3 41.36 13.79 0.32 0.813
## Residuals 3 129.20 43.07
##
## Error: block:canopy:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 4.371 4.371 13.370 0.0216 *
## canopy:variable 3 2.915 0.972 2.973 0.1600
## Residuals 4 1.308 0.327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 55 55.41 4.904 0.0273 *
## Residuals 447 5051 11.30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=NH4master,value~canopy*variable+Error(block/canopy/variable)+month))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 27.13 27.13
##
## Error: block:canopy
## Df Sum Sq Mean Sq F value Pr(>F)
## canopy 3 219.32 73.11 3.221 0.181
## Residuals 3 68.09 22.70
##
## Error: block:canopy:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 18.53 18.53 1.062 0.361
## canopy:variable 3 144.73 48.24 2.765 0.175
## Residuals 4 69.79 17.45
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 56 55.95 2.039 0.154
## Residuals 447 12265 27.44#ANOVAs of only Pre-adelgid-infestation data
summary(aov(data=NO3master[which(NO3master$invasion=="Pre-adelgid"),],
value~canopy*variable+Error(block/canopy/variable)+month))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 0.3914 0.3914
##
## Error: block:canopy
## Df Sum Sq Mean Sq F value Pr(>F)
## canopy 3 0.0521 0.0174 0.017 0.996
## Residuals 3 3.1481 1.0494
##
## Error: block:canopy:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 0.5723 0.5723 6.479 0.0636 .
## canopy:variable 3 1.6151 0.5384 6.095 0.0567 .
## Residuals 4 0.3533 0.0883
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 3.63 3.632 5.345 0.0217 *
## Residuals 223 151.53 0.679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(data=NH4master[which(NH4master$invasion=="Pre-adelgid"),],
value~canopy*variable+Error(block/canopy/variable)+month))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 5.829 5.829
##
## Error: block:canopy
## Df Sum Sq Mean Sq F value Pr(>F)
## canopy 3 46.53 15.51 0.312 0.818
## Residuals 3 149.17 49.72
##
## Error: block:canopy:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 4.43 4.43 0.191 0.685
## canopy:variable 3 105.32 35.11 1.511 0.341
## Residuals 4 92.96 23.24
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 10 9.964 0.634 0.427
## Residuals 223 3506 15.723#ANOVAs of only Post-adelgid-infestation data
#AME added this 10/21/14
summary(aov(data=NO3master[which(NO3master$invasion=="Post-adelgid"),],
value~canopy*variable+Error(block/canopy/variable)+month))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 0.2387 0.2387
##
## Error: block:canopy
## Df Sum Sq Mean Sq F value Pr(>F)
## canopy 3 86.89 28.96 0.302 0.824
## Residuals 3 287.93 95.98
##
## Error: block:canopy:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 4.955 4.955 9.244 0.0384 *
## canopy:variable 3 4.878 1.626 3.034 0.1558
## Residuals 4 2.144 0.536
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 4 3.507 0.156 0.693
## Residuals 207 4653 22.480summary(aov(data=NH4master[which(NH4master$invasion=="Post-adelgid"),],
value~canopy*variable+Error(block/canopy/variable)+month))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 99.92 99.92
##
## Error: block:canopy
## Df Sum Sq Mean Sq F value Pr(>F)
## canopy 3 662.5 220.8 1.252 0.429
## Residuals 3 529.2 176.4
##
## Error: block:canopy:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 16.12 16.12 5.624 0.0767 .
## canopy:variable 3 98.98 32.99 11.511 0.0195 *
## Residuals 4 11.47 2.87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 10 9.51 0.264 0.608
## Residuals 207 7448 35.98###############################################################################
#Carbon + Decomp
###############################################################################
summary(aov(data=carbon.rel,value~Treatment*variable+Error(DOY/Block/Treatment/variable))) #carbon ##
## Error: DOY
## Df Sum Sq Mean Sq
## Treatment 1 8.166 8.166
##
## Error: DOY:Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 37.94 37.94
##
## Error: DOY:Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 108.91 36.30 2.687 0.219
## Residuals 3 40.53 13.51
##
## Error: DOY:Block:Treatment:variable
## Df Sum Sq Mean Sq F value Pr(>F)
## variable 1 2.879 2.8792 2.155 0.216
## Treatment:variable 3 0.715 0.2382 0.178 0.906
## Residuals 4 5.345 1.3363
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 3.3 1.0975 0.561 0.641
## variable 1 0.0 0.0484 0.025 0.875
## Treatment:variable 3 0.4 0.1178 0.060 0.981
## Residuals 200 390.9 1.9547summary(aov(data = cell.coefs, m ~ Treatment * SubplotT + Error(Block/Treatment/SubplotT))) #cellulose##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 0.6982 0.6982
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 15.45 5.149 16.8 0.0222 *
## Residuals 3 0.92 0.307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 0.313 0.1564 0.421 0.670
## Treatment:SubplotT 6 3.673 0.6121 1.646 0.251
## Residuals 8 2.974 0.3718summary(aov(data = lig.coefs, m ~ Treatment * SubplotT + Error(Block/Treatment/SubplotT))) #lignin##
## Error: Block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 1.462 1.462
##
## Error: Block:Treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 3 6.439 2.146 1.633 0.348
## Residuals 3 3.943 1.314
##
## Error: Block:Treatment:SubplotT
## Df Sum Sq Mean Sq F value Pr(>F)
## SubplotT 2 1.573 0.7867 0.911 0.440
## Treatment:SubplotT 6 8.390 1.3983 1.619 0.258
## Residuals 8 6.909 0.8636Chunk 4: Set up summary dataframes and plot objects
###############################################################################
#Ants
###############################################################################
#Change factor on genus for faceting
ant.genus$Treatment<-factor(ant.genus$Treatment,levels(ant.genus$Treatment)[c(3,1,4,2)])
names <- c("HC", "G", "L", "HD")
levels(ant.genus$Treatment) <- names
#Summary dataframes of diversity
dsum <- ddply(ant.div, c("Treatment", "SubplotT"), summarise,
N = length(eH),e.H = mean(eH), sd = sd(eH), se = sd/sqrt(N))
dsum.pre <- ddply(ant.div[which(ant.div$invasion != "post"),],
c("Treatment", "SubplotT"), summarise,
N = length(eH),e.H = mean(eH), sd = sd(eH), se = sd/sqrt(N))
dsum.pre$invasion<-factor("Pre-adelgid")
dsum.post <- ddply(ant.div[which(ant.div$invasion != "pre"),],
c("Treatment", "SubplotT"), summarise,
N = length(eH),e.H = mean(eH), sd = sd(eH), se = sd/sqrt(N))
dsum.post$invasion <- factor("Post-adelgid")
dsum.invasion <- rbind(dsum.pre, dsum.post)
#Summary dataframes of Abundance
tsum<-ddply(ant.richness.count,c("Treatment","SubplotT"),summarise,
N=length(total),abundance=mean(total),sd=sd(total),se=sd/sqrt(N))
tsum.pre <- ddply(ant.richness.count[which(ant.richness.count$invasion != "post"),],
c("Treatment","SubplotT"),summarise,
N=length(total),abundance=mean(total),sd=sd(total),se=sd/sqrt(N))
tsum.pre$invasion <- as.factor("Pre-adelgid")
tsum.post <- ddply(ant.richness.count[which(ant.richness.count$invasion != "pre"),],
c("Treatment","SubplotT"),summarise,
N=length(total),abundance=mean(total),sd=sd(total),se=sd/sqrt(N))
tsum.post$invasion <- as.factor("Post-adelgid")
tsum.invasion <- rbind(tsum.pre, tsum.post)
#Base Plot Objects
int.abundance<-ggplot(tsum,aes(x=Treatment,y=abundance,group=SubplotT))+
geom_line(aes(linetype=SubplotT, color = SubplotT), alpha = 0.2)+
geom_point(aes(color = SubplotT), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax=abundance+se,ymin=abundance-se,color=SubplotT),
width=0,position=position_dodge(width = 0.1))
int.abundance.invasion <- ggplot(tsum.invasion,aes(x=Treatment,y=abundance,group=SubplotT))+
geom_line(aes(linetype=SubplotT, color = SubplotT), alpha = 0.2)+
geom_point(aes(color = SubplotT), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax=abundance+se,ymin=abundance-se,color=SubplotT),
width=0,position=position_dodge(width = 0.1))+facet_grid(.~invasion)
int.hill<-ggplot(dsum, aes(x = Treatment, y = e.H, group = SubplotT)) +
geom_line(aes(linetype = SubplotT, color = SubplotT), alpha = 0.2) +
geom_point(aes(color = SubplotT), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax = e.H + se, ymin = e.H - se,color = SubplotT),
width = 0, position = position_dodge(width = 0.1))
int.hill.invasion <- ggplot(dsum.invasion, aes(x = Treatment, y = e.H, group = SubplotT)) +
geom_line(aes(linetype = SubplotT, color = SubplotT), alpha = 0.2) +
geom_point(aes(color = SubplotT), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax = e.H + se, ymin = e.H - se,color = SubplotT),
width = 0, position = position_dodge(width = 0.1)) + facet_grid(.~invasion)
ts.genus <- ggplot(ant.genus, aes(x = Year, y = abundance, group = Genus)) +
geom_smooth(method = loess, se = FALSE, aes(color = Genus)) +
facet_grid(Treatment ~ SubplotT) +
geom_smooth(method = loess, se = FALSE, aes(x = Year, y = total), color = "black") +
theme(legend.title = element_blank())
ts.hill <- ggplot(ant.div, aes(x = Year, y = eH)) +
geom_smooth(method = loess, se = F, aes(color = Treatment, linetype = SubplotT))
###############################################################################
#Nitrogen
###############################################################################
#Format factors for proper plotting
NH4master$canopy <- factor(NH4master$canopy, levels(NH4master$canopy)[c(3, 1,
4, 2)])
NO3master$canopy <- factor(NO3master$canopy, levels(NO3master$canopy)[c(3, 1,
4, 2)])
NH4master$invasion <- factor(NH4master$invasion, levels(NH4master$invasion)[c(2, 1)])
NO3master$invasion <- factor(NO3master$invasion, levels(NO3master$invasion)[c(2, 1)])
levels(NH4master$canopy) <- c("HC", "G", "L", "HD")
levels(NO3master$canopy) <- c("HC", "G", "L", "HD")
#add before/after gap factor to break the line in 2012
NH4master$gap<-ifelse(NH4master$month<as.Date("2013-01-01"),"pre","post")
NO3master$gap<-ifelse(NO3master$month<as.Date("2013-01-01"),"pre","post")
#generate interaction dataframes of entire dataset
NH4sum <- ddply(NH4master, .(canopy, variable), summarise,
mean.NH4 = mean(value), N = length(value), sd = sd(value), se = sd/sqrt(N))
NO3sum <- ddply(NO3master, .(canopy, variable), summarise,
mean.NO3 = mean(value), N = length(value), sd = sd(value), se = sd/sqrt(N))
#generate interaction dataframes of pre/post invasion
NH4sum.invasion <- ddply(NH4master, .(canopy, variable, invasion), summarise,
mean.NH4 = mean(value), N = length(value), sd = sd(value), se = sd/sqrt(N))
NO3sum.invasion <- ddply(NO3master, .(canopy, variable, invasion), summarise,
mean.NO3 = mean(value), N = length(value), sd = sd(value), se = sd/sqrt(N))
#base timeseries plots
ts.NH4 <- ggplot(NH4master, aes(x = month, y = value, group = interaction(canopy, variable, gap))) +
geom_smooth(se = FALSE, method = loess, aes(color = canopy, linetype = variable))
ts.NO3 <- ggplot(NO3master, aes(x = month, y = value, group = interaction(canopy, variable, gap))) +
geom_smooth(se = FALSE, method = loess, aes(color = canopy, linetype = variable))
#generate base interaction plots for entire dataset
int.NH4 <- ggplot(NH4sum, aes(x = canopy, y = mean.NH4, group = variable)) +
geom_line(aes(linetype = variable, color = variable), alpha = 0.2) +
geom_point(aes(color = variable), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax = mean.NH4 + se, ymin = mean.NH4 - se, color = variable),
width = 0, position = position_dodge(width = 0.1))
int.NO3 <- ggplot(NO3sum, aes(x = canopy, y = mean.NO3, group = variable)) +
geom_line(aes(linetype = variable, color = variable), alpha = 0.2) +
geom_point(aes(color = variable), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax = mean.NO3 + se, ymin = mean.NO3 - se, color = variable),
width = 0, position = position_dodge(width = 0.1))
#generate base interaction plots for pre/post invasion
int.NH4.invasion <- ggplot(NH4sum.invasion, aes(x = canopy, y = mean.NH4, group = variable)) +
geom_line(aes(linetype = variable, color = variable), alpha = 0.2) +
geom_point(aes(color = variable), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax = mean.NH4 + se, ymin = mean.NH4 - se, color = variable),
width = 0, position = position_dodge(width = 0.1))+facet_grid(.~invasion)
int.NO3.invasion <- ggplot(NO3sum.invasion, aes(x = canopy, y = mean.NO3, group = variable)) +
geom_line(aes(linetype = variable, color = variable), alpha = 0.2) +
geom_point(aes(color = variable), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax = mean.NO3 + se, ymin = mean.NO3 - se, color = variable),
width = 0, position = position_dodge(width = 0.1))+facet_grid(.~invasion)+facet_grid(.~invasion)
###############################################################################
#Carbon and Decomp
###############################################################################
#reorder factors
cell.coefs$Treatment<-factor(cell.coefs$Treatment,levels(cell.coefs$Treatment)[c(3,1,4,2)])
lig.coefs$Treatment<-factor(lig.coefs$Treatment,levels(lig.coefs$Treatment)[c(3,1,4,2)])
carbon.rel$Treatment<-factor(carbon.rel$Treatment,levels(carbon.rel$Treatment)[c(3,1,4,2)])
#rename canopy treatments
levels(cell.coefs$Treatment) <- names
levels(lig.coefs$Treatment) <- names
levels(carbon.rel$Treatment) <- names
levels(cellmaster.extended$Treatment) <- names
levels(lignin.master$Treatment) <- names
#generate time series plot objects
ts.carbon<-ggplot(carbon.rel,aes(x=Date,y=value))+
geom_smooth(se=FALSE,method=loess,aes(color=Treatment,linetype=variable))
ts.lignin<-ggplot(lignin.master,aes(x=Date,y=Remaining))+
geom_smooth(se=FALSE,method=lm,aes(color=Treatment,linetype=SubplotT))
ts.cell<-ggplot(cellmaster.extended,aes(x=Date,y=Remaining))+
geom_smooth(se=FALSE,method=lm,aes(color=Treatment,linetype=SubplotT))
#Generate summary dataframes an dinteraction plots
carbonsum<-ddply(carbon.rel,.(Treatment,variable),summarise,mean.CO2=mean(value),N=length(value),sd=sd(value),se=sd/sqrt(N))
cellsum<-ddply(cell.coefs,.(Treatment,SubplotT),summarise,mean.cell.rate=mean(m),se=mean(slope.se))
ligsum<-ddply(lig.coefs,.(Treatment,SubplotT),summarise,mean.lig.rate=mean(m),se=mean(slope.se))
int.carbon<-ggplot(carbonsum,aes(x=Treatment,y=mean.CO2,group=variable))+
geom_line(aes(linetype=variable, color = variable), alpha = 0.2)+
geom_point(aes(color = variable), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax=mean.CO2+se,ymin=mean.CO2-se,color=variable),
width=0,position=position_dodge(width = 0.1))
int.cell<-ggplot(cellsum,aes(x=Treatment,y=mean.cell.rate,group=SubplotT))+
geom_line(aes(linetype=SubplotT, color = SubplotT), alpha = 0.2)+
geom_point(aes(color = SubplotT), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax=mean.cell.rate+se,ymin=mean.cell.rate-se,color=SubplotT),
width=0,position=position_dodge(width = 0.1))
int.lig<-ggplot(ligsum,aes(x=Treatment,y=mean.lig.rate,group=SubplotT))+
geom_line(aes(linetype=SubplotT, color = SubplotT), alpha = 0.2)+
geom_point(aes(color = SubplotT), size = 2,
position = position_dodge(width = 0.1)) +
geom_errorbar(aes(ymax=mean.lig.rate+se,ymin=mean.lig.rate-se,color=SubplotT),
width=0,position=position_dodge(width = 0.1))Chunk 5: Spiff up plots
#Common Elements
colors <- scale_color_manual(values = c("blue", "yellow3", "red", "plum3"),
name = "Canopy Treatment")
intcolors <- scale_color_manual(values = c("violetred3", "sienna3", "darkblue"),
name = "Ant Treatment")
intcolors.rel <- scale_color_manual(values = c("sienna3", "darkblue"),
name = "Ant Treatment")
lines <- scale_linetype_manual(values = c("dotted", "dashed", "solid"),
labels = c("Control", "Disturbance Control", "Ants Excluded"),
name = "Ant Treatment")
lines.rel <- scale_linetype_manual(values = c("dashed", "solid"),
labels = c ("Disturbance Control", "Ants Excluded"),
name = "Ant Treatment")
theme <- theme_bw() +
theme(text = element_text(size = 8),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.05),
plot.margin = unit(c(2, 2, 2, 2), "points"),
legend.margin = unit(-13, "points"))
nolegend <- theme(legend.position = "none")
noy <- theme(axis.title.y = element_blank(), axis.text.y = element_blank())
nox <- theme(axis.title.x = element_blank(), axis.text.x = element_blank())
timeframe = xlim(as.Date(c("2006-6-01","2014-11-01")))
#Recurring y-axis labels
yhill <- ylab(expression("Effective species " ~ (E^italic("H'"))))
ycarbon <- ylab(expression("Soil respiration"~(mg~C~m^"-2"~hr^"-1")))
yNH4 <- ylab(expression(Ammonium ~ (g ~ g ~ resin^"-1" ~ day^"-1")))
yNO3 <- ylab(expression(Nitrate ~ (g ~ g ~ resin^"-1" ~ day^"-1")))
#Add elements to panels
ts.genus <- ts.genus + theme + ylab("Abundance") + ggtitle("A") + xlab(NULL) +
scale_x_continuous(breaks = seq(2007, 2013, 2)) +
theme(legend.text = element_text(face = "italic"))
ts.hill <- ts.hill + theme + colors + lines + xlab(NULL) + yhill + ggtitle ("B")
int.abundance <- int.abundance + theme + lines + intcolors + nox + nolegend +
ylab("Abundance") + ggtitle("A")
int.hill <- int.hill + theme + lines + intcolors + nox + nolegend + yhill +
ggtitle("B")
int.abundance.invasion <- int.abundance.invasion + theme + lines + intcolors +
nox + nolegend + ylab("Abundance") + ggtitle("C")
int.hill.invasion <- int.hill.invasion + theme + lines + intcolors + yhill +
xlab ("Canopy Treatment") + ggtitle("C") + nolegend
int.legend = int.hill + ylim(5, 6) + nox + noy + theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line = element_blank(),
axis.text = element_blank(), axis.ticks = element_blank(), legend.position = c(.5,.8)) + ggtitle("")
ts.cell <- ts.cell + theme + lines + colors + nox + timeframe + nolegend +
ylab("% Cellulose Remaining") + ggtitle("A")
ts.lignin <- ts.lignin + theme + lines + colors + nox + timeframe + nolegend +
ylab("% Lignin Remaining") + ggtitle("B")
ts.carbon <- ts.carbon + theme + lines.rel + colors + xlab(NULL) + timeframe + nolegend +
ycarbon + ggtitle("C")
ts.NH4 <- ts.NH4 + theme + lines.rel + colors + nox + timeframe + nolegend +
yNH4 + ggtitle("D") + coord_cartesian(ylim = c(-4.25, 4.25))
ts.NO3 <- ts.NO3 + theme + lines.rel + colors + timeframe + nolegend +
yNO3 + ggtitle("E") + coord_cartesian(ylim = c(-2.5, 2.5)) + xlab(NULL)
ts.legend <- ts.hill + coord_cartesian(ylim = c(100,101)) + nox + noy +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(), axis.line = element_blank(),
axis.text = element_blank(), axis.ticks = element_blank(), legend.position = c(.5,.5),
legend.box = "horizontal") + ggtitle("")
int.cell <- int.cell + theme + lines + intcolors + nox + nolegend +
ylab(expression("Cellulose "~(~"%"~change~month^"-1"))) + ggtitle("A")
int.lig <- int.lig + theme + lines + intcolors + nox + nolegend +
ylab(expression("Lignin "~(~"%"~change~month^"-1"))) + ggtitle("B")
int.carbon <- int.carbon + theme + lines.rel + intcolors.rel + nolegend +
ycarbon + ggtitle("C") + xlab("Canopy Treatment")
int.NH4 <- int.NH4 + theme + lines.rel + intcolors.rel + nolegend + nox +
yNH4 + ggtitle("D")
int.NO3 <- int.NO3 + theme + lines.rel + intcolors.rel + nolegend +
yNO3 + ggtitle("E") + xlab("Canopy Treatment")
int.NH4.invasion <- int.NH4.invasion + theme + lines.rel + intcolors.rel +
nolegend + nox + yNH4 + ggtitle("A")
int.NO3.invasion <- int.NO3.invasion + theme + lines.rel + intcolors.rel +
nolegend + xlab("Canopy Treatment") + yNO3 + ggtitle("B")
int.legend.rel <- int.NH4 + coord_cartesian(ylim = c(100,101)) + nox + noy +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(), axis.line = element_blank(),
axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "top",
legend.box = "horizontal") + ggtitle("")Chunk 6: Generate Figures!
#multiplot function from R cookbook
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
require(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
pdf(file = "Figure 3- Ant Timeseries.pdf", width = 6.1, height = 6.8)
multiplot(ts.genus, ts.hill, cols=1)
dev.off()## pdf
## 2pdf(file = "Figure 4- Ant Interactions.pdf", width = 2.9, height = 6.8)
multiplot(int.abundance, int.hill, int.abundance.invasion, int.hill.invasion, int.legend, cols = 1)## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead## Warning: Removed 12 rows containing missing values (geom_point).dev.off()## pdf
## 2pdf(file = "Figure 6- Ecosystem Timeseries.pdf", width = 6.1, height = 5)
multiplot(ts.cell, ts.lignin, ts.carbon, ts.NH4, ts.NO3, ts.legend, cols = 2)## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 3 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).dev.off()## pdf
## 2pdf(file = "Figure 7- Ecosystem Interactions.pdf", width = 6.1, height = 5)
multiplot(int.cell, int.lig, int.carbon, int.NH4, int.NO3, int.legend, cols = 2)## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead## Warning: Removed 12 rows containing missing values (geom_point).dev.off()## pdf
## 2pdf(file = "Figure 8- Nitrogen Invasion Interactions.pdf", width = 2.9, height = 5)
multiplot(int.NH4.invasion, int.NO3.invasion, int.legend.rel, cols = 1)## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y insteaddev.off()## pdf
## 2Chunk 7: Path analysis
###############################################################################
#Cull all data and transform untransformed variables for consistency
###############################################################################
#decomp: put data into wide form with a seperate column for each subplot type
lig.rel<-cast(lig.coefs,Treatment+Block~SubplotT,value="m")
cell.rel<-cast(cell.coefs,Treatment+Block~SubplotT,value="m")
#decomp: subtract control from treatments
lig.rel$D.C<-lig.rel$D-lig.rel$C
lig.rel$X.C<-lig.rel$X-lig.rel$C
cell.rel$D.C<-cell.rel$D-cell.rel$C
cell.rel$X.C<-cell.rel$X-cell.rel$C
#decomp: put data back into long form
lig.rel<-as.data.frame(lig.rel)
lig.rel<-melt(lig.rel,id.vars=c("Treatment","Block"),measure.vars=c("D.C","X.C"))
cell.rel<-as.data.frame(cell.rel)
cell.rel<-melt(cell.rel,id.vars=c("Treatment","Block"),measure.vars=c("D.C","X.C"))
#Ants: build, clean, then dissassemble master df
ant.path<-ant.div
ant.path$total<-ant.richness.count$total
ant.path$richness<-ant.richness.count$richness
ant.path<-ant.path[(which(ant.path$Year==2009)),]
row.names(ant.path)<-NULL
hill.path<-ant.path[c(1,2,4,5,7)]
abundance.path<-ant.path[c(1,2,4,5,9)]
richness.path<-ant.path[c(1,2,4,5,10)]
#ants: transform
hill.path<-cast(hill.path, Plot+Block+Treatment~SubplotT,value="eH")
hill.path$D.C<-hill.path$"Disturbance Control"-hill.path$"Control"
hill.path$X.C<-hill.path$"Ants Excluded"-hill.path$"Control"
hill.path<-as.data.frame(hill.path)
hill.path<-melt(hill.path,id.vars=c("Plot","Block","Treatment"),measure.vars=c("D.C","X.C"))
abundance.path<-cast(abundance.path, Plot+Block+Treatment~SubplotT,value="total")
abundance.path$D.C<-abundance.path$"Disturbance Control"-abundance.path$"Control"
abundance.path$X.C<-abundance.path$"Ants Excluded"-abundance.path$"Control"
abundance.path<-as.data.frame(abundance.path)
abundance.path<-melt(abundance.path,id.vars=c("Plot","Block","Treatment"),measure.vars=c("D.C","X.C"))
richness.path<-cast(richness.path, Plot+Block+Treatment~SubplotT,value="richness")
richness.path$D.C<-richness.path$"Disturbance Control"-richness.path$"Control"
richness.path$X.C<-richness.path$"Ants Excluded"-richness.path$"Control"
richness.path<-as.data.frame(richness.path)
richness.path<-melt(richness.path,id.vars=c("Plot","Block","Treatment"),measure.vars=c("D.C","X.C"))
#Nitrogen: cull
NH4master$month<-as.Date(NH4master$month)
NO3master$month<-as.Date(NO3master$month)
NH4master<-NH4master[(which(NH4master$month>as.Date("2008-10-01")&NH4master$month<as.Date("2010-07-01"))),]
NO3master<-NO3master[(which(NO3master$month>as.Date("2008-10-01")&NO3master$month<as.Date("2010-07-01"))),]
NH4master<-ddply(NH4master,.(block,canopy,plot,variable),summarise,value=mean(value))
NO3master<-ddply(NO3master,.(block,canopy,plot,variable),summarise,value=mean(value))
###############################################################################
#Pool data for analysis
###############################################################################
#put everything in same order
cell.rel<-cell.rel[with(cell.rel,order(Treatment, Block,variable)),]
lig.rel<-lig.rel[with(lig.rel,order(Treatment, Block,variable)),]
richness.path<-richness.path[with(richness.path,order(Treatment, Block,variable)),]
abundance.path<-abundance.path[with(abundance.path,order(Treatment, Block,variable)),]
hill.path<-hill.path[with(hill.path,order(Treatment, Block,variable)),]
NH4master<-NH4master[with(NH4master,order(canopy, block,variable)),]
NO3master<-NO3master[with(NO3master,order(canopy, block,variable)),]
#create base of dataframe for path modeling
path.master<-hill.path[c(1,2,3,4)]
row.names(path.master)<-NULL
#Add canopy information as three binary variables telling whether plot is hardwood, girdle or log.
#These three treatments are then implicitly compared to hemlock control plots
path.master$hardcon<-ifelse(path.master$Plot>6,1,0)
path.master$girdle<-ifelse(path.master$Plot==1|path.master$Plot==5,1,0)
path.master$log<-ifelse(path.master$Plot==2|path.master$Plot==4,1,0)
#add a column for each variable
path.master$abundance<-abundance.path$value
path.master$richness<-richness.path$value
path.master$hill<-hill.path$value
path.master$lig<-lig.rel$value
path.master$cell<-cell.rel$value
path.master$NH4<-NH4master$value
path.master$NO3<-NO3master$value
###############################################################################
#Run model
###############################################################################
#construct inner model matrix. a value of 1 means the row variable affects the column variable
girdle<-c(0,0,0,1,1,1,1,1,1)
log<-c(0,0,0,1,1,1,1,1,1)
hardcon<-c(0,0,0,1,1,1,1,1,1)
abundance<-c(0,0,0,0,0,1,1,1,1)
diversity<-c(0,0,0,0,0,1,1,1,1)
lig<-c(0,0,0,0,0,0,0,1,1)
cell<-c(0,0,0,0,0,0,0,1,1)
NH4<-rep(0,9)
NO3<-rep(0,9)
inner.model<-cbind(hardcon,girdle,log,abundance,diversity,lig,cell,NH4,NO3)
rownames(inner.model)<-colnames(inner.model)
#identify which manifest variables define which latent variables in the inner model.
#only "diversity" is assigned more than one manifest variable (richness and hill)
outer.blocks<-list(5,6,7,8,9:10,11,12,13,14)
#run path model, show loadings and weights
pls<-plspm(path.master,inner.model,outer.blocks,modes=NULL)
plot(pls,what="loadings")
plot(pls,what="weights")
#output results
innerplot(pls)
summary(pls)## PARTIAL LEAST SQUARES PATH MODELING (PLS-PM)
##
## ----------------------------------------------------------
## MODEL SPECIFICATION
## 1 Number of Cases 16
## 2 Latent Variables 9
## 3 Manifest Variables 10
## 4 Scale of Data Standardized Data
## 5 Non-Metric PLS FALSE
## 6 Weighting Scheme centroid
## 7 Tolerance Crit 1e-06
## 8 Max Num Iters 100
## 9 Convergence Iters 3
## 10 Bootstrapping FALSE
## 11 Bootstrap samples NULL
##
## ----------------------------------------------------------
## BLOCKS DEFINITION
## Block Type Size Mode
## 1 hardcon Exogenous 1 A
## 2 girdle Exogenous 1 A
## 3 log Exogenous 1 A
## 4 abundance Endogenous 1 A
## 5 diversity Endogenous 2 A
## 6 lig Endogenous 1 A
## 7 cell Endogenous 1 A
## 8 NH4 Endogenous 1 A
## 9 NO3 Endogenous 1 A
##
## ----------------------------------------------------------
## BLOCKS UNIDIMENSIONALITY
## Mode MVs C.alpha DG.rho eig.1st eig.2nd
## hardcon A 1 1.00 1.000 1.00 0.000
## girdle A 1 1.00 1.000 1.00 0.000
## log A 1 1.00 1.000 1.00 0.000
## abundance A 1 1.00 1.000 1.00 0.000
## diversity A 2 0.83 0.922 1.71 0.291
## lig A 1 1.00 1.000 1.00 0.000
## cell A 1 1.00 1.000 1.00 0.000
## NH4 A 1 1.00 1.000 1.00 0.000
## NO3 A 1 1.00 1.000 1.00 0.000
##
## ----------------------------------------------------------
## OUTER MODEL
## weight loading communality redundancy
## hardcon
## 1 hardcon 1.000 1.000 1.000 0.000
## girdle
## 2 girdle 1.000 1.000 1.000 0.000
## log
## 3 log 1.000 1.000 1.000 0.000
## abundance
## 4 abundance 1.000 1.000 1.000 0.380
## diversity
## 5 richness 0.551 0.927 0.860 0.201
## 5 hill 0.530 0.921 0.849 0.199
## lig
## 6 lig 1.000 1.000 1.000 0.529
## cell
## 7 cell 1.000 1.000 1.000 0.451
## NH4
## 8 NH4 1.000 1.000 1.000 0.585
## NO3
## 9 NO3 1.000 1.000 1.000 0.664
##
## ----------------------------------------------------------
## CROSSLOADINGS
## hardcon girdle log abundance diversity lig
## hardcon
## 1 hardcon 1.000 -0.3333 -0.3333 -0.179 -0.4424 0.2993
## girdle
## 2 girdle -0.333 1.0000 -0.3333 0.068 0.0264 0.1689
## log
## 3 log -0.333 -0.3333 1.0000 -0.426 0.3283 -0.7212
## abundance
## 4 abundance -0.179 0.0680 -0.4262 1.000 -0.0695 0.3336
## diversity
## 5 richness -0.480 0.0685 0.2056 0.171 0.9275 -0.0415
## 5 hill -0.336 -0.0215 0.4054 -0.309 0.9213 -0.4689
## lig
## 6 lig 0.299 0.1689 -0.7212 0.334 -0.2715 1.0000
## cell
## 7 cell 0.583 -0.0919 -0.0381 -0.372 -0.1271 -0.3665
## NH4
## 8 NH4 0.372 -0.2077 -0.0696 -0.406 0.3200 0.0958
## NO3
## 9 NO3 0.226 -0.3095 -0.2011 0.268 0.0960 0.4933
## cell NH4 NO3
## hardcon
## 1 hardcon 0.5831 0.3721 0.2264
## girdle
## 2 girdle -0.0919 -0.2077 -0.3095
## log
## 3 log -0.0381 -0.0696 -0.2011
## abundance
## 4 abundance -0.3724 -0.4057 0.2681
## diversity
## 5 richness -0.3144 0.1621 0.2470
## 5 hill 0.0872 0.4350 -0.0758
## lig
## 6 lig -0.3665 0.0958 0.4933
## cell
## 7 cell 1.0000 0.3280 -0.0728
## NH4
## 8 NH4 0.3280 1.0000 0.1762
## NO3
## 9 NO3 -0.0728 0.1762 1.0000
##
## ----------------------------------------------------------
## INNER MODEL
## $abundance
## Estimate Std. Error t value Pr(>|t|)
## Intercept -2.09e-16 0.227 -9.19e-16 1.0000
## hardcon -5.37e-01 0.278 -1.93e+00 0.0774
## girdle -3.52e-01 0.278 -1.27e+00 0.2298
## log -7.23e-01 0.278 -2.60e+00 0.0233
##
## $diversity
## Estimate Std. Error t value Pr(>|t|)
## Intercept -1.31e-17 0.253 -5.18e-17 1.000
## hardcon -3.98e-01 0.309 -1.28e+00 0.223
## girdle -4.60e-02 0.309 -1.49e-01 0.884
## log 1.80e-01 0.309 5.83e-01 0.571
##
## $lig
## Estimate Std. Error t value Pr(>|t|)
## Intercept -1.09e-17 0.217 -5.04e-17 1.000
## hardcon 5.62e-02 0.324 1.73e-01 0.866
## girdle -4.58e-02 0.284 -1.61e-01 0.875
## log -6.91e-01 0.333 -2.07e+00 0.065
## abundance 5.12e-02 0.277 1.85e-01 0.857
## diversity -1.50e-02 0.249 -6.03e-02 0.953
##
## $cell
## Estimate Std. Error t value Pr(>|t|)
## Intercept 2.45e-16 0.234 1.04e-15 1.0000
## hardcon 7.31e-01 0.350 2.09e+00 0.0634
## girdle 2.13e-01 0.306 6.97e-01 0.5018
## log 1.59e-01 0.360 4.41e-01 0.6687
## abundance -1.80e-01 0.299 -6.01e-01 0.5610
## diversity 1.26e-01 0.269 4.69e-01 0.6494
##
## $NH4
## Estimate Std. Error t value Pr(>|t|)
## Intercept -1.37e-16 0.228 -6.02e-16 1.0000
## hardcon 1.70e-01 0.533 3.20e-01 0.7573
## girdle -2.54e-01 0.312 -8.16e-01 0.4383
## log -2.80e-01 0.479 -5.85e-01 0.5749
## abundance -4.54e-01 0.299 -1.52e+00 0.1682
## diversity 5.57e-01 0.268 2.08e+00 0.0711
## lig 2.58e-01 0.567 4.55e-01 0.6613
## cell 1.91e-01 0.526 3.63e-01 0.7259
##
## $NO3
## Estimate Std. Error t value Pr(>|t|)
## Intercept -2.11e-16 0.205 -1.03e-15 1.0000
## hardcon -5.47e-01 0.480 -1.14e+00 0.2875
## girdle -4.73e-01 0.281 -1.69e+00 0.1304
## log 6.33e-01 0.432 1.47e+00 0.1806
## abundance 3.37e-01 0.270 1.25e+00 0.2474
## diversity 1.98e-01 0.241 8.20e-01 0.4359
## lig 1.47e+00 0.511 2.88e+00 0.0205
## cell 9.17e-01 0.474 1.94e+00 0.0890
##
## ----------------------------------------------------------
## CORRELATIONS BETWEEN LVs
## hardcon girdle log abundance diversity lig
## hardcon 1.000 -0.3333 -0.3333 -0.1791 -0.4424 0.2993
## girdle -0.333 1.0000 -0.3333 0.0680 0.0264 0.1689
## log -0.333 -0.3333 1.0000 -0.4262 0.3283 -0.7212
## abundance -0.179 0.0680 -0.4262 1.0000 -0.0695 0.3336
## diversity -0.442 0.0264 0.3283 -0.0695 1.0000 -0.2715
## lig 0.299 0.1689 -0.7212 0.3336 -0.2715 1.0000
## cell 0.583 -0.0919 -0.0381 -0.3724 -0.1271 -0.3665
## NH4 0.372 -0.2077 -0.0696 -0.4057 0.3200 0.0958
## NO3 0.226 -0.3095 -0.2011 0.2681 0.0960 0.4933
## cell NH4 NO3
## hardcon 0.5831 0.3721 0.2264
## girdle -0.0919 -0.2077 -0.3095
## log -0.0381 -0.0696 -0.2011
## abundance -0.3724 -0.4057 0.2681
## diversity -0.1271 0.3200 0.0960
## lig -0.3665 0.0958 0.4933
## cell 1.0000 0.3280 -0.0728
## NH4 0.3280 1.0000 0.1762
## NO3 -0.0728 0.1762 1.0000
##
## ----------------------------------------------------------
## SUMMARY INNER MODEL
## Type R2 Block_Communality Mean_Redundancy AVE
## hardcon Exogenous 0.000 1.000 0.000 1.000
## girdle Exogenous 0.000 1.000 0.000 1.000
## log Exogenous 0.000 1.000 0.000 1.000
## abundance Endogenous 0.380 1.000 0.380 1.000
## diversity Endogenous 0.234 0.855 0.200 0.855
## lig Endogenous 0.529 1.000 0.529 1.000
## cell Endogenous 0.451 1.000 0.451 1.000
## NH4 Endogenous 0.585 1.000 0.585 1.000
## NO3 Endogenous 0.664 1.000 0.664 1.000
##
## ----------------------------------------------------------
## GOODNESS-OF-FIT
## [1] 0.6363
##
## ----------------------------------------------------------
## TOTAL EFFECTS
## relationships direct indirect total
## 1 hardcon -> girdle 0.0000 0.0000 0.0000
## 2 hardcon -> log 0.0000 0.0000 0.0000
## 3 hardcon -> abundance -0.5374 0.0000 -0.5374
## 4 hardcon -> diversity -0.3976 0.0000 -0.3976
## 5 hardcon -> lig 0.0562 -0.0215 0.0347
## 6 hardcon -> cell 0.7306 0.0465 0.7771
## 7 hardcon -> NH4 0.1705 0.1797 0.3502
## 8 hardcon -> NO3 -0.5472 0.5038 -0.0434
## 9 girdle -> log 0.0000 0.0000 0.0000
## 10 girdle -> abundance -0.3521 0.0000 -0.3521
## 11 girdle -> diversity -0.0460 0.0000 -0.0460
## 12 girdle -> lig -0.0458 -0.0173 -0.0631
## 13 girdle -> cell 0.2133 0.0575 0.2708
## 14 girdle -> NH4 -0.2542 0.1695 -0.0846
## 15 girdle -> NO3 -0.4730 0.0278 -0.4453
## 16 log -> abundance -0.7228 0.0000 -0.7228
## 17 log -> diversity 0.1805 0.0000 0.1805
## 18 log -> lig -0.6910 -0.0397 -0.7307
## 19 log -> cell 0.1586 0.1526 0.3112
## 20 log -> NH4 -0.2803 0.2992 0.0189
## 21 log -> NO3 0.6334 -0.9974 -0.3639
## 22 abundance -> diversity 0.0000 0.0000 0.0000
## 23 abundance -> lig 0.0512 0.0000 0.0512
## 24 abundance -> cell -0.1797 0.0000 -0.1797
## 25 abundance -> NH4 -0.4536 -0.0211 -0.4747
## 26 abundance -> NO3 0.3365 -0.0894 0.2471
## 27 diversity -> lig -0.0150 0.0000 -0.0150
## 28 diversity -> cell 0.1259 0.0000 0.1259
## 29 diversity -> NH4 0.5570 0.0202 0.5772
## 30 diversity -> NO3 0.1978 0.0933 0.2912
## 31 lig -> cell 0.0000 0.0000 0.0000
## 32 lig -> NH4 0.2581 0.0000 0.2581
## 33 lig -> NO3 1.4713 0.0000 1.4713
## 34 cell -> NH4 0.1910 0.0000 0.1910
## 35 cell -> NO3 0.9167 0.0000 0.9167
## 36 NH4 -> NO3 0.0000 0.0000 0.0000Chunk 8: additional analysis and figures: ant PCA, light, veg. Code by AME
###############################################################################
#Ant PCA
###############################################################################
colSums(ant.master.wide[,6:37])>5## aphful aphpic camher camnea camnov campen forarg forase forinc
## TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE FALSE
## forlas forneo1 forneo2 forper forsub1 forsub2 forsub3 lasali lasint
## FALSE TRUE FALSE FALSE TRUE TRUE TRUE TRUE FALSE
## lasnea lasumb myrame1 myrdet myrnea myrpun myrscu myrsmi stebre
## FALSE TRUE FALSE FALSE FALSE TRUE TRUE TRUE FALSE
## stedie steimp stesch tapses temlon
## FALSE TRUE FALSE TRUE TRUEant.wide.reduced <- ant.master.wide[,colSums(ant.master.wide[6:37])>5]
ant.wide.reduced2 <- cbind(ant.master.wide$Year, ant.wide.reduced)
names(ant.wide.reduced2)[1]<-"Year"
ants.pca <- prcomp(ant.wide.reduced2[,6:23], scale.=TRUE)
summary(ants.pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.22327 1.21668 1.13540 1.0941 1.04755 1.01416
## Proportion of Variance 0.08313 0.08224 0.07162 0.0665 0.06097 0.05714
## Cumulative Proportion 0.08313 0.16537 0.23699 0.3035 0.36446 0.42160
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 1.00827 1.00612 1.00463 1.00322 1.00267 1.00232
## Proportion of Variance 0.05648 0.05624 0.05607 0.05591 0.05585 0.05581
## Cumulative Proportion 0.47808 0.53432 0.59039 0.64630 0.70215 0.75797
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 0.95510 0.95328 0.8818 0.79655 0.76123 0.73756
## Proportion of Variance 0.05068 0.05049 0.0432 0.03525 0.03219 0.03022
## Cumulative Proportion 0.80865 0.85913 0.9023 0.93759 0.96978 1.00000plot(ants.pca)
ants.pca$rotation## PC1 PC2 PC3 PC4 PC5
## aphful 0.199343372 -0.07592543 0.2529752528 0.136745567 0.045515376
## camher -0.067191811 -0.06887589 -0.0611886070 0.078620143 -0.054922576
## campen 0.555113480 0.06015902 0.4243405456 -0.002903287 0.045518435
## forinc 0.341709460 0.02462506 0.4889798519 0.085856034 0.189700386
## forlas -0.063914799 0.02686698 -0.0371176160 0.103977867 0.080098551
## forneo1 0.464581230 -0.00741541 -0.5022060169 -0.075536812 0.141922655
## forneo2 -0.067913277 -0.02716756 -0.0363727425 0.132302315 0.030548591
## forsub2 0.493899651 0.05329816 -0.4479652032 -0.169414119 -0.052244536
## lasnea -0.139365989 0.27981832 -0.0145315137 -0.268119947 0.393671020
## lasumb 0.014813594 0.09135810 -0.1086566893 0.167159754 -0.136598762
## myrame1 -0.071326577 0.07372443 -0.0370406308 0.116777685 0.143673283
## myrpun -0.106997265 0.65429969 0.0006429475 0.072310431 0.173047817
## myrsmi 0.129182908 0.36597841 0.0455038888 -0.100645951 -0.536310726
## stebre -0.021946734 -0.10697247 0.0182136152 0.136296899 -0.057600079
## stedie -0.004861577 0.38969611 0.0604755342 0.084011858 -0.498999170
## steimp -0.080464506 0.07936647 0.1057974869 -0.677104809 0.002669575
## tapses 0.032778051 0.37393258 -0.0309664395 0.077329336 0.396932876
## temlon -0.037954949 -0.13172474 0.1609576339 -0.534888655 -0.069996770
## PC6 PC7 PC8 PC9 PC10
## aphful 0.105343089 0.694688575 0.207410143 0.318410300 0.2017304687
## camher -0.082689134 0.053168803 0.074269730 0.034802613 -0.1575699635
## campen 0.009389609 -0.012409185 -0.003867534 -0.003569784 0.0006991487
## forinc 0.119755754 -0.319751739 -0.107538486 -0.258059177 -0.1822480342
## forlas 0.038352956 0.056150530 0.076496187 0.004797892 -0.4886737712
## forneo1 -0.062019851 0.008407271 -0.022466135 -0.047843634 -0.0456347684
## forneo2 -0.027986955 0.088729930 0.247910978 -0.821057826 0.2769182854
## forsub2 0.078404567 0.084797649 -0.038870690 -0.067949094 -0.0643052484
## lasnea 0.438507724 0.293999608 -0.286815694 -0.132400893 0.0749188492
## lasumb 0.629333909 -0.456296155 0.110026797 0.251630893 0.2445358760
## myrame1 0.095720138 0.065961647 0.081674336 -0.004029815 -0.6573562046
## myrpun 0.040197485 0.003443942 0.001159544 -0.001201180 -0.0088665014
## myrsmi -0.062729887 0.008603854 -0.006247255 -0.001492342 0.0132702745
## stebre -0.164420746 0.018148121 -0.850630852 0.004655208 0.0631471273
## stedie -0.074169575 0.169805625 -0.052351118 -0.127631978 -0.1337235221
## steimp 0.023024910 0.016599188 -0.024440026 -0.003939807 0.0204052144
## tapses -0.550149687 -0.202070476 0.129747807 0.226615484 0.2186896584
## temlon -0.128818602 -0.156679295 0.160414701 0.034374499 -0.1080027704
## PC11 PC12 PC13 PC14 PC15
## aphful 0.0911326324 0.000000e+00 -0.05495802 -0.24190653 -0.086051106
## camher -0.8659038372 2.604691e-01 -0.12065373 -0.31707535 -0.028344754
## campen 0.0005571563 -1.526557e-15 -0.05533117 -0.10234514 0.029181007
## forinc -0.1687259717 1.183950e-02 0.06661269 0.13379314 -0.076036895
## forlas -0.0425841650 -8.050862e-01 -0.17799848 -0.18533297 -0.030971043
## forneo1 -0.0067074360 7.237266e-14 0.18978950 -0.02263616 -0.094411351
## forneo2 0.1104200411 -5.951073e-13 -0.18849499 -0.31634505 -0.044051328
## forsub2 -0.0057082965 1.180670e-13 -0.07845534 -0.07459712 0.025925806
## lasnea -0.1770354488 -1.004094e-12 0.10632774 0.06314306 0.311374545
## lasumb 0.0549683188 -3.142243e-13 -0.07267343 -0.36112738 -0.150325257
## myrame1 0.3627308260 5.327776e-01 -0.21146949 -0.12801728 -0.031253543
## myrpun 0.0083768248 5.914713e-14 0.09083459 -0.16137130 0.001882526
## myrsmi -0.0090080672 -7.334758e-14 -0.49287257 0.12074575 0.441745448
## stebre 0.1271272632 6.023029e-13 -0.14096676 -0.40902221 -0.055730551
## stedie -0.0143936207 2.404119e-13 0.53199696 -0.04779894 -0.348294766
## steimp -0.0121961835 -1.080455e-13 -0.34385913 0.01623416 -0.625457700
## tapses 0.0223460849 -3.915861e-13 -0.11149370 -0.13052123 -0.043318455
## temlon 0.1274692372 8.412993e-13 0.34567528 -0.54201269 0.373514585
## PC16 PC17 PC18
## aphful 0.069515000 0.321442264 0.07439658
## camher -0.026372961 -0.030301787 -0.03788070
## campen -0.244214615 -0.608474583 -0.24808127
## forinc 0.185646882 0.488784791 0.20520723
## forlas -0.078371705 -0.024728591 0.05971485
## forneo1 -0.504722304 0.368543261 -0.24738120
## forneo2 -0.058501207 -0.006405222 0.02743590
## forsub2 0.615687376 -0.192773831 0.26241283
## lasnea -0.233072934 -0.104345074 0.28683173
## lasumb -0.121410646 0.007987263 0.13198209
## myrame1 -0.109660393 -0.033415326 0.10419507
## myrpun 0.340008844 0.108427670 -0.60152450
## myrsmi -0.176569958 0.244054398 0.04304836
## stebre 0.008776798 0.063251450 -0.01968402
## stedie -0.127161453 -0.090677351 0.29314825
## steimp -0.036777155 0.045381610 -0.04959434
## tapses -0.099511671 -0.040218659 0.43787166
## temlon 0.040154977 0.116956506 0.01586088biplot(ants.pca)
antpca.out <- cbind(ant.wide.reduced2[,1:5],ants.pca$x[,1:2])
antpca.out$Year <- as.factor(antpca.out$Year)
str(antpca.out)## 'data.frame': 216 obs. of 7 variables:
## $ Year : Factor w/ 9 levels "2006","2007",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Block : Factor w/ 2 levels "Ridge","Valley": 2 2 2 2 2 2 2 2 2 1 ...
## $ Treatment: Factor w/ 4 levels "G","HC","HD",..: 1 1 1 4 4 4 2 2 2 4 ...
## $ Plot : num 1 1 1 2 2 2 3 3 3 4 ...
## $ SubplotT : Factor w/ 3 levels "Control","Disturbance Control",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ PC1 : num -0.223 -0.483 -0.483 -0.483 -0.483 ...
## $ PC2 : num -0.451 -0.479 -0.479 -0.479 -0.479 ...antpca.out$Treatment <- reorder(antpca.out$Treatment, new.order=c(2,1,4,3))
ggplot(antpca.out[antpca.out$Block=="Valley",], aes(x=PC1, y=PC2, colour=Year)) +
geom_path(aes(group=SubplotT), colour="black", arrow=arrow(angle=15, length=unit(5, "mm"))) +
geom_point(size=3) +
facet_grid(SubplotT~Treatment) +
theme_bw() +
ggtitle("Valley block") 
ggplot(antpca.out[antpca.out$Block=="Ridge",], aes(x=PC1, y=PC2, colour=Year)) +
geom_path(aes(group=SubplotT), colour="black", arrow=arrow(angle=15, length=unit(5, "mm"))) +
geom_point(size=3) +
facet_grid(SubplotT~Treatment) +
theme_bw() +
ggtitle("Ridge block")
#try to get both on the same plot
pdf("Figure 5- antpca.pdf", height=6, width=6.1)
ggplot(antpca.out[antpca.out$Block=="Valley",], aes(x=PC1, y=PC2, colour=Year)) +
geom_path(aes(group=SubplotT), colour="black", arrow=arrow(angle=15, length=unit(3, "mm"), type="closed")) +
# geom_point(size=3) +
facet_grid(Treatment~SubplotT) +
# geom_point(data=antpca.out[antpca.out$Block=="Ridge",],aes(x=PC1, y=PC2), size=5, shape=20, fill="White") +
geom_path(data=antpca.out[antpca.out$Block=="Ridge",],aes(group=SubplotT), colour="red", arrow=arrow(angle=15, length=unit(3, "mm"), type="closed")) +
theme +
#ggtitle("Both blocks")+
xlab("Principal Component 1") +
ylab("Principal Component 2")
dev.off()## pdf
## 2#pc loadings
ant.loadings <- ants.pca$rotation[,c(1,2)]
ant.1 <- data.frame(ant.loadings[,1][abs(ant.loadings[,1])>.2])
ant.1 <- cbind(rownames(ant.1), ant.1)
rownames(ant.1) <- NULL
names(ant.1) <- c("species", "PC1")
ant.2 <- data.frame(ant.loadings[,2][abs(ant.loadings[,2])>.2])
ant.2 <- cbind(rownames(ant.2), ant.2)
rownames(ant.2) <- NULL
names(ant.2) <- c("species", "PC2")
ant.big.loadings<- merge(ant.1, ant.2, all=TRUE)
ant.big.loadings## species PC1 PC2
## 1 campen 0.5551135 NA
## 2 forinc 0.3417095 NA
## 3 forneo1 0.4645812 NA
## 4 forsub2 0.4938997 NA
## 5 lasnea NA 0.2798183
## 6 myrpun NA 0.6542997
## 7 myrsmi NA 0.3659784
## 8 stedie NA 0.3896961
## 9 tapses NA 0.3739326###############################################################################
#Veg
###############################################################################
#Read in data, center and scale
veg <- read.csv("veg2014.csv")
str(veg)## 'data.frame': 48 obs. of 50 variables:
## $ block : Factor w/ 2 levels "Ridge","Valley": 2 2 2 2 2 2 1 1 1 1 ...
## $ canopy : Factor w/ 4 levels "Girdled","HardwoodControl",..: 1 1 4 4 3 3 4 4 1 1 ...
## $ plot : int 1 1 2 2 3 3 4 4 5 5 ...
## $ subplot : Factor w/ 2 levels "A","B": 2 1 1 2 2 1 1 2 1 2 ...
## $ sub_treatment : Factor w/ 3 levels "AntControl","AntDistControl",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Acer_rubrum : int 1 1 0 2 2 1 2 6 17 7 ...
## $ Acer_rubrum_BA : num 0 0 0 0 0 0 1.17 1.72 0 0 ...
## $ Aralia_nudicaulis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Betula_lenta : int 5 11 9 8 0 0 4 4 6 1 ...
## $ Betula_lenta_BA : num 5.24 6.08 8.81 6.53 0 ...
## $ Carex_spp_pct : int 0 0 1 1 0 0 10 20 5 0 ...
## $ Carya_sp : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Chimaphila_maculata : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Comptonia_peregrina : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Coptis_groenlandica : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Dendrolycopodium_dendroideum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Dendrolycopodium_obscurum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Epigaea_repens_pct : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Dennstaedtia_punctilobula_pct : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Erechtites_hieraciifolius : int 0 0 0 0 0 0 1 0 0 0 ...
## $ Fagus_grandifolia : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Fagus_grandifolia_BA : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Fraxinus_sp : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Gaultheria_procumbens : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Goodyera_repens : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Hamamelis_virginiana : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Lysimachia_borealis : int 0 5 0 0 0 0 0 0 0 0 ...
## $ Lysimachia_quadrifolia : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Maianthemum_canadense_pct : int 0 0 1 0 0 0 0 0 10 0 ...
## $ Medeola_virginiana : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Mitchella_repens_pct : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Monotropa_uniflora : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Oclemena_acuminata : int 0 0 0 0 0 0 0 0 0 1 ...
## $ Osmundastrum_cinnamomeum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Parathelypteris_noveboracensis: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Pinus_strobus : int 7 10 0 14 0 0 0 0 5 2 ...
## $ Pinus_strobus_BA : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Populus_sp : int 0 0 0 0 0 0 0 0 1 0 ...
## $ Pteridium.aquilinum_pct : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Prunus_sp : int 0 0 0 0 0 0 0 0 1 0 ...
## $ Quercus_alba : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Quercus_rubra : int 1 0 0 0 0 0 0 0 0 0 ...
## $ Rubus_flagellaris_pct : int 15 0 0 0 0 0 0 0 0 0 ...
## $ Rubus_sp : int 0 0 6 0 0 0 10 20 0 11 ...
## $ Spinulum_annotinum : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Tsuga_canadensis : int 0 1 1 0 0 0 0 0 2 0 ...
## $ Uvularia_sessilifolia : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Vaccinium_sp : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Viburnum_acerifolium : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Viburnum_sp : int 0 0 0 0 0 0 0 0 0 0 ...veg.scale <- scale(veg[,6:dim(veg)[[2]]])
use.veg<-data.frame(cbind(veg[,1:5],veg.scale))
uveg <- read.csv("understory-veg2014.csv")
uveg.scale <- scale(uveg[,6:dim(uveg)[[2]]])
use.uveg<-data.frame(cbind(uveg[,1:5],uveg.scale))
#pca
veg.pca <- prcomp(veg.scale[,5:45], scale=FALSE, center=FALSE)
veg.pca.out <- data.frame(cbind(veg[,1:5], veg.pca$x[,1:5]))
uveg.pca <-prcomp(uveg.scale[,5:26], scale=FALSE, center=FALSE)
uveg.pca.out <- data.frame(cbind(uveg[,1:5], uveg.pca$x[,1:5]))
#plots
#1. Screeplot suggests for uveg, PCs 1-2 are of interest
#2. biplots suggest 7BD is pulling PC1; and 1BX and 1BD are pulling out PC3.
par(mfrow=c(2,2))
screeplot(uveg.pca)
biplot(uveg.pca, choices=c(1,2), cex=c(1,.25))
biplot(veg.pca, choices=c(1,3), cex=c(1,.25))
biplot(veg.pca, choices=c(2,3), cex=c(1,.25))
#Anova suggests only canopy treatment affects PC2, no effects of treatments on PCs 1, 3, 4, and canopy*subtreatment interaction on PC5.
summary(aov(data=uveg.pca.out,PC1~treat*sub_treatment+ Error(block/treat/sub_treatment)))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 4.364 4.364
##
## Error: block:treat
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 3 37.34 12.446 2.764 0.213
## Residuals 3 13.51 4.503
##
## Error: block:treat:sub_treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## sub_treatment 2 0.826 0.4129 0.239 0.793
## treat:sub_treatment 6 12.933 2.1554 1.247 0.376
## Residuals 8 13.830 1.7287
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24 63.63 2.651summary(aov(data=uveg.pca.out,PC2~treat*sub_treatment+ Error(block/treat/sub_treatment)))##
## Error: block
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 5.568 5.568
##
## Error: block:treat
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 3 58.82 19.607 6.584 0.078 .
## Residuals 3 8.93 2.978
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: block:treat:sub_treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## sub_treatment 2 6.466 3.233 3.062 0.103
## treat:sub_treatment 6 6.027 1.004 0.951 0.510
## Residuals 8 8.448 1.056
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24 44.6 1.858#summary(aov(data=veg.pca.out,PC3~canopy*sub_treatment+ Error(block/canopy/sub_treatment)))
#summary(aov(data=veg.pca.out,PC4~canopy*sub_treatment+ Error(block/canopy/sub_treatment)))
#summary(aov(data=veg.pca.out,PC5~canopy*sub_treatment+ Error(block/canopy/sub_treatment)))
#biplot of PC2, PC5
par(mfrow=c(1,1))
biplot(veg.pca, choices=c(2,5), cex=c(1,.5))
#pretty plot
# visual elements to be used throughout
colors <- scale_color_manual(values = c("yellow","violet", "blue", "red" ))
shapes <- scale_shape_manual(values=c(0, 12, 15))
lines <- scale_linetype_manual(values = c("dashed", "dotted", "solid"), labels = c("Control", "Disturbance Control", "Ants Excluded"))
nolegend <- theme(legend.position = "none")
# corner<-theme(legend.position=c(.5,.9),legend.direction='horizontal')
legendtitle <- labs(color = "Canopy Treatment", linetype = "Ant Treatment")
veg.plot <- ggplot(veg.pca.out, aes(x=PC2, y=PC5)) +
geom_point(aes(colour = canopy, shape = sub_treatment),size=4) +
colors + shapes + theme + legendtitle
veg.plot
#what's loading?
veg.loadings <- veg.pca$rotation[,c(1,5)]
a.1 <- data.frame(veg.loadings[,1][abs(veg.loadings[,1])>.2])
a.1 <- cbind(rownames(a.1), a.1)
rownames(a.1) <- NULL
names(a.1) <- c("species", "PC2")
a.2 <- data.frame(veg.loadings[,2][abs(veg.loadings[,2])>.2])
a.2 <- cbind(rownames(a.2), a.2)
rownames(a.2) <- NULL
names(a.2) <- c("species", "PC5")
veg.big.loadings<- merge(a.1, a.2, all=TRUE)
veg.big.loadings## species PC2 PC5
## 1 Fagus_grandifolia 0.4164545 NA
## 2 Fagus_grandifolia_BA 0.4164545 NA
## 3 Quercus_alba 0.4164545 NA
## 4 Uvularia_sessilifolia 0.4164545 NA
## 5 Viburnum_acerifolium 0.3907133 NA
## 6 Carya_sp NA 0.2636838
## 7 Dendrolycopodium_dendroideum NA 0.2051748
## 8 Dendrolycopodium_obscurum NA 0.2636838
## 9 Dennstaedtia_punctilobula_pct NA -0.2153437
## 10 Gaultheria_procumbens NA 0.3170771
## 11 Lysimachia_borealis NA 0.2025959
## 12 Monotropa_uniflora NA -0.3898282
## 13 Pteridium.aquilinum_pct NA -0.2884622
## 14 Spinulum_annotinum NA -0.4089290
## 15 Vaccinium_sp NA 0.2752672
## 16 Viburnum_sp NA -0.2162726#Plot-level veg data
herb.shrub <- read.csv("hf106-03-shrub.csv")
x.hs <- aggregate(herb.shrub[,6:54], by=list(block=herb.shrub$block, treat=herb.shrub$treat, plot=herb.shrub$plot, year=herb.shrub$year), mean)
#pca
hs.pca <- prcomp(x.hs[,5:53], scale=TRUE, center=TRUE)
summary(hs.pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.9458 2.19850 2.13610 1.8783 1.55977 1.49750
## Proportion of Variance 0.1771 0.09864 0.09312 0.0720 0.04965 0.04577
## Cumulative Proportion 0.1771 0.27573 0.36886 0.4409 0.49050 0.53627
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 1.39081 1.32201 1.21723 1.17616 1.17354 1.11535
## Proportion of Variance 0.03948 0.03567 0.03024 0.02823 0.02811 0.02539
## Cumulative Proportion 0.57575 0.61141 0.64165 0.66988 0.69799 0.72338
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 1.08163 1.04934 1.03508 1.01775 0.98817 0.98374
## Proportion of Variance 0.02388 0.02247 0.02187 0.02114 0.01993 0.01975
## Cumulative Proportion 0.74725 0.76973 0.79159 0.81273 0.83266 0.85241
## PC19 PC20 PC21 PC22 PC23 PC24
## Standard deviation 0.91684 0.81904 0.80528 0.76521 0.75308 0.71117
## Proportion of Variance 0.01716 0.01369 0.01323 0.01195 0.01157 0.01032
## Cumulative Proportion 0.86956 0.88325 0.89649 0.90844 0.92001 0.93033
## PC25 PC26 PC27 PC28 PC29 PC30
## Standard deviation 0.6676 0.64097 0.57248 0.53171 0.51265 0.48561
## Proportion of Variance 0.0091 0.00838 0.00669 0.00577 0.00536 0.00481
## Cumulative Proportion 0.9394 0.94781 0.95450 0.96027 0.96563 0.97045
## PC31 PC32 PC33 PC34 PC35 PC36
## Standard deviation 0.46668 0.45795 0.43361 0.39298 0.35399 0.31912
## Proportion of Variance 0.00444 0.00428 0.00384 0.00315 0.00256 0.00208
## Cumulative Proportion 0.97489 0.97917 0.98301 0.98616 0.98872 0.99080
## PC37 PC38 PC39 PC40 PC41 PC42
## Standard deviation 0.2882 0.25894 0.24358 0.21828 0.20479 0.19480
## Proportion of Variance 0.0017 0.00137 0.00121 0.00097 0.00086 0.00077
## Cumulative Proportion 0.9925 0.99386 0.99507 0.99604 0.99690 0.99767
## PC43 PC44 PC45 PC46 PC47 PC48
## Standard deviation 0.17774 0.17037 0.15776 0.1204 0.09984 0.06348
## Proportion of Variance 0.00064 0.00059 0.00051 0.0003 0.00020 0.00008
## Cumulative Proportion 0.99832 0.99891 0.99942 0.9997 0.99992 1.00000
## PC49
## Standard deviation 1.308e-15
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00plot(hs.pca)
biplot(hs.pca)
hs.pca.out <- data.frame(cbind(x.hs[,1:4], hs.pca$x[,1:4]))
hs.pca.out$treat <- factor(hs.pca.out$treat, levels(hs.pca.out$treat)[c(3, 1, 4, 2)])
hs.loadings <- hs.pca$rotation[,c(1,2)]
a.1 <- data.frame(hs.loadings[,1][abs(hs.loadings[,1])>.2])
a.1 <- cbind(rownames(a.1), a.1)
rownames(a.1) <- NULL
names(a.1) <- c("species", "PC1")
a.2 <- data.frame(hs.loadings[,2][abs(hs.loadings[,2])>.2])
a.2 <- cbind(rownames(a.2), a.2)
rownames(a.2) <- NULL
names(a.2) <- c("species", "PC2")
hs.big.loadings<- merge(a.1, a.2, all=TRUE)
hs.big.loadings## species PC1 PC2
## 1 ARANUD -0.2872716 NA
## 2 ASTACU -0.2006033 NA
## 3 DENPUN -0.2514540 NA
## 4 LYCLUC -0.2178533 NA
## 5 LYCOBS -0.2970068 NA
## 6 MAICAN -0.2822933 NA
## 7 MITREP -0.2644930 NA
## 8 TRIBOR -0.2993972 NA
## 9 BERTHU NA -0.3162580
## 10 EPIREP NA -0.3366230
## 11 ILEVER NA -0.3233903
## 12 LONCAN NA -0.3165076
## 13 LYSQUA NA 0.2204926
## 14 OSMCIN NA -0.3520512
## 15 RHONUD NA -0.2124334#exploratory plots
ggplot(data=hs.pca.out, aes(x=PC1, y=PC2, colour=factor(year), group=block)) +
geom_point() +
scale_colour_grey(start=0.7, end=0.1) +
theme_bw() +
facet_grid(~treat)
ggplot(hs.pca.out[hs.pca.out$block=="Valley",], aes(x=PC1, y=PC2, colour=year)) +
geom_path(aes(group=block), colour="black", arrow=arrow(angle=15, length=unit(3, "mm"))) +
geom_path(data=hs.pca.out[hs.pca.out$block=="Ridge",], aes(group=block), colour="red", arrow=arrow(angle=15, length=unit(3, "mm"))) +
facet_grid(block~treat) +
theme_bw() +
#ggtitle("Both blocks")+
xlab("Principal Component 1") +
ylab("Principal Component 2")
#try to predict from 2014 ant plot data
antveg.2014 <- read.csv("veg2014-predict.csv")
x.veg2014 <- aggregate(antveg.2014[,6:54], by=list(block=antveg.2014$block, treat=antveg.2014$treat, plot=antveg.2014$plot, ants =antveg.2014$sub_treatment), mean)
veg2014.scale <- scale(x.veg2014[,5:dim(x.veg2014)[[2]]])
#a function for replacing NaNs with 0s
is.nan.data.frame <- function(x){
do.call(cbind, lapply(x, is.nan))
}
#run it
veg2014.scale[is.nan.data.frame(veg2014.scale)] <- 0
vs.p <- veg2014.scale
antveg.2014.p <- predict(hs.pca, vs.p)
av.p.out <- data.frame(cbind(x.veg2014[,1:4], antveg.2014.p[,1:4]))
av.p.out$treat <- factor(av.p.out$treat, levels(av.p.out$treat)[c(3, 1, 4, 2)])
#plots
hs.pca.out.old <- hs.pca.out
av.p.out.old <- av.p.out
#sqrt transform PC scores for clarity in plotting
hs.pca.out$PC1.t <- ifelse(hs.pca.out$PC1 >= 0, sqrt(hs.pca.out$PC1), -1*sqrt(abs(hs.pca.out$PC1)))## Warning in sqrt(hs.pca.out$PC1): NaNs producedhs.pca.out$PC2.t <- ifelse(hs.pca.out$PC2 >= 0, sqrt(hs.pca.out$PC2), -1*sqrt(abs(hs.pca.out$PC2)))## Warning in sqrt(hs.pca.out$PC2): NaNs producedav.p.out$PC1.t <- ifelse(av.p.out$PC1 >= 0, sqrt(av.p.out$PC1), -1*sqrt(abs(av.p.out$PC1)))## Warning in sqrt(av.p.out$PC1): NaNs producedav.p.out$PC2.t <- ifelse(av.p.out$PC2 >= 0, sqrt(av.p.out$PC2), -1*sqrt(abs(av.p.out$PC2)))## Warning in sqrt(av.p.out$PC2): NaNs producedlevels(hs.pca.out$treat) <- c("HC", "G", "L", "HD")
levels(av.p.out$treat) <- c("HC", "G", "L", "HD")
pdf("Figure 2- Veg PCA.pdf", height=5, width=6.1)
ggplot(hs.pca.out, aes(x=PC1.t, y=PC2.t, colour=factor(year))) +
#geom_path(arrow=arrow(angle=0,length=unit(5, "mm"), type="open")) +
geom_point(shape=19, size=2) +
scale_colour_grey(start=0.7, end=0.1) +
facet_grid(treat~block) +
geom_point(data=hs.pca.out[hs.pca.out$year==2003,], aes(x=PC1.t, y=PC2.t), colour="gray70", shape=19, size=3) +
geom_point(data=hs.pca.out[hs.pca.out$year==2014,], aes(x=PC1.t, y=PC2.t), colour="gray10", shape=19, size=3) +
geom_point(data=av.p.out, aes(x=PC1.t, y=PC2.t), shape=15, colour=c("violetred3","sienna3","darkblue"), size=3) +
theme +
theme(strip.text.y =element_text(size=8)) +
theme(legend.position="none") +
xlab(expression(paste("Understory vegetation (",sqrt(PC1),")"))) +
ylab(expression(paste("Understory vegetation (",sqrt(PC2),")")))
dev.off()## pdf
## 2#sapling data
main.sap <- read.csv("compare-all-saplings.csv")
main.sap.2013<- main.sap[main.sap$species=="BELE",]
main.sap.2013<-na.omit(main.sap.2013)
sap.lm <- aov(numperha~treat + location+plot/location, data=main.sap.2013)
summary(aov(sap.lm))## Df Sum Sq Mean Sq F value Pr(>F)
## treat 2 397877982 198938991 8.822 0.00619 **
## location 3 44207716 14735905 0.653 0.59879
## plot 1 223175932 223175932 9.897 0.01040 *
## location:plot 3 50157448 16719149 0.741 0.55128
## Residuals 10 225502059 22550206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov(sap.lm)## Call:
## aov(formula = sap.lm)
##
## Terms:
## treat location plot location:plot Residuals
## Sum of Squares 397877982 44207716 223175932 50157448 225502059
## Deg. of Freedom 2 3 1 3 10
##
## Residual standard error: 4748.706
## Estimated effects may be unbalancedxtabs(~numperha+plot+location, data=main.sap.2013)## , , location = antcontrol
##
## plot
## numperha 1 2 4 5 8
## 0 0 0 0 1 1
## 22 0 0 0 0 0
## 5000 0 0 0 0 0
## 7500 1 0 0 0 0
## 8750 0 0 0 0 0
## 8944 0 0 0 0 0
## 9011 0 0 0 0 0
## 9822 0 0 0 0 0
## 10000 0 1 1 0 0
## 15222 0 0 0 0 0
## 16250 0 0 0 0 0
## 20000 0 0 0 0 0
## 22500 0 0 0 0 0
##
## , , location = antDcontrol
##
## plot
## numperha 1 2 4 5 8
## 0 0 0 0 1 1
## 22 0 0 0 0 0
## 5000 1 0 0 0 0
## 7500 0 0 0 0 0
## 8750 0 0 1 0 0
## 8944 0 0 0 0 0
## 9011 0 0 0 0 0
## 9822 0 0 0 0 0
## 10000 0 0 0 0 0
## 15222 0 0 0 0 0
## 16250 0 0 0 0 0
## 20000 0 1 0 0 0
## 22500 0 0 0 0 0
##
## , , location = antExclosure
##
## plot
## numperha 1 2 4 5 8
## 0 0 0 0 1 1
## 22 0 0 0 0 0
## 5000 0 0 0 0 0
## 7500 0 0 1 0 0
## 8750 0 0 0 0 0
## 8944 0 0 0 0 0
## 9011 0 0 0 0 0
## 9822 0 0 0 0 0
## 10000 0 0 0 0 0
## 15222 0 0 0 0 0
## 16250 1 0 0 0 0
## 20000 0 0 0 0 0
## 22500 0 1 0 0 0
##
## , , location = center
##
## plot
## numperha 1 2 4 5 8
## 0 0 0 0 0 0
## 22 0 0 0 0 1
## 5000 0 0 0 0 0
## 7500 0 0 0 0 0
## 8750 0 0 0 0 0
## 8944 0 1 0 0 0
## 9011 0 0 1 0 0
## 9822 0 0 0 1 0
## 10000 0 0 0 0 0
## 15222 1 0 0 0 0
## 16250 0 0 0 0 0
## 20000 0 0 0 0 0
## 22500 0 0 0 0 0###############################################################################
#Light
###############################################################################
#Clean up the workspace to avoid strange problems with attach()
rm(list = ls())
#Read in file, get ready to go
simeslight.raw<-read.csv("hf107-01-light.csv")
simeslight.work <- simeslight.raw
simeslight.work$Year <- trunc(simeslight.work$Date/10000,0)
simeslight.agg <- aggregate(GSF~Year+Treatment+Block+Leaf, data=simeslight.work, mean)
simeslight.agg$adelgid <- ifelse(simeslight.agg$Year<2010, "pre-adelgid","post-adelgid")
simeslight.agg$adelgid <- as.factor(simeslight.agg$adelgid)
simeslight.agg$season <- ifelse(simeslight.agg$Leaf=="off", "Winter", "Summer")
simeslight.agg$season<-as.factor(simeslight.agg$season)
simeslight.agg$Treatment <- factor(simeslight.agg$Treatment, levels(simeslight.agg$Treatment)[c(3, 1, 4, 2)])
simeslight.agg$adelgid <- factor(simeslight.agg$adelgid, levels(simeslight.agg$adelgid)[c(2, 1)])
simeslight.agg$season <- factor(simeslight.agg$season, levels(simeslight.agg$season)[c(2, 1)])
simeslight.agg$block <- ifelse(simeslight.agg$Block==1, "Valley", "Ridge")
simeslight.agg$block <- as.factor(simeslight.agg$block)
head(simeslight.agg)## Year Treatment Block Leaf GSF adelgid season block
## 1 2004 girdle 1 off 0.0988 pre-adelgid Winter Valley
## 2 2005 girdle 1 off 0.0616 pre-adelgid Winter Valley
## 3 2006 girdle 1 off 0.1172 pre-adelgid Winter Valley
## 4 2007 girdle 1 off 0.4272 pre-adelgid Winter Valley
## 5 2008 girdle 1 off 0.2876 pre-adelgid Winter Valley
## 6 2009 girdle 1 off 0.3820 pre-adelgid Winter Valley#define some ggplot parameters
#Both of these are applied to each plot, eg: plot1 + colors + theme
colors <- scale_color_manual(values = c("blue", "yellow3", "red", "plum3"))
#initial Plot
#use data only from 2005 or later
simeslight.use <- simeslight.agg[simeslight.agg$Year>2004,]
#pdf("Simeslight.pdf", width=6)
ggplot(data=simeslight.use, aes(x=Year, y=GSF*100, group=adelgid, shape=block)) +
geom_point() +
facet_grid(Treatment~season) +
stat_smooth(aes(color=Treatment), method="lm", se=TRUE) +
colors +
ylab("Global site factor (%)") +
theme_bw()
#dev.off()
#Regression models
attach(simeslight.use)
#overall model:
lm.overall <- aov(GSF~Treatment+Error(block/Treatment)+Year)
#treatments differ; changes over years
#now let's look at each treatment/season, w/ and w/out breaks
#hemlock controls
simeslight.hc.winter <-subset(simeslight.use, Treatment=="hemlock.control" & season=="Winter")
lm.hc.winter <- lm(GSF~Year, data=simeslight.hc.winter)
summary(lm.hc.winter)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.hc.winter)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.083697 -0.016658 -0.007663 0.019690 0.080548
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -50.054729 6.173032 -8.109 2.02e-07 ***
## Year 0.025011 0.003072 8.142 1.91e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03946 on 18 degrees of freedom
## Multiple R-squared: 0.7865, Adjusted R-squared: 0.7746
## F-statistic: 66.29 on 1 and 18 DF, p-value: 1.906e-07hc.seg.winter <- segmented(lm.hc.winter, seg.Z = ~ Year, psi = 2010)
summary(hc.seg.winter)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.hc.winter, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2010.000 3.035
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -39.478532 18.250544 -2.163 0.0460 *
## Year 0.019742 0.009093 2.171 0.0453 *
## U1.Year 0.012341 0.012860 0.960 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04067 on 16 degrees of freedom
## Multiple R-Squared: 0.7984, Adjusted R-squared: 0.7606
##
## Convergence attained in 2 iterations with relative change -1.263614e-16slope(hc.seg.winter)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.01974 0.009093 2.171 0.0004646 0.03902
## slope2 0.03208 0.009093 3.528 0.0128100 0.05136simeslight.hc.summer <-subset(simeslight.use, Treatment=="hemlock.control" & season=="Summer")
lm.hc.summer <- lm(GSF~Year, data=simeslight.hc.summer)
summary(lm.hc.summer)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.hc.summer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.053209 -0.023141 -0.007906 0.011772 0.066059
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -26.01256 5.00434 -5.198 6.06e-05 ***
## Year 0.01301 0.00249 5.223 5.75e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03199 on 18 degrees of freedom
## Multiple R-squared: 0.6024, Adjusted R-squared: 0.5804
## F-statistic: 27.28 on 1 and 18 DF, p-value: 5.748e-05hc.seg.summer <- segmented(lm.hc.summer, seg.Z = ~ Year, psi = 2010)
summary(hc.seg.summer)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.hc.summer, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2013.000 1.305
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.389991 7.242927 -4.058 0.000914 ***
## Year 0.014689 0.003606 4.073 0.000885 ***
## U1.Year -0.024680 0.033247 -0.742 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03305 on 16 degrees of freedom
## Multiple R-Squared: 0.6228, Adjusted R-squared: 0.552
##
## Convergence attained in 2 iterations with relative change 0slope(hc.seg.summer)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.014690 0.003606 4.0730 0.007044 0.02233
## slope2 -0.009992 0.033050 -0.3023 -0.080060 0.06007#girdled plots
simeslight.g.winter <-subset(simeslight.use, Treatment=="girdle" & season=="Winter")
lm.g.winter <- lm(GSF~Year, data=simeslight.g.winter)
summary(lm.g.winter)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.g.winter)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11743 -0.06491 -0.01148 0.05035 0.17598
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -94.139089 13.354422 -7.049 1.41e-06 ***
## Year 0.047031 0.006646 7.077 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08536 on 18 degrees of freedom
## Multiple R-squared: 0.7356, Adjusted R-squared: 0.7209
## F-statistic: 50.08 on 1 and 18 DF, p-value: 1.341e-06g.seg.winter <- segmented(lm.g.winter, seg.Z = ~ Year, psi = 2010)
summary(g.seg.winter)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.g.winter, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2007.000 1.062
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -279.25335 144.63033 -1.931 0.0714 .
## Year 0.13930 0.07212 1.932 0.0713 .
## U1.Year -0.10796 0.07254 -1.488 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07212 on 16 degrees of freedom
## Multiple R-Squared: 0.8323, Adjusted R-squared: 0.8008
##
## Convergence attained in 3 iterations with relative change -5.695033e-16slope(g.seg.winter)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.13930 0.072120 1.932 -0.01358 0.29220
## slope2 0.03134 0.007869 3.983 0.01466 0.04802simeslight.g.summer <-subset(simeslight.use, Treatment=="girdle" & season=="Summer")
lm.g.summer <- lm(GSF~Year, data=simeslight.g.summer)
summary(lm.g.summer)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.g.summer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.07434 -0.03175 -0.01296 0.03503 0.10653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -24.670066 7.983361 -3.090 0.00631 **
## Year 0.012379 0.003973 3.116 0.00597 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05103 on 18 degrees of freedom
## Multiple R-squared: 0.3504, Adjusted R-squared: 0.3143
## F-statistic: 9.709 on 1 and 18 DF, p-value: 0.005967g.seg.summer <- segmented(lm.g.summer, seg.Z = ~ Year, psi = 2010)
summary(g.seg.summer)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.g.summer, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2008.0000 0.4325
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.135e+02 1.933e+01 -5.875 2.35e-05 ***
## Year 5.667e-02 9.632e-03 5.883 2.31e-05 ***
## U1.Year -6.285e-02 1.092e-02 -5.754 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03046 on 16 degrees of freedom
## Multiple R-Squared: 0.7943, Adjusted R-squared: 0.7557
##
## Convergence attained in 4 iterations with relative change 0slope(g.seg.summer)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.056670 0.009632 5.883 0.03625 0.077090
## slope2 -0.006183 0.005149 -1.201 -0.01710 0.004732#logged
simeslight.l.winter <-subset(simeslight.use, Treatment=="log" & season=="Winter")
lm.l.winter <- lm(GSF~Year, data=simeslight.l.winter)
summary(lm.l.winter)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.l.winter)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10151 -0.05011 -0.01205 0.04861 0.12698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.006931 11.191011 2.145 0.0458 *
## Year -0.011670 0.005569 -2.096 0.0505 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07154 on 18 degrees of freedom
## Multiple R-squared: 0.1961, Adjusted R-squared: 0.1515
## F-statistic: 4.391 on 1 and 18 DF, p-value: 0.05053l.seg.winter <- segmented(lm.l.winter, seg.Z = ~ Year, psi = 2010)
summary(l.seg.winter)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.l.winter, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2011.000 1.922
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.03947 24.43872 1.802 0.0904 .
## Year -0.02165 0.01217 -1.778 0.0943 .
## U1.Year 0.03349 0.02582 1.297 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07202 on 16 degrees of freedom
## Multiple R-Squared: 0.2757, Adjusted R-squared: 0.1399
##
## Convergence attained in 2 iterations with relative change 3.30457e-16slope(l.seg.winter)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -0.02165 0.01217 -1.7780 -0.04746 0.004157
## slope2 0.01184 0.02277 0.5198 -0.03644 0.060120simeslight.l.summer <-subset(simeslight.use, Treatment=="log" & season=="Summer")
lm.l.summer <- lm(GSF~Year, data=simeslight.l.summer)
summary(lm.l.summer)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.l.summer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.112605 -0.025201 0.002078 0.045549 0.064730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.393287 7.701360 9.14 3.50e-08 ***
## Year -0.034875 0.003832 -9.10 3.74e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04923 on 18 degrees of freedom
## Multiple R-squared: 0.8214, Adjusted R-squared: 0.8115
## F-statistic: 82.81 on 1 and 18 DF, p-value: 3.736e-08l.seg.summer <- segmented(lm.l.summer, seg.Z = ~ Year, psi = 2010)
summary(l.seg.summer)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.l.summer, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2008.000 1.481
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.822945 45.686450 0.390 0.702
## Year -0.008676 0.022775 -0.381 0.708
## U1.Year -0.036266 0.023574 -1.538 NA
##
## Residual standard error: 0.04555 on 16 degrees of freedom
## Multiple R-Squared: 0.8641, Adjusted R-squared: 0.8386
##
## Convergence attained in 18 iterations with relative change 0.0005927034slope(l.seg.summer)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -0.008676 0.022770 -0.3809 -0.05696 0.03960
## slope2 -0.044940 0.006087 -7.3830 -0.05785 -0.03204#hardwood control
simeslight.hd.winter <-subset(simeslight.use, Treatment=="hardwood.control" & season=="Winter")
lm.hd.winter <- lm(GSF~Year, data=simeslight.hd.winter)
summary(lm.hd.winter)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.hd.winter)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.22164 -0.07734 -0.01536 0.11700 0.20900
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.293725 19.478417 -0.682 0.504
## Year 0.006851 0.009693 0.707 0.489
##
## Residual standard error: 0.1245 on 18 degrees of freedom
## Multiple R-squared: 0.027, Adjusted R-squared: -0.02705
## F-statistic: 0.4996 on 1 and 18 DF, p-value: 0.4887hd.seg.winter <- segmented(lm.hd.winter, seg.Z = ~ Year, psi = 2010)
summary(hd.seg.winter)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.hd.winter, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2011.000 1.423
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.01918 41.02229 0.707 0.489
## Year -0.01423 0.02043 -0.696 0.496
## U1.Year 0.07525 0.04335 1.736 NA
##
## Residual standard error: 0.1209 on 16 degrees of freedom
## Multiple R-Squared: 0.1847, Adjusted R-squared: 0.03178
##
## Convergence attained in 2 iterations with relative change 1.1819e-16slope(hd.seg.winter)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -0.01423 0.02043 -0.6963 -0.05755 0.02909
## slope2 0.06102 0.03823 1.5960 -0.02002 0.14210simeslight.hd.summer <-subset(simeslight.use, Treatment=="hardwood.control" & season=="Summer")
lm.hd.summer <- lm(GSF~Year, data=simeslight.hd.summer)
summary(lm.hd.summer)##
## Call:
## lm(formula = GSF ~ Year, data = simeslight.hd.summer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.051130 -0.022700 -0.006482 0.013727 0.074937
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.31586 5.92795 -3.427 0.00301 **
## Year 0.01017 0.00295 3.448 0.00287 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03789 on 18 degrees of freedom
## Multiple R-squared: 0.3977, Adjusted R-squared: 0.3643
## F-statistic: 11.89 on 1 and 18 DF, p-value: 0.002871hd.seg.summer <- segmented(lm.hd.summer, seg.Z = ~ Year, psi = 2010)
summary(hd.seg.summer)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm.hd.summer, seg.Z = ~Year, psi = 2010)
##
## Estimated Break-Point(s):
## Est. St.Err
## 2012.0000 0.9891
##
## t value for the gap-variable(s) V: 0
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.807336 7.890892 -3.777 0.00165 **
## Year 0.014899 0.003929 3.792 0.00160 **
## U1.Year -0.045134 0.036221 -1.246 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03601 on 16 degrees of freedom
## Multiple R-Squared: 0.5166, Adjusted R-squared: 0.426
##
## Convergence attained in 3 iterations with relative change 0slope(hd.seg.summer)## $Year
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.01490 0.003929 3.7920 0.00657 0.02323
## slope2 -0.03023 0.036010 -0.8397 -0.10660 0.04610#write out graphing file for additional annotation
write.csv(simeslight.use, "simeslightuse.csv")
#Plot with identified breaks
#use data only from 2005 or later
simeslight.plot <- read.csv("simeslightuse-annotated.csv")
simeslight.plot$Treatment <- factor(simeslight.plot$Treatment, levels(simeslight.plot$Treatment)[c(3, 1, 4, 2)])
simeslight.plot$season <- factor(simeslight.plot$season, levels(simeslight.plot$season)[c(2, 1)])
levels(simeslight.plot$Treatment) <- c("HC", "G", "L", "HD")
pdf("Figure 1- Simeslight.pdf", height=5, width=6.1)
ggplot(data=simeslight.plot, aes(x=Year, y=GSF*100, group=seg.break, shape=block)) +
geom_point() +
facet_grid(Treatment~season) +
stat_smooth(aes(color=Treatment), method="lm", se=TRUE) +
colors +
ylab("Global site factor (%)") +
geom_vline(xintercept=2009.5, colour="black", linetype="longdash") +
scale_x_continuous(breaks=seq(2006,2014,2)) +
theme_bw() +
theme(legend.position="none",
text = element_text(size = 8),
legend.key = element_blank(),
plot.title = element_text(hjust = 0.05),
plot.margin = unit(c(2, 2, 2, 2), "points"),
legend.margin = unit(-13, "points"))
dev.off()## pdf
## 2