Spillman, Linear response with Plateau, and Qadraticwith Plateau

Spillman Vs Linear response with Plateau

# Spillman/ Assymptotic
#Niger, 2015
#Dosso
#Bengou
prate=c(0,7.5,15,22.5)
yield=c(826.5,1127.4, 1203.7, 1250)
#prate<- c(0,7.5,15,22.5, 30, 32.5)
#yield <- c(826.5,1035.4,1203.7,1250, 1250, 1250)
dat <- as.data.frame(cbind(prate,yield))
dat
mod9=nls(yield~a-b*(c^prate), start=list(a=3000, b=800,c=0.9),data=dat)
mod9
summary(mod9)
coef(mod9)

plot(yield~prate,ylab="Yield (kg/ha)", xlab="P input (kg/ha)",xlim=c(0,32.5), ylim=c(800,1300),las=1,cex.lab=1,cex.axis =0.9)
#points(8.7,1143,pch = 24)

M <-1256.2618098
A <-429.1096275
R <-0.8560028
lines(x<-c(0:32.5), M-A*(R^x),col='black ',lwd=2, lty=2)
library(nlstools)
mod10=yield~a-b*(c^prate)
preview(mod10,data=dat, 
        start=list(a=1500, b=800,c=0.9))
overview(mod9)
plotfit(mod9, smooth = TRUE)
residuals(mod9)
predict(mod9)
Resid=nlsResiduals(mod9)
plot(Resid)

#Linear response to plateau or also known as Linear stochastic response plateau 
#𝑦𝑖𝑗𝑡 =min(𝛽0+ 𝛽1𝑋𝑖𝑗𝑡,𝜇)+ 𝜀𝑖𝑗𝑡,𝑗𝑡,
#Important for convergence is the dataset to have a plateau
prate<- c(0,7.5,15,22.5, 30, 32.5)
yield <- c(826.5,1035.4,1203.7,1250, 1250, 1250)

dat1 <- as.data.frame(cbind(prate,yield))
dat1
f.lrp=function(prate,a,b,t.x){
  ifelse(prate>t.x,a+b*t.x,a+b*prate)
}
#t.x is maximum fertilizer that will give any result
m <- nls(yield ~ f.lrp(prate, a, b, t.x), data =dat1, start = list(a = 500, b = 100, t.x = 22.5), trace = T ) 
m
summary(m)
coefficients(m)
lines(fitted(m) ~ prate, lwd=2)

#Quadratic to plateau
prate<- c(0,7.5,15,22.5, 30, 32.5)
yield <- c(826.5,1035.4,1203.7,1250, 1250, 1250)
datf<- as.data.frame(cbind(prate,yield))
datf
fm<- nls(yield ~ (b0 + b1*prate + b2*I(prate^2))*(prate <= x0)
                + (b0 + b1*x0 + b2*I(x0^2))*(prate > x0),
                 data=datf,
                 start=list(b0=800, b1=1, b2=1, x0=30),
                 trace=T)
summary(fm)


lines(fitted(fm) ~ prate, lwd=2, lty=3)


#Adding legends 

legend(11, 1000, legend=c("Spillman (Assyptotic)", "Linear to Plateau", "Quadratic to Plateau"),
       col=c("black","black","black"),lty=c(2,1,3), cex=1.2)

#Goodness of fit tests

deviance(mod9)
logLik(mod9)
fitted(mod9)
plot(fitted(mod9), residuals(mod9))
abline(h=0)
standardRes <- residuals(mod9)/summary(mod9)$sigma
qqnorm(standardRes)
abline(a = 0, b = 1)