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