###################################################################################################################
# This R script contains functions for fitting power-regressions and linear-regression models
# Author: Logan Berner, OSU
# Date: Nov 4th, 2015
###################################################################################################################


fit.power <- function(x,y){
  outlist <- list()
  df <- data.frame(x=round(x),y=y)
  
  if (min(df$x)<0){
    offset <- min(df$x)*-1
    outlist$my.nlm <- nls(y~a*(x+offset)^b, start=list(a=1, b=1),df)
  } else {
    outlist$my.nlm <- nls(y~a*(x)^b, start=list(a=1, b=1),df)
  }
  
  my.nlm.summary <- summary(outlist$my.nlm)
  outlist$my.nlm.summary <- my.nlm.summary
  # model predictions
  outlist$preds <- data.frame(x=seq(min(df$x),max(df$x),1))
  outlist$preds$modeled <- predict(object=outlist$my.nlm, newdata=outlist$preds)
  resids <- outlist$my.nlm.summary$residuals
  resid.ss <- sum((df$y-predict(outlist$my.nlm))^2)
  total.ss <- (nrow(df)-1)*var(df$y)
  r2 <- round(1-resid.ss/total.ss,2)
  adj.r2 <- r2-(1-r2)*(1/(nrow(df)-1-1))
  rmse  <- sqrt(mean(resids^2))
  outlist$coefs <- data.frame(
    a = round(my.nlm.summary$coefficients[1],3), a.se = round(my.nlm.summary$coefficients[3],3),
    b = round(my.nlm.summary$coefficients[2],3), b.se = round(my.nlm.summary$coefficients[4],3),
    adj.r2 = adj.r2, p = round(my.nlm.summary$coefficients[2,4],3), 
    df = paste(my.nlm.summary$df[1], my.nlm.summary$df[2], sep=','),
    rmse = round(rmse,2))
  outlist
}


fit.lm <- function(x,y){
  outlist <- list()
  df <- data.frame(x=x,y=y)
  outlist$my.lm <- lm(y~x,df)
  my.lm.summary <- summary(outlist$my.lm)
  outlist$my.lm.summary <- my.lm.summary
  # model predictions
  outlist$preds <- data.frame(x=seq(min(df$x),max(df$x),0.01))
  outlist$preds$modeled <- predict(object=outlist$my.lm, newdata=outlist$preds)
  resids <- outlist$my.lm.summary$residuals
  resid.ss <- sum((df$y-predict(outlist$my.lm))^2)
  total.ss <- (nrow(df)-1)*var(df$y)
  r2 <- round(1-resid.ss/total.ss,3)
  adj.r2 <- r2-(1-r2)*(1/(nrow(df)-1-1))
  rmse  <- sqrt(mean(resids^2))
  outlist$coefs <- data.frame(
    a = round(my.lm.summary$coefficients[1],5), a.se = round(my.lm.summary$coefficients[3],5),
    b = round(my.lm.summary$coefficients[2],5), b.se = round(my.lm.summary$coefficients[4],5),
    adj.r2 = adj.r2, p = round(my.lm.summary$coefficients[2,4],4), 
    df = paste(my.lm.summary$df[1], my.lm.summary$df[2], sep=','),
    rmse = round(rmse,2))
  outlist
}
