USE OF HIGH-RESOLUTION IMAGE DATA OUTPERFORMS VEGETATION INDICES IN PREDICTION OF MAIZE YIELD
Supplementary methods
Number of LVs for PLS model
# Load pls library
library(pls)
# Loading data from heat and Drought Stress trials only
load('objectsDS.rdata')
# Create an object for scaled grain yieldy= scale(Y[,1])
# Initialize the list that will contain the number of LVsnLV.List=list()
#### The next will run for each of the single and multi-time pointsfor (win1:9){
if(w<6) {X= eval(parse(text=paste0("X",w)))}
if(w==6) {X= cbind(X4,X5)}
if(w==7) {X= cbind(X3,X4,X5)}
if(w==8) {X= cbind(X2,X3,X4,X5)}
if(w==9) {X= cbind(X1,X2,X3,X4,X5)}
nFolds= length(unique(Y[,2]))
nLV= seq(from=2,to=15,by=1)
folds= as.numeric(Y[,2])
COR.TST= rep(NA,nFolds)
nLV.BEST= rep(NA,nFolds)
for(iin unique(folds)){
trn= which(folds!=i)
COR=matrix(nrow=11,ncol=length(nLV))
for(jin1:11){
trn_trn= which(folds[trn]!=j)
for(kin1:length(nLV)){
COR[j,k] = cor(y[-trn_trn],predict(plsr(formula=y[trn_trn]~X[trn_trn,],ncomp=nLV[k]),ncomp=nLV[k],newdata=X[-trn_trn,]))
}
print(paste("subFold=",j,"Fold=",i,"Timepoint=",w))
}
COR= colMeans(COR)
nLV.BEST[i] =nLV[which(COR==min(COR))]
COR.TST[i] = cor(y[-trn],predict(plsr(formula=y[trn]~X[trn,],ncomp=nLV.BEST[i]),ncomp=nLV.BEST[i],newdata=X[-trn,]))
}
nLV.List[[w]] =nLV.BEST
}
The following script performs leave-one-trial out cross validation to obtain predictions for each model with data from heat and drought trials. Predictions of the irrigated trial were obtained separately with an adaptated version of this code.
# Load libraries
library(BGLR)
# Create an object for identification of trialstrials= as.integer(factor(Y[,2]))
# Obtain predictions from single time point datafor (kin1:5){
# initialize the object that will contain predictionsYHat=data.frame(matrix(nrow=length(y),ncol=9))
colnames(YHat)=c('NDVI','CWMI','mND','PRI','MTCI','MCARI2','OLS','PLS','BayesB')
X= eval(parse(text=paste0('X',k)))
X=X/sqrt(ncol(X))
# Cross Validation leave one trial outfor (iin unique(trials)){
print(i)
trn= which(trials!=i); tst= which(trials==i)
yTRN=y[trn] ; yTST=y[tst]
trialTRN=trials[trn] ; trialTST=trials[tst]
X.TRN=X[trn,] ; X.TST=X[tst,]
NDVI.TRN=NDVI[trn,k] ; NDVI.TST=NDVI[tst,k]
CWMI.TRN=CWMI[trn,k] ; CWMI.TST=CWMI[tst,k]
mND.TRN=mND[trn,k] ; mND.TST=mND[tst,k]
PRI.TRN=PRI[trn,k] ; PRI.TST=PRI[tst,k]
MTCI.TRN=MTCI[trn,k] ; MTCI.TST=MTCI[tst,k]
MCARI2.TRN=MCARI2[trn,k] ; MCARI2.TST=MCARI2[tst,k]
# NDVIModel1= lm(yTRN~NDVI.TRN)
pY1= as.matrix(NDVI.TST)%*%(coef(Model1)[-1]) + coef(Model1)[1]
YHat[tst,'NDVI'] =pY1# CWMIModel2= lm(yTRN~CWMI.TRN)
pY2= as.matrix(CWMI.TST)%*%(coef(Model2)[-1]) + coef(Model2)[1]
YHat[tst,'CWMI'] =pY2# mNDModel3= lm(yTRN~mND.TRN)
pY3= as.matrix(mND.TST)%*%(coef(Model3)[-1]) + coef(Model3)[1]
YHat[tst,'mND'] =pY3# PRIModel4= lm(yTRN~PRI.TRN)
pY4= as.matrix(PRI.TST)%*%(coef(Model4)[-1]) + coef(Model4)[1]
YHat[tst,'PRI'] =pY4# MTCIModel5= lm(yTRN~MTCI.TRN)
pY5= as.matrix(MTCI.TST)%*%(coef(Model5)[-1]) + coef(Model5)[1]
YHat[tst,'MTCI'] =pY5# MCARI2Model6= lm(yTRN~MCARI2.TRN)
pY6= as.matrix(MCARI2.TST)%*%(coef(Model6)[-1]) + coef(Model6)[1]
YHat[tst,'MCARI2'] =pY6# OLSModel7= lm(yTRN~X.TRN)
pY7= as.matrix(X.TST)%*%(coef(Model7)[-1]) + coef(Model7)[1]
YHat[tst,'OLS'] =pY7# PLSModel8= plsr(formula=yTRN~X.TRN, ncomp=nLV.List[[k]][i])
pY8= predict(Model8,ncomp=nLV.List[[k]][i],newdata=X.TST)
YHat[tst,'PLS'] =pY8# BayesBModel9= BGLR(y=yTRN,ETA=list(list(X=X.TRN,model='BayesB')),nIter=120000,burnIn=5000, verbose=F,rmExistingFiles=T,saveAt="~/tmp_uavs",groups=trials[trn])
pY9= as.matrix(X.TST)%*%Model9$ETA[[1]]$b+Model9$muYHat[tst,'BayesB'] =pY9
}
# Save the predictions in files separated by time point
save(y,YHat,trials,file=paste0('YHATDSCV',timePoint,'.RData'))
}