# 2013_limpet_mass_analysis.R
# 
# Author: Luke Miller Jan 14, 2015
###############################################################################

fname = '2013_limpet_mass_master.csv'
wt2 = read.csv(fname)
# Remove all rows with chipped shell
wt3 = wt2[-(which(wt2$notes == 'chipped shell')),]
# Convert tissue mass in grams to milligrams
wt3$dry.tissue.mg = wt3$net.dry.tissue.mass.g * 1000
# Reorder the species 
wt3$Species = factor(as.character(wt3$Species),
		levels = c("L. austrodigitalis","L. scabra","L. limatula","L. pelta"))

# Fit linear regression to log-transformed values, extract back-transformed
# coefficients for use in a power-law equation to predict dry tissue mass from
# shell area. 
wt3 = wt3[-(is.na(wt3$shell.area.mm2)),] # remove rows with no shell area value
sppvec = levels(wt3$Species)
wt.fits = data.frame(Species = character(4), alpha = numeric(4),
		beta = numeric(4), MSE = numeric(4), R2 = numeric(4))
wt.fits$Species = factor(sppvec) # copy over species names
# Reorder the species
wt.fits$Species = factor(as.character(wt.fits$Species),
		levels = c("L. austrodigitalis","L. scabra","L. limatula","L. pelta"))
for (i in 1:length(sppvec)){
	# fit linear model to ln-transformed data
	mod = lm(log(dry.tissue.mg)~log(shell.area.mm2), 
			data = wt3[wt3$Species == sppvec[i],])
	summod = summary(mod)
	# Calculate mean square error for the model
	wt.fits[i,'MSE'] = deviance(mod)/df.residual(mod)
	# extract slope and intercept from linear model
	wt.fits[i,'alpha'] = mod$coefficients[1]
	wt.fits[i,'beta'] = mod$coefficients[2]
	wt.fits[i,'R2'] = summod$r.squared
}
# The contents of wt.fits are still in log-transformed (geometric) space. To 
# predict a mass from them, you need to back-transform to arithmetic space.
# The exp(MSE/2) factor in the equation below is based on the correction 
# proposed in Mascaro et al 2014, Biol. J. Linn. Soc. v111 p230, based on
# Baskerville 1972, Can. J. Forest Research, v2 p49.
# mass in mg = (exp(alpha) * (shell.areas ^ beta)) * exp(MSE/2)

###############################################################
# Now calculate the same kind of fit for empty shell mass
# Create a column of shell mass in milligrams
wt3$empty.shell.mass.mg = wt3$empty.shell.mass.g * 1000
# Add extra columns to wt.fits data frame
wt.fits$alphaShell = 0
wt.fits$betaShell = 0 
wt.fits$MSEShell = 0
wt.fits$R2Shell = 0
for (i in 1:length(sppvec)){
	# fit linear model to ln-transformed data
	mod = lm(log(empty.shell.mass.mg)~log(shell.area.mm2), 
			data = wt3[wt3$Species == sppvec[i],])
	summod = summary(mod)
	# Calculate mean square error for the model
	wt.fits[i,'MSEShell'] = deviance(mod)/df.residual(mod)
	# extract slope and intercept from linear model
	wt.fits[i,'alphaShell'] = mod$coefficients[1]
	wt.fits[i,'betaShell'] = mod$coefficients[2]
	wt.fits[i,'R2Shell'] = summod$r.squared
}



########################################################################
# Create a 4-panel figure showing the relationship between shell projected area
# and dry tissue mass (mg)
par(mfrow = c(2,2)) # declare 4-panel layout
for (i in 1:length(sppvec)){ # cycle through each species contained in the
							 # vector 'sppvec' created earlier
	plot(dry.tissue.mg~shell.area.mm2, 
			data = wt3[wt3$Species == sppvec[i],],
			main = sppvec[i],
			las = 1,
			ylab = 'Dry tissue mass, mg',
			xlab = expression(Shell~projected~list(area,mm^2)), 
			pch = 20,
			col = rgb(0,0,0,0.6), # add alpha transparency to plotted points
			bg = rgb(0,0,0,0.6))
	# get sample size
	n = length(wt3[wt3$Species == sppvec[i],'shell.area.mm2'])
	# get min and max shell areas
	lohi = c(min(wt3[wt3$Species == sppvec[i],'shell.area.mm2'], na.rm=TRUE),
			max(wt3[wt3$Species == sppvec[i],'shell.area.mm2'], na.rm=TRUE))
	# make sequence of new x values
	shell.areas = seq(lohi[1],lohi[2],0.1)
	# Get the power law coefficients from wt.fits data frame
	alpha = wt.fits[wt.fits$Species == sppvec[i],'alpha']
	beta = wt.fits[wt.fits$Species == sppvec[i],'beta']
	MSE = wt.fits[wt.fits$Species == sppvec[i], 'MSE']
	# calculate back-transformed masses
	wts = (exp(alpha) * (shell.areas ^ beta)) * exp(MSE/2) 
	# Plot back-transformed power law fit
	lines(x = shell.areas, y = wts, lty = 1)
	# Plot the R^2 and power law equation
	rp = vector('expression',3)
	rp[1] = substitute(expression(italic(R)^2 == MYVALUE),
			list(MYVALUE = format(wt.fits[i,'R2'],digits = 2)))[2]
	rp[2] = substitute(expression(italic(Y) == MYVALUE2~X^MYVALUE3),
			list(MYVALUE2 = format(exp(wt.fits[i,'alpha']), digits = 3),
					MYVALUE3 = format(wt.fits[i,'beta'],digits = 4)))[2]
	rp[3] = substitute(expression(plain(n)==MYVALUE4),
			list(MYVALUE4 = n))[2]
	legend('topleft', legend = rp, bty = 'n')
}

#############################################################
# Same plot, using shell mass instead of tissue mass
windows() # open a new plot window, on Mac use quartz() instead
par(mfrow = c(2,2))
for (i in 1:length(sppvec)){
	plot(empty.shell.mass.mg~shell.area.mm2, 
			data = wt3[wt3$Species == sppvec[i],],
			main = sppvec[i],
			las = 1,
			ylab = 'Shell mass, mg',
			xlab = expression(Shell~projected~list(area,mm^2)),
			pch = 20,
			col = rgb(0,0,0,0.6),
			bg = rgb(0,0,0,0.6)) 
	# get sample size
	n = length(wt3[wt3$Species == sppvec[i],'shell.area.mm2'])		
	# get min and max shell areas
	lohi = c(min(wt3[wt3$Species == sppvec[i],'shell.area.mm2'], na.rm=TRUE),
			max(wt3[wt3$Species == sppvec[i],'shell.area.mm2'], na.rm=TRUE))
	# make sequence of new x values
	shell.areas = seq(lohi[1],lohi[2],0.1)
	# Get the power law coefficients from wt.fits data frame
	alpha = wt.fits[wt.fits$Species == sppvec[i],'alphaShell']
	beta = wt.fits[wt.fits$Species == sppvec[i],'betaShell']
	MSE = wt.fits[wt.fits$Species == sppvec[i], 'MSEShell']
	# calculate back-transformed masses
	wts = (exp(alpha) * (shell.areas ^ beta)) * exp(MSE/2) 
	# Plot back-transformed power law fit
	lines(x = shell.areas, y = wts, lty = 1)
	
	rp = vector('expression',3)
	rp[1] = substitute(expression(italic(R)^2 == MYVALUE),
			list(MYVALUE = format(wt.fits[i,'R2Shell'],digits = 2)))[2]
	rp[2] = substitute(expression(italic(Y) == MYVALUE2~X^MYVALUE3),
			list(MYVALUE2 = format(exp(wt.fits[i,'alphaShell']), digits = 3),
					MYVALUE3 = format(wt.fits[i,'betaShell'],digits = 4)))[2]
	rp[3] = substitute(expression(plain(n)==MYVALUE4),
			list(MYVALUE4 = n))[2]
	legend('topleft', legend = rp, bty = 'n')
	
}