setwd("H:/SERC census")
library(MuMIn); library(ggplot2)
options(na.action = "na.fail") #needed for MuMIn

# read in and prepare allometry dataset
dat.allom <- read.delim("Phrag allometry data 2011-2015.txt")
dat.allom$yr <- as.factor(dat.allom$yr)
dat.allom <- dat.allom[!is.na(dat.allom$diam.mm),]
dat.allom$bms.ln <- log(dat.allom$bms.g)


# plot all data
ggplot(dat.allom) + aes(x=ht.cm, y=bms.g, col=yr) + geom_point()  #options: bms.ln, bms.g


# evaluate linear models 
mod.1 <- lm(sqrt(bms.g) ~ ht.cm, data=dat.allom)
mod.2 <- lm(sqrt(bms.g) ~ ht.cm * diam.mm, data=dat.allom)
mod.3 <- lm(log(bms.g) ~ ht.cm, data=dat.allom)
mod.4 <- lm(log(bms.g) ~ ht.cm * diam.mm, data=dat.allom)

# evaluate normality of residuals
par(mfrow=c(2,2))
qqnorm(resid(mod.1)) # linear
qqnorm(resid(mod.2)) # very linear
qqnorm(resid(mod.3)) # slightly curved
qqnorm(resid(mod.4)) # slightly curved

# compare all models
dredge(mod.2, extra=c("R^2"))

# choose models for use with function
mod.best.w.diam <- mod.2
mod.best.no.diam <- mod.1 


# function to predict biomass based on a diameter cutoff
bms.diam.split <- function(data.set=dat.chamb, cutoff=200, trans=NULL) {
	# options for trans are c("sqrt", "log")
	
	dat.over <- data.set[data.set$ht.cm >= cutoff, ]
	dat.under <- data.set[data.set$ht.cm < cutoff, ]
	
	bms.over <- predict(mod.best.w.diam, newdata=dat.over, interval=c("prediction"), level=0.95, type="response")[,"fit"]
	bms.under <- predict(mod.best.no.diam, newdata=dat.under, interval=c("prediction"), level=0.95, type="response")[,"fit"]
	# generate output data table
	out <- rbind(data.frame(yr=dat.over$yr, chamber=dat.over$chamber, height=dat.over$ht.cm, width=dat.over$diam.mm, bms=bms.over),
			data.frame(yr=dat.under$yr, chamber=dat.under$chamber, height=dat.under$ht.cm, width=dat.under$diam.mm, bms=bms.under))
	# apply transformations
	if (trans=="log") out$bms <- exp(out$bms)
	if (trans=="sqrt") out$bms <- (out$bms)^2
	
	return(out)
}

# read in and prepare data from censuses
dat.chamb <- read.delim("Phrag CO2xN Master Phrag Census 2011-2014 (08-29-2015).txt", na.strings=-99)
dat.chamb$yr <- as.factor(dat.chamb$Year) #convert to a categorical variable
# create columns with names that match what function looks for
dat.chamb$diam.mm <- dat.chamb$Width
dat.chamb$ht.cm <- dat.chamb$Total_Height
dat.chamb$chamber <- dat.chamb$Chamber
# exclude rows with no plant or no height for plant 
dat.chamb <- dat.chamb[!is.na(dat.chamb$Total_Height),] # plant but height missing
dat.chamb <- dat.chamb[dat.chamb$Total_Height!=0,] # no plant in sub-plot

# define treatment status of the chambers in the Phragmites experiment
treats <- data.frame(
		chamber=c(1,   2,    3,   4,   5,    6,   7,    8,   9,  10,   11,   12),
		co2=c(FALSE,TRUE,FALSE,TRUE,TRUE,FALSE,TRUE,FALSE,TRUE,TRUE,FALSE,FALSE),
		nit=c(FALSE,FALSE,TRUE,TRUE,FALSE,TRUE,TRUE,FALSE,TRUE,FALSE,FALSE,TRUE),
		treat=c("A",  "E", "N","EN",  "E", "N","EN",  "A","EN",  "E",  "A", "N")
)
treats.treat <- factor(treats$treat, levels=levels(treats$treat)[c(1,2,4,3)]) #changer order


# estimate biomass for Phrag in chambers
bms.chamb <- bms.diam.split(dat.chamb, trans="sqrt", cutoff=200)
# tack on treatment information
bms.chamb <- merge(bms.chamb, treats)

# output table to file
write.table(bms.chamb, "phrag morphometry and bms pred (2011-14) 18Sept2015.txt", sep="\t", row.names=F)