R/loess_function.R

In AgronomiaR/seedreg: Regression Analysis for Seed Germination as a Function of Temperature

#' Analysis: loess regression
#'
#' Fit a polynomial surface determined by one or more numerical predictors, using local fitting.
#' @param trat Numerical or complex vector with treatments
#' @param resp Numerical vector containing the response of the experiment.
#' @param ylab Variable response name (Accepts the \emph{expression}() function)
#' @param xlab treatments name (Accepts the \emph{expression}() function)
#' @param theme ggplot2 theme (\emph{default} is theme_bw())
#' @param error Error bar (It can be SE - \emph{default}, SD or FALSE)
#' @param legend.position legend position (\emph{default} is c(0.3,0.8))
#' @param cardinal defines the value of y considered extreme (\emph{default} considers 0 germination)
#' @param scale Sets x scale (\emph{default} is none, can be "log")
#' @return The function returns the loess regression; optimal and cardinal temperatures and graph using ggplot2.
#' @seealso \link{loess}
#' @export
#' @note if the maximum predicted value is equal to the maximum x, the curve does not have a maximum point within the studied range. If the minimum value is less than the lowest point studied, disregard the value.
#' @author Gabriel Danilo Shimizu
#' @author Leandro Simoes Azeredo Goncalves
#' @examples
#' library(seedreg)
#' data("aristolochia")
#' attach(aristolochia)
#' loess_model(trat,resp)

loess_model=function(trat,
                     resp,
                     ylab="Germination (%)",
                     xlab=expression("Temperature ("^"o"*"C)"),
                     theme=theme_classic(),
                     error="SE",
                     cardinal=0,
                     legend.position="top",
                     scale="none"){
  requireNamespace("ggplot2")
  requireNamespace("crayon")
  ymean=tapply(resp,trat,mean)
  if(error=="SE"){ysd=tapply(resp,trat,sd)/sqrt(tapply(resp,trat,length))}
  if(error=="SD"){ysd=tapply(resp,trat,sd)}
  if(error=="FALSE"){ysd=0}
  desvio=ysd
  xmean=tapply(trat,trat,mean)
  mod=loess(resp~trat)
  xp=seq(min(trat),max(trat),length.out = 1000)
  preditos=data.frame(x=xp,
                      y=predict(mod,newdata = data.frame(trat=xp)))
  x=preditos$x
  y=preditos$y
  data=data.frame(xmean,ymean)
  data1=data.frame(trat=xmean,resp=ymean)
  s="~~~ Loess~regression"
  graph=ggplot(data,aes(x=xmean,y=ymean))
  if(error!="FALSE"){graph=graph+geom_errorbar(aes(ymin=ymean-ysd,ymax=ymean+ysd),
                                               width=0.5)}
  graph=graph+geom_point(aes(color="black"),size=4.5,shape=21,fill="gray")+
    theme+
    geom_line(data=preditos,aes(x=preditos$x,
                                y=preditos$y,color="black"),size=0.8)+
    scale_color_manual(name="",values=1,label="Loess regression")+
    theme(axis.text = element_text(size=12,color="black"),
          legend.position = legend.position,
          legend.text = element_text(size=12),
          legend.direction = "vertical",
          legend.text.align = 0,
          legend.justification = 0)+
    ylab(ylab)+xlab(xlab)
  if(scale=="log"){graph=graph+scale_x_log10()}
  temp1=seq(min(trat),max(trat),length.out=10000)
  result=predict(mod,newdata = data.frame(trat=temp1),type="response")
  maximo=temp1[which.max(result)]
  respmax=result[which.max(result)]
  result=round(result,0)
  fa=temp1[result<=cardinal & temp1>maximo]
  if(length(fa)>0){maxl=max(temp1[result<=cardinal & temp1>maximo])}else{maxl=NA}
  fb=temp1[result<=cardinal & temp1<maximo]
  if(length(fb)>0){minimo=max(temp1[result<=cardinal & temp1<maximo])}else{minimo=NA}
  graphs=data.frame("Parameter"=c("optimum temperature",
                                  "Maximum response",
                                  "Predicted maximum value",
                                  "Predicted minimum value"),
                    "values"=c(maximo,respmax,maxl,minimo))
  graficos=list("test"=graphs,graph)
  print(graficos)
}