Bayes Emulation of simulated Crop Yield

Hasan, M.M., Cumming, J.A. (2021). Bayes Linear Emulation of Simulated Crop Yield. In: Chaubey, Y.P., Lahmiri, S., Nebebe, F., Sen, A. (eds) Applied Statistics and Data Science. CCAS 2021. Springer Proceedings in Mathematics & Statistics, vol 375. Springer, Cham. https://doi.org/10.1007/978-3-030-86133-9_7

df=read.csv("Data/MAIZ1-yld.csv")

attach(df)

x1=as.factor(df$St) #Steepness

x2=as.factor(df$So) #Soil

x3=as.factor(df$WYear) #Weather

x4=as.factor(df$SYear) #Simulation Year

Data=data.frame(x1,x2,x3,x4)

a=as.factor(paste(Data$x1,Data$x2,Data$x3,Data$x4,sep='-'))

Data = cbind(Data,Unique=a,GYLD,N,P)

my.levels = levels(a)

#one unique level

d4=subset(Data,Data$Unique==my.levels[1])

dfmb$<-$d4[d4$N$>$0 & d4$P$>$0 ,]

# Training Data

s=c(t1 =.2,t2=.2,t3=.2,t4=.2,test=.2)

g = sample(cut(seq(nrow(dfmb)), nrow(dfmb)*cumsum(c(0,s)), labels = names(s)))

res = split(dfmb, g)

df=rbind(res$t1,res$t2,res$t3,res$t4)

y1=df$GYLD #yield

xn=df$N #Nitrogen

zp=df$P #Phosphorus

D=y1

pt=lm(y1~xn+zp) # Linear Modelling

ED=pt$coefficients[1]+pt$coefficients[2]*xn +pt$coefficients[3]*zp

#Testing Data

ytest=res$test[,6]

xt=res$test[,7]

zt=res$test[,8]

p=expand.grid(xt,zt)

v=matrix(vcov(pt),nrow =3,ncol=3)

p1=data.frame(1,xn,zp)

px=as.matrix(p1,nrow=length(xn),ncol=3)

VSigma=px%*%v%*%t(px) # Equation 3 basis of variance (First part)

p=expand.grid(xt,zt)

k1=data.frame(1,p[,1],p[,2])

k=as.matrix(k1,nrow=length(p[,1]),ncol=3)

#variance of general basis of Equ3

VS=k%*%v%*%t(k) # Equation 3 basis of variance (First part)

VSn=k%*%v%*%t(px) # Equation 3 basis of variance (First part)

EB=pt$coefficients[1]+pt$coefficients[2]*p[,1] +pt$coefficients[3]*p[,2]

rs=sum((y1-ED)^2)/(length(xn)-2)#Error SS

Snew=VSn+rs*exp(-.015*outer(p[,1],xn,'-')^2) *exp(-0.015*outer(p[,2], zp,'-')^2) #Equation 4

Sstar=VS+rs*exp(-.015*outer(p[,1],p[,1],'-')^2) *exp(-0.015 *outer(p[,2],p[,2],'-')^2) #Equation 4

Sigma=VSigma+rs*exp(-.015* outer(xn,xn,'-')^2) *exp(-0.015*outer(zp,zp,'-')^2) # Equation 4

# Equation 1

ED_B=EB+Snew%*%solve(Sigma)%*%(D - ED)

#Equation 2

VarD_B=Sstar-Snew%*%solve(Sigma)%*%t(Snew)

# Diagnostics

library(plotly)

Emean=c(ED_B) # Expectation

df1 <-data.frame(x=x1,y=x2,z=Emean)

p <-plot_ly(data=df1,x=~x1,y=~x2, z=~Emean,

type = "contour", colorscale='Jet') #Fig. 2

sd=c(sqrt(diag(VarD_B))) # Variance

df2 <-data.frame(x=x1,y=x2,z=sd)

p <-plot_ly(data=df2,x=~x1,y=~x2, z=~sd,

type ="contour",colorscale='Jet') #Fig. 2

RD_B=1-(diag(VarD_B)/diag(Sigmastar)) #Equation 5

x1=p[,1]; x2=p[,2]

df3 <-data.frame(x=x1,y=x2,z=RD_B)

p<-plot_ly(data=df3, x=~x1,y=~x2, z=~RD_B,

type = "contour", colorscale='Jet')#Fig3

ui=(ytest-ED_B)/sqrt(diag(VarD_B))#Equ 6

plot(xt,c(ui),xlab="N points",

ylab="Standardized Prediction Errors") #Fig3

abline(h=3,col="blue")

abline(h=-3,col="blue")

plot(zt,c(ui),xlab="P points",

ylab="Standardized Prediction Errors") #Fig3

abline(h=3,col="blue")

abline(h=-3,col="blue")