###################################################################################################################
# This R script fits models AGB, ANPP, and LAI as a function of CMIgy using power-regression models 
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(reshape2)
require(plotrix)
setwd("")
source("fun_fit_models.R")

plot.dat <- read.csv('data/plot_summaries/plot_master_data.csv', header = T)
voi <- c('tree.kgCm2.2014.avg','tree.anpp.gCm2yr.avg','lai.hc.avg','cmi.gy.avg','vpd.gy.avg','freeze.gy.avg')
plot.dat <- plot.dat[,which(colnames(plot.dat)%in%voi)]

# AGB
agb.cmi.nls <- with(plot.dat, fit.power(cmi.gy.avg,tree.kgCm2.2014.avg))
plot(tree.kgCm2.2014.avg~cmi.gy.avg, plot.dat, xlim=c(-1300,100),ylim=c(0,25))
lines(agb.cmi.nls$preds$x, agb.cmi.nls$preds$modeled)

# NPP
npp.cmi.nls <- with(plot.dat, fit.power(cmi.gy.avg,tree.anpp.gCm2yr.avg))
plot(tree.anpp.gCm2yr.avg~cmi.gy.avg, plot.dat, xlim=c(-1300,100),ylim=c(0,0.3))
lines(npp.cmi.nls$preds$x, npp.cmi.nls$preds$modeled)
npp.cmi.nls$coefs

# LAI
summary(with(plot.dat, lm(lai.hc.avg~cmi.gy.avg)))
lai.cmi.nls <- with(plot.dat, fit.power(cmi.gy.avg,lai.hc.avg))
plot(lai.hc.avg~cmi.gy.avg, plot.dat, xlim=c(-1300,100),ylim=c(0,10))
lines(lai.cmi.nls$preds$x, lai.cmi.nls$preds$modeled)
lai.cmi.nls$coefs

# summary across models
model.coefs <- rbind(agb.cmi.nls$coefs,npp.cmi.nls$coefs,lai.cmi.nls$coefs)
model.coefs$y.var <- c('agb','npp','lai')
model.coefs
write.table(model.coefs, 'outputs/agb_npp_lai_by_cmi_nls.csv', sep = ',', col.names = T, row.names = F)

save.image('workspaces/forest_by_cmi_nls.RDATA')
