LINEAR PLATEAU MODEL

# https://agronomy4future.org/?p=17137

pacman::p_load(tidyverse, nlraa, minpack.lm)
github <- "https://raw.githubusercontent.com/agronomy4future/raw_data_practice/main/sulphur%20application.csv"
dataA <- readr::read_csv(github)
  
# QUADRATIC PLATEAU
# to find reasonable initial values for parameters
fit.lm<-lm(yield ~ poly(sulphur,2, raw=TRUE),data=dataA)

a_parameter<-fit.lm$coefficients[1]
b_parameter<-fit.lm$coefficients[2]
c_parameter<-fit.lm$coefficients[3]
x_mean<-mean(dataA$sulphur)

# to define quadratic  plateau function
# a = intercept
# b = slope
# c = quadratic term (curvy bit)
# jp = join point = break point = critical concentration

qp <- function(x, a, b, jp) {
  c <- -0.5 * b / jp
  if_else(condition = x < jp,
          true  = a + (b * x) + (c * x * x),
          false = a + (b * jp) + (c * jp * jp))
}

# to find the best fit parameters
model <- nls(formula=yield ~ qp(sulphur, a, b, jp),
             data=dataA,
             start=list(a=a_parameter, b=b_parameter, jp=x_mean))

summary(model)

dataA %>% 
  ggplot(aes(sulphur, yield)) +
  geom_point(size=4, alpha = 0.5) +
  geom_line(stat="smooth",
            method="nls",
            formula=y~SSquadp3xs(x,a,b,jp),
            se=FALSE,
            color="Dark red") +
  geom_vline(xintercept=25.0656, linetype="solid", color="grey") +
  annotate("text", label=paste("sulphur=","25.1","kg/ha"),
           x=25.2, y=1000, angle=90, hjust=0, vjust=1.5, alpha=0.5)+
  labs(x="Sulphur application (kg/ha)", y="Yield (kg/ha)") 

# to find reasonable initial values for parameters
fit.lm<-lm(yield~sulphur,data=dataA)
a_parameter<-fit.lm$coefficients[1]
b_parameter<-fit.lm$coefficients[2]
x_mean<-mean(dataA$sulphur)

# to define quadratic  plateau function
#a = intercept
#b = slope
#jp = join point or break point 

linplat<- function(x, a, b, jp){
  ifelse(condition = x<jp, 
         true = a+b*x,
         false = a+b*jp)
}

model <- nlsLM(formula = yield ~ linplat(sulphur, a, b, jp),
               data = dataA,
               start=list(a=a_parameter,
                          b=b_parameter,
                          jp=x_mean))

summary(model)

dataA %>% 
  ggplot(aes(sulphur, yield)) +
  geom_point(size=4, alpha = 0.5) +
  geom_line(stat="smooth",
            method="nlsLM",
            formula=y~SSlinp(x,a,b,jp),
            se=FALSE,
            color="Dark red") +
  geom_vline(xintercept=23.2722, linetype="solid", color="grey") +
  annotate("text", label=paste("sulphur=","23.3","kg/ha"),
           x=23.3, y=1000, angle=90, hjust=0, vjust=1.5, alpha=0.5)+
  labs(x="Sulphur application (kg/ha)", y="Yield (kg/ha)")