Fecundity.R

###############################################################################
# 22/01/2013
# Author: Alastair V Harry - alastair.harry@gmail.com
# Centre For Sustainable Tropical Fisheries & Aquaculture
# James Cook University, Townsville
# 
# From Harry et al (2013) Age, growth, and reproductive biology of the spot-tail
# shark, Carcharhinus sorrah, and the Australian blacktip shark, C. tilstoni, 
# from the Great Barrier Reef World Heritage Area, north-eastern Australia. 
# Marine and Freshwater Research
#
# Runs a linear regression of fecundity and length and plots the data in the 
# style of Walker (2005) Reproduction in fisheries science
###############################################################################

library(RCurl)

# Load data from github
url <- "https://raw.githubusercontent.com/alharry/spot-tail/master/SSS.csv"
data <- getURL(url)                

# If this step doesnt work use this, but read this first
# http://ademar.name/blog/2006/04/curl-ssl-certificate-problem-v.html
#data <- getURL(url,ssl.verifypeer = FALSE)

data <- read.csv(textConnection(data))

# Tabulate records of total embryo counts
table(data$Emb)

# Run a linear regression of number of embryos (depedent variable) against 
# stretched total length (STL) in mm. 
Fecundity.Model<-lm(Emb~STL,data)
summary(Fecundity.Model)

# Plot and add confidence and prediction intervals

# Define minimum and maximum range of pregnant females with embryos
min.x<-min(data$STL[data$Emb>0],na.rm=T)
max.x<-max(data$STL[data$Emb>0],na.rm=T)

# Specify range of lengths over which to predict fecundit
new <- data.frame(STL = seq(min.x,max.x,0.1))
# Calculate confidence intervals
pred.w.plim <- predict(Fecundity.Model, new, interval="prediction")
# Calculate prediction intervals
pred.w.clim <- predict(Fecundity.Model, new, interval="confidence")
# Plot the linear model with confidence and prediction intervals
matplot(new$STL,cbind(pred.w.clim, pred.w.plim[,-1]), lty=c(1,2,2,3,3),
        	col=1, type="l",las=1,bty="l",xlab="STL (mm)",
		ylab="",
		xlim=c(min.x,max.x),ylim=c(min(data$Emb,na.rm=T),max(data$Emb,na.rm=T)))
mtext(bquote(paste("No. ",italic("in"), " ", italic(utero), " embryos")),side=2,padj=-30)
points(data$Emb~data$STL,pch=19)