Skip to content

Instantly share code, notes, and snippets.

@MartinMacharia
Last active December 23, 2015 10:50
Show Gist options
  • Save MartinMacharia/c6cc27b929488c11754b to your computer and use it in GitHub Desktop.
Save MartinMacharia/c6cc27b929488c11754b to your computer and use it in GitHub Desktop.
Kenya Legacy analysis
Kenya <- read.csv("C:/Users/machariam/Desktop/Spatial_layers/Kenya.csv", na.strings=c("NA"))
names(Kenya)
#create a new dataset with less variable of interst
str(Kenya)
Kenya_1=Kenya[ , c("Country","District","Center","Site","Longitudes","Latitudes","Elevation","Year","Crop_1","Variety_1",
"Crop_2","Crop_2_Variety","Crop_3" ,"Crop_3_Variety","PlantMon","PrevCrop","AppliMethod",
"Manure","Inorganic", "Crop_Fallow","pH_water","PH_Other","Soil_type","Nrate","Prate","Krate",
"Srate","Znrate","Curate","MgRate","Brate","Other_Nutrients","FYM","Compost","Grain_Yield",
"Stover_yield", "Additional_comments")]
summary(Kenya_1)
str(Kenya_1)
#Subset the data by district
Nakuru=Kenya_1[Kenya_1$District== "Nakuru", ]
summary(Nakuru)
#Printing the data in excel
write.csv(Nakuru, "/Users/machariam/Desktop/Spatial_layers/Kenya/Nakuru.csv")
#Subset Nakuru data into sites
OlNgarua = Nakuru[Nakuru$Center=="OlNgarua",]
summary(OlNgarua)
Grain_Yield=as.numeric(Grain_Yield)
#write.csv(OlNgarua, "/Users/machariam/Desktop/Spatial_layers/Kenya/OlNgarua.csv")
#OlNgarua site/center has the same treatment structure of 5 N levels across seasons and years
# Are season and year effects significant?
plot(OlNgarua$Nrate,OlNgarua$Grain_Yield)
#Too much variability within an N level, hence investigate by year across crops
OlNgarua_year=OlNgarua[OlNgarua$Year=="1990",]
summary(OlNgarua_year)
plot(OlNgarua_year$Nrate,OlNgarua_year$Grain_Yield)
# Per single crop
OlNgarua_maize=OlNgarua[OlNgarua$Crop_1=="Maize",]
plot(OlNgarua_maize$Nrate,OlNgarua_maize$Grain_Yield)
# subset by crop and year
OlNgarua_maize_1990=subset(OlNgarua, Crop_1=="Maize" & Year=="1990" & PlantMon=="S1M1" & Variety_1 == "H625" & Prate==50 )
OlNgarua_maize_1990
write.csv(OlNgarua_maize_1990, "/Users/machariam/Desktop/Spatial_layers/Kenya/OlNgarua_maize_1990.csv")
plot(OlNgarua_maize_1990$Nrate,OlNgarua_maize_1990$Grain_Yield)
#Overlying the scatter plots genearted by country and crop and nutrient basis
#add the argument type = "l" to get a line
plot=plot(OlNgarua_maize_1990$Nrate,OlNgarua_maize_1990$Grain_Yield, xlim=c(0,100), ylim=c(1000,2000)
,xlab="Nitrogen Rate", ylab="Yield(t/ha)", pch= 19,col="Red",cex=0.9, cex.lab=1.2)
Grain_Yield_num=as.numeric(OlNgarua_maize_1990$Grain_Yield)
str(OlNgarua_maize_1990$Grain_Yield)
str(Grain_Yield_num)
#Adding a straight line
plot1=lm(Grain_Yield_num ~ OlNgarua_maize_1990$Nrate)
abline(plot1, col = 'red')
#Adding a smooth line
smoothplot1=lowess(Grain_Yield_num ~ OlNgarua_maize_1990$Nrate)
lines(smoothplot1, col = 'red')
library(locfit)
smoothplot1=locfit(Grain_Yield_num ~lp(OlNgarua_maize_1990$Nrate, nn = 0.80, deg = 1))
lines(smoothplot1, col="red", lwd=1.4)
#nn=specifies smoothness of parameter
#deg=specifies the degree (1,2...) of the polynomial
#overlying the scatter plots
str(OlNgarua)
par(mfrow=c(2,2))
Dataplot1=subset(OlNgarua, Crop_1=="Maize" & Year=="1990" & PlantMon=="S1M1"
& Variety_1 == "H625" & Prate==0)
Dataplot1
plot1=plot(Dataplot1$Nrate,Dataplot1$Grain_Yield)
Dataplot2=subset(OlNgarua, Crop_1=="Maize" & Year=="1990" &
PlantMon=="S1M1" & Variety_1 == "H625" & Prate==25 )
Dataplot2
plot2=plot(Dataplot2$Nrate,Dataplot2$Grain_Yield)
Dataplot3=subset(OlNgarua, Crop_1=="Maize" & Year=="1990" &
PlantMon=="S1M1" & Variety_1 == "H625" & Prate==50)
Dataplot3
plot3=plot(Dataplot3$Nrate,Dataplot3$Grain_Yield,xlim=c(0,100), ylim=c(1000,2000))
Dataplot4=subset(OlNgarua, Crop_1=="Maize" & Year=="1990" &
PlantMon=="S1M1" & Variety_1 == "H625" & Prate==75)
Dataplot4
plot4=plot(Dataplot4$Nrate,Dataplot4$Grain_Yield, xlim=c(0,100), ylim=c(1000,2000))
#combining the above individual scatter plots
plot1=plot(Dataplot1$Nrate,Dataplot1$Grain_Yield,pch= 19,xlim=c(0,85),
#axis(1, at = seq(0, 85, by = 10), las=2),
ylim=c(1000,1600), xlab="N Rate (Kg/ha)"
,ylab="Grain yield (t/ha)",main="Maize yield response to N & P fertlizers"
,cex.main=0.9,cex.axis=0.7, abline(v=25,lty=5))
plot2=points(Dataplot2$Nrate,Dataplot2$Grain_Yield,pch= 19, col="red")
plot3=points(Dataplot3$Nrate,Dataplot3$Grain_Yield, pch= 19,col="purple")
plot4=points(Dataplot4$Nrate,Dataplot4$Grain_Yield, pch= 19,col="Orange")
# Fitting lowess curves on the above points and draving lines
#0P
smoothplot2=locfit(Dataplot1$Grain_Yield ~lp(Dataplot1$Nrate, nn = 0.80, deg = 1))
lines(smoothplot2,lwd=1.7)
#25P
smoothplot3=locfit(Dataplot2$Grain_Yield ~lp(Dataplot2$Nrate, nn = 0.80, deg = 1))
lines(smoothplot3, col="red", lwd=1.7)
#50P
smoothplot4=locfit(Dataplot3$Grain_Yield ~lp(Dataplot3$Nrate, nn = 0.80, deg = 1))
lines(smoothplot4, col="purple", lwd=1.7)
#75P
smoothplot5=locfit(Dataplot2$Grain_Yield ~lp(Dataplot2$Nrate, nn = 0.80, deg = 2))
lines(smoothplot5, col="Orange", lwd=1.7)
#Adding legends
legend(x=40,y=1200, 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","red","purple","orange"))
#Pch- plot symbols as in scatter plot
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment