MartinMacharia · October 31, 2017 05:02
diff --git a/OFRA b/OFRA
 #1988

 #Import data
 OlNgarua=read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgarua.csv")
 OlNgarua
 names(OlNgarua)
 str(OlNgarua)
 OlNgarua_1=OlNgarua[ ,c("District","Center","Longitudes","Latitudes","Year","Crop",
                        "Variety","PlantMon","AppliMethod","Nrate","Prate","Grain_Yield")]
 summary(OlNgarua_1)
 str(OlNgarua_1)
 #Subset data to SIMI and consider year 1988 with the different P levels
 #P=0
 OlNgarua_2=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1"  & Prate==0)
 View(OlNgarua_2)
 #P=25
 OlNgarua_3=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1"  & Prate==25)
 View(OlNgarua_3)
 #P=50
 OlNgarua_4=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1"  & Prate==50)
 View(OlNgarua_4)
 #P=75
 OlNgarua_5=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1"  & Prate==75)
 View(OlNgarua_5)

 #Basic Plots

 plot=plot(OlNgarua_2$Nrate,OlNgarua_2$Grain_Yield, xlim=c(0,100), ylim=c(8,12)
     ,xlab="Nitrogen Rate", ylab="Yield(t/ha)", pch= 19,col="Red",cex=0.9, cex.lab=1.2)

 #Joining points with a line
 plot1=lowess(OlNgarua_2$Grain_Yield ~ OlNgarua_2$Nrate)
 lines(smoothplot1, col = 'red')
 #Making the line smooth
 library(locfit)
 #Locfit is a software package performing local regression,likelihood
 #and related smoothing procedures.
 plot2=locfit(OlNgarua_2$Grain_Yield ~lp(OlNgarua_2$Nrate, nn = 0.80, deg = 1))
 lines(plot2, col="red", lwd=1.4)
 #nn=specifies smoothness of parameter
 #deg=specifies the degree (1,2...) of the polynomial

 #Investigate the individual plots side by side
 par(mfrow=c(2,2))
 plot3=plot(OlNgarua_2$Nrate,OlNgarua_2$Grain_Yield)#P=0
 plot4=plot(OlNgarua_3$Nrate,OlNgarua_3$Grain_Yield)#P=25
 plot5=plot(OlNgarua_4$Nrate,OlNgarua_4$Grain_Yield)#P=50
 plot6=plot(OlNgarua_5$Nrate,OlNgarua_5$Grain_Yield)#P=75

 #combining the above individual scatter plots
 plot7=plot(OlNgarua_2$Nrate,OlNgarua_2$Grain_Yield,pch= 19,xlim=c(0,85),
           ylim=c(6.5,12), xlab="N Rate (Kg/ha)"
           ,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
           ,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")
 abline(v=25,lty=5)
           
 axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.7)
 axis(2, at = seq(6.5,12, by = 0.5), las=1,cex.axis=0.7)
            
 #las use las to rotate the labels 
 ##yaxt gives the to overwrite the axis labels

 plot8=points(OlNgarua_3$Nrate,OlNgarua_3$Grain_Yield,pch= 19, col="red")#P=25
 plot9=points(OlNgarua_4$Nrate,OlNgarua_4$Grain_Yield, pch= 19,col="purple")#P=50
 plot10=points(OlNgarua_5$Nrate,OlNgarua_5$Grain_Yield, pch= 19,col="Orange")#P=75

 # Fitting lowess curves on the above points
 #0P
 plot11=locfit(OlNgarua_2$Grain_Yield ~lp(OlNgarua_2$Nrate, nn = 0.80, deg = 1))
 lines(plot11,lwd=1.7)
 #25P
 plot12=locfit(OlNgarua_3$Grain_Yield ~lp(OlNgarua_3$Nrate, nn = 0.80, deg = 1))
 lines(plot12, col="red", lwd=1.7)
 #50P
 plot13=locfit(OlNgarua_4$Grain_Yield ~lp(OlNgarua_4$Nrate, nn = 0.80, deg = 1))
 lines(plot13, col="purple", lwd=1.7)
 #75P
 plot14=locfit(OlNgarua_5$Grain_Yield ~lp(OlNgarua_5$Nrate, nn = 0.80, deg = 1))
 lines(plot14, col="Orange", lwd=1.7)



 #Adding legends 
 #clip(0, 5, 0, 10) 
 par(xpd=TRUE)# To turn off clip
 legend(x=50,y=8.4, c("0P with 0,25,50,75N","25P with 0,25,50,75N",
                      "50P,with 0,25,50,75N","75P with 0,25,50,75N"), 
                        text.font=0.5,cex=0.8,
                      pch=c(19,19,19,19), lty=c(1,1,1,1),
                      col=c("black","purple","red","orange"))

 #Pch- plot symbols as in scatter plot
 #Other legend symbols
 # xpd = TRUE tells R that it is OK to plot outside the region horiz = TRUE
 # tells R that I want a horizontal legend inset = c(x,y) tells R how to move
 # the legend relative to the 'bottom' location bty = 'n' means that 'no' box
 # will be drawn around it pch and col are the types and colors of points cex
 # = 2 makes the legend twice as large as other fonts

 #Anova
 OlNgarua_6=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1")
 View(OlNgarua_6)
 names(OlNgarua_6)
 mod=aov(OlNgarua_6$Grain_Yield~OlNgarua_6$Nrate*OlNgarua_6$Prate)
 summary(mod)
 plot(mod)
 #drop the following combinations 0N,50P;0N,75P;25N,50P;25N,75P
 write.csv(OlNgarua_6, "/Users/machariam/Desktop/OlNgarua_7.csv")


 OlNgarua_7 <- read.csv("C:/Users/machariam/Desktop/OlNgarua_7.csv")
 mod4=aov(OlNgarua_7$Grain_Yield~OlNgarua_7$Nrate*OlNgarua_7$Prate)
 summary(mod4)

 ###################################################################################
 #Combined season s1M1 and S2M2
 OlNgaruaS1_2 <- read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgaruaS1_2.csv")
 View(OlNgaruaS1_2)
 names(OlNgaruaS1_2 )
 #Subset data to SIMI and consider year 1988 with the different P levels
 OlNgaruaA=subset(OlNgaruaS1_2,Year=="1988"& Prate==0)
 OlNgaruaB=subset(OlNgaruaS1_2,Year=="1988"  & Prate==25)
 OlNgaruaC=subset(OlNgaruaS1_2,Year=="1988" & Prate==50)
 OlNgaruaD=subset(OlNgaruaS1_2,Year=="1988" & Prate==75)
 #Investigate the individual plots side by side
 par(mfrow=c(2,2))
 plotA=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield)#P=0
 plotB=plot(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield)#P=25
 plotC=plot(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield)#P=50
 plotD=plot(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield)#P=75

 #combining the above individual scatter plots
 plotE=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield,pch= 19,xlim=c(0,85),
           ylim=c(5.5,12), xlab="N Rate (Kg/ha)"
           ,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
           ,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")
 abline(v=25,lty=5)

 axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.6)
 axis(2, at = seq(5.5,12, by = 0.5), las=1,cex.axis=0.6)

 #las use las to rotate the labels 
 ##yaxt gives the to overwrite the axis labels

 plotF=points(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield,pch= 19, col="red")#P=25
 plotG=points(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield, pch= 19,col="purple")#P=50
 plotH=points(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield, pch= 19,col="Orange")#P=75

 # Fitting lowess curves on the above points
 plotI=locfit(OlNgaruaA$Grain_Yield ~lp(OlNgaruaA$Nrate, nn = 0.80, deg = 1))#0P
 lines(plotI,lwd=1.7)
 plotJ=locfit(OlNgaruaB$Grain_Yield ~lp(OlNgaruaB$Nrate, nn = 0.80, deg = 1))#25P
 lines(plotJ, col="red", lwd=1.7)
 plotK=locfit(OlNgaruaC$Grain_Yield ~lp(OlNgaruaC$Nrate, nn = 0.80, deg = 1))#50P
 lines(plotK, col="purple", lwd=1.7)
 plotL=locfit(OlNgaruaD$Grain_Yield ~lp(OlNgaruaD$Nrate, nn = 0.80, deg = 1))#75P
 lines(plotL, col="Orange", lwd=1.7)



 #Adding legends 
 #clip(0, 5, 0, 10) 
 par(xpd=TRUE)# To turn off clip
 legend(x=50,y=12, c("0P with 0,25,50,75N","25P with 0,25,50,75N",
                     "50P,with 0,25,50,75N","75P with 0,25,50,75N"), 
       text.font=0.5,cex=0.8,
       pch=c(19,19,19,19), lty=c(1,1,1,1),
       col=c("black","purple","red","orange"))

 #Pch- plot symbols as in scatter plot
 #Other legend symbols
 # xpd = TRUE tells R that it is OK to plot outside the region horiz = TRUE
 # tells R that I want a horizontal legend inset = c(x,y) tells R how to move
 # the legend relative to the 'bottom' location bty = 'n' means that 'no' box
 # will be drawn around it pch and col are the types and colors of points cex
 # = 2 makes the legend twice as large as other fonts

 #Anova
 View(OlNgaruaS1_2)
 mod=aov(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate*OlNgaruaS1_2$Prate)
 summary(mod)
 plot(mod)
 #drop the following combinations 0N,50P;0N,75P;25N,50P;25N,75P
 write.csv(OlNgarua_6, "/Users/machariam/Desktop/OlNgarua_7.csv")


 OlNgarua_7 <- read.csv("C:/Users/machariam/Desktop/OlNgarua_7.csv")
 mod4=aov(OlNgarua_7$Grain_Yield~OlNgarua_7$Nrate*OlNgarua_7$Prate)
 summary(mod4)

 ####################################################################################
 ################################1989################################################
 ####################################################################################
 #Combined season s1M1 and S2M2
 OlNgaruaS1_2 <- read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgaruaS1_2.csv")
 View(OlNgaruaS1_2)
 names(OlNgaruaS1_2)
 #Subset data to SIMI and consider year 1988 with the different P levels
 OlNgaruaA=subset(OlNgaruaS1_2,Year=="1989"& Prate==0)
 View(OlNgaruaA)
 OlNgaruaB=subset(OlNgaruaS1_2,Year=="1989"  & Prate==25)
 OlNgaruaC=subset(OlNgaruaS1_2,Year=="1989" & Prate==50)
 OlNgaruaD=subset(OlNgaruaS1_2,Year=="1989" & Prate==75)
 View(OlNgaruaD)
 #Investigate the individual plots side by side
 par(mfrow=c(2,2))
 plotA=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield)#P=0
 plotB=plot(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield)#P=25
 plotC=plot(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield)#P=50
 plotD=plot(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield)#P=75

 #combining the above individual scatter plots
 plotE=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield,pch= 19,xlim=c(0,85),
           ylim=c(3,8), xlab="N Rate (Kg/ha)"
           ,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
           ,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")
 #abline(v=25,lty=5)

 axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.6)
 axis(2, at = seq(3,8, by = 0.5), las=1,cex.axis=0.6)

 #las use las to rotate the labels 
 ##yaxt gives the to overwrite the axis labels

 plotF=points(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield,pch= 19, col="red")#P=25
 plotG=points(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield, pch= 19,col="purple")#P=50
 plotH=points(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield, pch= 19,col="Orange")#P=75

 # Fitting lowess curves on the above points
 plotI=locfit(OlNgaruaA$Grain_Yield ~lp(OlNgaruaA$Nrate, nn = 0.80, deg = 1))#0P
 lines(plotI,lwd=1.7)
 plotJ=locfit(OlNgaruaB$Grain_Yield ~lp(OlNgaruaB$Nrate, nn = 0.80, deg = 1))#25P
 lines(plotJ, col="red", lwd=1.7)
 plotK=locfit(OlNgaruaC$Grain_Yield ~lp(OlNgaruaC$Nrate, nn = 0.80, deg = 1))#50P
 lines(plotK, col="purple", lwd=1.7)
 plotL=locfit(OlNgaruaD$Grain_Yield ~lp(OlNgaruaD$Nrate, nn = 0.80, deg = 1))#75P
 lines(plotL, col="Orange", lwd=1.7)



 #Adding legends 
 #clip(0, 5, 0, 10) 
 par(xpd=TRUE)# To turn off clip
 legend(x=34,y=4.5, c("0P with 0,25,50,75N","25P with 0,25,50,75N",
                    "50P,with 0,25,50,75N","75P with 0,25,50,75N"), 
       text.font=0.5,cex=0.8,
       pch=c(19,19,19,19), lty=c(1,1,1,1),
       col=c("black","purple","red","orange"))

 #Pch- plot symbols as in scatter plot
 #Other legend symbols
 # xpd = TRUE tells R that it is OK to plot outside the region horiz = TRUE
 # tells R that I want a horizontal legend inset = c(x,y) tells R how to move
 # the legend relative to the 'bottom' location bty = 'n' means that 'no' box
 # will be drawn around it pch and col are the types and colors of points cex
 # = 2 makes the legend twice as large as other fonts

 #Anova
 View(OlNgaruaS1_2)
 mod=aov(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate*OlNgaruaS1_2$Prate)
 summary(mod)
 plot(mod)

 #######################################################################
 OlNgaruaS1_2 <- read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgaruaS1_2.csv")
 View(OlNgaruaS1_2)
 plotM=plot(OlNgaruaS1_2$Grain_Yield ~ OlNgaruaS1_2$Nrate)
 plotN=plot(OlNgaruaS1_2$Grain_Yield ~ OlNgaruaS1_2$Nrate,pch= 19,xlim=c(0,85),
           ylim=c(4,9), xlab="N Rate (Kg/ha)"
           ,ylab="Grain yield (t/ha)",main="Maize yield response to N"
           ,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")


 axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.6)
 axis(2, at = seq(4,9, by = 0.5), las=1,cex.axis=0.6)
 plotN=locfit(OlNgaruaS1_2$Grain_Yield ~lp(OlNgaruaS1_2$Nrate, nn = 0.80, deg = 1))#0P
 lines(plotN,lwd=1.7)
 summary(OlNgaruaS1_2$Grain_Yield)
 #Anova and linear model
 names(OlNgaruaS1_2)
 mod1=aov(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate*OlNgaruaS1_2$Prate)
 summary(mod1)
 mod2=lm(mod1)
 summary(mod2)
 str(OlNgaruaS1_2)
 mod3=aov(mod2)
 summary(mod3)
 mod4=lm(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate)
 summary(mod4)
 ############################################################################
 #######fitting a spillmans curve###########################
 #y=M-AR^X
 write.csv(OlNgaruaS1_2, "/Users/machariam/Desktop/OlNgaruaSP2.csv")

 #fit an NlS of the form;y=M-AR^X
 #A = increase in output due to X (i.e. max. Y – min. Y) 
 #Theoretical maximum increase in yield obtainable by increasing X indefinitely
 #R = ratio of successive increments in output to total output in which a zero
 #level of X produces no Y, M= A and the spillman devolve to; 
 #The ratio of a decreasing geometric series the terms of which are the respective
 #increments in yield due to successive unit increments in X 
 #M = maximum total production obtained by the use of resource X (It is the limit 
 #approached by Y as X increases indefinitely or the theoretical maximum yield 
 #possible with any number of units of the growth factor)
 mod5=nls(Yield~M-A*(R^Nrate), start=list(M=5.5, A=3,R=0.1),data=OlNgaruaSP2)
 mod5
 #Make linear
 mod6=nls(Yield~(-A)*(R^Nrate)+ M, start=list(M=5.5, A=3,R=0.1),data=OlNgaruaSP2)
 mod6
 #Use nls2
 # The challage is trying to fit 3 unknown parameters from a dataset with 5 points
 #The results make no sense if the top plateau of a sigmoid curve is far larger 
 #than the highest data point, the same applies to the theoretical maximum increase in yield,
 #due the applied input(s)
 require(nlmrt)
 mod7=nlxb(Yield~M-A*(R^Nrate), start=list(M=5.5, A=3,R=1),data=OlNgaruaSP2)
 mod7

 ####################################### Adding a line of best based on a model###############
 #par(mfrow=c(1,2))
 mod8=nls(Yld~M-A*(R^Prate), start=list(M=800, A=190,R=0.8),data=Book2)
 mod8
 summary(mod8)
 coef(mod8)
 plot(Yld~Prate, ylab="Yield (kg)", xlab="P input (kg)", Book2)
 M <-1243.9270758
 A <-584.6197832
 R <-0.9865372
 lines(x<-c(0:50),M-A*(R^x),col='red')
 ####################################An example of a dataset that fits a spillmans function#####
 ###################################
 x <- c(0,10,20,30,40,50)
 y <- c(2748,3162,3307,3571,3576,3544)
 dat <- as.data.frame(cbind(x,y))
 dat
 plot(dat)

 mod9=nls(y~M-A*(R^x), start=list(M=800, A=190,R=0.8),data=dat)
 mod9
 summary(mod9)
 coef(mod9)

 plot(y~x, ylab="Yield (kg)", xlab="P input (kg)", dat)
 M <-3628.7228167
 A <-883.4230643
 R <-0.9398198
 lines(x<-c(0:50),M-A*(R^x),col='red')
	#1988

	#Import data
	OlNgarua=read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgarua.csv")
	OlNgarua
	names(OlNgarua)
	str(OlNgarua)
	OlNgarua_1=OlNgarua[ ,c("District","Center","Longitudes","Latitudes","Year","Crop",
	"Variety","PlantMon","AppliMethod","Nrate","Prate","Grain_Yield")]
	summary(OlNgarua_1)
	str(OlNgarua_1)
	#Subset data to SIMI and consider year 1988 with the different P levels
	#P=0
	OlNgarua_2=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1" & Prate==0)
	View(OlNgarua_2)
	#P=25
	OlNgarua_3=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1" & Prate==25)
	View(OlNgarua_3)
	#P=50
	OlNgarua_4=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1" & Prate==50)
	View(OlNgarua_4)
	#P=75
	OlNgarua_5=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1" & Prate==75)
	View(OlNgarua_5)

	#Basic Plots

	plot=plot(OlNgarua_2$Nrate,OlNgarua_2$Grain_Yield, xlim=c(0,100), ylim=c(8,12)
	,xlab="Nitrogen Rate", ylab="Yield(t/ha)", pch= 19,col="Red",cex=0.9, cex.lab=1.2)

	#Joining points with a line
	plot1=lowess(OlNgarua_2$Grain_Yield ~ OlNgarua_2$Nrate)
	lines(smoothplot1, col = 'red')
	#Making the line smooth
	library(locfit)
	#Locfit is a software package performing local regression,likelihood
	#and related smoothing procedures.
	plot2=locfit(OlNgarua_2$Grain_Yield ~lp(OlNgarua_2$Nrate, nn = 0.80, deg = 1))
	lines(plot2, col="red", lwd=1.4)
	#nn=specifies smoothness of parameter
	#deg=specifies the degree (1,2...) of the polynomial

	#Investigate the individual plots side by side
	par(mfrow=c(2,2))
	plot3=plot(OlNgarua_2$Nrate,OlNgarua_2$Grain_Yield)#P=0
	plot4=plot(OlNgarua_3$Nrate,OlNgarua_3$Grain_Yield)#P=25
	plot5=plot(OlNgarua_4$Nrate,OlNgarua_4$Grain_Yield)#P=50
	plot6=plot(OlNgarua_5$Nrate,OlNgarua_5$Grain_Yield)#P=75

	#combining the above individual scatter plots
	plot7=plot(OlNgarua_2$Nrate,OlNgarua_2$Grain_Yield,pch= 19,xlim=c(0,85),
	ylim=c(6.5,12), xlab="N Rate (Kg/ha)"
	,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
	,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")
	abline(v=25,lty=5)

	axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.7)
	axis(2, at = seq(6.5,12, by = 0.5), las=1,cex.axis=0.7)

	#las use las to rotate the labels
	##yaxt gives the to overwrite the axis labels

	plot8=points(OlNgarua_3$Nrate,OlNgarua_3$Grain_Yield,pch= 19, col="red")#P=25
	plot9=points(OlNgarua_4$Nrate,OlNgarua_4$Grain_Yield, pch= 19,col="purple")#P=50
	plot10=points(OlNgarua_5$Nrate,OlNgarua_5$Grain_Yield, pch= 19,col="Orange")#P=75

	# Fitting lowess curves on the above points
	#0P
	plot11=locfit(OlNgarua_2$Grain_Yield ~lp(OlNgarua_2$Nrate, nn = 0.80, deg = 1))
	lines(plot11,lwd=1.7)
	#25P
	plot12=locfit(OlNgarua_3$Grain_Yield ~lp(OlNgarua_3$Nrate, nn = 0.80, deg = 1))
	lines(plot12, col="red", lwd=1.7)
	#50P
	plot13=locfit(OlNgarua_4$Grain_Yield ~lp(OlNgarua_4$Nrate, nn = 0.80, deg = 1))
	lines(plot13, col="purple", lwd=1.7)
	#75P
	plot14=locfit(OlNgarua_5$Grain_Yield ~lp(OlNgarua_5$Nrate, nn = 0.80, deg = 1))
	lines(plot14, col="Orange", lwd=1.7)



	#Adding legends
	#clip(0, 5, 0, 10)
	par(xpd=TRUE)# To turn off clip
	legend(x=50,y=8.4, c("0P with 0,25,50,75N","25P with 0,25,50,75N",
	"50P,with 0,25,50,75N","75P with 0,25,50,75N"),
	text.font=0.5,cex=0.8,
	pch=c(19,19,19,19), lty=c(1,1,1,1),
	col=c("black","purple","red","orange"))

	#Pch- plot symbols as in scatter plot
	#Other legend symbols
	# xpd = TRUE tells R that it is OK to plot outside the region horiz = TRUE
	# tells R that I want a horizontal legend inset = c(x,y) tells R how to move
	# the legend relative to the 'bottom' location bty = 'n' means that 'no' box
	# will be drawn around it pch and col are the types and colors of points cex
	# = 2 makes the legend twice as large as other fonts

	#Anova
	OlNgarua_6=subset(OlNgarua_1,Year=="1988" & PlantMon=="S1M1")
	View(OlNgarua_6)
	names(OlNgarua_6)
	mod=aov(OlNgarua_6$Grain_Yield~OlNgarua_6$Nrate*OlNgarua_6$Prate)
	summary(mod)
	plot(mod)
	#drop the following combinations 0N,50P;0N,75P;25N,50P;25N,75P
	write.csv(OlNgarua_6, "/Users/machariam/Desktop/OlNgarua_7.csv")


	OlNgarua_7 <- read.csv("C:/Users/machariam/Desktop/OlNgarua_7.csv")
	mod4=aov(OlNgarua_7$Grain_Yield~OlNgarua_7$Nrate*OlNgarua_7$Prate)
	summary(mod4)

	###################################################################################
	#Combined season s1M1 and S2M2
	OlNgaruaS1_2 <- read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgaruaS1_2.csv")
	View(OlNgaruaS1_2)
	names(OlNgaruaS1_2 )
	#Subset data to SIMI and consider year 1988 with the different P levels
	OlNgaruaA=subset(OlNgaruaS1_2,Year=="1988"& Prate==0)
	OlNgaruaB=subset(OlNgaruaS1_2,Year=="1988" & Prate==25)
	OlNgaruaC=subset(OlNgaruaS1_2,Year=="1988" & Prate==50)
	OlNgaruaD=subset(OlNgaruaS1_2,Year=="1988" & Prate==75)
	#Investigate the individual plots side by side
	par(mfrow=c(2,2))
	plotA=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield)#P=0
	plotB=plot(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield)#P=25
	plotC=plot(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield)#P=50
	plotD=plot(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield)#P=75

	#combining the above individual scatter plots
	plotE=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield,pch= 19,xlim=c(0,85),
	ylim=c(5.5,12), xlab="N Rate (Kg/ha)"
	,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
	,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")
	abline(v=25,lty=5)

	axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.6)
	axis(2, at = seq(5.5,12, by = 0.5), las=1,cex.axis=0.6)

	#las use las to rotate the labels
	##yaxt gives the to overwrite the axis labels

	plotF=points(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield,pch= 19, col="red")#P=25
	plotG=points(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield, pch= 19,col="purple")#P=50
	plotH=points(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield, pch= 19,col="Orange")#P=75

	# Fitting lowess curves on the above points
	plotI=locfit(OlNgaruaA$Grain_Yield ~lp(OlNgaruaA$Nrate, nn = 0.80, deg = 1))#0P
	lines(plotI,lwd=1.7)
	plotJ=locfit(OlNgaruaB$Grain_Yield ~lp(OlNgaruaB$Nrate, nn = 0.80, deg = 1))#25P
	lines(plotJ, col="red", lwd=1.7)
	plotK=locfit(OlNgaruaC$Grain_Yield ~lp(OlNgaruaC$Nrate, nn = 0.80, deg = 1))#50P
	lines(plotK, col="purple", lwd=1.7)
	plotL=locfit(OlNgaruaD$Grain_Yield ~lp(OlNgaruaD$Nrate, nn = 0.80, deg = 1))#75P
	lines(plotL, col="Orange", lwd=1.7)



	#Adding legends
	#clip(0, 5, 0, 10)
	par(xpd=TRUE)# To turn off clip
	legend(x=50,y=12, c("0P with 0,25,50,75N","25P with 0,25,50,75N",
	"50P,with 0,25,50,75N","75P with 0,25,50,75N"),
	text.font=0.5,cex=0.8,
	pch=c(19,19,19,19), lty=c(1,1,1,1),
	col=c("black","purple","red","orange"))

	#Pch- plot symbols as in scatter plot
	#Other legend symbols
	# xpd = TRUE tells R that it is OK to plot outside the region horiz = TRUE
	# tells R that I want a horizontal legend inset = c(x,y) tells R how to move
	# the legend relative to the 'bottom' location bty = 'n' means that 'no' box
	# will be drawn around it pch and col are the types and colors of points cex
	# = 2 makes the legend twice as large as other fonts

	#Anova
	View(OlNgaruaS1_2)
	mod=aov(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate*OlNgaruaS1_2$Prate)
	summary(mod)
	plot(mod)
	#drop the following combinations 0N,50P;0N,75P;25N,50P;25N,75P
	write.csv(OlNgarua_6, "/Users/machariam/Desktop/OlNgarua_7.csv")


	OlNgarua_7 <- read.csv("C:/Users/machariam/Desktop/OlNgarua_7.csv")
	mod4=aov(OlNgarua_7$Grain_Yield~OlNgarua_7$Nrate*OlNgarua_7$Prate)
	summary(mod4)

	####################################################################################
	################################1989################################################
	####################################################################################
	#Combined season s1M1 and S2M2
	OlNgaruaS1_2 <- read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgaruaS1_2.csv")
	View(OlNgaruaS1_2)
	names(OlNgaruaS1_2)
	#Subset data to SIMI and consider year 1988 with the different P levels
	OlNgaruaA=subset(OlNgaruaS1_2,Year=="1989"& Prate==0)
	View(OlNgaruaA)
	OlNgaruaB=subset(OlNgaruaS1_2,Year=="1989" & Prate==25)
	OlNgaruaC=subset(OlNgaruaS1_2,Year=="1989" & Prate==50)
	OlNgaruaD=subset(OlNgaruaS1_2,Year=="1989" & Prate==75)
	View(OlNgaruaD)
	#Investigate the individual plots side by side
	par(mfrow=c(2,2))
	plotA=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield)#P=0
	plotB=plot(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield)#P=25
	plotC=plot(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield)#P=50
	plotD=plot(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield)#P=75

	#combining the above individual scatter plots
	plotE=plot(OlNgaruaA$Nrate,OlNgaruaA$Grain_Yield,pch= 19,xlim=c(0,85),
	ylim=c(3,8), xlab="N Rate (Kg/ha)"
	,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
	,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")
	#abline(v=25,lty=5)

	axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.6)
	axis(2, at = seq(3,8, by = 0.5), las=1,cex.axis=0.6)

	#las use las to rotate the labels
	##yaxt gives the to overwrite the axis labels

	plotF=points(OlNgaruaB$Nrate,OlNgaruaB$Grain_Yield,pch= 19, col="red")#P=25
	plotG=points(OlNgaruaC$Nrate,OlNgaruaC$Grain_Yield, pch= 19,col="purple")#P=50
	plotH=points(OlNgaruaD$Nrate,OlNgaruaD$Grain_Yield, pch= 19,col="Orange")#P=75

	# Fitting lowess curves on the above points
	plotI=locfit(OlNgaruaA$Grain_Yield ~lp(OlNgaruaA$Nrate, nn = 0.80, deg = 1))#0P
	lines(plotI,lwd=1.7)
	plotJ=locfit(OlNgaruaB$Grain_Yield ~lp(OlNgaruaB$Nrate, nn = 0.80, deg = 1))#25P
	lines(plotJ, col="red", lwd=1.7)
	plotK=locfit(OlNgaruaC$Grain_Yield ~lp(OlNgaruaC$Nrate, nn = 0.80, deg = 1))#50P
	lines(plotK, col="purple", lwd=1.7)
	plotL=locfit(OlNgaruaD$Grain_Yield ~lp(OlNgaruaD$Nrate, nn = 0.80, deg = 1))#75P
	lines(plotL, col="Orange", lwd=1.7)



	#Adding legends
	#clip(0, 5, 0, 10)
	par(xpd=TRUE)# To turn off clip
	legend(x=34,y=4.5, c("0P with 0,25,50,75N","25P with 0,25,50,75N",
	"50P,with 0,25,50,75N","75P with 0,25,50,75N"),
	text.font=0.5,cex=0.8,
	pch=c(19,19,19,19), lty=c(1,1,1,1),
	col=c("black","purple","red","orange"))

	#Pch- plot symbols as in scatter plot
	#Other legend symbols
	# xpd = TRUE tells R that it is OK to plot outside the region horiz = TRUE
	# tells R that I want a horizontal legend inset = c(x,y) tells R how to move
	# the legend relative to the 'bottom' location bty = 'n' means that 'no' box
	# will be drawn around it pch and col are the types and colors of points cex
	# = 2 makes the legend twice as large as other fonts

	#Anova
	View(OlNgaruaS1_2)
	mod=aov(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate*OlNgaruaS1_2$Prate)
	summary(mod)
	plot(mod)

	#######################################################################
	OlNgaruaS1_2 <- read.csv("C:/Users/machariam/Desktop/Legacy_data_Analysis/OlNgaruaS1_2.csv")
	View(OlNgaruaS1_2)
	plotM=plot(OlNgaruaS1_2$Grain_Yield ~ OlNgaruaS1_2$Nrate)
	plotN=plot(OlNgaruaS1_2$Grain_Yield ~ OlNgaruaS1_2$Nrate,pch= 19,xlim=c(0,85),
	ylim=c(4,9), xlab="N Rate (Kg/ha)"
	,ylab="Grain yield (t/ha)",main="Maize yield response to N"
	,cex.main=0.9,cex.lab=0.8,xaxt="n", yaxt="n")


	axis(1, at = seq(0, 85, by = 5), las=1,cex.axis=0.6)
	axis(2, at = seq(4,9, by = 0.5), las=1,cex.axis=0.6)
	plotN=locfit(OlNgaruaS1_2$Grain_Yield ~lp(OlNgaruaS1_2$Nrate, nn = 0.80, deg = 1))#0P
	lines(plotN,lwd=1.7)
	summary(OlNgaruaS1_2$Grain_Yield)
	#Anova and linear model
	names(OlNgaruaS1_2)
	mod1=aov(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate*OlNgaruaS1_2$Prate)
	summary(mod1)
	mod2=lm(mod1)
	summary(mod2)
	str(OlNgaruaS1_2)
	mod3=aov(mod2)
	summary(mod3)
	mod4=lm(OlNgaruaS1_2$Grain_Yield~OlNgaruaS1_2$Nrate)
	summary(mod4)
	############################################################################
	#######fitting a spillmans curve###########################
	#y=M-AR^X
	write.csv(OlNgaruaS1_2, "/Users/machariam/Desktop/OlNgaruaSP2.csv")

	#fit an NlS of the form;y=M-AR^X
	#A = increase in output due to X (i.e. max. Y – min. Y)
	#Theoretical maximum increase in yield obtainable by increasing X indefinitely
	#R = ratio of successive increments in output to total output in which a zero
	#level of X produces no Y, M= A and the spillman devolve to;
	#The ratio of a decreasing geometric series the terms of which are the respective
	#increments in yield due to successive unit increments in X
	#M = maximum total production obtained by the use of resource X (It is the limit
	#approached by Y as X increases indefinitely or the theoretical maximum yield
	#possible with any number of units of the growth factor)
	mod5=nls(Yield~M-A*(R^Nrate), start=list(M=5.5, A=3,R=0.1),data=OlNgaruaSP2)
	mod5
	#Make linear
	mod6=nls(Yield~(-A)*(R^Nrate)+ M, start=list(M=5.5, A=3,R=0.1),data=OlNgaruaSP2)
	mod6
	#Use nls2
	# The challage is trying to fit 3 unknown parameters from a dataset with 5 points
	#The results make no sense if the top plateau of a sigmoid curve is far larger
	#than the highest data point, the same applies to the theoretical maximum increase in yield,
	#due the applied input(s)
	require(nlmrt)
	mod7=nlxb(Yield~M-A*(R^Nrate), start=list(M=5.5, A=3,R=1),data=OlNgaruaSP2)
	mod7

	####################################### Adding a line of best based on a model###############
	#par(mfrow=c(1,2))
	mod8=nls(Yld~M-A*(R^Prate), start=list(M=800, A=190,R=0.8),data=Book2)
	mod8
	summary(mod8)
	coef(mod8)
	plot(Yld~Prate, ylab="Yield (kg)", xlab="P input (kg)", Book2)
	M <-1243.9270758
	A <-584.6197832
	R <-0.9865372
	lines(x<-c(0:50),M-A*(R^x),col='red')
	####################################An example of a dataset that fits a spillmans function#####
	###################################
	x <- c(0,10,20,30,40,50)
	y <- c(2748,3162,3307,3571,3576,3544)
	dat <- as.data.frame(cbind(x,y))
	dat
	plot(dat)

	mod9=nls(y~M-A*(R^x), start=list(M=800, A=190,R=0.8),data=dat)
	mod9
	summary(mod9)
	coef(mod9)

	plot(y~x, ylab="Yield (kg)", xlab="P input (kg)", dat)
	M <-3628.7228167
	A <-883.4230643
	R <-0.9398198
	lines(x<-c(0:50),M-A*(R^x),col='red')