R Markdown
## d13C by site #, stdev
p1 <- ggplot(data=data, mapping=aes(x=factor(data$Location, levels=c("Bauza\nIsland", "First\nArm", "Espinosa\nPoint", "Crooked\nArm", "Campbell's\nKingdom", "Hall\nArm", "Deep\nCove")), y=data$d13C, fill=factor(data$Depth), group=factor(data$Depth)))
p1 <- p1 + geom_bar(stat="identity", position=position_dodge(0.9))
#p1 <- p1 + geom_errorbar(aes(ymin=d13C-stdevC, ymax=d13C+stdevC), width=.25, position=position_dodge(.9))
p1 <- p1 + theme_few() + scale_fill_brewer(palette="Dark2", name="Depth (m)")
p1 <- p1 + xlab("Site Name") + ylab("d13C (per mil)")
p1

# d13C line graph
p1.1 <- ggplot(data=data, mapping=aes(x=factor(data$Location, levels=c("Bauza\nIsland", "First\nArm", "Espinosa\nPoint", "Crooked\nArm", "Campbell's\nKingdom", "Hall\nArm", "Deep\nCove")), y=data$d13C, color=factor(data$Depth), group=factor(data$Depth)))
p1.1 <- p1.1 + geom_line(stat="identity", size=1)
p1.1 <- p1.1 + theme_few() + scale_color_brewer(palette="Dark2", name="Depth (m)")
p1.1 <- p1.1 + xlab("Site Name") + ylab("d13C (per mil)")
p1.1

# d d13c bargraph
datadC <- read.csv("Carbon deltas.csv")
levels(datadC$Location) <- gsub(" ", "\n", levels(datadC$Location))
p1.5 <- ggplot(data=datadC, mapping=aes(x=factor(datadC$Location, levels=c("Bauza\nIsland", "First\nArm", "Espinosa\nPoint", "Crooked\nArm", "Campbell's\nKingdom", "Hall\nArm", "Deep\nCove")), y=datadC$dd13C, fill=datadC$Depth, group=datadC$Depth))
p1.5 <- p1.5 + geom_bar(stat="identity", size=1)
p1.5 <- p1.5 + theme_few() + scale_fill_manual(values=c("dodgerblue3", "deepskyblue1"), name="Water layer") + scale_y_continuous(breaks = seq(-6,1,by = 1), limits=c(-6,0.25))
p1.5 <- p1.5 + xlab("Site Name") + ylab("Change in d13C (per mil)") + geom_abline(aes(slope=0, intercept=0), lty=2)
p1.5

# stats
datadCstats <- read.csv("Carbon Deltas stats.csv")
supermodel <- lm(datadCstats$dd13C ~ datadCstats$variable)
summary(supermodel)
##
## Call:
## lm(formula = datadCstats$dd13C ~ datadCstats$variable)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.8442 -0.5256 0.2675 -0.5289 -0.0725 -1.0392 1.0544
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3544 0.4856 -2.789 0.03848 *
## datadCstats$variableOuter -3.2164 0.6424 -5.007 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8411 on 5 degrees of freedom
## Multiple R-squared: 0.8337, Adjusted R-squared: 0.8005
## F-statistic: 25.07 on 1 and 5 DF, p-value: 0.00408
d13C interpolated
## d13C Interpolated
# 5m
data5 <- subset(data, Depth==5)
data5$x <- data5$Long
data5$y <- data5$Lat
coordinates(data5) = ~x+y
#plot(data5)
x.range = as.numeric(c(166.83333333333333, 167.25))
y.range = as.numeric(c(-45.56, -45.21666666666667))
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by = 0.01), y=seq(from=y.range[1], to=y.range[2], by=0.01))
coordinates(grd) <- ~x+y
gridded(grd) <- TRUE
plot(grd, cex=1.5, col="black")
points(data5, pch=1, col="red", cex=1)

idw5 <- idw(formula = d13C ~ 1, locations = data5, newdata=grd)
## [inverse distance weighted interpolation]
idw5.output = as.data.frame(idw5)
names(idw5.output)[1:3] <- c("Long", "Lat", "d13C")
nz <- map("nzHires", fill=TRUE)

p1.2 <- ggplot() + geom_tile(data=idw5.output, aes(x=Long, y=Lat, fill=d13C)) + geom_polygon(data=nz, aes(x=long, y=lat, group=group)) + coord_fixed(1.3) + coord_map(xlim=x.range, ylim=y.range) + geom_point(data=as.data.frame(data5), aes(x=Long, y=Lat),color="red")
p1.2 <- p1.2 + scale_fill_continuous(breaks=seq(-28,-17,by=1.5), low="steelblue4", high="papayawhip", name="d13C at 5m") # steelblue4, papayawhip
p1.2 <- p1.2 + geom_text(data=as.data.frame(data5), aes(x=Long, y=Lat, label=Location, color="red"), show.legend=FALSE)
p1.2 <- p1.2 + ggtitle("5m") + theme(plot.title=element_text(hjust=0.5))
p1.2

# 10m
data10 <- subset(data, Depth==10)
data10$x <- data10$Long
data10$y <- data10$Lat
coordinates(data10) = ~x+y
idw10 <- idw(formula = d13C ~ 1, locations = data10, newdata=grd)
## [inverse distance weighted interpolation]
idw10.output = as.data.frame(idw10)
names(idw10.output)[1:3] <- c("Long", "Lat", "d13C")
nz <- map("nzHires", fill=TRUE)

p1.3 <- ggplot() + geom_tile(data=idw10.output, aes(x=Long, y=Lat, fill=d13C)) + geom_polygon(data=nz, aes(x=long, y=lat, group=group)) + coord_fixed(1.3) + coord_map(xlim=x.range, ylim=y.range) + geom_point(data=as.data.frame(data10), aes(x=Long, y=Lat),color="red")
p1.3 <- p1.3 + scale_fill_continuous(breaks=seq(-28,-17,by=1.5), low="steelblue4", high="papayawhip", name="d13C at 10m")
p1.3 <- p1.3 + geom_text(data=as.data.frame(data10), aes(x=Long, y=Lat, label=Location, color="red"), show.legend=FALSE)
p1.3 <- p1.3 + ggtitle("10m") + theme(plot.title=element_text(hjust=0.5))
p1.3

# 20m
data20 <- subset(data, Depth==20)
data20$x <- data20$Long
data20$y <- data20$Lat
coordinates(data20) = ~x+y
idw20 <- idw(formula = d13C ~ 1, locations = data20, newdata=grd)
## [inverse distance weighted interpolation]
idw20.output = as.data.frame(idw20)
names(idw20.output)[1:3] <- c("Long", "Lat", "d13C")
nz <- map("nzHires", fill=TRUE)

p1.4 <- ggplot() + geom_tile(data=idw20.output, aes(x=Long, y=Lat, fill=d13C)) + geom_polygon(data=nz, aes(x=long, y=lat, group=group)) + coord_fixed(1.3) + coord_map(xlim=x.range, ylim=y.range) + geom_point(data=as.data.frame(data20), aes(x=Long, y=Lat), color="red")
p1.4 <- p1.4 + scale_fill_continuous(breaks=seq(-28,-17,by=1.5), low="steelblue4", high="papayawhip", name="d13C at 20m")
p1.4 <- p1.4 + geom_text(data=as.data.frame(data20), aes(x=Long, y=Lat, label=Location, color="red"), show.legend=FALSE)
p1.4 <- p1.4 + ggtitle("20m") + theme(plot.title=element_text(hjust=0.5))
p1.4

# arrange together
get_legend<-function(myggplot){
tmp <- ggplot_gtable(ggplot_build(myggplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)
}
plegend <- p1.2 + scale_fill_continuous(breaks=seq(-28.5,-16,by=1.5), low="steelblue4", high="papayawhip", name="d13C")
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
legend <- get_legend(plegend)
p1.2f <- p1.2 + theme(legend.position="none")
p1.3f <- p1.3 + theme(legend.position="none")
p1.4f <- p1.4 + theme(legend.position="none")
grid.arrange(p1.2f, p1.3f, p1.4f, legend, ncol=4, widths=c(4,4,4,1.3))

# Source mixing
dataCSM <- read.csv("d13C source mixing 20m.csv")
dataCSM <- melt(dataCSM)
## Using Location as id variables
levels(dataCSM$Location) <- gsub(" ", "\n", levels(dataCSM$Location))
dataCSM$sd_min <- dataCSM$value-0.0515
dataCSM$sd_max <- dataCSM$value+0.0515
dataCSM$sd_min[dataCSM$variable=="Terrestrial"] <- with(dataCSM,sd_min[variable=="Marine"] + sd_min[variable=="Terrestrial"] - value[variable=="Terrestrial"] + .0515)
dataCSM$sd_max[dataCSM$variable=="Terrestrial"] <- with(dataCSM,sd_max[variable=="Marine"] + sd_max[variable=="Terrestrial"] - value[variable=="Terrestrial"] - .0515)
p1.6 <- ggplot(data=dataCSM, mapping=aes(x=factor(Location, levels=c("Bauza\nIsland", "First\nArm", "Espinosa\nPoint", "Crooked\nArm", "Campbell's\nKingdom", "Hall\nArm", "Deep\nCove")), y=value, fill=variable))
p1.6 <- p1.6 + geom_bar(stat="identity", size=1)
p1.6 <- p1.6 + geom_errorbar(aes(ymin=sd_min, ymax=sd_max), width=.25, position="identity")
p1.6 <- p1.6 + theme_few() + scale_fill_manual(values=c("burlywood", "deepskyblue1"), name="Source of\norganic matter\nat 20m") + scale_y_continuous(labels = scales::percent)
p1.6 <- p1.6 + xlab("Site Name") + ylab("Composition of d13C signature at 20m")
p1.6

# stats
dataCSMstat <- read.csv("d13CSM stats.csv")
leveneTest(dataCSMstat$value ~ dataCSMstat$variable)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.7396 0.4291
## 5
supermodel <- lm(dataCSMstat$value ~ dataCSMstat$variable)
summary(supermodel)
##
## Call:
## lm(formula = dataCSMstat$value ~ dataCSMstat$variable)
##
## Residuals:
## 1 2 3 4 5 6 7
## -0.01017 0.07717 0.06888 -0.10189 -0.21984 0.15283 0.03300
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.37895 0.06854 5.529 0.00265 **
## dataCSMstat$variableTerrestrial 0.44884 0.10470 4.287 0.00781 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1371 on 5 degrees of freedom
## Multiple R-squared: 0.7861, Adjusted R-squared: 0.7433
## F-statistic: 18.38 on 1 and 5 DF, p-value: 0.007812
## C:N by site
p2 <- ggplot(data=data, mapping=aes(x=factor(data$Location, levels=c("Bauza\nIsland", "First\nArm", "Espinosa\nPoint", "Crooked\nArm", "Campbell's\nKingdom", "Hall\nArm", "Deep\nCove")), y=data$CNr, fill=factor(data$Depth), group=factor(data$Depth)))
p2 <- p2 + geom_bar(stat="identity", position=position_dodge(0.9))
p2 <- p2 + theme_few() + scale_fill_brewer(palette="Dark2", name="Depth (m)")
p2 <- p2 + xlab("Site Name") + ylab("C:N Ratio")
p2
