#required libraries library(nlme) library(multcomp) library(ggplot2) library(grid) library(gridExtra) library(lattice) library(plyr) library(MASS) library(car) #required files: #"Initial_Urchin_data_raw.csv" #"Final_urchin_data_raw.csv" #"Feeding_Data_raw_7Sep2014.csv" #"urchinControls.csv" #"CHN.csv" #read in raw urchin data and perform calculations, put in new dataframe initial<-read.csv("Initial_Urchin_data_raw.csv") final<-read.csv("Final_urchin_data_raw.csv") expData<-data.frame( BOX=initial$URCHIN, TREATMENT=rep(c("MAPY","PTCA","R","CC","POLY"),times=7), #mean of 5 initial TD measurements INITIAL_MEAN_TD=apply(initial[,3:7],MARGIN=1,FUN=mean), #mean of 5 final TD measurements FINAL_MEAN_TD=apply(final[,2:6],MARGIN=1,FUN=mean)) #test change expData$TEST_CHANGE<-round(with(expData,FINAL_MEAN_TD - INITIAL_MEAN_TD),digits=1) #wet wt change expData$WET_WT_CHANGE<-round(final$Wet_urchin_wt_no_boat - initial$WET_WT,digits=2) #gonad dry wt expData$FINAL_GONAD_DRY_WT<-round(final$Dry_gonad_and_boat_wt - final$Gonad_boat_wt,digits=2) #change in jaw length expData$JAW_CHANGE<-final$Fluorescent_Band_Length #read consumption data into dataframe consumptionData<-read.csv("Feeding_Data_raw_7Sep2014.csv", na.strings=".") #create function to calculate the sample size within each treatment (does not count NA's) calcN<-function(x) {sum(!is.na(x))} #reshape consumption data consumptionData$added<-as.POSIXct(paste(consumptionData$ADD_DATE, consumptionData$ADD_TIME, sep=" "), format="%m/%d/%y %H:%M") consumptionData$removed<-as.POSIXct(paste(consumptionData$REM_DATE, consumptionData$REM_TIME, sep=" "), format="%m/%d/%y %H:%M") consumptionData$days<-as.numeric(with(consumptionData, removed-added)) consData<-ddply(consumptionData, .(TRIAL, BOX), function(adf){ trt<-adf$SP_CODE[1] if(nrow(adf)>1) trt<-"POLY" data.frame(TRT=trt, days=adf$days[1], delta=sum(adf$ADD_WT-adf$REM_WT), ADD_WT=sum(adf$ADD_WT), RM_WT=sum(adf$REM_WT)) }) consData$rate<-consData$delta/consData$days #calculate sum consumed and mean cons rate, append to expData dataframe sumCons<-ddply(consData,"BOX",function(x){ data.frame( SUM_CONSUMED=sum(x$delta,na.rm=TRUE), CONSUMPTION_RATE=mean(x$rate,na.rm=TRUE)) }) expData$SUM_CONSUMED<-sumCons$SUM_CONSUMED expData$CONSUMPTION_RATE<-sumCons$CONSUMPTION_RATE #subset consumptionData into mixed diet treatments only (mixed diet BOX #'s are multiples of 5) mixedData<-consumptionData[consumptionData$BOX %% 5 == 0,] mixedData$delta<-mixedData$ADD_WT - mixedData$REM_WT mixedData$rate<-mixedData$delta/mixedData$days #reshape to get species-specific consumption rate table mixed.summary<-ddply(mixedData,.(BOX,SP_CODE),function(x){ data.frame(CONSUMPTION_RATE=mean(x$rate,na.rm=TRUE)) }) #############fit linear models (not including consumption--see below) LMtestChange <- lm(formula = TEST_CHANGE ~ TREATMENT, data = expData, na.action = na.omit) LMwetWtChange <- lm(formula = WET_WT_CHANGE ~ TREATMENT, data = expData, na.action = na.omit) LMgonad <- lm(formula = FINAL_GONAD_DRY_WT ~ TREATMENT, data = expData, na.action = na.omit) LMjaw <- lm(formula = JAW_CHANGE ~ TREATMENT, data = expData, na.action = na.omit) #fit linear mixed models with consumption data LMconsumption <- lm(formula = CONSUMPTION_RATE ~ TREATMENT, data = expData, na.action = na.omit) mixedLME<-lme(CONSUMPTION_RATE ~ SP_CODE, random=~1|BOX, data=mixed.summary,na.action=na.omit) anova(mixedLME) #check model assumptions mixedLME$data[["residuals"]]<-mixedLME$residuals[1:28] #plot histogram of residuals of all models h1<-as.list( hist(LMtestChange$residuals,main="Test Change"), hist(LMwetWtChange$residuals,main="Wet Weight Change"), hist(LMgonad$residuals,main="Gonad Weight"), hist(LMjaw$residuals,main="Change in Jaw Length"), hist(LMconsumption$residuals,main="Consumption Rate"), hist(mixedLME$data[["residuals"]]) ) #plot normal q-q plot for residuals of all models q1<-as.list(qqnorm(studres(LMtestChange),main="Test Change"), qqnorm(studres(LMwetWtChange),main="Wet Weight Change"), qqnorm(studres(LMgonad),main="Gonad Weight"), qqnorm(studres(LMjaw),main="Change in Jaw Length"), qqnorm(studres(LMconsumption),main="Consumption Rate") ) s1<-as.list( spreadLevelPlot(LMtestChange,main="Test Change"), spreadLevelPlot(LMwetWtChange,main="Wet Weight Change"), spreadLevelPlot(LMgonad,main="Gonad Weight"), qplot(LMjaw$model[["TREATMENT"]],LMjaw$residuals,main="Jaw Growth"), spreadLevelPlot(LMconsumption,main="Consumption Rate"), qplot(mixedLME$data[["SP_CODE"]],mixedLME$data[["residuals"]],main="Consumption w/in Mixture") ) #shapiro-wilk test on residuals of all models stest<-shapiro.test(LMtestChange$residuals) sWW<-shapiro.test(LMwetWtChange$residuals) sgonad<-shapiro.test(LMgonad$residuals) sjaw<-shapiro.test(LMjaw$residuals) scons<-shapiro.test(LMconsumption$residuals) smixedLME<-shapiro.test(mixedLME$residuals) #barlett test for all models LMtestChange$model[[3]]<-as.numeric(LMtestChange$residuals) LMwetWtChange$model[[3]]<-as.numeric(LMwetWtChange$residuals) LMgonad$model[[3]]<-as.numeric(LMgonad$residuals) LMjaw$model[[3]]<-as.numeric(LMjaw$residuals) LMconsumption$model[[3]]<-as.numeric(LMconsumption$residuals) mixedLME$data["residuals"]<-as.numeric(mixedLME$residuals[1:28]) btest<-bartlett.test(LMtestChange$model[[3]]~LMtestChange$model[[2]]) bWW<-bartlett.test(LMwetWtChange$model[[3]]~LMwetWtChange$model[[2]]) bgonad<-bartlett.test(LMgonad$model[[3]]~LMgonad$model[[2]]) bjaw<-bartlett.test(LMjaw$model[[3]]~LMjaw$model[[2]]) bconsumption<-bartlett.test(LMconsumption$model[[3]]~LMconsumption$model[[2]]) #log10 transformation for wet weight, test change, and mixed models expData["LOG_WW"]<-log(expData[[6]],base=10) expData["LOG_CONS_RATE"]<-log(expData["CONSUMPTION_RATE"],base=10) mixed.summary[["LOG_CONS"]]<-log(mixed.summary$CONSUMPTION_RATE,base=10) expData$TREATMENT<-as.factor(expData$TREATMENT) LMwwLog<-lm(LOG_WW~TREATMENT,data=expData,na.action=na.omit) LMconsLog<-lm(LOG_CONS_RATE~TREATMENT,data=expData,na.action=na.omit) mixedLMElog<-lme(LOG_CONS ~ SP_CODE, random=~1|BOX, data=mixed.summary,na.action=na.omit) summary(glht(LMwwLog,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) summary(glht(LMconsLog,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) summary(glht(mixedLMElog,linfct=mcp('SP_CODE'="Tukey"),test=adjusted("fdr"))) #now test assumptions again LMwwLog$model[[3]]<-as.numeric(LMwwLog$residuals) LMconsLog$model[[3]]<-as.numeric(LMconsLog$residuals) mixedLMElog$data[["residuals"]]<-mixedLMElog$residuals[1:28] swwLog<-shapiro.test(LMwwLog$residuals) sconsLog<-shapiro.test(LMconsLog$residuals) smixedLog<-shapiro.test(mixedLMElog$data[["residuals"]]) bwwLog<-bartlett.test(LMwwLog$model[[3]]~LMwwLog$model[[2]]) bconsLog<-bartlett.test(LMconsLog$model[[3]]~LMconsLog$model[[2]]) bmixedLog<-bartlett.test(mixedLMElog$data[["residuals"]]~mixedLMElog$data[["SP_CODE"]]) #plots for model assumptions h2<-as.list( hist(LMwwLog$residuals,main="Wet Weight Change"), hist(LMconsLog$residuals,main="Consumption Rate"), hist(mixedLMElog$data[["residuals"]]) ) q2<-as.list( qqnorm(studres(LMwwLog),main="Wet Weight Change"), qqnorm(studres(LMconsLog),main="Consumption Rate"), qqnorm(mixedLMElog,main="Consumption w/in Mixture") ) s2<-as.list( spreadLevelPlot(LMwwLog), spreadLevelPlot(LMconsLog), qplot(mixedLMElog$data[["SP_CODE"]],mixedLMElog$data[["residuals"]]) ) #0 intercept for anova table purposes LMwwLog0<-lm(LOG_WW~TREATMENT +0,data=expData,na.action=na.omit) LMconsLog0<-lm(LOG_CONS_RATE~TREATMENT + 0,data=expData,na.action=na.omit) mixedLMElog0<-lme(LOG_CONS ~ SP_CODE + 0, random=~1|BOX, data=mixed.summary,na.action=na.omit) LMtestChange0<-lm(TEST_CHANGE~TREATMENT + 0,data=expData,na.action=na.omit) LMgonad0<-lm(FINAL_GONAD_DRY_WT~TREATMENT + 0,data=expData,na.action=na.omit) LMjaw0<-lm(JAW_CHANGE~TREATMENT + 0,data=expData,na.action=na.omit) LMconsLog0<-lm(LOG_CONS_RATE~TREATMENT + 0,data=expData,na.action=na.omit) #write ANOVA tables for one-way ANOVAs write.table(rbind (anova(LMtestChange), anova(LMwwLog), anova(LMgonad), anova(LMjaw), anova(LMconsLog) ),"Table1.csv", sep=",", row.names=T) #write ANOVA table for mixed LME write.table(anova(mixedLMElog0), "Table2.csv", sep=",", row.names=T) #post hoc Tukey tests summary(glht(LMtestChange,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) summary(glht(LMwwLog,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) summary(glht(LMgonad,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) summary(glht(LMjaw,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) summary(glht(mixedLMElog,linfct=mcp('SP_CODE'="Tukey"),test=adjusted("fdr"))) summary(glht(LMconsLog,linfct=mcp('TREATMENT'="Tukey"),test=univariate())) ##########figures mBarPlotLog<-function(x) { y<-as.data.frame(summary(x)$coefficients) rownames(y)<-NULL z<-data.frame( UL=as.numeric(10^(y[,"Estimate"]+y[,"Std. Error"])), LL=as.numeric(10^(y[,"Estimate"]-y[,"Std. Error"])), mean=as.numeric(10^(y[,"Estimate"])), treatment=c("Chondracanthus","Macrocystis","Mixture","Pterygophora","Rhodymenia") ) return(z) } mBarPlot<-function(x) { y<-summary(x)$coefficients rownames(y)<-NULL z<-data.frame( UL=as.numeric(y[,"Estimate"]+y[,"Std. Error"]), LL=as.numeric(y[,"Estimate"]-y[,"Std. Error"]), mean=as.numeric(y[,"Estimate"]), treatment=c("Chondracanthus","Macrocystis","Mixture","Pterygophora","Rhodymenia") ) return(z) } #rename treatment levels expData$TREATMENT<-revalue(expData$TREATMENT, c("CC"="Chondracanthus", "MAPY"="Macrocystis","POLY"="Mixture","PTCA"="Pterygophora","R"="Rhodymenia")) #create summary tables for response variables testChange.summary<-mBarPlot(LMtestChange0) wetWeight.summary<-mBarPlotLog(LMwwLog0) gonad.summary<-mBarPlot(LMgonad0) jaw.summary<-mBarPlot(LMjaw0) consumption.summary<-mBarPlotLog(LMconsLog0) ##create individual plots #test change p1<-ggplot(testChange.summary, aes(x = treatment, y = mean)) +  geom_bar(position = position_dodge(), stat="identity", fill="gray50",width=0.5,colour="black",size=0.3) + geom_linerange(aes(ymin=LL, ymax=UL),size=0.3) + theme_bw() + theme( panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("Change in Test Diameter (mm)\n") + geom_text(aes(x=1:5, y=8.25, label=c("A", "A", "AB", "BC", "BC")),size=3.5) + ggtitle("(a)") #wet weight change p2<-ggplot(wetWeight.summary, aes(x = treatment, y = mean)) +  geom_bar(position = position_dodge(), stat="identity", fill="gray50",width=0.5,colour="black",size=0.3) + geom_linerange(aes(ymin=LL, ymax=UL),size=0.3) + theme_bw() + theme(panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("Change in Wet Weight (g)\n") + geom_text(aes(x=1:5, y=14.1, label=c("A", "A", "A", "B", "B")),size=3.5)+ ggtitle("(b)") #gonad mass p3<-ggplot(gonad.summary, aes(x = treatment, y = mean)) +  geom_bar(position = position_dodge(), stat="identity", fill="gray50",width=0.5,colour="black",size=0.3) + geom_linerange(aes(ymin=LL, ymax=UL),size=0.3) + theme_bw() + theme(panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("Gonad Dry Weight (g)\n") + geom_text(aes(x=1:5, y=1.1, label=c("A", "A", "A", "A", "B")),size=3.5)+ ggtitle("(c)") #jaw growth p4<-ggplot(jaw.summary, aes(x = treatment, y = mean)) +  geom_bar(position = position_dodge(), stat="identity", fill="gray50",width=0.5,colour="black",size=0.3) + geom_linerange(aes(ymin=LL, ymax=UL),size=0.3) + theme_bw() + theme(panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("Change in Jaw Length (mm)\n") + geom_text(aes(x=1:5, y=0.6, label=c("A", "AB", "A", "AB", "B")),size=3.5) + ggtitle("(d)") #consumption rate p5<-ggplot(consumption.summary, aes(x = treatment, y = mean)) +  geom_bar(position = position_dodge(), stat="identity", fill="gray71",width=0.5,colour="black",size=0.3) + geom_linerange(aes(ymin=LL, ymax=UL),size=0.3) + theme_bw() + theme(panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("Consumption Rate\n (g per 24 h)\n") + geom_text(aes(x=1:5, y=2.75, label=c("A", "A", "A", "B", "B")),size=3.5) + ggtitle("(b)") #put 4 "performance measures" plots in one figure postscript('Figure1.eps',paper="special",horizontal=FALSE,onefile=FALSE,height=11,width=11) grid.arrange(p1,p2,p3,p4,widths=unit(0.5,"npc"),heights=unit(0.5,"npc")) dev.off() m<-summary(mixedLMElog0)$tTable ##########################plot consumption within mixed diet #create summary data table for mixed diet consumption mixedPlot.summary<-data.frame(species=c("CC"="Chondracanthus", "MAPY"="Macrocystis","PTCA"="Pterygophora","R"="Rhodymenia"), mean=10^(m[,"Value"]), UL=10^(m[,"Value"]+m[,"Std.Error"]), LL=10^(m[,"Value"]-m[,"Std.Error"]) ) #rename species/treatment levels # mixedPlot.summary$species<-c("CC"="Chondracanthus", "MAPY"="Macrocystis","PTCA"="Pterygophora","R"="Rhodymenia") mixPlot<-ggplot(mixedPlot.summary, aes(x=species, y=mean)) + geom_bar(position = position_dodge(), stat="identity", fill="gray71",width=0.5,weight=1,colour="black",size=0.3) + geom_linerange(aes(ymin=LL, ymax=UL),size=0.3) + theme_bw() + theme(panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("Consumption Rate\n (g per 24 h)\n") + geom_text(aes(x=1:4, y=1.35, label=c("A", "B", "C", "C")),size=3.5) + ggtitle("(a)") #put 2 consumption plots in one figure postscript('Figure2.eps',paper="special",horizontal=FALSE, onefile=FALSE, height=11,width=11) grid.arrange(mixPlot,p5,nrow=1,widths=unit(0.5,"npc"),heights=unit(0.5,"npc")) dev.off() ############data from dissected urchins from collection site (control) #read into data frame controlDataRaw<-read.csv("urchinControls.csv") #calculations controlData<-data.frame( TREATMENT=controlDataRaw$Treatment, FINAL_MEAN_TD=apply(controlDataRaw[,5:9],MARGIN=1,FUN=mean), FINAL_GONAD_DRY_WT=round(with(controlDataRaw,Gonad_and_boat_dry_wt - Gonad_boat_wt),digits=2) ) #find range of initial test diameter for experimental urchins expTDrange<-range(expData$INITIAL_MEAN_TD,na.rm=TRUE) #subset control urchins to restrict within this range controlDataSub<-subset(controlData,controlData$FINAL_MEAN_TD > expTDrange[1] & controlData$FINAL_MEAN_TD < expTDrange[2]) #get mean and se of test diameter for these urchins mean(controlDataSub$FINAL_MEAN_TD,na.rm=TRUE) sd(controlDataSub$FINAL_MEAN_TD,na.rm=TRUE)/sqrt(sum(!is.na(controlDataSub$FINAL_MEAN_TD))) #set negative gonad weights to zero u<-controlDataSub$FINAL_GONAD_DRY_WT < 0 controlDataSub$FINAL_GONAD_DRY_WT[which(u)] <- 0 #get mean and se of gonad dry weights mean(controlDataSub$FINAL_GONAD_DRY_WT,na.rm=TRUE) sd(controlDataSub$FINAL_GONAD_DRY_WT,na.rm=TRUE)/sqrt(sum(!is.na(controlDataSub$FINAL_GONAD_DRY_WT))) #convert from factor to character, otherwise can't append() controlDataSub$TREATMENT<-as.character(controlDataSub$TREATMENT) expData$TREATMENT<-as.character(expData$TREATMENT) #create new dataframe with control urchin and Rhodymenia urchin gonad weights controlDF<-data.frame( Type=append(controlDataSub$TREATMENT,expData$TREATMENT[expData$TREATMENT=="Rhodymenia"]), FINAL_MEAN_TD=append(controlDataSub$FINAL_MEAN_TD,expData$FINAL_MEAN_TD[expData$TREATMENT=="Rhodymenia"]), FINAL_GONAD_DRY_WT=append(controlDataSub$FINAL_GONAD_DRY_WT,expData$FINAL_GONAD_DRY_WT[expData$TREATMENT=="Rhodymenia"]) ) #remove NA's controlDF1<-controlDF[!is.na(controlDF$FINAL_GONAD_DRY_WT),] #look at distributions, both normal shapiro.test(controlDF1$FINAL_GONAD_DRY_WT[controlDF1$Type=="Control"]) shapiro.test(controlDF1$FINAL_GONAD_DRY_WT[controlDF1$Type=="Rhodymenia"]) #Welch's t-test t.test(controlDF1$FINAL_GONAD_DRY_WT~controlDF1$Type) ##########CN Data #read into data frame chnData<-read.csv("CHN.csv") #convert to date format and order chnData$Date<-as.POSIXct(chnData$Date,format="%m/%d/%y") cnDateSp<-ddply(chnData,.(Date,Species),summarize,meanCN = mean(C.N)) chnData$Date<-strftime(chnData$Date,"%d-%b-%Y") cnDateSp$Date<-strftime(cnDateSp$Date,"%d-%b-%Y") chnData$Date<-as.factor(chnData$Date) cnDateSp$Date<-as.factor(cnDateSp$Date) chnData$Date<-relevel(chnData$Date,"12-Dec-2010","10-Jan-2011","24-Jan-2011") cnDateSp$Date<-relevel(cnDateSp$Date,"12-Dec-2010","10-Jan-2011","24-Jan-2011") #create function to make summary data frames (mean,n,sd,sem); calcN function defined above chnSummary<-function(param) { chnData.summary<-data.frame( mean=tapply(param,chnData$Species,mean,na.rm=TRUE), n=tapply(param,chnData$Species,calcN), sd=tapply(param,chnData$Species,sd,na.rm=TRUE), min=tapply(param,chnData$Species,min,na.rm=TRUE), max=tapply(param,chnData$Species,max,na.rm=TRUE)) chnData.summary$sem<-chnData.summary$sd/sqrt(chnData.summary$n) return(chnData.summary) } #create summary tables C.NSummary<-chnSummary(chnData$C.N) CSummary<-chnSummary(chnData$Wt_Percent_C) NSummary<-chnSummary(chnData$Wt_Percent_N) #test whether C:N ratio different across treatment LMchn1<-lme(C.N ~ Species,random=~1|Date,data=chnData) #post hoc LMchn1GLHT<-glht(LMchn1,linfct=mcp('Species'="Tukey"),test=adjusted("fdr")) #remove random effect term LMchn2<-lm(C.N ~ Species,data=chnData) #post hoc LMchn2GLHT<-glht(LMchn2,linfct=mcp('Species'="Tukey"),test=univariate()) #add mean C:N value for each species to expData dataframe--Mixture treatment is NA tempCN<-as.numeric(rep(c(C.NSummary$mean[[2]],C.NSummary$mean[[3]],C.NSummary$mean[[4]],C.NSummary$mean[[1]],NA),times=7)) expData$C.N<-round(tempCN,2) #test whether C:N correlated to diet treatment-- LMcnTest<-lm(TEST_CHANGE ~ C.N,data=expData) LMcnWW<-lm(WET_WT_CHANGE ~ C.N,data=expData) LMcnJaw<-lm(JAW_CHANGE ~ C.N,data=expData) LMcnGonad<-lm(FINAL_GONAD_DRY_WT ~ C.N,data=expData) #significant ##plot C:N vs. Date*Species postscript("Figure3.eps",height=4,width=5.5,horizontal=FALSE,onefile=FALSE,paper="special",pointsize=9) print(chn1<-ggplot(chnData,aes(Date,C.N)) + geom_point(data=chnData,aes(Date,C.N),size=1) + geom_path(data=cnDateSp,aes(x=Date,y=meanCN,group=Species))+ theme_bw() + theme( panel.grid.major = element_blank(), axis.title.x = element_text(family="Helvetica",size=10), axis.title.y = element_text(family="Helvetica",size=10)) + xlab("\nDate of Collection") + ylab("C : N\n") + geom_text(aes(x=3.3,y=c(C.NSummary$mean[[3]],C.NSummary$mean[[2]],C.NSummary$mean[[1]],C.NSummary$mean[[4]]), label=c("Pterygophora", "Macrocystis", "Chondracanthus", "Rhodymenia"),fontface="italic"),size=2.8) + theme(legend.position="none") ) dev.off() #####plot C:N vs. Species #create variable w/ treatment names C.NSummary$species<-rownames(C.NSummary) #rename treatment levels C.NSummary$species<-revalue(C.NSummary$species, c("CC"="Chondracanthus", "MAPY"="Macrocystis","PTCA"="Pterygophora","R"="Rhodymenia")) pdf('cn.pdf') print(CNPlot<-ggplot(C.NSummary, aes(x = species, y = mean)) +  geom_bar(position = position_dodge(), stat="identity", fill="black",width=0.5,colour="black",size=0.3) + geom_linerange(aes(ymin=mean-sem, ymax=mean+sem),size=0.3) + theme_bw() + theme( panel.grid.major = element_blank(), axis.text.x = element_text(face="italic",angle=45,hjust=.99,size=11), axis.title.y = element_text(family="Helvetica",size=11)) + xlab("") + ylab("C : N\n") + geom_text(aes(x=1:4, y=16.75, label=c("A", "B", "C", "D")),size=3.5) ) dev.off() #########calculate mean +/- se of initial test diameter of experimental urchins mean(expData$INITIAL_MEAN_TD,na.rm=TRUE) sd(expData$INITIAL_MEAN_TD,na.rm=TRUE)/sqrt(sum(!is.na(expData$INITIAL_MEAN_TD))) ######calculate the proportion of time Rhodymenia absent timeSummary<-tapply(consData$days,consData$TRIAL,mean,na.rm=TRUE) timeSummary[6]/sum(timeSummary[1:9]) #proportion of time Rhodymenia absent #calculate performance by treatment as a proportion of performance in Rhodymenia trt testChange.summary$RDiff<-with(testChange.summary, mean/mean[treatment=="Rhodymenia"]) wetWeight.summary$RDiff<-with(wetWeight.summary, mean/mean[treatment=="Rhodymenia"]) gonad.summary$RDiff<-with(gonad.summary, mean/mean[treatment=="Rhodymenia"]) jaw.summary$RDiff<-with(jaw.summary, mean/mean[treatment=="Rhodymenia"]) ##########calculate mean algal weight supplied to urchins #mean weight supplied to monospecific diets mean(consData$ADD_WT,na.rm=TRUE) #mean weight of each individual algal sp given in mixed diet mean(mixedData$ADD_WT,na.rm=TRUE) #calculate time urchins exposed to algae sum(timeSummary) ##calculate duration of experiment consumptionData$removed[length(consumptionData$removed)]-consumptionData$added[1]