# If you copy and paste this code to your R console, you can obtain the same graph above.

if(!require(remotes)) install.packages("readr")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(minpack.lm)) install.packages("minpack.lm")
if(!require(nlraa)) install.packages("nlraa")
library (readr)
library(ggplot2)
library(minpack.lm)
library(nlraa)

github="https://raw.githubusercontent.com/agronomy4future/raw_data_practice/main/sulphur%20application.csv"
df=data.frame(read_csv(url(github),show_col_types = FALSE))

fit.lm= lm(yield ~ poly(sulphur,2, raw=TRUE),data=df)
a_parameter= fit.lm$coefficients[1]
b_parameter= fit.lm$coefficients[2]
c_parameter= fit.lm$coefficients[3]
x_mean= mean(df$sulphur)

qp= function(x, a, b, jp) {
  c= -0.5 * b / jp
  ifelse(x < jp,
          a + (b * x) + (c * x * x),
          a + (b * jp) + (c * jp * jp))
}

quadratic= lm(yield~poly(sulphur,2,raw= TRUE), data= df)
parameter= coef(quadratic)

model2= nls(formula=yield~qp(sulphur, a, b, jp),
            data=df,
            start=list(a=parameter[1],
                       b=parameter[2],
                       jp=mean(df$sulphur)))

df %>%
  ggplot(aes(sulphur, yield)) +
  geom_point(size=4, alpha = 0.5, color="grey15") +
  geom_line(stat="smooth",
            method="nls",
            formula=y ~ SSquadp3xs(x, a, b, jp),
            se=FALSE,
            color="darkred") +
          
  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, size=5)+  
                            
  scale_x_continuous(breaks = seq(18,28,2), limits = c(18,28)) +
  scale_y_continuous(breaks = seq(1000,1500,100), limits = c(1000,1500)) +
  labs(x="Sulphur application (kg/ha)", y="Yield (kg/ha)") +
  theme_classic (base_size=18, base_family="serif")+
  theme(legend.position="none",
        axis.line=element_line(linewidth=0.5, colour="black"))

Previous post
□ How to analyze linear plateau model in R Studio?

In my previous post, I explained how to analyze the linear plateau model. I simulated yield data for five different crop varieties with varying sulphur applications and suggested that the optimum sulphur application would be 23.3 kg/ha based on the linear plateau model. In this post, I’ll explain how to analyze the quadratic plateau model using the same data in R Studio.

1) Data upload

If you run the below code, the data which I simulated will be uploaded to R.

if(!require(remotes)) install.packages("readr")
library (readr)

github="https://raw.githubusercontent.com/agronomy4future/raw_data_practice/main/sulphur%20application.csv"
df=data.frame(read_csv(url(github),show_col_types = FALSE))

print(head(df,3))
  variety sulphur yield
1     CV1   20.67  1286
2     CV1   19.90  1165
3     CV1   19.73  1176
.
.
.

Or, you can directly generate data like below.

variety= rep(c("CV1","CV2","CV3","CV4","CV5"), each=9)
sulphur= c(20.67, 19.9, 19.73,21.16, 21.61, 20.9, 22.36, 21.97, 21.9, 21.34, 21.28, 
           21.9, 20.7, 22.65, 22.16, 22.03, 20.4, 22.47, 22.65, 21.82, 24.41, 22.34, 
           24.03, 24.38, 23.6, 23.44, 23.51, 23.89, 24.7, 22.8, 23.08, 25.05, 24.09, 
           24.3, 25.23, 23.14, 24.06, 26.98, 25.71, 23.93, 25.16, 23.79, 26.34, 24.64, 
           24.77)
yield= c(1286, 1165, 1176, 1262, 1271, 1221, 1374, 1288, 1285, 1261, 1281, 1302, 1199, 
         1368, 1340, 1272, 1276, 1329, 1361, 1267, 1380, 1386, 1408, 1380, 1378, 1413, 
         1400, 1412, 1422, 1381, 1410, 1410, 1381, 1422, 1389, 1403, 1420, 1403, 1400, 
         1420, 1379, 1393, 1417, 1415, 1383)
df= data.frame(variety,sulphur,yield)

Let’s assume that this data is about yield data for five different crop varieties with different sulphur applications.

2) General quadratic regression model

First of all, I’ll analyze general quadratic regression model.

if(!require(remotes)) install.packages("readr")
library (ggplot2)

ggplot(data=df, aes(x=sulphur, y=yield)) +
  geom_point(aes(shape=variety, fill=variety),col="black", size=5, stroke = 0.5) +
  geom_smooth(method=lm, level=0.95, se=FALSE, linetype=1, size=0.5, color="Red", 
              formula=y~poly(x,2, raw=TRUE)) +
  scale_shape_manual(values=rep(c(21),5)) +
  scale_fill_manual(values=rep(c("grey55"),5))+
  scale_x_continuous(breaks=seq(15,30,5), limits=c(15,30)) +
  scale_y_continuous(breaks=seq(500,1800,300), limits=c(500,1800)) +
  labs(x="sulphur application (kg/ha)", y="Yield (kg/ha)") +
  theme_classic(base_size=15, base_family="serif")+
  theme(legend.position="none",
        axis.line=element_line(linewidth=0.5, colour="black"))

I obtained the graph like above. Now I’ll do statistical analysis for the quadratic regression model.

regression= lm(yield ~ poly(sulphur,2, raw=TRUE), data=df)
summary(regression)

Call:
lm(formula = yield ~ poly(sulphur, 2, raw = TRUE), data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-51.536 -23.671   1.901  17.072  60.441 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   -3707.735    643.594  -5.761 8.76e-07 ***
poly(sulphur, 2, raw = TRUE)1   403.477     55.829   7.227 6.90e-09 ***
poly(sulphur, 2, raw = TRUE)2    -7.948      1.206  -6.588 5.67e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 27.53 on 42 degrees of freedom
Multiple R-squared:  0.8637,	Adjusted R-squared:  0.8572 
F-statistic: 133.1 on 2 and 42 DF,  p-value: < 2.2e-16

It shows the quadratic term is significant, and also regression variance is significant. Therefore, I can suggest the model equation; y= -3707.735 + 403.4772x - 7.948x2.

Now, I’d like to find the optimum sulphur amount to bring about the greatest yield. Simply we can ‘predict’ according to the model equation, but I’d like to know the critical value as a break point (the x-value above which there is no further increase in y), indicating the ‘plateau value’ (the statistically highest value that y reaches) as quadratic regression model. This is called ‘Quadratic plateau model.’

3) Quadratic plateau model

Now, I’ll introduce how to analyze quadratic plateau model in R step by step.

3.1) to find reasonable initial values for parameters

fit.lm= lm(yield ~ poly(sulphur,2, raw=TRUE),data=df)
a_parameter= fit.lm$coefficients[1]
b_parameter= fit.lm$coefficients[2]
c_parameter= fit.lm$coefficients[3]
x_mean= mean(dataA$sulphur)

Now, the general values for quadratic regression were calculated again.

a_parameter = intercept, 
b_parameter = slope x, 
c_parameter = slope x2 
x_mean = average of sulphur. 

These values are the same as model equation, y= -3707.735 + 403.4772x - 7.948x2

3.2) to define quadratic plateau function

Quadratic plateau model is defined as below.

# a = intercept
# b = slope
# c = quadratic term (curvy bit)
# jp = join point = break point

qp= function(x, a, b, jp) {
  c= -0.5 * b / jp
  ifelse(x < jp,
          a + (b * x) + (c * x * x),
          a + (b * jp) + (c * jp * jp))
}

FYI, this is the definition of linear plateau model.

#a = intercept
#b = slope
#jp = join point or break point 

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

3.3) to find the best fit parameters

if(!require(dplyr)) install.packages("dplyr")
library(dplyr)

model= nls(formula=yield ~ qp(sulphur, a, b, jp),
           data=df,
           start=list(a=a_parameter, b=b_parameter, jp=x_mean))

summary(model)
Formula: yield ~ qp(sulphur, a, b, jp)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
a  -4174.2608   923.9016  -4.518 5.01e-05 ***
b    445.3539    81.5181   5.463 2.34e-06 ***
jp    25.0656     0.4925  50.897  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 27.33 on 42 degrees of freedom

Number of iterations to convergence: 10 
Achieved convergence tolerance: 5.685e-06

I obtained the quadratic plateau model equation.

y≈ -4174.26 + 445.35x + __x².

This code does not provide quadratic parameter, but simply we can calculate by hand.

c= -0.5 * b / jp

c= -0.5 * 445.3539 / 25.0656 ≈ -8.88

and critical value is 25.0656. Then, quadratic plateau model would be written as

y≈ -4174.26 + 445.35x -8.88x2 (if x < 25.1)

This is another code to find the best fit parameters

quadratic= lm(yield~poly(sulphur,2,raw= TRUE), data= df)
parameter= coef(quadratic)

model= nls(formula=yield~qp(sulphur, a, b, jp),
            data=df,
            start=list(a=parameter[1],
                       b=parameter[2],
                       jp=mean(df$sulphur)))

summary(model)

----------- or ----------- 

if(!require(minpack.lm)) install.packages("minpack.lm")
library(minpack.lm)

model= nlsLM(formula=yield~qp(sulphur, a, b, jp),
             data=df,
             start=list(a=a_parameter,
                         b=b_parameter,
                         jp=x_mean))

summary(model)

---

If you use nlraa() package, you can obtain the model equation at once without extra coding. However, it would be better to follow what I suggested step by step to understand the concept of quadratic plateau model.

if(!require(nlraa)) install.packages("nlraa")
library(nlraa)

model= nlsLM(formula=yield~SSquadp3xs(sulphur,a,b,jp),
       data=df)

summary(model)
Formula: yield ~ SSquadp3xs(sulphur, a, b, jp)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
a  -4174.1127   923.8811  -4.518 5.01e-05 ***
b    445.3407    81.5162   5.463 2.34e-06 ***
jp    25.0657     0.4925  50.896  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 27.33 on 42 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 1.49e-08

3.4) Visualization

if(!require(ggplot2)) install.packages("ggplot2")
if(!require(minpack.lm)) install.packages("minpack.lm")
if(!require(nlraa)) install.packages("nlraa")
library(ggplot2)
library(minpack.lm)
library(nlraa)

df %>%
  ggplot(aes(sulphur, yield)) +
  geom_point(size=4, alpha = 0.5, color="grey15") +
  geom_line(stat="smooth",
            method="nls",
            formula=y ~ SSquadp3xs(x, a, b, jp),
            se=FALSE,
            color="darkred") +
  scale_x_continuous(breaks = seq(18,28,2), limits = c(18,28)) +
  scale_y_continuous(breaks = seq(1000,1500,100), limits = c(1000,1500)) +
  labs(x="Sulphur application (kg/ha)", y="Yield (kg/ha)") +
  theme_classic (base_size=18, base_family="serif")+
  theme(legend.position="none",
        axis.line=element_line(linewidth=0.5, colour="black"))

The breaking point is 25.1 and at this point, yield would be 1,410 according to the model equation, y≈ -4174.26 + 445.35x -8.88x2 (if x < 25.1).

To visualize more clearly, more code will be added.

if(!require(ggplot2)) install.packages("ggplot2")
if(!require(minpack.lm)) install.packages("minpack.lm")
if(!require(nlraa)) install.packages("nlraa")
library(ggplot2)
library(minpack.lm)
library(nlraa)

df %>%
  ggplot(aes(sulphur, yield)) +
  geom_point(size=4, alpha = 0.5, color="grey15") +
  geom_line(stat="smooth",
            method="nls",
            formula=y ~ SSquadp3xs(x, a, b, jp),
            se=FALSE,
            color="darkred") +
          
  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, size=5)+    
             
  scale_x_continuous(breaks = seq(18,28,2), limits = c(18,28)) +
  scale_y_continuous(breaks = seq(1000,1500,100), limits = c(1000,1500)) +
  labs(x="Sulphur application (kg/ha)", y="Yield (kg/ha)") +
  theme_classic (base_size=18, base_family="serif")+
  theme(legend.position="none",
        axis.line=element_line(linewidth=0.5, colour="black"))

3.5) How to interpret quadratic plateau model?

y≈ -4174.26 + 445.35x -8.88x2 (if x < 25.1)

In the model equation; y≈ -4174.26 + 445.35x -8.88x2 (if x < 25.1), the important thing is critical value, 25.1. Statistically, at this point, yield does not increase although sulphur application increases. Therefore, critical value is regarded as the highest value that yield reaches.

and when sulphur application is 25.1, yield would be

y≈ -4174.26 + 445.35*25.0656 – 8.88*25.0656^2 ≈ 1409.5

That is, when 25.1 S kg/ha was applied, yield would be 1409.5 kg/ha and no more yield increases, indicating the optimum sulphur application would be 25.1 kg/ha.

4) p-value and pseudo R-squared

nullfunct= function(x, m){m}
m.ini= mean(dataA$yield)
null= nls(yield~nullfunct(sulphur, m),
         data=df,
         start=list(m=m.ini),
         trace=FALSE,
         nls.control(maxiter= 1000))

library(rcompanion)
nagelkerke(model, null)

■ Code summary: https://github.com/agronomy4future/r_code/blob/main/How_to_analyze_quadratic_plateau_model_in_R_Studio.ipynb

---

Reference
https://gradcylinder.org/post/linear-plateau/
https://rcompanion.org/handbook/I_11.html

---