require("lme4") require("lmerTest") require("ggplot2") require("plyr") require("dplyr") require(tidyverse) # setwd("/home/akilimo/projects/SoilNPK") #reading in and preparing the data dsY <- read.csv("NOT_Sep2020.csv") ## this data comes from NOT_DataAnalysis_July2020.R in D:\ACAI_Wrapper\cloud_compute\LINTUL dsY <- dsY[rowSums(is.na(dsY)) <= ncol(dsY),] names(dsY)[names(dsY)=="rootY_tha"] <- "rootY" dsY$treat <- factor(dsY$treat, levels = c("NPK_micro", "NPK", "PK", "NK", "NP", "half_NPK", "CON")) head(dsY) ## count number treatplot <- unique(dsY[, c("trialID", "plotID","treat")]) Q <- as.data.frame(table(treatplot$treat, treatplot$trialID)) QNPK <- Q[Q$Var1 == "NPK", ] QNPKnot <- Q[!Q$Var1 == "NPK", ] QNPK[QNPK$Freq > 2, ] QNPKnot[QNPKnot$Freq > 1, ] ###################################################################################################################################### #### plots with double yield data, slect the original plotID when compared to replaced or from plotconfirm. keep tials with 4 plots/treat ## script used is in setwd("D:/ACAI_Wrapper/Rwanda") ###################################################################################################################################### ### getting submission id for root yield setwd("C:/Users/Meklit/CGIAR/ACAI - Documents/ODK_Data") rootYield1 <- read.csv("Yield Assessment Cassava _Roots_.csv") rootYield2 <- read.csv("Yield Assessment Cassava _Roots__yieldAssessment.csv") rootYield <- merge(rootYield1, rootYield2, by.x="KEY", by.y="PARENT_KEY") rootYield <- rootYield[, c("SubmissionDate", "today", "plantID", "plotID", "rootFW.tuberizedRootsFW", "rootFW.tuberizedDiseasedRootsFW", "rootFW.tuberizedSmallRootsFW", "rootFW.tuberizedMarketableRootsFW", "tuberizedRootsFWss", "plantSampleID_tuberizedRoots", "tuberizedSmallRootsFWss", "plantSampleID_tuberizedSmallRoots", "tuberizedMarketableRootsFWss", "plantSampleID_tuberizedMarketableRoots", "username")] setwd("C:/Users/Meklit/CGIAR/AKILIMO_serverDownload") ONA_RY1 <- read.csv("Assess_Root_Yield_Cassava_AC.csv") ONA_RY2 <- read.csv("Assess_Root_Yield_Cassava_AC-yieldAssessment.csv") ONA_RY1 <- subset(ONA_RY1, select = -c(uuid,start,deviceid, subscriberid, email, simserial,phonenumber, banner, intro, project, country, login, firstName, surName, geopoint.Altitude, geopoint.Accuracy, generated_table_list_label_23, reserved_name_for_field_list_labels_24, SET.OF.yieldAssessment, end, instanceID)) names(ONA_RY1) <- c("SubmissionDate","today","username","geopoint.Latitude.ONA","geopoint.Longitude.ONA","entity","diseaseScoring","rootQuality", "sampling","fixedSize", "method","densityFixed","nrRowsFixed","nrPlantsRowFixed", "plotDimNote","L1Fixed","W1Fixed","L2Fixed","W2Fixed", "betweenRowFixed", "withinRowFixed", "densityFixedCalc", "maxStandFixed","plotSizeFixed","noteFixed","repeat.","KEY") ONA_RY2 <- ONA_RY2[, c("plantID","plotID","tuberizedRootsFW","tuberizedDiseasedRootsFW", "tuberizedSmallRootsFW","tuberizedMarketableRootsFW", "tuberizedRootsFWss","plantSampleID_tuberizedRoots","tuberizedSmallRootsFWss", "plantSampleID_tuberizedSmallRoots", "tuberizedMarketableRootsFWss","plantSampleID_tuberizedMarketableRoots","PARENT_KEY")] colnames(ONA_RY2) <- c( "plantID","plotID", "rootFW.tuberizedRootsFW","rootFW.tuberizedDiseasedRootsFW", "rootFW.tuberizedSmallRootsFW","rootFW.tuberizedMarketableRootsFW","tuberizedRootsFWss", "plantSampleID_tuberizedRoots", "tuberizedSmallRootsFWss","plantSampleID_tuberizedSmallRoots", "tuberizedMarketableRootsFWss", "plantSampleID_tuberizedMarketableRoots","PARENT_KEY") ONA_RY <- merge(ONA_RY1, ONA_RY2, by.x="KEY", by.y="PARENT_KEY") ONA_RY <- subset(ONA_RY, select=-c(KEY)) length(unique(ONA_RY$plotID))##3462 ONA_rootYield <- ONA_RY[, colnames(rootYield)] Root_Yield_Cassava_ONA_ODK <- unique(rbind(ONA_rootYield, rootYield)) head(Root_Yield_Cassava_ONA_ODK) names(Root_Yield_Cassava_ONA_ODK) <- gsub("rootFW.", "", colnames(Root_Yield_Cassava_ONA_ODK)) Root_Yield_Cassava_ONA_ODK$Yieldsum <- rowSums(subset(Root_Yield_Cassava_ONA_ODK,select=c(tuberizedRootsFW, tuberizedDiseasedRootsFW, tuberizedSmallRootsFW, tuberizedMarketableRootsFW)),na.rm=TRUE) Root_Yield_Cassava_ONA_ODK <- Root_Yield_Cassava_ONA_ODK[, c("SubmissionDate" ,"today" ,"plotID", "Yieldsum")] Root_Yield_Cassava_ONA_ODK <- Root_Yield_Cassava_ONA_ODK[!Root_Yield_Cassava_ONA_ODK$plotID == "", ] setwd("D:/ACAI_sample_processing/Version3/linkedForms") Plot_ODK_ONA <- read.csv("Plot_ODK_ONA.csv") plotconfirm_ODK_ONA <- read.csv("plotconfirm_ODK_ONA.csv") ReplacePlot_ODK_ONA <- read.csv("ReplacePlot_ODK_ONA.csv") ReplaceTrial_ODK_ONA <- read.csv("ReplaceTrial_ODK_ONA.csv") ## adding on rootyield the detsils of plot on the origianlly assigned plotIDs OR_plotYield <- unique(merge(Root_Yield_Cassava_ONA_ODK, Plot_ODK_ONA, by="plotID", all.x = TRUE)) head(OR_plotYield) OR_plotYieldGood <- OR_plotYield[!is.na(OR_plotYield$trialID), ] OR_plotYieldNot <- OR_plotYield[is.na(OR_plotYield$trialID), c("SubmissionDate","today" ,"plotID", "Yieldsum")] ## get plot detail from plotconfirm_ODK_ONA confirm_plotYield <- unique(merge(OR_plotYieldNot, plotconfirm_ODK_ONA, by="plotID", all.x = TRUE)) head(confirm_plotYield) confirm_plotYieldGood <- confirm_plotYield[!is.na(confirm_plotYield$trialID), ] confirm_plotYieldNot <- confirm_plotYield[is.na(confirm_plotYield$trialID), c("SubmissionDate","today" ,"plotID", "Yieldsum")] ## get plot detail from ReplacePlot_ODK_ONA Replace_plotYield <- unique(merge(confirm_plotYieldNot, ReplacePlot_ODK_ONA, by.x="plotID", by.y="newPlotID",all.x = TRUE)) head(Replace_plotYield) Replace_plotYieldGood <- Replace_plotYield[!is.na(Replace_plotYield$trialID),] Replace_plotYieldNot <- Replace_plotYield[is.na(Replace_plotYield$trialID), c("SubmissionDate","today" ,"plotID", "Yieldsum")] names(Replace_plotYieldGood) <- c("plotID","SubmissionDate", "today","Yieldsum","plantLabel", "plantID","trialID","trialCode", "treatCode_label") ### merge the three confirm_plotYieldGood <- confirm_plotYieldGood[, colnames(OR_plotYieldGood)] OR_plotYieldGood$plotID_status <- "Original" confirm_plotYieldGood$plotID_status <- "plotConfirm" Y1 <- rbind(confirm_plotYieldGood, OR_plotYieldGood) Y1 <- subset(Y1, select=-c(treatCode)) Replace_plotYieldGood$plotID_status <- "Replaced" Replace_plotYieldGood <- Replace_plotYieldGood[, colnames(Y1)] PlotYield <- rbind(Y1, Replace_plotYieldGood) PlotYield <- unique(PlotYield) head(PlotYield) PlotYield$treatCode_label <- gsub("NOT-1_", "", PlotYield$treatCode_label) PlotYield$treatCode_label <- gsub("NOT-2_", "", PlotYield$treatCode_label) PlotYield$treatCode_label <- gsub("NOT-3_", "", PlotYield$treatCode_label) unique(PlotYield$treatCode_label) ## does every treatment has only one value ina plot (for NPK two) treatplot <- unique(dsY[, c("trialID", "plotID","treat")]) Q <- as.data.frame(table(treatplot$treat, treatplot$trialID)) QNPK <- Q[Q$Var1 == "NPK", ] QNPKnot <- Q[!Q$Var1 == "NPK", ] length(unique(QNPK[QNPK$Freq > 2, ]$Var2)) ### solving NPK setwd("D:/ACAI_sample_processing/Version3/data") ODK_ONA_ID_GPS <- read.csv("ODK_ONA_ID_GPS.csv") ODK_ONA_ID_GPS <- unique(ODK_ONA_ID_GPS[, c("plotID", "newPlotID","treatCode", "treatCode_label","trialID", "newTrialID", "trialCode", "fieldID", "newFieldID")]) QNPK[QNPK$Freq > 2, ] ## NPK treat 15 dsY[dsY$trialID == "ACTLNG000794" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG000794" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG001113" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG001113" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG001449" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG001449" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG001449",] dsY[dsY$trialID == "ACTLNG001471" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG001471" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG003021" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG003021" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG004514" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG004514" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG005415" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG005415" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLNG007526" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLNG007526" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLTZ000536" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLTZ000536" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLTZ000536",] dsY[dsY$trialID == "ACTLTZ000574" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLTZ000574" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLTZ001096" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLTZ001096" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLTZ001982" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLTZ001982" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLTZ002053" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLTZ002053" & treatCode_label %in% c("NPK rep1", "NPK rep2")) dsY[dsY$trialID == "ACTLTZ003783" & dsY$treat == "NPK", ] filter(PlotYield, trialID == "ACTLTZ003783" & treatCode_label %in% c("NPK rep1", "NPK rep2")) QNPK[QNPK$Freq > 2, ] FourTrials <- c("ACTLNG001449", "ACTLNG001471", "ACTLNG003021", "ACTLNG007526", "ACTLTZ003783") ## 4 plots per treatment assigned and all have yield: the four trials set at one location dropPlot <- c("ACPONG013121", "ACPONG013786", "ACPONG012388", "ACPONG034643", "ACPONG034667", "ACPOTZ017376","ACPOTZ017388", "ACPOTZ019233", "ACPOTZ019246", "ACPOTZ007557", "ACPOTZ007559", "ACPOTZ017375", "ACPOTZ017387", "ACPOTZ017400", "ACPOTZ017865", "ACPOTZ017873", "ACPOTZ017917", "ACPOTZ017841", "ACPOTZ017881", "ACPOTZ017912") ### the other treatments should have only one value QNPKnot <- QNPKnot[!QNPKnot$Var2 %in% FourTrials, ] QNPKnot[QNPKnot$Freq > 1, ] length(unique(QNPKnot[QNPKnot$Freq > 1, ]$Var2))##12 trials dsY[dsY$trialID == "ACTLNG000794" & dsY$treat %in% c("PK", "half_NPK"), ] filter(PlotYield, trialID == "ACTLNG000794" & treatCode_label %in% c("PK", "half_NPK")) dsY[dsY$trialID == "ACTLNG001093" & dsY$treat %in% c("NK", "CON"), ] filter(PlotYield, trialID == "ACTLNG001093" & treatCode_label %in% c("NK", "CON")) filter(PlotYield, trialID == "ACTLNG001093") filter(ODK_ONA_ID_GPS, plotID == "ACPONG014086") dsY[dsY$trialID == "ACTLNG005404" & dsY$treat %in% c("NP"), ] filter(PlotYield, trialID == "ACTLNG005404" & treatCode_label %in% c("NP")) filter(PlotYield, trialID == "ACTLNG005404") filter(ODK_ONA_ID_GPS, plotID == "ACTLNG005404") dsY[dsY$trialID == "ACTLNG005415" & dsY$treat %in% c("PK", "NK","NP", "half_NPK"), ] filter(PlotYield, trialID == "ACTLNG005415" & treatCode_label %in% c("PK", "NK", "NP", "half_NPK")) filter(PlotYield, trialID == "ACTLNG005415") filter(ODK_ONA_ID_GPS, plotID == "ACTLNG005415") dsY[dsY$trialID == "ACTLTZ000419" & dsY$treat %in% c("NP"), ] filter(PlotYield, trialID == "ACTLTZ000419" & treatCode_label %in% c("NP")) dsY[dsY$trialID == "ACTLTZ000574" & dsY$treat %in% c("PK", "NP", "CON"), ] filter(PlotYield, trialID == "ACTLTZ000574" & treatCode_label %in% c("PK", "NP", "CON")) filter(PlotYield, trialID == "ACTLTZ000574") filter(ODK_ONA_ID_GPS, plotID == "ACTLTZ000574") dsY[dsY$trialID == "ACTLTZ000580" & dsY$treat %in% c("PK"), ] filter(PlotYield, trialID == "ACTLTZ000580" & treatCode_label %in% c("PK")) dsY[dsY$trialID == "ACTLTZ001096" & dsY$treat %in% c("PK", "NK", "NP", "half_NPK", "CON"), ] filter(PlotYield, trialID == "ACTLTZ001096" & treatCode_label %in% c("PK", "NK", "NP", "half_NPK", "CON")) filter(PlotYield, trialID == "ACTLTZ001096") dsY[dsY$trialID == "ACTLTZ001982" & dsY$treat %in% c("PK", "NK", "NP", "half_NPK", "CON"), ] filter(PlotYield, trialID == "ACTLTZ001982" & treatCode_label %in% c("PK", "NK", "NP", "half_NPK", "CON")) dsY[dsY$trialID == "ACTLTZ002053" & dsY$treat %in% c("PK", "NK", "NP", "half_NPK", "CON"), ] filter(PlotYield, trialID == "ACTLTZ002053" & treatCode_label %in% c("PK", "NK", "NP", "half_NPK", "CON")) dsY[dsY$trialID == "ACTLTZ002527" & dsY$treat %in% c("CON"), ] filter(PlotYield, trialID == "ACTLTZ002527" & treatCode_label %in% c("CON")) dsY[dsY$trialID == "ACTLTZ002667" & dsY$treat %in% c("NP"), ] filter(PlotYield, trialID == "ACTLTZ002667" & treatCode_label %in% c("NP")) QNPKnot[QNPKnot$Freq > 1, ] FourTrials <- c("ACTLNG001449", "ACTLNG001471", "ACTLNG003021", "ACTLNG007526", "ACTLTZ003783") ## 4 plots per treatment assigned and all have yield: the four trials set at one location dropPlot <- c("ACPONG013121", "ACPONG013786", "ACPONG012388", "ACPONG034643", "ACPONG034667", "ACPOTZ017376","ACPOTZ017388", "ACPOTZ019233", "ACPOTZ019246", "ACPOTZ007557", "ACPOTZ007559", "ACPOTZ017375", "ACPOTZ017387", "ACPOTZ017400", "ACPOTZ017865", "ACPOTZ017873", "ACPOTZ017917", "ACPOTZ017841", "ACPOTZ017881", "ACPOTZ017912", "ACPONG014083", "ACPONG003591", "ACPONG034900", "ACPONG012340", "ACPONG034687", "ACPONG034682", "ACPONG034688","ACPOTZ019272", "ACPOTZ003961", "ACPOTZ003962", "ACPOTZ019245", "ACPOTZ019793", "ACPOTZ007501", "ACPOTZ007584", "ACPOTZ007589", "ACPOTZ007590", "ACPOTZ007586", "ACPOTZ017377", "ACPOTZ017389", "ACPOTZ017370", "ACPOTZ017382", "ACPOTZ017694", "ACPOTZ017697", "ACPOTZ017373", "ACPOTZ017385", "ACPOTZ017396", "ACPOTZ017384", "ACPOTZ017372", "ACPOTZ017384", "ACPOTZ017396", "ACPOTZ017374", "ACPOTZ017386", "ACPOTZ017398", "ACPOTZ017843", "ACPOTZ017879", "ACPOTZ017935", "ACPOTZ017888", "ACPOTZ017845", "ACPOTZ017903", "ACPOTZ017889", "ACPOTZ017840", "ACPOTZ017836", "ACPOTZ017887", "ACPOTZ017913", "ACPOTZ017844", "ACPOTZ017880", "ACPOTZ017905", "ACPOTZ015719", "ACPOTZ015968") dsY2 <- filter(dsY, (plotID %in% dropPlot) == FALSE) nrow(dsY) nrow(dsY2) length(unique(dsY$trialID)) length(unique(dsY2$trialID)) ### in order to see what is causing the higher yeilds we need to see how many plants are harvested wihtin a plot ## if a trial has plant level yield, then it should be done across all treatments and teoretically there will be only 25 plants/plot ## if a trial does not have plant level yield then teoretically there are 30 plants/plot setwd("/home/akilimo/projects/SoilNPK/Rwanda") PA_ODK_ONA <- read.csv("Plant_ODK_ONA.csv") ONA_ODK_rootYield <- read.csv("Root_Yield_Cassava_ONA_ODK.csv") ReplacePlot_ODK_ONA <- read.csv("ReplacePlot_ODK_ONA.csv") plantYield <- droplevels(ONA_ODK_rootYield[!ONA_ODK_rootYield$plantID == "", ]) head(plantYield) plantYield$plantYeild <- rowSums(subset(plantYield,select=c(tuberizedRootsFW, tuberizedDiseasedRootsFW, tuberizedSmallRootsFW, tuberizedMarketableRootsFW)),na.rm=TRUE) plantYield <- unique(plantYield[, c("today","plantID", "plantYeild")]) plantYieldLined1 <- unique(merge(PA_ODK_ONA, plantYield, by="plantID")) ## not replaced plantids head(plantYieldLined1) plantYield[!plantYield$plantID %in% plantYieldLined1$plantID, ] ## replaced plantids repPlant <- ReplacePlot_ODK_ONA[ReplacePlot_ODK_ONA$plantID %in% plantYield[!plantYield$plantID %in% plantYieldLined1$plantID, ]$plantID, ] repPlant <- repPlant[, c("plantID", "newPlotID")] colnames(repPlant) <- c("plantID", "plotID") plantYieldLined2 <- unique(merge(repPlant, plantYield, by="plantID")) head(plantYieldLined2) plantYieldLinked <- rbind(plantYieldLined1, plantYieldLined2) plantYield[!plantYield$plantID %in% plantYieldLinked$plantID, ] ## no plot link head(dsY) head(plantYieldLinked) dsY$plantYield <- ifelse(dsY$plotID %in% unique(plantYieldLinked$plotID), "Yes", "No") ### checking if plant yield was taken in a trial, is it then done across all treats? trialplantyield <- unique(dsY[, c("trialID", "plantYield")]) Q <- as.data.frame(table(trialplantyield$trialID)) ppyield <- dsY[dsY$trialID %in% Q[Q$Freq > 1, ]$Var1, c("trialID", "plotID","treat", "plantYield", "nrPlantsNP","rootY")] ppyield <- ppyield[order(ppyield$trialID, ppyield$treat), ] length(unique(ppyield$trialID))##74 length(unique(dsY$trialID))##649 head(ppyield, 10) plantYieldLinked[plantYieldLinked$plotID == "ACPONG014251", ] ppyield[!is.na(ppyield$nrPlantsNP), ] ReplacePlot_ODK_ONA[ReplacePlot_ODK_ONA$trialID == "ACTLTZ001969" & ReplacePlot_ODK_ONA$treats == ""] ggplot(ppyield, aes(treat, rootY, col=plantYield))+ geom_point()+ facet_wrap(~treat) ####################################################################### ###################################################################################################### ### data get corrected for nr of plants per plot NOT_Oct2020 and tehn teh pldates were searched manually NOT_Nov2020 ###################################################################################################### require(lsmeans) require(emmeans) require(lme4) require(lmerTest) require(tidyr) require(ggplot2) require(party) require(plyr) require(gtools) library(raster) library(sf) require(MASS) #require(glmnet) require(caret) require(mlbench) require(psych) library(rgdal) require(randomForest) require(rgeos) require(rpart) setwd("/home/akilimo/projects/SoilNPK") # dsY <- read.csv("NOT_Oct2020.csv") dsY <- read.csv("NOT_Nov2020.csv") dsY <- dsY[rowSums(is.na(dsY)) <= ncol(dsY),] names(dsY)[names(dsY)=="rootY_tha"] <- "rootY" dsY$treat <- factor(dsY$treat, levels = c("NPK_micro", "NPK", "PK", "NK", "NP", "half_NPK", "CON")) dropPlot <- c("ACPONG013121", "ACPONG013786", "ACPONG012388", "ACPONG034643", "ACPONG034667", "ACPOTZ017376","ACPOTZ017388", "ACPOTZ019233", "ACPOTZ019246", "ACPOTZ007557", "ACPOTZ007559", "ACPOTZ017375", "ACPOTZ017387", "ACPOTZ017400", "ACPOTZ017865", "ACPOTZ017873", "ACPOTZ017917", "ACPOTZ017841", "ACPOTZ017881", "ACPOTZ017912", "ACPONG014083", "ACPONG003591", "ACPONG034900", "ACPONG012340", "ACPONG034687", "ACPONG034682", "ACPONG034688","ACPOTZ019272", "ACPOTZ003961", "ACPOTZ003962", "ACPOTZ019245", "ACPOTZ019793", "ACPOTZ007501", "ACPOTZ007584", "ACPOTZ007589", "ACPOTZ007590", "ACPOTZ007586", "ACPOTZ017377", "ACPOTZ017389", "ACPOTZ017370", "ACPOTZ017382", "ACPOTZ017694", "ACPOTZ017697", "ACPOTZ017373", "ACPOTZ017385", "ACPOTZ017396", "ACPOTZ017384", "ACPOTZ017372", "ACPOTZ017384", "ACPOTZ017396", "ACPOTZ017374", "ACPOTZ017386", "ACPOTZ017398", "ACPOTZ017843", "ACPOTZ017879", "ACPOTZ017935", "ACPOTZ017888", "ACPOTZ017845", "ACPOTZ017903", "ACPOTZ017889", "ACPOTZ017840", "ACPOTZ017836", "ACPOTZ017887", "ACPOTZ017913", "ACPOTZ017844", "ACPOTZ017880", "ACPOTZ017905", "ACPOTZ015719", "ACPOTZ015968") dsY <- filter(dsY, (plotID %in% dropPlot) == FALSE) head(dsY) boxplot(dsY$rootY ~ dsY$country) #need to check for outliers - especially values above 150 t/ha dsY <- dsY[dsY$rootY < 150,] #removing impossible values ggplot(dsY, aes(y=rootY, x=trialID)) + geom_boxplot() + facet_wrap(~country, scales="free_x") #likely still some individual outliers, and some suspect trials with overall high yield values ggplot(dsY, aes(y=rootY, x=treat)) + geom_boxplot() + facet_wrap(~country) #simple model to detect and remove outliers mod1 <- lmer(data=dsY, rootY ~ treat + country + (1|trialID)) plot(mod1) #some issues with outliers and heteroscedasticity - transformation needed mod2 <- lmer(data=dsY, log10(rootY+1) ~ treat + country + (1|trialID)) plot(mod2) #looks much better dsY <- dsY[which(resid(mod2)<0.5),] #outliers removed dsY <- dsY[which(resid(mod2)>-0.6),] #outliers removed dsY$Region <- as.character(dsY$Region) dsY[is.na(dsY$lat),]$Region <- as.character(dsY[is.na(dsY$lat),]$country) #trials without GPS coordinates lack Region - identify the region based on neighbouring sequence numbers for trialID for (i in na.omit(unique(dsY[is.na(dsY$lat),]$trialID))){ nr <- as.numeric(substr(i, 7, 12)) rg <- unique(dsY[dsY$trialID %in% paste0(substr(i, 1, 6), ifelse(nr<1000, "000", "00"), (nr-100):(nr+100)),]$Region) rg <- rg[!rg %in% c("Nigeria", "Tanzania")] if(length(rg)==1){ dsY[dsY$trialID==i & !is.na(dsY$trialID),]$Region <- rg print(paste0("Region for trialID ", i, " corrected to ", rg)) }else{ if(length(rg)>1){ rg <- unique(dsY[dsY$trialID %in% paste0(substr(i, 1, 6), ifelse(nr<1000, "000", "00"), (nr-50):(nr+50)),]$Region) rg <- rg[!rg %in% c("Nigeria", "Tanzania")] dsY[dsY$trialID==i & !is.na(dsY$trialID),]$Region <- rg print(paste0("Region for trialID ", i, " corrected to ", rg)) }else{ print(paste0("No region found for trialID", i)) } } } dsY$Region <- droplevels(as.factor(revalue(dsY$Region, c("NG_east" = "SE Nigeria", "NG_west" = "SW Nigeria", "TZ_lake" = "Tanzania Lake Region", "TZ_south" = "Tanzania Coastal Region")))) dsY <- dsY[!is.na(dsY$Region),] acaibox <- ggplot(dsY, aes(treat, rootY, fill=treat)) + geom_boxplot() + facet_grid(harvestYear~Region) + ylab("Root yield (t/ha)") + theme_bw() + theme(axis.title.x = element_blank(), axis.title.y = element_text(size=12, face="bold"), axis.text.y = element_text(size=12), axis.text.x = element_text(size=11, angle=35, hjust=0.9, vjust=0.9), strip.text = element_text(size = 13, face="bold"), strip.background = element_blank(), legend.title = element_blank(), legend.text = element_text(size=14), legend.position = "none") ggsave("acai_box.png", acaibox, width = 10, height = 6) ggvioline <- ggplot(dsY, aes(y=rootY, x=treat)) + geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), fill = "grey85", ) + geom_jitter(aes(shape = as.factor(harvestYear)), height=0, width=.3, size=2, alpha=.2) + scale_shape_manual(values = 15:18) + facet_wrap(~Region) + ylab("Cassava fresh root yield (t/ha)") + theme_bw()+ theme(axis.title.x = element_blank(), axis.title.y = element_text(size=12, face="bold"), axis.text.y = element_text(size=12), axis.text.x = element_text(size=12, angle=20, hjust=0.8), strip.text = element_text(size = 13, face="bold"), strip.background = element_blank(), legend.title = element_blank(), legend.text = element_text(size=14), legend.position = "none" ) ggsave("acai_ggvioline.png", ggvioline, width = 10, height = 6) head(dsY) summary(dsY$rootY ) trialmeans <- ddply(dsY, .(trialID, trialCode, harvestYear, Region), summarize, meanT = mean(rootY)) head(trialmeans) table(trialmeans$Region, trialmeans$harvestYear) trialmean <- ggplot(trialmeans, aes(y=meanT, x=trialID)) + geom_boxplot() + facet_wrap(~Region, scales="free_x")+ xlab("Trials") + ylab("Average yield (t/ha) at trial level aross all treatments") + theme_bw()+ theme(axis.text.x = element_blank(), axis.title = element_text(size=14)) ggsave("trialmean.png", trialmean, width=10, height = 7) highYieldTrials <- droplevels(filter(trialmeans, meanT > 60)) table(highYieldTrials$Region, highYieldTrials$harvestYear) FR_forCKAN <- droplevels(dsY[!dsY$trialID %in% trialmeans[trialmeans$meanT > 60, ]$trialID, ]) FR_forCKAN <- FR_forCKAN[FR_forCKAN$rootY <= 75, ] ggplot(FR_forCKAN, aes(y=rootY, x=trialID)) + geom_boxplot() + facet_wrap(~country, scales="free_x") FR_forCKAN$treat <- factor(FR_forCKAN$treat, levels = c("NPK", "half_NPK", "NPK_micro", "NP", "NK", "PK", "CON")) ggplot(FR_forCKAN, aes(treat, rootY, fill=treat)) + geom_boxplot()+ facet_grid(harvestYear~Region)+ ylab("Root yield (t/ha)")+ theme_bw()+ theme(axis.title.y = element_text(size=16, face="bold"), axis.title.x = element_blank(), axis.text = element_text(size=14), axis.text.x = element_text(angle=-90, hjust=0, vjust=.5), strip.text = element_text(size=16, face="bold")) FR_forCKAN <- FR_forCKAN[,c("country", "Region","Admin1","Admin2","lon", "lat","trialID", "trialCode", "plotID", "treat", "harvestDate", "harvestMonth", "harvestYear", "plantSeason", "rootY")] FR_forCKAN <- unique(FR_forCKAN) names(FR_forCKAN) <- c("country", "Region","Admin1","Admin2","lon", "lat","trialID", "trialCode", "plotID", "treat", "harvestDate", "harvestMonth", "harvestYear", "plantSeason", "rootYield (t/ha)") head(FR_forCKAN) setwd("D:/IITADocs/CKAN_data") write.csv(FR_forCKAN, "ACAI_FR_forCKAN.csv", row.names = FALSE) ds_formodel <- droplevels(dsY[!dsY$trialID %in% trialmeans[trialmeans$meanT > 60, ]$trialID, ]) ds_formodel <- unique(subset(ds_formodel, select=-c(X, index))) ggplot(ds_formodel, aes(treat, rootY, fill=treat)) + geom_boxplot()+ facet_grid(harvestYear~Region)+ ylab("Root yield (t/ha)")+ theme_bw()+ theme(axis.title.y = element_text(size=16, face="bold"), axis.title.x = element_blank(), axis.text = element_text(size=14), axis.text.x = element_text(angle=-90, hjust=0, vjust=.5), strip.text = element_text(size=16, face="bold")) str(ds_formodel) ggfin <- ggplot(ds_formodel, aes(treat, rootY, fill=treat)) + geom_boxplot()+ facet_grid(harvestYear~Region)+ ylab("Root yield (t/ha)")+ theme_bw()+ theme(axis.title.x = element_blank(), axis.title.y = element_text(size=12, face="bold"), axis.text.y = element_text(size=12), axis.text.x = element_text(size=11, angle=35, hjust=0.9, vjust=0.9), strip.text = element_text(size = 13, face="bold"), strip.background = element_blank(), legend.title = element_blank(), legend.text = element_text(size=14), legend.position = "none") ggsave("ACAI_box.png", ggfin, width=10, height = 6) ds_formodel <- droplevels(dsY[!dsY$trialID %in% trialmeans[trialmeans$meanT > 60, ]$trialID, ]) ds_formodel <- unique(subset(ds_formodel, select=-c(X, index))) ds_formodel[ds_formodel$harvestYear == 2018 & ds_formodel$Region == "Tanzania Coastal Region" & ds_formodel$rootY>65, ] ds_formodel <- ds_formodel[!ds_formodel$plotID == "ACPOTZ002908", ] unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/22-Aug-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/22-Aug-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "08/22/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 08 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("22/Aug/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-234)+278 filter(ds_formodel, plantingDate == "NA/22-Jul-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/22-Jul-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "07/22/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 07 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("22/Jul/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-203)+264 filter(ds_formodel, plantingDate == "NA/12-Aug-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/12-Aug-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "08/12/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 08 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("12/Aug/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-224)+268 filter(ds_formodel, plantingDate == "NA/23-Jun-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/23-Jun-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "06/23/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 06 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("23/Jun/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-174)+262 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/15-Jun-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/15-Jun-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "06/15/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 06 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("15/Jun/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-166)+241 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/21-Jun-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/21-Jun-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "06/21/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 06 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("21/Jun/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-172)+242 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/27-Jun-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/27-Jun-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "06/21/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 06 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("21/Jun/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-172)+255 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/10-Jul-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/10-Jul-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "07/10/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 07 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("10/Jul/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-191)+206 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/11-Jul-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/11-Jul-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "07/11/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 07 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("11/Jul/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-191)+206 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/26-Jul-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/26-Jul-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "07/26/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 07 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("26/Jul/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-207)+219 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/31-Jul-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/31-Jul-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "07/31/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 07 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("31/Jul/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-212)+230 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/28-Dec-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/28-Dec-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "12/28/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 07 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("28/Dec/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-362)+276 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/02-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/02-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/02/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("02/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-2)+281 ds_formodel$plantingDate <- as.character(ds_formodel$plantingDate) unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/27-Dec-17/NA") plpld <- c("ACPOTZ002303", "ACPOTZ002310", "ACPOTZ002307", "ACPOTZ002308", "ACPOTZ002304", "ACPOTZ002305", "ACPOTZ002309", "ACPOTZ002306") ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "11/27/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 11 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("27/Dec/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-361)+277 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/19-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/19-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/19/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("19/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-19)+307 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/02-Dec-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/02-Dec-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "11/02/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 11 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("02/Dec/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-336)+275 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/04-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/04-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/04/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("01/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-1)+287 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/03-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/03-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/03/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("03/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-3)+285 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/20-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/20-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/20/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("20/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-20)+312 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/08-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/08-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/08/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("08/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-8)+313 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/10-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/10-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/10/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("10/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-10)+309 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/15-Jan-17/NA") plpld <- filter(ds_formodel, plantingDate == "NA/15-Jan-17/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "01/15/2017" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2017 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 01 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("15/Jan/2017", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-15)+308 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/10-Feb-18/NA") plpld <- c("ACPOTZ015865","ACPOTZ015862","ACPOTZ015863","ACPOTZ015860","ACPOTZ015864","ACPOTZ015861") ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "02/10/2018" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2018 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 02 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("10/Feb/2018", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-41)+35 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/20-Feb-18/NA") plpld <- filter(ds_formodel, plantingDate == "NA/20-Feb-18/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "02/20/2018" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2018 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 02 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("20/Feb/2018", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-51)+31 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/19-Feb-18/NA") plpld <- filter(ds_formodel, plantingDate == "NA/19-Feb-18/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "02/19/2018" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2018 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 02 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("19/Feb/2018", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-50)+78 unique(ds_formodel$plantingDate) filter(ds_formodel, plantingDate == "NA/13-Feb-18/NA") plpld <- filter(ds_formodel, plantingDate == "NA/13-Feb-18/NA")$plotID ds_formodel[ds_formodel$plotID %in% plpld, ]$plantingDate <- "02/13/2018" ds_formodel[ds_formodel$plotID %in% plpld, ]$plantYear <- 2018 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantMonth <- 02 ds_formodel[ds_formodel$plotID %in% plpld, ]$plantDOY <- lubridate::yday(as.Date("13/Feb/2018", format = '%d/%b/%Y')) ds_formodel[ds_formodel$plotID %in% plpld, ] ds_formodel[ds_formodel$plotID %in% plpld, ]$DAP <- (365-44)+73 saveRDS(ds_formodel, "NOT_formodel.RDS") head(ds_formodel) #################################################################################### ## linear mixed model to get BLUP :: this is made on my POC because the server does not have space #################################################################################### setwd("/home/akilimo/projects/SoilNPK") dsY <- readRDS("NOT_formodel.RDS") ## made in teh server 64 dsY$harvestYear <- as.factor(dsY$harvestYear) str(dsY) hist(dsY$DAP) summary(dsY$DAP) head(dsY) ddply(dsY, .(treat), summarise, treatmean = mean(rootY)) dsY$DAPclass <- ifelse(dsY$DAP < 300, "short", ifelse(dsY$DAP > 500, "long", "10to15")) dap_trial <- unique(dsY[, c("DAPclass", "trialID", "country")]) table(dap_trial$DAPclass) table(dap_trial$DAPclass, dap_trial$country) ####################################################################################### ### data filterd for planting dates ####################################################################################### dsY_DAP <- filter(dsY, DAP <= 500) ## trials harvested below 15 MAP. Trials with out planting dates and >500 DAP are dropped summary(dsY_DAP$DAP) dsY_DAP$DAPclass <- as.factor(ifelse(dsY_DAP$DAP <= 300, "10MAP", ifelse(dsY_DAP$DAP > 300 & dsY_DAP$DAP <=335, "11MAP", ifelse(dsY_DAP$DAP > 335 & dsY_DAP$DAP <= 365, "12MAP", ifelse(dsY_DAP$DAP > 365 & dsY_DAP$DAP < 395, "13MAP", ifelse(dsY_DAP$DAP >=395 & dsY_DAP$DAP < 426, "14MAP", ifelse(dsY_DAP$DAP >=426 & dsY_DAP$DAP < 455, "15MAP", ifelse(dsY_DAP$DAP >=455, "16MAP", "check")))))))) summary(dsY_DAP[dsY_DAP$DAPclass == "check", ]$DAP) dsY_DAP$DAPclass <- as.factor(dsY_DAP$DAPclass) ggplot(dsY_DAP, aes(DAPclass, rootY))+ geom_boxplot() modAC1_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat*Region*harvestYear + DAPclass+ (0+treat|trialID), REML = FALSE) modAC2_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat + treat:Region + DAPclass+(0+treat|trialID), REML = FALSE) modAC3_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat + treat:harvestYear + DAPclass+(0+treat|trialID), REML = FALSE) modAC4_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat + harvestYear:Region + DAPclass+(0+treat|trialID), REML = FALSE) modAC5_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat + treat:Region:harvestYear + DAPclass+(0+treat|trialID), REML = FALSE) modAC6_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat:Region:harvestYear + DAPclass+(0+treat|trialID), REML = FALSE) modAC7_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat:Region:harvestYear + (0+treat|trialID), REML = FALSE) anova(modAC1_DAP, modAC2_DAP) anova(modAC1_DAP, modAC3_DAP) anova(modAC1_DAP, modAC4_DAP) anova(modAC1_DAP, modAC5_DAP) anova(modAC1_DAP, modAC6_DAP) anova(modAC5_DAP, modAC6_DAP) anova(modAC6_DAP, modAC7_DAP) modAC_DAP <- lmer(data=dsY_DAP, sqrt(rootY) ~ treat:Region:harvestYear + DAPclass + (0+treat|trialID)) plot(modAC_DAP) summary(modAC_DAP) lsm3_DAP <- lsmeans(modAC_DAP, c("treat", "Region", "harvestYear")) lsmDF3_DAP <- as.data.frame(lsm3_DAP) lsmDF3_DAP <- lsmDF3_DAP[,1:4] head(lsmDF3_DAP) saveRDS(lsmDF3_DAP, "lsmDF3_ACAI_NOT_DAP15.RDS") names(lsmDF3_DAP)[1] <- "treatment" row.names(lsmDF3_DAP) <- NULL ran3_DAP <- ranef(modAC_DAP)$trialID ran3_DAP$trialID <- row.names(ran3_DAP) names(ran3_DAP) <- gsub("treat", "", names(ran3_DAP)) ran3_DAP <- melt(ran3_DAP, measure.vars=c("NPK", "PK", "NK", "NP", "half_NPK", "CON"), variable.name="treatment", id.vars="trialID", value.name="blup") ran3_DAP <- merge(ran3_DAP, unique(dsY_DAP[, c("trialID", "harvestYear", "Region")]), by="trialID") ran3_DAP <- merge(ran3_DAP, lsmDF3_DAP, by=c("treatment", "harvestYear", "Region")) ran3_DAP$blup <- ran3_DAP$blup + ran3_DAP$lsmean ran3_DAP <- subset(ran3_DAP, select=-c(lsmean)) ran3_DAP$blup <- (ran3_DAP$blup)**2 head(ran3_DAP) saveRDS(ran3_DAP, "ran3AC_DAP.RDS") ####################################################################################### ### data without filtering for planting dates info ####################################################################################### modAC1 <- lmer(data=dsY, sqrt(rootY) ~ treat*Region*harvestYear + (0+treat|trialID), REML = FALSE) modAC2 <- lmer(data=dsY, sqrt(rootY) ~ treat + treat:Region + (0+treat|trialID), REML = FALSE) modAC3 <- lmer(data=dsY, sqrt(rootY) ~ treat + treat:harvestYear + (0+treat|trialID), REML = FALSE) modAC4 <- lmer(data=dsY, sqrt(rootY) ~ treat + harvestYear:Region + (0+treat|trialID), REML = FALSE) modAC5 <- lmer(data=dsY, sqrt(rootY) ~ treat + treat:Region:harvestYear + (0+treat|trialID), REML = FALSE) modAC6 <- lmer(data=dsY, sqrt(rootY) ~ treat:Region:harvestYear + (0+treat|trialID), REML = FALSE) anova(modAC1, modAC2) anova(modAC1, modAC3) anova(modAC1, modAC4) anova(modAC1, modAC5) anova(modAC1, modAC6) anova(modAC5, modAC6) ## mod 5 is the best model modAC <- lmer(data=dsY, sqrt(rootY) ~ treat:Region:harvestYear + (0+treat|trialID)) plot(modAC) summary(modAC) # dsY <- droplevels(dsY[which(resid(modAC) < 2),] ) # modAC <- lmer(data=dsY, sqrt(rootY_tha) ~ treat + treat:Region:harvestYear + (0+treat|trialID)) # plot(modAC) dsY$predicted <- (predict(modAC)**2) # Save the predicted values dsY$residuals <- (residuals(modAC)**2) # Save the residual values ggplot(dsY, aes(rootY, predicted, col=Region)) + geom_point() + geom_abline(slope=1, intercept = 0)+ facet_grid(harvestYear~Region, scale="free")+ ylab("Predicted root yield")+ xlab("Observed root yield") + theme_bw()+theme(legend.position = "none", strip.text = element_text(size=12), axis.text = element_text(size=10), axis.title = element_text(size=12)) ref.grid(modAC, emm_options(pbkrtest.limit = 4866)) lsm3 <- lsmeans(modAC, c("treat", "Region", "harvestYear")) lsmDF3 <- as.data.frame(lsm3) lsmDF3 <- lsmDF3[,1:4] head(lsmDF3) saveRDS(lsmDF3, "lsmDF3_ACAI_NOT_allDAPs.RDS") # lsmDF3 <- readRDS("lsmDF3_ACAI_NOT_cleaned.RDS") names(lsmDF3)[1] <- "treatment" row.names(lsmDF3) <- NULL ran3 <- ranef(modAC)$trialID ran3$trialID <- row.names(ran3) names(ran3) <- gsub("treat", "", names(ran3)) ran3 <- melt(ran3, measure.vars=c("NPK", "PK", "NK", "NP", "half_NPK", "CON"), variable.name="treatment", id.vars="trialID", value.name="blup") ran3 <- merge(ran3, unique(dsY[, c("trialID", "harvestYear", "Region")]), by="trialID") ran3 <- merge(ran3, lsmDF3, by=c("treatment", "harvestYear", "Region")) ran3$blup <- ran3$blup + ran3$lsmean ran3 <- subset(ran3, select=-c(lsmean)) ran3$blup <- (ran3$blup)**2 head(ran3) saveRDS(ran3, "ran3AC_allDAPs.RDS") ################################################################################### ## plot model output done on PC ################################################################################### setwd("D:/ACAI_Wrapper/Rwanda") setwd("/home/akilimo/projects/SoilNPK/Rwanda") lsmDF3_DAP_15 <- readRDS("lsmDF3_ACAI_NOT_DAP15.RDS") lsmDF3_DAP <- readRDS("lsmDF3_ACAI_NOT_allDAPs.RDS") ran3_DAP_15 <- readRDS("ran3AC_DAP.RDS") ran3_DAP <- readRDS("ran3AC_allDAPs.RDS") nrow(lsmDF3_DAP_15) nrow(lsmDF3_DAP) nrow(ran3_DAP_15) nrow(ran3_DAP) ### plotinthte BLUPs ran_AC <- ran3_DAP ## all without filtering for DAP rm(ran_AC) ran_AC <- ran3_DAP_15 ## DAP upto 15 MAP ran_NG_2017E <- ran_AC[ran_AC$harvestYear == 2017 & ran_AC$Region == "SE Nigeria", ] ran_NG_2018E <- ran_AC[ran_AC$harvestYear == 2018 & ran_AC$Region == "SE Nigeria", ] ran_NG_2019E <- ran_AC[ran_AC$harvestYear == 2019 & ran_AC$Region == "SE Nigeria", ] ran_NG_2017W <- ran_AC[ran_AC$harvestYear == 2017 & ran_AC$Region == "SW Nigeria", ] ran_NG_2018W <- ran_AC[ran_AC$harvestYear == 2018 & ran_AC$Region == "SW Nigeria", ] ran_NG_2019W <- ran_AC[ran_AC$harvestYear == 2019 & ran_AC$Region == "SW Nigeria", ] ran_TZ_2017L <- ran_AC[ran_AC$harvestYear == 2017 & ran_AC$Region == "Tanzania Lake Region", ] ran_TZ_2018L <- ran_AC[ran_AC$harvestYear == 2018 & ran_AC$Region == "Tanzania Lake Region", ] ran_TZ_2017S <- ran_AC[ran_AC$harvestYear == 2017 & ran_AC$Region == "Tanzania Coastal Region", ] ran_TZ_2018S <- ran_AC[ran_AC$harvestYear == 2018 & ran_AC$Region == "Tanzania Coastal Region", ] ran_TZ_2019S <- ran_AC[ran_AC$harvestYear == 2019 & ran_AC$Region == "Tanzania Coastal Region", ] ranCY <- rbind(ran_NG_2017E, ran_NG_2018E, ran_NG_2019E,ran_NG_2017W, ran_NG_2018W, ran_NG_2019W, ran_TZ_2017L, ran_TZ_2018L, ran_TZ_2017S, ran_TZ_2018S, ran_TZ_2019S) head(ranCY) ddply(ranCY, .(Region, treatment), summarize, avyeartreat = mean(blup)) ran.w <- spread(ran_AC, treatment, blup) head(ran.w) hist(ran.w$NPK - ran.w$CON) ran.w$harvestYear <- as.numeric(as.character(ran.w$harvestYear)) saveRDS(ran.w, "ranw_allDAP.RDS") # saveRDS(ran.w, "ranw_DAP15.RDS") ran.w <- readRDS("ranw_allDAP.RDS") head(ran.w) ggplot(ran.w, aes(PK, NPK, col=harvestYear)) + geom_point() + facet_grid(harvestYear~Region, scales="free") + geom_abline(slope=1, intercept = 0) + theme_bw() longData <- gather(ran.w, treat, value, PK:CON) head(longData) longData$YieldDiff <- longData$NPK - longData$value longData$NutrientMissing <- ifelse(longData$treat == "NP", "K", ifelse(longData$treat == "NK", "P", ifelse(longData$treat == "PK", "N", "NPK"))) longData$TreatDiff <- ifelse(longData$treat == "NP", "NPK - NP (K Eff.)", ifelse(longData$treat == "NK", "NPK - NK (P Eff.)", ifelse(longData$treat == "PK", "NPK - PK (N Eff.)", "NPK - CON (NPK Eff.)"))) longData$harvestYear <- as.factor(longData$harvestYear) head(longData) gg17 <- ggplot(data=longData[longData$harvestYear == 2017, ], aes(value, y=YieldDiff)) + facet_grid(harvestYear+Region ~ TreatDiff) + geom_point(size=1.2, col="grey4") + stat_density2d(col = "orangered") + geom_abline(intercept=0,slope=0, col="darkgrey") + ##, chartreuse3, col = "deepskyblue1", orangered xlab("") + ylab("Yield difference [t/ha]") + theme_bw() + theme(axis.text = element_text(size=13), axis.title = element_text(size=15), strip.text = element_text(size=14), legend.text =element_text(size=15), legend.title = element_blank()) # ggsave("NOT_BLUPS_AC_2017_DAPall.png", gg17, width=12, height = 6) ggsave("NOT_BLUPS_AC_2017_DAP15.png", gg17, width=12, height = 6) gg18 <- ggplot(data=longData[longData$harvestYear == 2018, ], aes(value, y=YieldDiff)) + facet_grid(harvestYear+Region ~ TreatDiff) + geom_point(size=1.2, col="grey4") + stat_density2d(col = "chartreuse3") + geom_abline(intercept=0,slope=0, col="darkgrey") + ##, chartreuse3, col = "deepskyblue1", orangered xlab("") + ylab("Yield difference [t/ha]") + theme_bw() + theme(axis.text = element_text(size=13), axis.title = element_text(size=15), strip.text = element_text(size=14), legend.text =element_text(size=15), legend.title = element_blank()) # ggsave("NOT_BLUPS_AC_2018_DAPall.png", gg18, width=12, height = 6) ggsave("NOT_BLUPS_AC_2018_DAP15.png", gg18, width=12, height = 6) gg19 <- ggplot(data=longData[longData$harvestYear == 2019, ], aes(value, y=YieldDiff)) + facet_grid(harvestYear+Region ~ TreatDiff) + geom_point(size=1.2, col="grey4") + stat_density2d(col = "deepskyblue1") + geom_abline(intercept=0,slope=0, col="darkgrey") + ##, chartreuse3, col = "deepskyblue1", orangered xlab("") + ylab("Yield difference [t/ha]") + theme_bw() + theme(axis.text = element_text(size=13), axis.title = element_text(size=15), strip.text = element_text(size=14), legend.text =element_text(size=15), legend.title = element_blank()) # ggsave("NOT_BLUPS_AC_2019_DAPall.png", gg19, width=12, height = 6) ggsave("NOT_BLUPS_AC_2019_DAP15.png", gg19, width=12, height = 6) ##############################################################################3 ## get ISIRIC data ############################################################################# setwd("/home/akilimo/lintul") source("sourceISRIC.R") head(dsY) ACAI_loc <- unique(dsY[, c("lon","lat", "trialID")]) ACAI_loc <- ACAI_loc[!is.na(ACAI_loc$lon), ] str(ACAI_loc) ##522 data points ISRIC_AC_NOT <- NULL for(i in 1:nrow(ACAI_loc)){ Sdata <- getISRICData(lat = ACAI_loc$lat[i], lon = ACAI_loc$lon[i]) Sdata$trialID <- ACAI_loc$trialID[i] ISRIC_AC_NOT <- rbind(ISRIC_AC_NOT, Sdata) } head(ISRIC_AC_NOT) setwd("/home/akilimo/projects/SoilNPK/Rwanda") saveRDS(ISRIC_AC_NOT, "ISRIC_ACAI_Cassava.RDS") ###########################################################################################333 ## prepare data for reverse QUEFTS usinf isda soils ############################################################################################### setwd("/home/akilimo/projects/SoilNPK/Coordinates") ACAI_landCropCover <- readRDS("ACAI_isda_cover.RDS" ) ACAI_isda <- readRDS("ACAI_isda.RDS") nrow(ACAI_isda) nrow(ACAI_landCropCover) ACAI_landCropCover$land_cover_2015 <- revalue(ACAI_landCropCover$land_cover_2015, c("Closed forest, unknown" = "losedForest", "Cultivated and managed vegetation/agriculture (cropland)" = "Cropland", "Open forest, unknown"= " OpenForest", "Urban / built up" = "Urban", "Shrubs" = "Shrubs", "Herbaceous vegetation" = "herbaceousVegetation")) ACAI_landCropCover$land_cover_2016 <- revalue(ACAI_landCropCover$land_cover_2016, c("Closed forest, unknown" = "losedForest", "Cultivated and managed vegetation/agriculture (cropland)" = "Cropland", "Open forest, unknown"= " OpenForest", "Urban / built up" = "Urban", "Shrubs" = "Shrubs", "Herbaceous vegetation" = "herbaceousVegetation")) ACAI_landCropCover$land_cover_2017 <- revalue(ACAI_landCropCover$land_cover_2017, c("Closed forest, unknown" = "losedForest", "Cultivated and managed vegetation/agriculture (cropland)" = "Cropland", "Open forest, unknown"= " OpenForest", "Urban / built up" = "Urban", "Shrubs" = "Shrubs", "Herbaceous vegetation" = "herbaceousVegetation")) ACAI_landCropCover$land_cover_2018 <- revalue(ACAI_landCropCover$land_cover_2018, c("Closed forest, unknown" = "losedForest", "Cultivated and managed vegetation/agriculture (cropland)" = "Cropland", "Open forest, unknown"= " OpenForest", "Urban / built up" = "Urban", "Shrubs" = "Shrubs", "Herbaceous vegetation" = "herbaceousVegetation")) ACAI_landCropCover$land_cover_2019 <- revalue(ACAI_landCropCover$land_cover_2019, c("Closed forest, unknown" = "losedForest", "Cultivated and managed vegetation/agriculture (cropland)" = "Cropland", "Open forest, unknown"= " OpenForest", "Urban / built up" = "Urban", "Shrubs" = "Shrubs", "Herbaceous vegetation" = "herbaceousVegetation")) ACAI_landCropCover <- ACAI_landCropCover[, c("lat","lon", "slope_angle", "fcc", "land_cover_2015", "land_cover_2016", "land_cover_2017", "land_cover_2018", "land_cover_2019")] ACAI_isda_main <- ACAI_isda[!ACAI_isda$Variable == "bedrock_depth", ] ACAI_bedRock_200cm <- ACAI_isda[ACAI_isda$Variable == "bedrock_depth", ] ACAI_bedRock_200cm$value <- as.numeric(as.character(ACAI_bedRock_200cm$value)) bedrockdepth1 <- ACAI_bedRock_200cm[seq(1,180, by=2), ] bedrockdepth2 <- ACAI_bedRock_200cm[seq(2,180, by=2), ] summary(bedrockdepth1$value) summary(bedrockdepth2$value) ACAI_isda <- rbind(ACAI_isda_main, bedrockdepth1) head(ACAI_isda) ACAI_isda$Variable <- paste(ACAI_isda$Variable, ACAI_isda$depth, sep="_") ACAI_isda_wide <- spread(ACAI_isda[, c("lat","lon","Variable", "value")], Variable, value, ) colnames(ACAI_isda_wide) <- gsub("-", "_", colnames(ACAI_isda_wide)) head(ACAI_isda_wide) ACAI_isda_wide <- unique(merge(ACAI_isda_wide, ACAI_landCropCover, by=c("lat", "lon"))) head(ACAI_isda_wide) str(ACAI_isda_wide) ## covert to numeric require(tidyverse) ACAI_isda_wide %>% str() ACAI_isda_wide$texture_class_0_20 <- as.factor(ACAI_isda_wide$texture_class_0_20) ACAI_isda_wide$texture_class_20_50 <- as.factor(ACAI_isda_wide$texture_class_20_50) ACAI_isda_wide <- ACAI_isda_wide %>% mutate_if(is.character, as.numeric) ACAI_isda_wide %>% str() ### get orgnaic matter ACAI_isda_wide$soilOM_0_20 <- ACAI_isda_wide$carbon_organic_0_20 * 2 ACAI_isda_wide$soilOM_20_50 <- ACAI_isda_wide$carbon_organic_20_50 * 2 ACAI_isda_wide$percentSOM_0_20 <- ACAI_isda_wide$soilOM_0_20/10 ACAI_isda_wide$percentSOM_20_50 <- ACAI_isda_wide$soilOM_20_50/10 ## use pedo-transfer equations and get FC WP and SC for the two depth apart ##### WP ###### ACAI_isda_wide$wp_0_20A <- (-0.024*ACAI_isda_wide$sand_content_0_20/100) + 0.487*ACAI_isda_wide$clay_content_0_20/100 + 0.006*ACAI_isda_wide$percentSOM_0_20 + 0.005*(ACAI_isda_wide$sand_content_0_20/100 *ACAI_isda_wide$percentSOM_0_20) - 0.013*(ACAI_isda_wide$clay_content_0_20/100 * ACAI_isda_wide$percentSOM_0_20) + 0.068*(ACAI_isda_wide$sand_content_0_20/100 * ACAI_isda_wide$clay_content_0_20/100 ) + 0.031 ACAI_isda_wide$wp_0_20 <- (ACAI_isda_wide$wp_0_20A + (0.14 * ACAI_isda_wide$wp_0_20A - 0.02))*100 ACAI_isda_wide$wp_20_50A <- (-0.024*ACAI_isda_wide$sand_content_20_50/100) + 0.487*ACAI_isda_wide$clay_content_20_50/100 + 0.006*ACAI_isda_wide$percentSOM_20_50 + 0.005*(ACAI_isda_wide$sand_content_20_50/100 *ACAI_isda_wide$percentSOM_20_50) - 0.013*(ACAI_isda_wide$clay_content_20_50/100 * ACAI_isda_wide$percentSOM_20_50) + 0.068*(ACAI_isda_wide$sand_content_20_50/100 * ACAI_isda_wide$clay_content_20_50/100 ) + 0.031 ACAI_isda_wide$wp_20_50 <- (ACAI_isda_wide$wp_20_50A + (0.14 * ACAI_isda_wide$wp_20_50A - 0.02))*100 ACAI_isda_wide <- subset(ACAI_isda_wide, select=-c(wp_0_20A, wp_20_50A)) ##### FC ###### ACAI_isda_wide$FC_0_20_A <- -0.251*ACAI_isda_wide$sand_content_0_20/100 + 0.195*ACAI_isda_wide$clay_content_0_20/100 + 0.011*ACAI_isda_wide$percentSOM_0_20 + 0.006*(ACAI_isda_wide$sand_content_0_20/100 * ACAI_isda_wide$percentSOM_0_20) - 0.027*(ACAI_isda_wide$clay_content_0_20/100 * ACAI_isda_wide$percentSOM_0_20) + 0.452*(ACAI_isda_wide$sand_content_0_20/100 * ACAI_isda_wide$clay_content_0_20/100) + 0.299 ACAI_isda_wide$FC_0_20 <- (ACAI_isda_wide$FC_0_20_A + (1.283 * ACAI_isda_wide$FC_0_20_A^2 - 0.374 * ACAI_isda_wide$FC_0_20_A - 0.015))*100 ACAI_isda_wide$FC_20_50_A <- -0.251*ACAI_isda_wide$sand_content_20_50/100 + 0.195*ACAI_isda_wide$clay_content_20_50/100 + 0.011*ACAI_isda_wide$percentSOM_20_50 + 0.006*(ACAI_isda_wide$sand_content_20_50/100 * ACAI_isda_wide$percentSOM_20_50) - 0.027*(ACAI_isda_wide$clay_content_20_50/100 * ACAI_isda_wide$percentSOM_20_50) + 0.452*(ACAI_isda_wide$sand_content_20_50/100 * ACAI_isda_wide$clay_content_20_50/100) + 0.299 ACAI_isda_wide$FC_20_50 <- (ACAI_isda_wide$FC_20_50_A + (1.283 * ACAI_isda_wide$FC_20_50_A^2 - 0.374 * ACAI_isda_wide$FC_20_50_A - 0.015))*100 ACAI_isda_wide <- subset(ACAI_isda_wide, select=-c(FC_0_20_A, FC_20_50_A)) ##### soil water at ###### ACAI_isda_wide$sws_0_20_A <- 0.278*(ACAI_isda_wide$sand_content_0_20/100)+0.034*(ACAI_isda_wide$clay_content_0_20/100)+0.022*ACAI_isda_wide$percentSOM_0_20- 0.018*(ACAI_isda_wide$sand_content_0_20/100*ACAI_isda_wide$percentSOM_0_20)- 0.027*(ACAI_isda_wide$clay_content_0_20/100*ACAI_isda_wide$percentSOM_0_20)-0.584*(ACAI_isda_wide$sand_content_0_20/100*ACAI_isda_wide$clay_content_0_20/100)+0.078 ACAI_isda_wide$sws_0_20_B <- (ACAI_isda_wide$sws_0_20_A +(0.636*ACAI_isda_wide$sws_0_20_A-0.107))*100 ACAI_isda_wide$sws_0_20 <- (ACAI_isda_wide$FC_0_20/100+ACAI_isda_wide$sws_0_20_B/100-(0.097*ACAI_isda_wide$sand_content_0_20/100)+0.043)*100 ACAI_isda_wide$sws_20_50_A <- 0.278*(ACAI_isda_wide$sand_content_20_50/100)+0.034*(ACAI_isda_wide$clay_content_20_50/100)+0.022*ACAI_isda_wide$percentSOM_20_50- 0.018*(ACAI_isda_wide$sand_content_20_50/100*ACAI_isda_wide$percentSOM_20_50)- 0.027*(ACAI_isda_wide$clay_content_20_50/100*ACAI_isda_wide$percentSOM_20_50)-0.584*(ACAI_isda_wide$sand_content_20_50/100*ACAI_isda_wide$clay_content_20_50/100)+0.078 ACAI_isda_wide$sws_20_50_B <- (ACAI_isda_wide$sws_20_50_A +(0.636*ACAI_isda_wide$sws_20_50_A-0.107))*100 ACAI_isda_wide$sws_20_50 <- (ACAI_isda_wide$FC_20_50/100+ACAI_isda_wide$sws_20_50_B/100-(0.097*ACAI_isda_wide$sand_content_20_50/100)+0.043)*100 ACAI_isda_wide <- subset(ACAI_isda_wide, select=-c(sws_0_20_A, sws_0_20_B, sws_20_50_A, sws_20_50_B)) str(ACAI_isda_wide) ACAI_isda_wide$location <- paste(ACAI_isda_wide$lon, ACAI_isda_wide$lat, sep="_") setwd("/home/akilimo/projects/SoilNPK") dsY <- readRDS("NOT_formodel.RDS") ### get trial id and GPS head(dsY) setwd("/home/akilimo/projects/SoilNPK/Rwanda") ### add gps to blup by trilaiD ranall <- readRDS("ranw_allDAP.RDS") ranDAP15 <- readRDS("ranw_DAP15.RDS") kk <- unique(dsY[, c("trialID", "lon", "lat", "harvestDate_DOY")]) kk <- ddply(kk, .(trialID), summarize, lon=mean(lon), lat=mean(lat), harvestDate_DOY= mean(harvestDate_DOY)) tt <- as.data.frame(table(kk$trialID)) tt[tt$Freq>1, ] ranall <- unique(merge(ranall, kk, by="trialID")) ranDAP15 <- unique(merge(ranDAP15, kk, by="trialID")) ranall$location <- paste(ranall$lon, ranall$lat, sep="_") ranDAP15$location <- paste(ranDAP15$lon, ranDAP15$lat, sep="_") ACAI_isda_wide_allDAP <- unique(merge(ACAI_isda_wide, unique(subset(ranall, select=-c(lon, lat))), by="location" )) ## get trialID ACAI_isda_wide_DAP15 <- unique(merge(ACAI_isda_wide, unique(subset(ranDAP15, select=-c(lon, lat))), by="location" )) ## get trialID ACAI_isda_wide_allDAP$country <- "undefined" ACAI_isda_wide_allDAP$country[grep("NG", ACAI_isda_wide_allDAP$trialID)] <- "NG" ACAI_isda_wide_allDAP$country[grep("TZ", ACAI_isda_wide_allDAP$trialID)] <- "TZ" unique(ACAI_isda_wide_allDAP$country ) ACAI_isda_wide_DAP15$country <- "undefined" ACAI_isda_wide_DAP15$country[grep("NG", ACAI_isda_wide_DAP15$trialID)] <- "NG" ACAI_isda_wide_DAP15$country[grep("TZ", ACAI_isda_wide_DAP15$trialID)] <- "TZ" unique(ACAI_isda_wide_DAP15$country ) head(ACAI_isda_wide_allDAP) saveRDS(ACAI_isda_wide_allDAP, "ACAI_isda_GPS_allDAP.RDS") saveRDS(ACAI_isda_wide_DAP15, "ACAI_isda_GPS_DAP15.RDS") nrow(ACAI_isda_wide_allDAP) nrow(ACAI_isda_wide_DAP15) #################################################################################################################################################### ## prepare data for reverse QUEFTS with isric soil #################################################################################################################################################### setwd("/home/akilimo/projects/SoilNPK/Rwanda") ISRIC_ACAI <- readRDS("ISRIC_ACAI_Cassava.RDS") head(ISRIC_ACAI) ds_GPS_AC <- unique(merge(ISRIC_ACAI, kk[, c("trialID","harvestDate_DOY")]), by="trialID") head(ds_GPS_AC) setwd("/home/akilimo/projects/SoilNPK/Rwanda") ### add gps to blup by trilaiD ranall <- readRDS("ranw_allDAP.RDS") ranDAP15 <- readRDS("ranw_DAP15.RDS") ds_GPS_AC_DAPall <- unique(merge(ds_GPS_AC, ranall, by="trialID")) ds_GPS_AC_DAP15 <- unique(merge(ds_GPS_AC, ranDAP15, by="trialID")) ds_GPS_AC_DAPall$country <- "undefined" ds_GPS_AC_DAPall$country[grep("NG", ds_GPS_AC_DAPall$trialID)] <- "NG" ds_GPS_AC_DAPall$country[grep("TZ", ds_GPS_AC_DAPall$trialID)] <- "TZ" unique(ds_GPS_AC_DAPall$country ) ds_GPS_AC_DAP15$country <- "undefined" ds_GPS_AC_DAP15$country[grep("NG", ds_GPS_AC_DAP15$trialID)] <- "NG" ds_GPS_AC_DAP15$country[grep("TZ", ds_GPS_AC_DAP15$trialID)] <- "TZ" unique(ds_GPS_AC_DAP15$country ) saveRDS(ds_GPS_AC_DAPall, "ACAI_isric_GPS_allDAP.RDS") saveRDS(ds_GPS_AC_DAP15, "ACAI_isric_GPS_DAP15.RDS") ######################################################################################### ### change fresh matter to dry matter ######################################################################################### fd <- read.csv("/home/akilimo/lintul/lintul/dataSources/fd2.csv") getRFY <- function(HD, RDY, country = c("NG", "TZ"), fd){ d <- HD # fd <- read.csv("/home/akilimo/lintul/dataSources/fd2.csv") DC <- merge(data.frame(dayNr=d), fd[fd$country==country,], sort=FALSE)$DMCont RFY <- RDY / DC * 100 return(RFY) } setwd("/home/akilimo/projects/SoilNPK/Rwanda") ACAI_isric_GPS_allDAP <- readRDS("ACAI_isric_GPS_allDAP.RDS") ACAI_isric_GPS_DAP15 <- readRDS("ACAI_isric_GPS_DAP15.RDS") ACAI_isda_wide_allDAP <- readRDS("ACAI_isda_GPS_allDAP.RDS") ACAI_isda_wide_DAP15 <- readRDS("ACAI_isda_GPS_DAP15.RDS") getRDY <- function(HD, RFY, country, fd){ #SHORT DEF: Function to convert root FM yield into root dry matter yield (RDY): user define CY in FM in ton/ha, QUEFTS require Cy in DM kg/ha #RETURNS: RDY: root dry yield in the same units as root FM yield input #INPUT: HD: harvest date (Date format) # RFY: root fresh matter yield (user's units) # country = c("NG", "TZ") d <- HD # fd <- read.csv("/home/akilimo/lintul/dataSources/fd2.csv") DC <- merge(data.frame(dayNr=d), fd[fd$country==country,], sort=FALSE)$DMCont RDY <- (RFY * DC)/100 return(RDY) } allTreat_DM <- function(ds_GPS_AC){ ds_GPS_AC$harvestDate_DOY <- round(ds_GPS_AC$harvestDate_DOY, digits = 0) ds_GPS_AC$DM_NPK <- NULL ds_GPS_AC$DM_PK <- NULL ds_GPS_AC$DM_NK <- NULL ds_GPS_AC$DM_NP <- NULL ds_GPS_AC$DM_half_NPK <- NULL ds_GPS_AC$DM_CON <- NULL for(k in 1:nrow(ds_GPS_AC)){ print(k) ds_GPS_AC$DM_NPK[k] <- getRDY(HD=ds_GPS_AC$harvestDate_DOY[k], RFY = ds_GPS_AC$NPK[k], country= ds_GPS_AC$country[k], fd=fd)*1000 ds_GPS_AC$DM_PK[k] <- getRDY(HD=ds_GPS_AC$harvestDate_DOY[k], RFY = ds_GPS_AC$PK[k], country = ds_GPS_AC$country[k], fd=fd)*1000 ds_GPS_AC$DM_NK[k] <- getRDY(HD=ds_GPS_AC$harvestDate_DOY[k], RFY = ds_GPS_AC$NK[k], country = ds_GPS_AC$country[k], fd=fd)*1000 ds_GPS_AC$DM_NP[k] <- getRDY(HD=ds_GPS_AC$harvestDate_DOY[k], RFY = ds_GPS_AC$NP[k], country = ds_GPS_AC$country[k], fd=fd)*1000 ds_GPS_AC$DM_half_NPK[k] <- getRDY(HD=ds_GPS_AC$harvestDate_DOY[k], RFY = ds_GPS_AC$half_NPK[k], country = ds_GPS_AC$country[k], fd=fd)*1000 ds_GPS_AC$DM_CON[k] <- getRDY(HD=ds_GPS_AC$harvestDate_DOY[k], RFY = ds_GPS_AC$CON[k], country = ds_GPS_AC$country[k], fd=fd)*1000 } return(ds_GPS_AC) } ACAI_isric_GPS_allDAP <- allTreat_DM (ds_GPS_AC=ACAI_isric_GPS_allDAP) ACAI_isric_GPS_DAP15 <- allTreat_DM (ds_GPS_AC=ACAI_isric_GPS_DAP15) ACAI_isda_wide_allDAP <- allTreat_DM (ds_GPS_AC=ACAI_isda_wide_allDAP) ACAI_isda_wide_DAP15 <- allTreat_DM (ds_GPS_AC=ACAI_isda_wide_DAP15) ACAI_isda_wide_allDAP$fcc <- gsub(",", "_", ACAI_isda_wide_allDAP$fcc) ACAI_isda_wide_allDAP$fcc <- gsub(" ", "", ACAI_isda_wide_allDAP$fcc) ACAI_isda_wide_DAP15$fcc <- gsub(",", "_", ACAI_isda_wide_DAP15$fcc) ACAI_isda_wide_DAP15$fcc <- gsub(" ", "", ACAI_isda_wide_DAP15$fcc) ACAI_isda_wide_DAP15$fcc <- as.factor(ACAI_isda_wide_DAP15$fcc) ACAI_isda_wide_allDAP$fcc <- as.factor(ACAI_isda_wide_allDAP$fcc) ######################################################################################### ### Revers QUEFTS ######################################################################################### RecoveryFraction <- data.frame(rec_N=0.5, rec_P=0.21, rec_K=0.49) source("QUEFTS_function.R") head(ds_GPS_AC) getsoilNPK_valresult <- function(ds_GPS_AC){ soilNPK_AC <- NULL validation_result_AC <- NULL for(i in 1:nrow(ds_GPS_AC)){ oneTrial <- ds_GPS_AC[i,] WLY <- oneTrial$DM_NPK# is in kg/ha dry root addedFertilizer <- as.data.frame(matrix(nrow=5,ncol=3)) colnames(addedFertilizer) <- c("FN", "FP", "FK") addedFertilizer[,1] <- c(0, 0, 150, 150, 75) addedFertilizer[,2] <- c(0, 40, 0, 40, 20) addedFertilizer[,3] <- c(0, 180, 180, 0, 90) YieldMatrix <- as.data.frame(matrix(nrow=5,ncol=2)) colnames(YieldMatrix) <- c("treatment", "NOT_yield") YieldMatrix[,1] <- c("Control", "PK", "NK", "NP", "halfNPK") YieldMatrix[,2] <- c(oneTrial$DM_CON, oneTrial$DM_PK, oneTrial$DM_NK, oneTrial$DM_NP, oneTrial$DM_half_NPK) ## optimizing by min TSS. soil_NPK <- optim(par=c(0,0,0), optim_INS, lower=c(0.1, 0.1, 0.1), method = "L-BFGS-B", addedFertilizer=addedFertilizer, YieldMatrix=YieldMatrix, WLY=WLY)$par Yield_halfnpk <- Yield_S(SN = addedFertilizer[5,1] + soil_NPK[1] , SP = addedFertilizer[5,2] + soil_NPK[2], SK = addedFertilizer[5,3] + soil_NPK[3], WLY = WLY) oneTrial$halfNPK_observed <- oneTrial$DM_half_NPK oneTrial$halfNPK_estimated <- Yield_halfnpk Yield_control <- Yield_S(SN = addedFertilizer[1,1] + soil_NPK[1] , SP = addedFertilizer[1,2] + soil_NPK[2], SK = addedFertilizer[1,3] + soil_NPK[3], WLY = WLY) oneTrial$Control_observed <- oneTrial$DM_CON oneTrial$Control_estimated <- Yield_control validation_result_AC <- rbind(validation_result_AC, oneTrial) snpk <- oneTrial snpk$soilN <- soil_NPK[1] snpk$soilP <- soil_NPK[2] snpk$soilK <- soil_NPK[3] soilNPK_AC <- rbind(soilNPK_AC, snpk) } return(list(soilNPK_AC, validation_result_AC)) } ACAI_isric_GPS_allDAP_soilval <- getsoilNPK_valresult(ACAI_isric_GPS_allDAP) ACAI_isric_GPS_DAP15_soilval <- getsoilNPK_valresult(ACAI_isric_GPS_DAP15) ACAI_isda_wide_allDAP_soilval <- getsoilNPK_valresult(ACAI_isda_wide_allDAP) ACAI_isda_wide_DAP15_soilval <- getsoilNPK_valresult(ACAI_isda_wide_DAP15) ACAI_isric_GPS_allDAP_soilNPK <- ACAI_isric_GPS_allDAP_soilval[[1]] ACAI_isric_GPS_allDAP_VAl <- ACAI_isric_GPS_allDAP_soilval[[2]] ACAI_isric_GPS_DAP15_soilNPK <- ACAI_isric_GPS_DAP15_soilval[[1]] ACAI_isric_GPS_DAP15_VAl <- ACAI_isric_GPS_DAP15_soilval[[2]] ACAI_isda_wide_allDAP_soilNPK <- ACAI_isda_wide_allDAP_soilval[[1]] ACAI_isda_wide_allDAP_VAl <- ACAI_isda_wide_allDAP_soilval[[2]] ACAI_isda_wide_DAP15_soilNPK <- ACAI_isda_wide_DAP15_soilval[[1]] ACAI_isda_wide_DAP15_VAl <- ACAI_isda_wide_DAP15_soilval[[2]] head(ACAI_isric_GPS_allDAP_VAl) head(ACAI_isric_GPS_allDAP_soilNPK) saveRDS(ACAI_isric_GPS_allDAP_VAl, "validation_AC_isric_allDAP.RDS") saveRDS(ACAI_isric_GPS_allDAP_soilNPK, "soilNPK_AC_isric_allDAP.rds") saveRDS(ACAI_isric_GPS_DAP15_VAl, "validation_AC_isric_DAP15.RDS") saveRDS(ACAI_isric_GPS_DAP15_soilNPK, "soilNPK_AC_isric_DAP15.rds") saveRDS(ACAI_isda_wide_allDAP_VAl, "validation_AC_isda_allDAP.RDS") saveRDS(ACAI_isda_wide_allDAP_soilNPK, "soilNPK_AC_isda_allDAP.rds") saveRDS(ACAI_isda_wide_DAP15_VAl, "validation_AC_isda_DAP15.RDS") saveRDS(ACAI_isda_wide_DAP15_soilNPK, "soilNPK_AC_isda_DAP15.rds") ## translate dry matter estimates to fresh matter setwd("/home/akilimo/projects/SoilNPK/Rwanda") ACAI_isric_GPS_allDAP_VAl <- readRDS("validation_AC_isric_allDAP.RDS") ACAI_isric_GPS_allDAP_soilNPK <- readRDS("soilNPK_AC_isric_allDAP.rds") ACAI_isric_GPS_DAP15_VAl <- readRDS("validation_AC_isric_DAP15.RDS") ACAI_isric_GPS_DAP15_soilNPK <- readRDS("soilNPK_AC_isric_DAP15.rds") ACAI_isda_wide_allDAP_VAl <- readRDS("validation_AC_isda_allDAP.RDS") ACAI_isda_wide_allDAP_soilNPK <- readRDS("soilNPK_AC_isda_allDAP.rds") ACAI_isda_wide_DAP15_VAl <- readRDS("validation_AC_isda_DAP15.RDS") ACAI_isda_wide_DAP15_soilNPK <- readRDS("soilNPK_AC_isda_DAP15.rds") ACAI_isda_data <- subset(ACAI_isda_wide_allDAP_soilNPK, select = -c(NPK, PK,NK,NP,half_NPK,CON, DM_NPK,DM_PK,DM_NK,DM_NP,DM_half_NPK, DM_CON,halfNPK_observed,halfNPK_estimated, Control_observed,Control_estimated,soilN, soilP,soilK)) saveRDS(ACAI_isda_data, "ACAI_isda_data.RDS") ### converting QUEFTS output (dry matter) to fresh matter getVAlresult_FW <- function(validation_result_AC, fd){ validation_result_AC$CON_estFW <- NULL for(k in 1:nrow(validation_result_AC)){ validation_result_AC$CON_estFW[k] <- getRFY(HD=validation_result_AC$harvestDate_DOY[k], RDY = validation_result_AC$Control_estimated[k], country = validation_result_AC$country[k], fd=fd) } validation_result_AC$halfNPK_estFW <- NULL for(k in 1:nrow(validation_result_AC)){ validation_result_AC$halfNPK_estFW[k] <- getRFY(HD=validation_result_AC$harvestDate_DOY[k], RDY = validation_result_AC$halfNPK_estimated[k], country = validation_result_AC$country[k], fd=fd) } validation_result_AC$cc <- ifelse(validation_result_AC$country == "NG", "Nigeria", "Tanzania") return(validation_result_AC) } ACAI_isric_GPS_allDAP_VAl <- getVAlresult_FW(validation_result_AC=ACAI_isric_GPS_allDAP_VAl, fd=fd) ACAI_isric_GPS_DAP15_VAl <- getVAlresult_FW(validation_result_AC=ACAI_isric_GPS_DAP15_VAl, fd=fd) ACAI_isda_GPS_allDAP_VAl <- getVAlresult_FW(validation_result_AC=ACAI_isda_wide_allDAP_VAl, fd=fd) ACAI_isda_GPS_DAP15_VAl <- getVAlresult_FW(validation_result_AC=ACAI_isda_wide_DAP15_VAl, fd=fd) head(ACAI_isric_GPS_allDAP_VAl) # ACAI_isric_GPS_allDAP_VAl$index <- c(1:nrow(ACAI_isric_GPS_allDAP_VAl)) # ACAI_isric_GPS_allDAP_VAl[ACAI_isric_GPS_allDAP_VAl$CON_estFW > 80000, ]##265 # ACAI_isric_GPS_allDAP_VAl[ACAI_isric_GPS_allDAP_VAl$CON_estFW > 40000 & ACAI_isric_GPS_allDAP_VAl$country=="TZ", ]##342 # validation_result_AC <- validation_result_AC[!validation_result_AC$index %in% c(265, 342), ] ggcontrol <- ggplot(ACAI_isda_GPS_allDAP_VAl, aes(CON, CON_estFW/1000)) + geom_point() + geom_abline(intercept=0, slope=1, col="grey") + geom_smooth(method=lm) + facet_wrap(~cc, scales="free") + stat_density2d() + xlab("Observed control yield (t/ha)") + ylab("Estimated control yield using reverse QUEFTS (kg/ha)") + ggtitle("ACAI all") + theme_bw()+ theme(axis.text = element_text(size=16), strip.text.x = element_text(size=16)) gghalf <- ggplot(ACAI_isric_GPS_allDAP_VAl, aes(half_NPK, halfNPK_estFW/1000)) + geom_point() + geom_abline(intercept=0, slope=1, col="grey") + geom_smooth(method=lm) + facet_wrap(~cc, scales="free") + stat_density2d() + xlab("Observed half NPK (t/ha)") + ylab("Estimated half NPK using reverse QUEFTS (kg/ha)") + ggtitle("ACAI all") + theme_bw()+ theme(axis.text = element_text(size=16), strip.text.x = element_text(size=16)) ggsave("Control_Observed_QUEFTS_isda_All.pdf", ggcontrol, width = 6, height = 6) ggsave("halfNPK_Observed_QUEFTS_isda_all.pdf", gghalf, width = 6, height = 6) ggsave("Control_Observed_QUEFTS_isda_DAP15.pdf", ggcontrol, width = 6, height = 6) ggsave("halfNPK_Observed_QUEFTS_isda_DAP15.pdf", gghalf, width = 6, height = 6) ################################################################################# ## random forest model for AC data: trying both isric and iSDA data ################################################################################# ACAI_isric_GPS_allDAP_soilNPK ACAI_isric_GPS_DAP15_soilNPK ACAI_isda_wide_allDAP_soilNPK ACAI_isda_wide_DAP15_soilNPK GIS_soilINS_modData_AC <- ACAI_isric_GPS_DAP15_soilNPK[, c("soilN", "soilP", "soilK", "exchK", "P_Mehlich", "Clay_5","Clay_15","Clay_30","percentSOM_5","percentSOM_15","percentSOM_30", "pH_5","pH_15","pH_30", "silt_5", "silt_15", "silt_30","BD_5", "BD_15", "BD_30", "CEC_5","CEC_15","CEC_30", "percentSOC_5", "percentSOC_15", "percentSOC_30", "FC_5", "FC_15", "FC_30", "wp_5", "wp_15", "wp_30", "sws_5", "sws_15", "sws_30", "TotalN", "Mn","B", "Ca","Fe", "Cu","Al", "Mg", "Na","ncluster", "country", "CON", "trialID")] GIS_soilINS_modData_AC$Clay_5 <- as.numeric(GIS_soilINS_modData_AC$Clay_5) GIS_soilINS_modData_AC$Clay_15 <- as.numeric(GIS_soilINS_modData_AC$Clay_15) GIS_soilINS_modData_AC$Clay_30 <- as.numeric(GIS_soilINS_modData_AC$Clay_30) GIS_soilINS_modData_AC$silt_5 <- as.numeric(GIS_soilINS_modData_AC$silt_5) GIS_soilINS_modData_AC$silt_15 <- as.numeric(GIS_soilINS_modData_AC$silt_15) GIS_soilINS_modData_AC$silt_30 <- as.numeric(GIS_soilINS_modData_AC$silt_30) GIS_soilINS_modData_AC$BD_5 <- as.numeric(GIS_soilINS_modData_AC$BD_5) GIS_soilINS_modData_AC$BD_15 <- as.numeric(GIS_soilINS_modData_AC$BD_15) GIS_soilINS_modData_AC$BD_30 <- as.numeric(GIS_soilINS_modData_AC$BD_30) GIS_soilINS_modData_AC$CEC_5 <- as.numeric(GIS_soilINS_modData_AC$CEC_5) GIS_soilINS_modData_AC$CEC_15 <- as.numeric(GIS_soilINS_modData_AC$CEC_15) GIS_soilINS_modData_AC$CEC_30 <- as.numeric(GIS_soilINS_modData_AC$CEC_30) GIS_soilINS_modData_AC$TotalN <- as.numeric(GIS_soilINS_modData_AC$TotalN) GIS_soilINS_modData_AC$Mn <- as.numeric(GIS_soilINS_modData_AC$Mn) GIS_soilINS_modData_AC$B <- as.numeric(GIS_soilINS_modData_AC$B) GIS_soilINS_modData_AC$Ca <- as.numeric(GIS_soilINS_modData_AC$Ca) GIS_soilINS_modData_AC$Fe <- as.numeric(GIS_soilINS_modData_AC$Fe) GIS_soilINS_modData_AC$Cu <- as.numeric(GIS_soilINS_modData_AC$Cu) GIS_soilINS_modData_AC$Al <- as.numeric(GIS_soilINS_modData_AC$Al) GIS_soilINS_modData_AC$Mg <- as.numeric(GIS_soilINS_modData_AC$Mg) GIS_soilINS_modData_AC$Na <- as.numeric(GIS_soilINS_modData_AC$Na) GIS_soilINS_modData_AC$ncluster <- as.factor(GIS_soilINS_modData_AC$ncluster) GIS_soilINS_modData_AC$CON <- as.numeric(GIS_soilINS_modData_AC$CON) GIS_soilINS_modData_AC$country <- as.factor(GIS_soilINS_modData_AC$country ) ### create classes of control GIS_soilINS_modData_AC$CONclass <- ifelse(GIS_soilINS_modData_AC$CON < 7.5, "class1", ifelse(GIS_soilINS_modData_AC$CON >= 7.5 & GIS_soilINS_modData_AC$CON < 15, "class2", ifelse(GIS_soilINS_modData_AC$CON >= 15 & GIS_soilINS_modData_AC$CON < 22.5, "class3", ifelse(GIS_soilINS_modData_AC$CON >= 22.5 & GIS_soilINS_modData_AC$CON < 30, "class4", "class5")))) GIS_soilINS_modData_AC$CONclass <- as.factor(GIS_soilINS_modData_AC$CONclass) ACAI_isric_GPS_allDAP_soilNPK <- GIS_soilINS_modData_AC ACAI_isric_GPS_DAP15_soilNPK <- GIS_soilINS_modData_AC saveRDS(ACAI_isric_GPS_allDAP_soilNPK, "ACAI_isric_GPS_allDAP_soilNPK.RDS") saveRDS(ACAI_isric_GPS_DAP15_soilNPK, "ACAI_isric_GPS_DAP15_soilNPK.RDS") nrow(ACAI_isric_GPS_allDAP_soilNPK) nrow(ACAI_isric_GPS_DAP15_soilNPK) ACAI_isric_GPS_allDAP_soilNPK <- readRDS("ACAI_isric_GPS_allDAP_soilNPK.RDS") ACAI_isric_GPS_DAP15_soilNPK <- readRDS("ACAI_isric_GPS_DAP15_soilNPK.RDS") ## "bedrock_depth_0_200", has many missing data so it is taken out any way the values are always 198 or 200 ACAI_isda_wide_allDAP_soilNPK <- ACAI_isda_wide_allDAP_soilNPK[, c("soilN", "soilP", "soilK","aluminium_extractable_0_20", "aluminium_extractable_20_50", "bulk_density_0_20","bulk_density_20_50", "calcium_extractable_0_20", "calcium_extractable_20_50", "carbon_organic_0_20", "carbon_organic_20_50", "carbon_total_0_20", "carbon_total_20_50", "cation_exchange_capacity_0_20", "cation_exchange_capacity_20_50","clay_content_0_20","clay_content_20_50","iron_extractable_0_20","iron_extractable_20_50", "magnesium_extractable_0_20", "magnesium_extractable_20_50", "nitrogen_total_0_20", "nitrogen_total_20_50", "ph_0_20", "ph_20_50", "phosphorous_extractable_0_20", "phosphorous_extractable_20_50","potassium_extractable_0_20", "potassium_extractable_20_50", "sand_content_0_20","sand_content_20_50", "silt_content_0_20", "silt_content_20_50", "stone_content_0_20", "stone_content_20_50", "sulphur_extractable_0_20", "sulphur_extractable_20_50", "texture_class_0_20", "texture_class_20_50", "zinc_extractable_0_20", "zinc_extractable_20_50", "slope_angle", "fcc", "land_cover_2015", "land_cover_2016","land_cover_2017", "land_cover_2018","land_cover_2019", "soilOM_0_20", "soilOM_20_50", "percentSOM_0_20", "percentSOM_20_50", "wp_0_20", "wp_20_50", "FC_0_20", "FC_20_50", "sws_0_20", "sws_20_50", "country", "CON", "trialID")] ### create classes of control ACAI_isda_wide_allDAP_soilNPK$CONclass <- ifelse(ACAI_isda_wide_allDAP_soilNPK$CON < 7.5, "class1", ifelse(ACAI_isda_wide_allDAP_soilNPK$CON >= 7.5 & ACAI_isda_wide_allDAP_soilNPK$CON < 15, "class2", ifelse(ACAI_isda_wide_allDAP_soilNPK$CON >= 15 & ACAI_isda_wide_allDAP_soilNPK$CON < 22.5, "class3", ifelse(ACAI_isda_wide_allDAP_soilNPK$CON >= 22.5 & ACAI_isda_wide_allDAP_soilNPK$CON < 30, "class4", "class5")))) ACAI_isda_wide_allDAP_soilNPK$CONclass <- as.factor(ACAI_isda_wide_allDAP_soilNPK$CONclass) ACAI_isda_wide_DAP15_soilNPK <- ACAI_isda_wide_DAP15_soilNPK[, c("soilN", "soilP", "soilK","aluminium_extractable_0_20", "aluminium_extractable_20_50", "bulk_density_0_20","bulk_density_20_50", "calcium_extractable_0_20", "calcium_extractable_20_50", "carbon_organic_0_20", "carbon_organic_20_50", "carbon_total_0_20", "carbon_total_20_50", "cation_exchange_capacity_0_20", "cation_exchange_capacity_20_50","clay_content_0_20","clay_content_20_50","iron_extractable_0_20","iron_extractable_20_50", "magnesium_extractable_0_20", "magnesium_extractable_20_50", "nitrogen_total_0_20", "nitrogen_total_20_50", "ph_0_20", "ph_20_50", "phosphorous_extractable_0_20", "phosphorous_extractable_20_50","potassium_extractable_0_20", "potassium_extractable_20_50", "sand_content_0_20","sand_content_20_50", "silt_content_0_20", "silt_content_20_50", "stone_content_0_20", "stone_content_20_50", "sulphur_extractable_0_20", "sulphur_extractable_20_50", "texture_class_0_20", "texture_class_20_50", "zinc_extractable_0_20", "zinc_extractable_20_50", "slope_angle", "fcc", "land_cover_2015", "land_cover_2016","land_cover_2017", "land_cover_2018","land_cover_2019", "soilOM_0_20", "soilOM_20_50", "percentSOM_0_20", "percentSOM_20_50", "wp_0_20", "wp_20_50", "FC_0_20", "FC_20_50", "sws_0_20", "sws_20_50", "country", "CON", "trialID")] ### create classes of control ACAI_isda_wide_DAP15_soilNPK$CONclass <- ifelse(ACAI_isda_wide_DAP15_soilNPK$CON < 7.5, "class1", ifelse(ACAI_isda_wide_DAP15_soilNPK$CON >= 7.5 & ACAI_isda_wide_DAP15_soilNPK$CON < 15, "class2", ifelse(ACAI_isda_wide_DAP15_soilNPK$CON >= 15 & ACAI_isda_wide_DAP15_soilNPK$CON < 22.5, "class3", ifelse(ACAI_isda_wide_DAP15_soilNPK$CON >= 22.5 & ACAI_isda_wide_DAP15_soilNPK$CON < 30, "class4", "class5")))) ACAI_isda_wide_DAP15_soilNPK$CONclass <- as.factor(ACAI_isda_wide_DAP15_soilNPK$CONclass) setwd("/home/akilimo/projects/SoilNPK/Rwanda") saveRDS(ACAI_isda_wide_allDAP_soilNPK, "ACAI_isda_wide_allDAP_soilNPK2.RDS") saveRDS(ACAI_isda_wide_DAP15_soilNPK, "ACAI_isda_wide_DAP15_soilNPK2.RDS") ACAI_isda_wide_allDAP_soilNPK <- readRDS("ACAI_isda_wide_allDAP_soilNPK2.RDS") head(ACAI_isda_wide_allDAP_soilNPK) nrow(ACAI_isda_wide_allDAP_soilNPK) # # ### create classes of control # ## if the data is allowed to define cuts for discritization # require(caret) # require(rpart) # Q <- rpart(soilN ~ CON, data=GIS_soilINS_modData_AC, control=rpart.control(maxdepth = 3)) # plot(Q) # text(Q, use.n=TRUE, all=TRUE, cex = 0.8) # Q$frame # # tree_count <- ctree(soilN ~ CON, data=GIS_soilINS_modData_AC, controls=ctree_control(maxdepth = 2)) # plot(tree_count) # # # GIS_soilINS_modData_AC$Yield_con <- ifelse(GIS_soilINS_modData_AC$CON < 0, 0, as.numeric(as.character(GIS_soilINS_modData_AC$CON))) ## because of BLUPS control yield can be negative # # control_classes <- data.frame( # ConClass_1 = c(round_any(quantile(GIS_soilINS_modData_AC$Yield_con, probs=c(0.25, 0.75)), 10)), # ConClass_2 = c(1643, 4763)) # # # if(discrete == TRUE & C_class == "DataQuant"){ # GIS_soilINS_modData_AC$CONclass <- ifelse(GIS_soilINS_modData_AC$Yield_con < control_classes$ConClass_1[1], "low", # ifelse(GIS_soilINS_modData_AC$Yield_con >= control_classes$ConClass_1[1] & # GIS_soilINS_modData_AC$Yield_con < control_classes$ConClass_1[2], "medium", "high")) # # }else if (discrete == TRUE & C_class == "dataClass") { # GIS_soilINS_modData_AC$CONclass <- ifelse(GIS_soilINS_modData_AC$Yield_con < control_classes$ConClass_2[1], "low", # ifelse(GIS_soilINS_modData_AC$Yield_con >= control_classes$ConClass_2[1] & # GIS_soilINS_modData_AC$Yield_con < control_classes$ConClass_2[2], "medium", "high")) # # } harvestWeekVar100 <- readRDS("TestHarvestDate_100.RDS") harvestWeekVar92 <- readRDS("TestHarvestDate_92.RDS") harvestWeekVar127 <- readRDS("TestHarvestDate_127.RDS") harvestWeekVar159 <- readRDS("TestHarvestDate_159.RDS") harvestWeekVar547 <- readRDS("TestHarvestDate_547.RDS") harvestWeekVar587 <- readRDS("TestHarvestDate_587.RDS") harvestWeekVar1039 <- readRDS("TestHarvestDate_1039.RDS") harvestWeekVar1087 <- readRDS("TestHarvestDate_1087.RDS") harvestWeekVar174 <- readRDS("TestHarvestDate_174.RDS") harvestWeekVar609 <- readRDS("TestHarvestDate_609.RDS") harvestWeekVar1111 <- readRDS("TestHarvestDate_1111.RDS") harvestWeekVar183 <- readRDS("TestHarvestDate_183.RDS") harvestWeekVar621 <- readRDS("TestHarvestDate_621.RDS") harvestWeekVar1123 <- readRDS("TestHarvestDate_1123.RDS") harvestWeekVar <- rbind(harvestWeekVar100, harvestWeekVar92, harvestWeekVar127, harvestWeekVar159, harvestWeekVar547, harvestWeekVar587, harvestWeekVar1039, harvestWeekVar1087, harvestWeekVar174, harvestWeekVar609, harvestWeekVar1111, harvestWeekVar183, harvestWeekVar621, harvestWeekVar1123) ### Data partioning take 70:30 by country ACAI_isric_GPS_allDAP_soilNPK ACAI_isric_GPS_DAP15_soilNPK ACAI_isda_wide_allDAP_soilNPK ACAI_isda_wide_DAP15_soilNPK ACAI_isda_wide_allDAP_soilNPK <- unique(droplevels(ACAI_isda_wide_allDAP_soilNPK[complete.cases(ACAI_isda_wide_allDAP_soilNPK), ])) ACAI_isda_wide_DAP15_soilNPK <- unique(droplevels(ACAI_isda_wide_DAP15_soilNPK[complete.cases(ACAI_isda_wide_DAP15_soilNPK), ])) ACAI_isric_GPS_allDAP_soilNPK <- droplevels(ACAI_isric_GPS_allDAP_soilNPK[ACAI_isric_GPS_allDAP_soilNPK$trialID %in% unique(ACAI_isda_wide_allDAP_soilNPK$trialID), ]) ACAI_isric_GPS_DAP15_soilNPK <- droplevels(ACAI_isric_GPS_DAP15_soilNPK[ACAI_isric_GPS_DAP15_soilNPK$trialID %in% unique(ACAI_isda_wide_DAP15_soilNPK$trialID), ]) ACAI_isric_GPS_allDAP_soilNPK$country <- as.factor(ACAI_isric_GPS_allDAP_soilNPK$country) ACAI_isric_GPS_DAP15_soilNPK$country <- as.factor(ACAI_isric_GPS_DAP15_soilNPK$country) ACAI_isda_wide_allDAP_soilNPK$country <- as.factor(ACAI_isda_wide_allDAP_soilNPK$country) ACAI_isda_wide_DAP15_soilNPK$country <- as.factor(ACAI_isda_wide_DAP15_soilNPK$country) splitdata <- function(probs=c(0.7,0.3), soil_BLUPSData){ set.seed(444) trianData <- NULL testDatA <- NULL for(country in unique(soil_BLUPSData$country)){ Country_data <- droplevels(soil_BLUPSData[soil_BLUPSData$country == country,]) indexCountry <- sample(2, nrow(Country_data), replace=TRUE, prob=probs)## where conrtol yield is used as a covariate trainDataCountry <- Country_data[indexCountry == 1, ] testDataCountry <- Country_data[indexCountry == 2, ] trianData <- rbind(trianData, trainDataCountry) testDatA <- rbind(testDatA, testDataCountry) } return(list(trianData, testDatA)) } traintestData_isric_AC_all <- splitdata(soil_BLUPSData = ACAI_isric_GPS_allDAP_soilNPK) trainData_AC_isric_all <- traintestData_isric_AC_all[[1]] testData_AC_isric_all <- traintestData_isric_AC_all[[2]] traintestData_isric_AC_DAP15 <- splitdata(soil_BLUPSData = ACAI_isric_GPS_DAP15_soilNPK) trainData_AC_isric_DAP15 <- traintestData_isric_AC_DAP15[[1]] testData_AC_isric_DAP15 <- traintestData_isric_AC_DAP15[[2]] traintestData_isda_AC_all <- splitdata(soil_BLUPSData = ACAI_isda_wide_allDAP_soilNPK) trainData_AC_isda_all <- traintestData_isda_AC_all[[1]] testData_AC_isda_all <- traintestData_isda_AC_all[[2]] traintestData_isda_AC_DAP15 <- splitdata(soil_BLUPSData = ACAI_isda_wide_DAP15_soilNPK) trainData_AC_isda_DAP15 <- traintestData_isda_AC_DAP15[[1]] testData_AC_isda_DAP15 <- traintestData_isda_AC_DAP15[[2]] # trainData_AC_isda_all <- droplevels(ACAI_isda_wide_allDAP_soilNPK[ACAI_isda_wide_allDAP_soilNPK$trialID %in% trainData_AC_isric_all$trialID, ]) # testData_AC_isda_all <- droplevels(ACAI_isda_wide_allDAP_soilNPK[ACAI_isda_wide_allDAP_soilNPK$trialID %in% unique(testData_AC_isric_all$trialID), ]) # # trainData_AC_isda_DAP15 <- droplevels(ACAI_isda_wide_DAP15_soilNPK[ACAI_isda_wide_DAP15_soilNPK$trialID %in% trainData_AC_isric_DAP15$trialID, ]) # testData_AC_isda_DAP15 <- droplevels(ACAI_isda_wide_DAP15_soilNPK[ACAI_isda_wide_DAP15_soilNPK$trialID %in% unique(testData_AC_isric_DAP15$trialID), ]) nrow(unique(trainData_AC_isda_all)) nrow(trainData_AC_isric_all) trainData_AC_isric_all <-unique(subset(trainData_AC_isric_all, select=-c(CON, trialID))) testData_AC_isric_all <- unique(subset(testData_AC_isric_all, select=-c(CON, trialID))) trainData_AC_isda_all <- unique(subset(trainData_AC_isda_all, select=-c(CON, trialID))) testData_AC_isda_all <- unique(subset(testData_AC_isda_all, select=-c(CON, trialID))) trainData_AC_isric_DAP15 <- unique(subset(trainData_AC_isric_DAP15, select=-c(CON, trialID))) testData_AC_isric_DAP15 <- unique(subset(testData_AC_isric_DAP15, select=-c(CON, trialID))) trainData_AC_isda_DAP15 <- unique(subset(trainData_AC_isda_DAP15, select=-c(CON, trialID))) testData_AC_isda_DAP15 <- unique(subset(testData_AC_isda_DAP15, select=-c(CON, trialID))) ## isric all Ndata_Train_isric_all <- subset(trainData_AC_isric_all, select=-c(soilP, soilK)) Pdata_Train_isric_all <- subset(trainData_AC_isric_all, select=-c(soilN, soilK)) Kdata_Train_isric_all <- subset(trainData_AC_isric_all, select=-c(soilN, soilP)) Ndata_Valid_isric_all <- subset(testData_AC_isric_all, select=-c(soilP, soilK)) Pdata_Valid_isric_all <- subset(testData_AC_isric_all, select=-c(soilN, soilK)) Kdata_Valid_isric_all <- subset(testData_AC_isric_all, select=-c(soilN, soilP)) ## isda all Ndata_Train_isda_all <- subset(trainData_AC_isda_all, select=-c(soilP, soilK)) Pdata_Train_isda_all <- subset(trainData_AC_isda_all, select=-c(soilN, soilK)) Kdata_Train_isda_all <- subset(trainData_AC_isda_all, select=-c(soilN, soilP)) Ndata_Valid_isda_all <- subset(testData_AC_isda_all, select=-c(soilP, soilK)) Pdata_Valid_isda_all <- subset(testData_AC_isda_all, select=-c(soilN, soilK)) Kdata_Valid_isda_all <- subset(testData_AC_isda_all, select=-c(soilN, soilP)) ## isric DAP15 Ndata_Train_isric_DAP15 <- subset(trainData_AC_isric_DAP15, select=-c(soilP, soilK)) Pdata_Train_isric_DAP15 <- subset(trainData_AC_isric_DAP15, select=-c(soilN, soilK)) Kdata_Train_isric_DAP15 <- subset(trainData_AC_isric_DAP15, select=-c(soilN, soilP)) Ndata_Valid_isric_DAP15 <- subset(testData_AC_isric_DAP15, select=-c(soilP, soilK)) Pdata_Valid_isric_DAP15 <- subset(testData_AC_isric_DAP15, select=-c(soilN, soilK)) Kdata_Valid_isric_DAP15 <- subset(testData_AC_isric_DAP15, select=-c(soilN, soilP)) ## isda DAP15 Ndata_Train_isda_DAP15 <- subset(trainData_AC_isda_DAP15, select=-c(soilP, soilK)) Pdata_Train_isda_DAP15 <- subset(trainData_AC_isda_DAP15, select=-c(soilN, soilK)) Kdata_Train_isda_DAP15 <- subset(trainData_AC_isda_DAP15, select=-c(soilN, soilP)) Ndata_Valid_isda_DAP15 <- subset(testData_AC_isda_DAP15, select=-c(soilP, soilK)) Pdata_Valid_isda_DAP15 <- subset(testData_AC_isda_DAP15, select=-c(soilN, soilK)) Kdata_Valid_isda_DAP15 <- subset(testData_AC_isda_DAP15, select=-c(soilN, soilP)) require(gtools) ## Coustome control parameter #custom <- trainControl(method="repeatedcv", number=10, repeats=5, verboseIter=TRUE) custom <- trainControl(method="oob", number=10) ########################################################################## ## Random Forest soilN: isric all data = 0.70, isda all data = 0.80; isricd DAP15 = 0.84; isda DAP15 = 0.83 ########################################################################## set.seed(773) Ndata_Train <- Ndata_Train_isric_all##0.71 Ndata_Valid <- Ndata_Valid_isric_all Ndata_Train <- Ndata_Train_isda_all##0.79 Ndata_Valid <- Ndata_Valid_isda_all Ndata_Train <- Ndata_Train_isric_DAP15##0.78 Ndata_Valid <- Ndata_Valid_isric_DAP15 Ndata_Train <- Ndata_Train_isda_DAP15##0.82 Ndata_Valid <- Ndata_Valid_isda_DAP15 rm(RF_N1) RF_N1 <- randomForest(soilN ~ ., Ndata_Train, importance=TRUE, ntree=1000) Ndata_Valid$predN <- predict(RF_N1, Ndata_Valid) sst <- sum((Ndata_Valid$soilN - mean(Ndata_Valid$soilN))^2) sse <- sum((Ndata_Valid$soilN - Ndata_Valid$predN)^2) rsq <- 1 - sse / sst rsq varImpPlot(RF_N1) varImpPlot(RF_N1) varimportance <- data.frame(RF_N1$importance) varimportance$vars <- rownames(varimportance) varimportance <- varimportance[order(varimportance$X.IncMSE, decreasing = TRUE), ] rownames(varimportance) <- NULL varimportance$Importance <- c(1:nrow(varimportance)) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilN_isric_all.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilN_isda_all.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilN_isric_DAP15.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilN_isda_DAP15.csv", row.names = FALSE) ggn <- ggplot(Ndata_Valid, aes(soilN, predN, col=country)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + xlim(0,200) + ylim(0,200)+ xlab("soil N from QUEFTS optimization")+ ylab("perdicted soil N from GIS data")+ theme_bw() ggsave("AC_isric_all_QUEFTS_GIS_RF_soilN.pdf", ggn, width=5, height=5) ggsave("AC_isda_all_QUEFTS_GIS_RF_soilN.pdf", ggn, width=5, height=5) ggsave("AC_isric DAP15_QUEFTS_GIS_RF_soilN.pdf", ggn, width=5, height=5) ggsave("AC_isda DAP15_QUEFTS_GIS_RF_soilN.pdf", ggn, width=5, height=5) ########################################################################## ## Random Forest soilP: isric all =0.79; isad all = 0.81; isrics DAP15= 0.79; isda DAP15 = 0.72 ########################################################################## set.seed(773) Pdata_Train <- Pdata_Train_isric_all##0.69 Pdata_Valid <- Pdata_Valid_isric_all Pdata_Train <- Pdata_Train_isda_all##0.75 Pdata_Valid <- Pdata_Valid_isda_all Pdata_Train <- Pdata_Train_isric_DAP15##0.67 Pdata_Valid <- Pdata_Valid_isric_DAP15 Pdata_Train <- Pdata_Train_isda_DAP15##0.76 Pdata_Valid <- Pdata_Valid_isda_DAP15 rm(RF_P1) RF_P1 <- randomForest(soilP ~ ., Pdata_Train, importance=TRUE, ntree=1000) Pdata_Valid$predP <- predict(RF_P1, Pdata_Valid) sst <- sum((Pdata_Valid$soilP - mean(Pdata_Valid$soilP))^2) sse <- sum((Pdata_Valid$soilP - Pdata_Valid$predP)^2) rsq <- 1 - sse / sst rsq varImpPlot(RF_P1) varimportance <- data.frame(RF_P1$importance) varimportance$vars <- rownames(varimportance) varimportance <- varimportance[order(varimportance$X.IncMSE, decreasing = TRUE), ] rownames(varimportance) <- NULL varimportance$Importance <- c(1:nrow(varimportance)) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilP_isric_all.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilP_isda_all.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilP_isric_DAP15.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilP_isda_DAP15.csv", row.names = FALSE) ggp <- ggplot(Pdata_Valid, aes(soilP, predP, col=country)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + xlim(0,150) + ylim(0,150)+ xlab("soil P from QUEFTS optimization")+ ylab("perdicted soil P from GIS data")+ theme_bw() ggsave("AC_isric_all_QUEFTS_GIS_RF_soilP.pdf", ggp, width=5, height=5) ggsave("AC_isda_all_QUEFTS_GIS_RF_soilP.pdf", ggp, width=5, height=5) ggsave("AC_isric_DAP15_QUEFTS_GIS_RF_soilP.pdf", ggp, width=5, height=5) ggsave("AC_isda_DAP15_QUEFTS_GIS_RF_soilP.pdf", ggp, width=5, height=5) ########################################################################## ## Random Forest soilK" isric all = 0.54, isda all= 0.54; isirc DAP15 = 0.54; isda DAP15 = 0.37 ########################################################################## set.seed(773) Kdata_Train <- Kdata_Train_isric_all##0.40 Kdata_Valid <- Kdata_Valid_isric_all Kdata_Train <- Kdata_Train_isda_all##0.51 Kdata_Valid <- Kdata_Valid_isda_all Kdata_Train <- Kdata_Train_isric_DAP15##0.36 Kdata_Valid <- Kdata_Valid_isric_DAP15 Kdata_Train <- Kdata_Train_isda_DAP15##0.38 Kdata_Valid <- Kdata_Valid_isda_DAP15 rm(RF_K1) RF_K1 <- randomForest(soilK ~ ., Kdata_Train, importance=TRUE, ntree=1000) Kdata_Valid$predK <- predict(RF_K1, Kdata_Valid) sst <- sum((Kdata_Valid$soilK - mean(Kdata_Valid$soilK))^2) sse <- sum((Kdata_Valid$soilK - Kdata_Valid$predK)^2) rsq <- 1 - sse / sst rsq varImpPlot(RF_K1) varimportance <- data.frame(RF_P1$importance) varimportance$vars <- rownames(varimportance) varimportance <- varimportance[order(varimportance$X.IncMSE, decreasing = TRUE), ] rownames(varimportance) <- NULL varimportance$Importance <- c(1:nrow(varimportance)) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilK_isric_all.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilK_isda_all.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilK_isric_DAP15.csv", row.names = FALSE) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilK_isda_DAP15.csv", row.names = FALSE) ggk <- ggplot(Kdata_Valid, aes(soilK, predK, col=country)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + xlim(0,350) + ylim(0,350)+ xlab("soil K from QUEFTS optimization")+ ylab("perdicted soil K from GIS data")+ theme_bw() ggsave("AC_isric_all_QUEFTS_GIS_RF_soilK.pdf", ggk, width=5, height=5) ggsave("AC_isda_all_QUEFTS_GIS_RF_soilK.pdf", ggk, width=5, height=5) ggsave("AC_isric_DAP15_QUEFTS_GIS_RF_soilK.pdf", ggk, width=5, height=5) ggsave("AC_isda_DAP15_QUEFTS_GIS_RF_soilK.pdf", ggk, width=5, height=5) ######################################################################################################################################################## ## using ACAI data to train the modle and CIALCA data for validation ######################################################################################################################################################## ACAI_isric_GPS_allDAP_soilNPK ACAI_isric_GPS_DAP15_soilNPK ACAI_isda_wide_allDAP_soilNPK ACAI_isda_wide_DAP15_soilNPK setwd("/home/akilimo/projects/SoilNPK/Rwanda") soilNPK <- readRDS("soilNPK_isrics.rds") Cialca_isric_soilnpk <- soilNPK[, c("soilN", "soilP", "soilK", "exchK", "P_Mehlich", "Clay_5","Clay_15","Clay_30","percentSOM_5","percentSOM_15","percentSOM_30", "pH_5","pH_15","pH_30", "silt_5", "silt_15", "silt_30","BD_5", "BD_15", "BD_30", "CEC_5","CEC_15","CEC_30", "percentSOC_5", "percentSOC_15", "percentSOC_30", "FC_5", "FC_15", "FC_30", "wp_5", "wp_15", "wp_30", "sws_5", "sws_15", "sws_30", "TotalN", "Mn","B", "Ca","Fe", "Cu","Al", "Mg", "Na","ncluster", "Country", "Control")] cialca_soilNPK_isad <- readRDS("soilNPK_isda.rds") str(cialca_soilNPK_isad) Cialca_isda_soilnpk <- cialca_soilNPK_isad[, c("soilN", "soilP", "soilK","aluminium_extractable_0_20", "aluminium_extractable_20_50", "bedrock_depth_0_200", "bulk_density_0_20","bulk_density_20_50", "calcium_extractable_0_20", "calcium_extractable_20_50", "carbon_organic_0_20", "carbon_organic_20_50", "carbon_total_0_20", "carbon_total_20_50", "cation_exchange_capacity_0_20", "cation_exchange_capacity_20_50","clay_content_0_20","clay_content_20_50","iron_extractable_0_20","iron_extractable_20_50", "magnesium_extractable_0_20", "magnesium_extractable_20_50", "nitrogen_total_0_20", "nitrogen_total_20_50", "ph_0_20", "ph_20_50", "phosphorous_extractable_0_20", "phosphorous_extractable_20_50","potassium_extractable_0_20", "potassium_extractable_20_50", "sand_content_0_20","sand_content_20_50", "silt_content_0_20", "silt_content_20_50", "stone_content_0_20", "stone_content_20_50", "sulphur_extractable_0_20", "sulphur_extractable_20_50", "texture_class_0_20", "texture_class_20_50", "zinc_extractable_0_20", "zinc_extractable_20_50", "slope_angle", "fcc", "land_cover_2015", "land_cover_2016","land_cover_2017", "land_cover_2018","land_cover_2019", "soilOM_0_20", "soilOM_20_50", "percentSOM_0_20", "percentSOM_20_50", "wp_0_20", "wp_20_50", "FC_0_20", "FC_20_50", "sws_0_20", "sws_20_50", "Country", "Control")] Cialca_isric_soilnpk$Clay_5 <- as.numeric(Cialca_isric_soilnpk$Clay_5) Cialca_isric_soilnpk$Clay_15 <- as.numeric(Cialca_isric_soilnpk$Clay_15) Cialca_isric_soilnpk$Clay_30 <- as.numeric(Cialca_isric_soilnpk$Clay_30) Cialca_isric_soilnpk$silt_5 <- as.numeric(Cialca_isric_soilnpk$silt_5) Cialca_isric_soilnpk$silt_15 <- as.numeric(Cialca_isric_soilnpk$silt_15) Cialca_isric_soilnpk$silt_30 <- as.numeric(Cialca_isric_soilnpk$silt_30) Cialca_isric_soilnpk$BD_5 <- as.numeric(Cialca_isric_soilnpk$BD_5) Cialca_isric_soilnpk$BD_15 <- as.numeric(Cialca_isric_soilnpk$BD_15) Cialca_isric_soilnpk$BD_30 <- as.numeric(Cialca_isric_soilnpk$BD_30) Cialca_isric_soilnpk$CEC_5 <- as.numeric(Cialca_isric_soilnpk$CEC_5) Cialca_isric_soilnpk$CEC_15 <- as.numeric(Cialca_isric_soilnpk$CEC_15) Cialca_isric_soilnpk$CEC_30 <- as.numeric(Cialca_isric_soilnpk$CEC_30) Cialca_isric_soilnpk$TotalN <- as.numeric(Cialca_isric_soilnpk$TotalN) Cialca_isric_soilnpk$Mn <- as.numeric(Cialca_isric_soilnpk$Mn) Cialca_isric_soilnpk$B <- as.numeric(Cialca_isric_soilnpk$B) Cialca_isric_soilnpk$Ca <- as.numeric(Cialca_isric_soilnpk$Ca) Cialca_isric_soilnpk$Fe <- as.numeric(Cialca_isric_soilnpk$Fe) Cialca_isric_soilnpk$Cu <- as.numeric(Cialca_isric_soilnpk$Cu) Cialca_isric_soilnpk$Al <- as.numeric(Cialca_isric_soilnpk$Al) Cialca_isric_soilnpk$Mg <- as.numeric(Cialca_isric_soilnpk$Mg) Cialca_isric_soilnpk$Na <- as.numeric(Cialca_isric_soilnpk$Na) Cialca_isric_soilnpk$ncluster <- as.factor(Cialca_isric_soilnpk$ncluster) Cialca_isric_soilnpk$Control <- as.numeric(Cialca_isric_soilnpk$Control) Cialca_isric_soilnpk$Country <- as.factor(Cialca_isric_soilnpk$Country ) ### create classes of control Cialca_isric_soilnpk$Control <- Cialca_isric_soilnpk$Control/1000 Cialca_isric_soilnpk$CONclass <- ifelse(Cialca_isric_soilnpk$Control < 7.5, "class1", ifelse(Cialca_isric_soilnpk$Control >= 7.5 & Cialca_isric_soilnpk$Control < 15, "class2", ifelse(Cialca_isric_soilnpk$Control >= 15 & Cialca_isric_soilnpk$Control < 22.5, "class3", ifelse(Cialca_isric_soilnpk$Control >= 22.5 & Cialca_isric_soilnpk$Control < 30, "class4", "class5")))) Cialca_isric_soilnpk$CONclass <- as.factor(Cialca_isric_soilnpk$CONclass) Cialca_isric_soilnpk$country <- Cialca_isric_soilnpk$Country Cialca_isric_soilnpk <- subset(Cialca_isric_soilnpk, select = -c(Country)) ### create classes of control Cialca_isda_soilnpk$Control <- Cialca_isda_soilnpk$Control/1000 Cialca_isda_soilnpk$CONclass <- ifelse(Cialca_isda_soilnpk$Control < 7.5, "class1", ifelse(Cialca_isda_soilnpk$Control >= 7.5 & Cialca_isda_soilnpk$Control < 15, "class2", ifelse(Cialca_isda_soilnpk$Control >= 15 & Cialca_isda_soilnpk$Control < 22.5, "class3", ifelse(Cialca_isda_soilnpk$Control >= 22.5 & Cialca_isda_soilnpk$Control < 30, "class4", "class5")))) Cialca_isda_soilnpk$CONclass <- as.factor(Cialca_isda_soilnpk$CONclass) Cialca_isda_soilnpk$Country <- as.factor(Cialca_isda_soilnpk$Country) Cialca_isda_soilnpk <- Cialca_isda_soilnpk[complete.cases(Cialca_isda_soilnpk), ] Cialca_isda_soilnpk$country <- Cialca_isda_soilnpk$Country Cialca_isda_soilnpk <- subset(Cialca_isda_soilnpk, select = -c(Country)) #### merging ACAI and CILACA data Cialca_isric_soilnpk$CON <- Cialca_isric_soilnpk$Control Cialca_isric_soilnpk <- subset(Cialca_isric_soilnpk, select = -c(Control)) ACAI_isric_GPS_allDAP_soilNPK <- ACAI_isric_GPS_allDAP_soilNPK[, colnames(Cialca_isric_soilnpk)] AC_CI_isric <- rbind(ACAI_isric_GPS_allDAP_soilNPK, Cialca_isric_soilnpk) unique(AC_CI_isric$country) saveRDS(AC_CI_isric, "AC_CIALCA_isric.RDS") Cialca_isda_soilnpk$CON <- Cialca_isda_soilnpk$Control Cialca_isda_soilnpk <- subset(Cialca_isda_soilnpk, select = -c(Control, bedrock_depth_0_200)) ACAI_isda_wide_allDAP_soilNPK <- ACAI_isda_wide_allDAP_soilNPK[, colnames(Cialca_isda_soilnpk)] AC_CI_isda <- rbind(ACAI_isda_wide_allDAP_soilNPK, Cialca_isda_soilnpk) unique(AC_CI_isda$country) saveRDS(AC_CI_isda, "AC_CIALCA_isda.RDS") ############################### ## not by country but just random split 70 30 AC_CI_isric_DRC <- droplevels(AC_CI_isric[!AC_CI_isric$country =="DRC", ]) AC_CI_isda_DRC <- droplevels(AC_CI_isda[!AC_CI_isda$country =="DRC", ]) traintestData_AC_CI_isric <- splitdata(soil_BLUPSData = AC_CI_isric_DRC) trainData_AC_CI_isric <- traintestData_AC_CI_isric[[1]] testData_AC_CI_isric <- traintestData_AC_CI_isric[[2]] traintestData_AC_CI_isda <- splitdata(soil_BLUPSData = AC_CI_isda_DRC) trainData_AC_CI_isda <- traintestData_AC_CI_isda[[1]] testData_AC_CI_isda <- traintestData_AC_CI_isda[[2]] trainData_AC_CI_isda$texture_class_0_20 <- factor(trainData_AC_CI_isda$texture_class_0_20, levels=levels(AC_CI_isda$texture_class_0_20)) trainData_AC_CI_isda$texture_class_20_50 <- factor(trainData_AC_CI_isda$texture_class_20_50 , levels=levels(AC_CI_isda$texture_class_20_50 )) trainData_AC_CI_isda$fcc <- factor(trainData_AC_CI_isda$fcc, levels=levels(AC_CI_isda$fcc)) trainData_AC_CI_isda$land_cover_2015 <- factor(trainData_AC_CI_isda$land_cover_2015, levels=levels(AC_CI_isda$land_cover_2015)) trainData_AC_CI_isda$land_cover_2016 <- factor(trainData_AC_CI_isda$land_cover_2016, levels=levels(AC_CI_isda$land_cover_2016)) trainData_AC_CI_isda$land_cover_2017 <- factor(trainData_AC_CI_isda$land_cover_2017, levels=levels(AC_CI_isda$land_cover_2017)) trainData_AC_CI_isda$land_cover_2018 <- factor(trainData_AC_CI_isda$land_cover_2018, levels=levels(AC_CI_isda$land_cover_2018)) trainData_AC_CI_isda$land_cover_2019 <- factor(trainData_AC_CI_isda$land_cover_2019, levels=levels(AC_CI_isda$land_cover_2019)) trainData_AC_CI_isda$country <- factor(trainData_AC_CI_isda$country, levels=levels(AC_CI_isda$country)) testData_AC_CI_isda$texture_class_0_20 <- factor(testData_AC_CI_isda$texture_class_0_20, levels=levels(AC_CI_isda$texture_class_0_20)) testData_AC_CI_isda$texture_class_20_50 <- factor(testData_AC_CI_isda$texture_class_20_50 , levels=levels(AC_CI_isda$texture_class_20_50 )) testData_AC_CI_isda$fcc <- factor(testData_AC_CI_isda$fcc, levels=levels(AC_CI_isda$fcc)) testData_AC_CI_isda$land_cover_2015 <- factor(testData_AC_CI_isda$land_cover_2015, levels=levels(AC_CI_isda$land_cover_2015)) testData_AC_CI_isda$land_cover_2016 <- factor(testData_AC_CI_isda$land_cover_2016, levels=levels(AC_CI_isda$land_cover_2016)) testData_AC_CI_isda$land_cover_2017 <- factor(testData_AC_CI_isda$land_cover_2017, levels=levels(AC_CI_isda$land_cover_2017)) testData_AC_CI_isda$land_cover_2018 <- factor(testData_AC_CI_isda$land_cover_2018, levels=levels(AC_CI_isda$land_cover_2018)) testData_AC_CI_isda$land_cover_2019 <- factor(testData_AC_CI_isda$land_cover_2019, levels=levels(AC_CI_isda$land_cover_2019)) testData_AC_CI_isda$country <- factor(testData_AC_CI_isda$country, levels=levels(AC_CI_isda$country)) Ndata_Train_ACCI_isric <- subset(trainData_AC_CI_isric, select=-c(soilP, soilK)) Pdata_Train_ACCI_isric <- subset(trainData_AC_CI_isric, select=-c(soilN, soilK)) Kdata_Train_ACCI_isric <- subset(trainData_AC_CI_isric, select=-c(soilN, soilP)) Ndata_Valid_ACCI_isric <- subset(testData_AC_CI_isric, select=-c(soilP, soilK)) Pdata_Valid_ACCI_isric <- subset(testData_AC_CI_isric, select=-c(soilN, soilK)) Kdata_Valid_ACCI_isric <- subset(testData_AC_CI_isric, select=-c(soilN, soilP)) Ndata_Train_ACCI_isda <- subset(trainData_AC_CI_isda, select=-c(soilP, soilK)) Pdata_Train_ACCI_isda <- subset(trainData_AC_CI_isda, select=-c(soilN, soilK)) Kdata_Train_ACCI_isda <- subset(trainData_AC_CI_isda, select=-c(soilN, soilP)) Ndata_Valid_ACCI_isda <- subset(testData_AC_CI_isda, select=-c(soilP, soilK)) Pdata_Valid_ACCI_isda <- subset(testData_AC_CI_isda, select=-c(soilN, soilK)) Kdata_Valid_ACCI_isda <- subset(testData_AC_CI_isda, select=-c(soilN, soilP)) #######################################3 # split by country using ACAI to train the model and CI to validate trainData_isric <- droplevels(AC_CI_isric[AC_CI_isric$country %in% c("NG", "TZ"),]) testData_isric <- droplevels(AC_CI_isric[AC_CI_isric$country %in% c("BI", "Rw"),]) testData_isric$ncluster <- factor(testData_isric$ncluster, levels = levels(AC_CI_isric$ncluster)) trainData_isric$country <- factor(trainData_isric$country, levels = levels(AC_CI_isric$country)) testData_isric$country <- factor(testData_isric$country, levels = levels(AC_CI_isric$country)) trainData_isda <- droplevels(AC_CI_isda[AC_CI_isda$country %in% c("NG", "TZ"),]) testData_isda <- droplevels(AC_CI_isda[AC_CI_isda$country %in% c("BI", "Rw"),]) trainData_isda$texture_class_0_20 <- factor(trainData_isda$texture_class_0_20, levels=levels(AC_CI_isda$texture_class_0_20)) trainData_isda$texture_class_20_50 <- factor(trainData_isda$texture_class_20_50 , levels=levels(AC_CI_isda$texture_class_20_50 )) trainData_isda$fcc <- factor(trainData_isda$fcc, levels=levels(AC_CI_isda$fcc)) trainData_isda$land_cover_2015 <- factor(trainData_isda$land_cover_2015, levels=levels(AC_CI_isda$land_cover_2015)) trainData_isda$land_cover_2016 <- factor(trainData_isda$land_cover_2016, levels=levels(AC_CI_isda$land_cover_2016)) trainData_isda$land_cover_2017 <- factor(trainData_isda$land_cover_2017, levels=levels(AC_CI_isda$land_cover_2017)) trainData_isda$land_cover_2018 <- factor(trainData_isda$land_cover_2018, levels=levels(AC_CI_isda$land_cover_2018)) trainData_isda$land_cover_2019 <- factor(trainData_isda$land_cover_2019, levels=levels(AC_CI_isda$land_cover_2019)) trainData_isda$country <- factor(trainData_isda$country, levels=levels(AC_CI_isda$country)) testData_isda$texture_class_0_20 <- factor(testData_isda$texture_class_0_20, levels=levels(AC_CI_isda$texture_class_0_20)) testData_isda$texture_class_20_50 <- factor(testData_isda$texture_class_20_50 , levels=levels(AC_CI_isda$texture_class_20_50 )) testData_isda$fcc <- factor(testData_isda$fcc, levels=levels(AC_CI_isda$fcc)) testData_isda$land_cover_2015 <- factor(testData_isda$land_cover_2015, levels=levels(AC_CI_isda$land_cover_2015)) testData_isda$land_cover_2016 <- factor(testData_isda$land_cover_2016, levels=levels(AC_CI_isda$land_cover_2016)) testData_isda$land_cover_2017 <- factor(testData_isda$land_cover_2017, levels=levels(AC_CI_isda$land_cover_2017)) testData_isda$land_cover_2018 <- factor(testData_isda$land_cover_2018, levels=levels(AC_CI_isda$land_cover_2018)) testData_isda$land_cover_2019 <- factor(testData_isda$land_cover_2019, levels=levels(AC_CI_isda$land_cover_2019)) testData_isda$country <- factor(testData_isda$country, levels=levels(AC_CI_isda$country)) Ndata_Train_isric <- subset(trainData_isric, select=-c(soilP, soilK)) Pdata_Train_isric <- subset(trainData_isric, select=-c(soilN, soilK)) Kdata_Train_isric <- subset(trainData_isric, select=-c(soilN, soilP)) Ndata_Valid_isric <- subset(testData_isric, select=-c(soilP, soilK)) Pdata_Valid_isric <- subset(testData_isric, select=-c(soilN, soilK)) Kdata_Valid_isric <- subset(testData_isric, select=-c(soilN, soilP)) Ndata_Train_isda <- subset(trainData_isda, select=-c(soilP, soilK)) Pdata_Train_isda <- subset(trainData_isda, select=-c(soilN, soilK)) Kdata_Train_isda <- subset(trainData_isda, select=-c(soilN, soilP)) Ndata_Valid_isda <- subset(testData_isda, select=-c(soilP, soilK)) Pdata_Valid_isda <- subset(testData_isda, select=-c(soilN, soilK)) Kdata_Valid_isda <- subset(testData_isda, select=-c(soilN, soilP)) custom <- trainControl(method="oob", number=10) ########################################################################## ## Random Forest soilN: with isric = 0.52; isda = 0.57 ########################################################################## set.seed(773) Ndata_Train <- subset(Ndata_Train_isric, select=-c(CON))## country based isric 0.57 Ndata_Valid <- Ndata_Valid_isric summary(Ndata_Train_isric$soilN) summary(Ndata_Valid$soilN) Ndata_Train <- subset(Ndata_Train_isda, select=-c(CON))## country based isda 0.61 Ndata_Valid <- Ndata_Valid_isda summary(Ndata_Train$soilN) summary(Ndata_Valid$soilN) Ndata_Train <- subset(Ndata_Train_ACCI_isric, select=-c(CON)) ## not country based isric rsq = 0.71 Ndata_Valid <- Ndata_Valid_ACCI_isric Ndata_Train <- subset(Ndata_Train_ACCI_isda, select=-c(CON)) ## not country based isda rsq = 0.78 Ndata_Valid <- Ndata_Valid_ACCI_isda rm(RF_N1) RF_N1 <- randomForest(soilN ~ ., Ndata_Train , importance=TRUE, ntree=1000) Ndata_Valid$predN <- predict(RF_N1, Ndata_Valid) sst <- sum((Ndata_Valid$soilN - mean(Ndata_Valid$soilN))^2) sse <- sum((Ndata_Valid$soilN - Ndata_Valid$predN)^2) rsq <- 1 - sse / sst rsq varImpPlot(RF_N1) varimportance <- data.frame(RF_N1$importance) varimportance$vars <- rownames(varimportance) varimportance <- varimportance[order(varimportance$X.IncMSE, decreasing = TRUE), ] rownames(varimportance) <- NULL varimportance$Importance <- c(1:nrow(varimportance)) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilN_ACCI_isric_together.csv", row.names = FALSE) ggnn <- ggplot(Ndata_Valid, aes(soilN, predN, col=country)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + xlim(0,250) + ylim(0,250)+ xlab("soil N from QUEFTS optimization, ACAI")+ ylab("perdicted soil N from GIS data, CIALCA")+ ggtitle("Soil Grids + control class, AC and CI data together")+ theme_bw() ggsave("QUEFTS_GIS_RF_soilN_isric_AC_CI.pdf", ggnn, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilN_isda_AC_CI.pdf", ggnn, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilN_isric_ACCI_together.pdf", ggnn, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilN_isda_ACCI_together.pdf", ggnn, width=5, height=5) ########################################################################## ## Random Forest soilP: with isric = 0.51; isda=0.54 ########################################################################## set.seed(773) Pdata_Train <- subset(Pdata_Train_isric, select=-c(CON))##0.69 Pdata_Valid <- Pdata_Valid_isric summary(Pdata_Train_isric$soilP) summary(Pdata_Valid_isric$soilP) Pdata_Train <- subset(Pdata_Train_isda, select=-c(CON))##0.73 Pdata_Valid <- Pdata_Valid_isda summary(Pdata_Train_isric$soilP) summary(Pdata_Valid_isric$soilP) Pdata_Train <- subset(Pdata_Train_ACCI_isric, select=-c(CON)) ## rsq = 0.70 Pdata_Valid <- Pdata_Valid_ACCI_isric Pdata_Train <- subset(Pdata_Train_ACCI_isda, select=-c(CON)) ## rsq = 0.74 Pdata_Valid <- Pdata_Valid_ACCI_isda rm(RF_P1) RF_P1 <- randomForest(soilP ~ ., Pdata_Train , importance=TRUE, ntree=1000) Pdata_Valid$predP <- predict(RF_P1, Pdata_Valid) sst <- sum((Pdata_Valid$soilP - mean(Pdata_Valid$soilP))^2) sse <- sum((Pdata_Valid$soilP - Pdata_Valid$predP)^2) rsq <- 1 - sse / sst rsq varImpPlot(RF_P1) varimportance <- data.frame(RF_P1$importance) varimportance$vars <- rownames(varimportance) varimportance <- varimportance[order(varimportance$X.IncMSE, decreasing = TRUE), ] rownames(varimportance) <- NULL varimportance$Importance <- c(1:nrow(varimportance)) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilP_ACCI_isric_together.csv", row.names = FALSE) ggpp <- ggplot(Pdata_Valid, aes(soilP, predP, col=country)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + xlim(0,150) + ylim(0,150) + ggtitle("Soil Grids + control class, AC and CI data together")+ xlab("soil P from QUEFTS optimization, ACAI") + ylab("perdicted soil P from GIS data, CIALCA") + theme_bw() ggsave("QUEFTS_GIS_RF_soilP_isric_AC_CI.pdf", ggpp, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilP_isda_AC_CI.pdf", ggpp, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilP_isda_ACCI_together.pdf", ggpp, width=5, height=5) ########################################################################## ## Random Forest soilK: with isric = 0.44, isda=0.16 ########################################################################## Kdata_Train <- subset(Kdata_Train_isric, select=-c(CON))##.48 Kdata_Valid <- Kdata_Valid_isric summary(Kdata_Train$soilK) summary(Kdata_Valid$soilK) Kdata_Train <- subset(Kdata_Train_isda, select=-c(CON)) ##0.12 Kdata_Valid <- Kdata_Valid_isda summary(Kdata_Train$soilK) summary(Kdata_Valid$soilK) Kdata_Train <- subset(Kdata_Train_ACCI_isric, select=-c(CON))##0.44 Kdata_Valid <- Kdata_Valid_ACCI_isric Kdata_Train <- subset(Kdata_Train_ACCI_isda , select=-c(CON))##0.55 Kdata_Valid <- Kdata_Valid_ACCI_isda rm(RF_K1) set.seed(773) RF_K1 <- randomForest(soilK ~ ., Kdata_Train , importance=TRUE, ntree=1000) Kdata_Valid$predK <- predict(RF_K1, Kdata_Valid) sst <- sum((Kdata_Valid$soilK - mean(Kdata_Valid$soilK))^2) sse <- sum((Kdata_Valid$soilK - Kdata_Valid$predK)^2) rsq <- 1 - sse / sst rsq varImpPlot(RF_K1) varimportance <- data.frame(RF_K1$importance) varimportance$vars <- rownames(varimportance) varimportance <- varimportance[order(varimportance$X.IncMSE, decreasing = TRUE), ] rownames(varimportance) <- NULL varimportance$Importance <- c(1:nrow(varimportance)) write.csv(varimportance[, c('vars', 'Importance')], "AC_RFvars_soilK_ACCI_isric_together.csv", row.names = FALSE) ggkk <- ggplot(Kdata_Valid, aes(soilK, predK, col=country)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + xlim(0,225) + ylim(0,225)+ xlab("soil K from QUEFTS optimization, ACAI")+ ggtitle("Soil Grids + control class, AC and CI data together")+ ylab("perdicted soil K from GIS data, CIALCA")+ theme_bw() ggsave("QUEFTS_GIS_RF_soilK_isric_AC_CI.pdf", ggkk, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilK_isda_AC_CI.pdf", ggkk, width=5, height=5) ggsave("QUEFTS_GIS_RF_soilK_isda_ACCI_together.pdf", ggkk, width=5, height=5) ggplot(Kdata_Valid, aes(soilK, predK)) + geom_point()+ geom_abline(slope = 1, intercept = 0) + geom_smooth(method=lm) + xlim(0,225) + ylim(0,225)+ xlab("soil K from QUEFTS optimization, ACAI")+ ylab("perdicted soil K from GIS data, CIALCA")+ theme_bw()