###################################################################################################################
# This R script fits linear and power-regression models relating traits, stand-level variable, and climate 
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################
rm(list=ls())
require(lattice)
require(reshape2)
require(plotrix)
setwd("")
source('fun_fit_models.R')

#---------------------------------------------------------------------------------------------------
# load plot data and identify predictor and response vaiables
#---------------------------------------------------------------------------------------------------
plot.dat <- read.csv('data/plot_summaries/plot_master_data.csv', header = T)
colnames(plot.dat)
# resp.var <- c('lai.hc.avg', 'tree.anpp.gCm2yr.avg', 'tree.kgCm2.2014.avg',
#               'rwc.pcnt.wavg','ht.m.wavg','lai.2.sap.wavg','sg.g.cm3.wavg',
#               'sla.cm2.gC.wavg','leaf.C.g.g.wavg','leaf.N.mg.g.wavg')
# pred.var <- c('cmi.gy.avg','ppt.gy.avg','tmean.gy.avg')
# resp.var <- c('lai.hc.avg', 'tree.anpp.gCm2yr.avg', 'tree.kgCm2.2014.avg',
#               'ht.m.wavg','lai.2.sap.wavg','rw.cmi.cor.avg','sg.g.cm3.wavg')
# pred.var <- c('cmi.gy.avg','lai.hc.avg', 'tree.anpp.gCm2yr.avg')
resp.var <- c('ht.m.wavg','lai.2.sap.wavg','sg.g.cm3.wavg')
pred.var <- c('ht.m.wavg','lai.2.sap.wavg','sg.g.cm3.wavg')

n.resp <- length(resp.var)
n.pred <- length(pred.var)

#---------------------------------------------------------------------------------------------------
# create output, summary data frame
#---------------------------------------------------------------------------------------------------
model.fits.df <- data.frame(matrix(NA, nrow = 0, ncol = 11))
colnames(model.fits.df) <- c('response','predictor','form','a','a.se','b','b.se','adj.r2','p','df','rmse')

#---------------------------------------------------------------------------------------------------
# fit models
#---------------------------------------------------------------------------------------------------
for (i in 1:n.resp){
  for (j in 1:n.pred){
    # grab data
    resp.dat <- plot.dat[,names(plot.dat)%in%resp.var[i]]
    pred.dat <- plot.dat[,names(plot.dat)%in%pred.var[j]]
    # fit models
    #power <- fit.power(x = pred.dat, y = resp.dat)
    linear <- fit.lm(x = pred.dat, y = resp.dat)
    # put models in workspace
    #assign(paste(resp.var[i],'by',pred.var[j], 'power',sep = '_'), power)
    assign(paste(resp.var[i],'by',pred.var[j], 'linear',sep = '_'), linear)
    # save model summary stats 
    #model.fits.df <- rbind(model.fits.df, cbind(data.frame(response=resp.var[i],predictor=pred.var[j],form='power',power$coefs)))
    model.fits.df <- rbind(model.fits.df, cbind(data.frame(response=resp.var[i],predictor=pred.var[j],form='linear',linear$coefs)))
  }
}

model.fits.df

# 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')
