#============================================================================================================
#
#    Licensing:
#
#    Copyright Maria Grünewald and Ola Hössjer (2010)
#
#    This code is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version. See <http://www.gnu.org/licenses/>.
#
#    This code is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#
#    Citation:
#
#    The methodology used in this library is described in Grünewald & Hössjer (2011) A General Statistical Framework 
#	  for Multistage Designs, ...
#    Please refer to this paper when using the library.
#============================================================================================================
#
#    Contents:
#
#    This library computes and plots efficiencies for sequential design in Three-stage sampling design with two Gaussian outcomes
#    For description of the model see Example 5 in
#    Grünewald & Hössjer (2011) A General Statistical Framework 
#	  for Multistage Designs, ...

#
#
#   Output:
#
#   File    "EfficiencyTreestage.txt" saved to R working directory
#           containing settings and numerical results from simulation run
#
#   Figures shown in R Graphics devise
#           *Plots illustrating cost adjusted efficiency
#
#
#============================================================================================================
#
#    The code is intended for use in the software R (version 2.10.0):
#    R Development Core Team (2008). R: A language and environment for
#    statistical computing. R Foundation for Statistical Computing,
#    Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.
#
#============================================================================================================

#============================================================================================================
# Setting parameters for your calculations
#============================================================================================================

# Model
# x2->x1, x1,x2->y  x1,y->A
# X2~Bin, X1~Bin, Y~N

#------------------------------------------------------------------------------------------------------------
# Parameters specified by user
#------------------------------------------------------------------------------------------------------------


# Number of simulated individuals. This number does not refer to the number of subjects in the actual study,
# but only affects how well the Monte Calo method approximates theoretical results.

N<-100

#
B0x2<-0
#px2= exp(B0x2)/(1+exp(B0x2))
x2groups<-2
#Here programmed only to work with x1groups=2. If more groups, the ascertainment
#scheme has to be updated.
x1groups<-2
x1s<-c(0,1)


B0x1<--3
Bx2x1<-2
B0y<-0
Bx2y<-1
Bx1y<-1
SDy<-1



















# for sampling scheme with only two regions per variable "cuty" and "cutx1" determines these.
# See function "weights" for further specification

cuty<-3
cutx1<-0.5




# Specifying sampling probabilities
# (Vector of length >=1)

#lambda_2 for y<cuty
PAyl<-c(0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
#lambda_2 for y>=cuty
PAyh<-c(rep(1,(length(PAyl)-1)),1)

#lambda_3 for y<cuty,x1<cutx1)
pAylx1l<-c(0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
#lambda_3 for y>=cuty,x1<cutx1
pAyhx1l<-c(rep(1,(length(pAylx1l)-1)),1)
#lambda_3 for y>=cuty,x1>=cutx1
pAyhx1h<-c(rep(1,(length(pAylx1l)-1)),1)
#lambda_3 for y<cuty,x1>=cutx1
pAylx1h<-c(rep(1,(length(pAylx1l)-1)),1)


#------------------------------------------------------------------------------------------------------------
# The following parameters follow from settings above
#------------------------------------------------------------------------------------------------------------


#number of parameters in the model
parnum<-7
# number of parameters in explanatory variable
x2parnum<-1

parnumfirsty<-x2parnum+1


#============================================================================================================
# Functions to simulate data and create weights for Monte carlo calculations
#============================================================================================================



# Simulating data


simdata<-function(N,B0x2,x2groups,B0x1,Bx2x1,x1groups,B0y,Bx2y,Bx1y,SDy){

sim<-matrix(NA,N,3)

#x2
x2<-rbinom(N,x2groups-1,exp(B0x2)/(1+exp(B0x2)))

#x1
x1<-rbinom(N,x1groups-1,exp(B0x1+Bx2x1*x2)/(1+exp(B0x1+Bx2x1*x2)))

#y
y<-rnorm(N,B0y+Bx2y*x2+Bx1y*x1,SDy)

sim[,1]<-x2
sim[,2]<-x1
sim[,3]<-y

return(sim)
}



# The weights below corresponds to the sampling scheme
# used in Example 5
# Change the weight function to alter the sampling scheme



weightsX2<-function(x1,y,cutx1,cuty,pAyhx1h,pAyhx1l,pAylx1h,pAylx1l){
ascprob<-matrix(NA,length(x1),1)
for (i in 1:length(x1)){
    if (y[i]> cuty) {
        if (x1[i]> cutx1) {
            ascprob[i,]<-pAyhx1h}
        if (x1[i]< cutx1) {
            ascprob[i,]<-pAyhx1l}
    }
    if (y[i] < cuty){
        if (x1[i]> cutx1) {
            ascprob[i,]<-pAylx1h}
        if (x1[i]< cutx1) {
            ascprob[i,]<-pAylx1l}
    }
}
return(ascprob)
}


weightsX1<-function(y,cuty,PAyl,PAyh){
ascprob<-matrix(NA,length(y),1)
for (i in 1:length(y)){
    if (y[i]> cuty) {
        ascprob[i,]<-PAyh

        }
    if (y[i] < cuty){
            ascprob[i,]<-PAyl
        }
    }
return(ascprob)
}




#============================================================================================================
# Functions to calculate efficiency
#============================================================================================================



#------------------------------------------------------------------------------------------------------------
# Function to calculate \psi
#------------------------------------------------------------------------------------------------------------

psi<-function(N,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,x2,x1,y){

    dB0x2<-(exp(B0x2)*(x2-(x2groups-1))+x2)/(1+exp(B0x2))
    dB0x1<-(exp(B0x1+Bx2x1*x2)*(x1-(x1groups-1))+x1)/(1+exp(B0x1+Bx2x1*x2))
    dBx2x1<-x2*(exp(B0x1+Bx2x1*x2)*(x1-(x1groups-1))+x1)/(1+exp(B0x1+Bx2x1*x2))
    dB0y<--(B0y+Bx2y*x2+Bx1y*x1-y)/SDy^2
    dBx2y<-x2*(-(B0y+Bx2y*x2+Bx1y*x1-y))/SDy^2
    dBx1y<-x1*(-(B0y+Bx2y*x2+Bx1y*x1-y))/SDy^2
    dSDy<--(SDy^2-(B0y+Bx2y*x2+Bx1y*x1-y)^2)/SDy^3


    dpar<-matrix(c(dB0x2,dB0x1,dBx2x1,dB0y,dBx2y,dBx1y,dSDy),N,parnum)

    return(dpar)

}


#------------------------------------------------------------------------------------------------------------
# Function to calculate P(z)
#------------------------------------------------------------------------------------------------------------


Pz<-function(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,x2,x1,y){
    dnorm(y, mean=(B0y+Bx2y*x2+Bx1y*x1), sd=SDy)*dbinom(x1, size=(x1groups-1), prob=exp(B0x1+Bx2x1*x2)/(1+exp(B0x1+Bx2x1*x2)))*dbinom(x2, size=(x2groups-1), prob=exp(B0x2)/(1+exp(B0x2)))
}

#------------------------------------------------------------------------------------------------------------
# Function to calculate P(y)
#------------------------------------------------------------------------------------------------------------


Py<-function(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,y){
    py<-0
    for(j in 1:x2groups){
        for(k in 1:x1groups){
        py<-Pz(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(j-1),(k-1),y)+py
        }
    }
    return(py)
}



#------------------------------------------------------------------------------------------------------------
# Function to calculate P(yx1)
#------------------------------------------------------------------------------------------------------------



Pyx1<-function(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,x1,y){
    pyx1<-0
    for(j in 1:x2groups){
        pyx1<-Pz (B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(j-1),x1,y)+pyx1
    }
    return(pyx1)
}

#------------------------------------------------------------------------------------------------------------
# Function to calculate P(x2x1)
#------------------------------------------------------------------------------------------------------------



Px2x1<-function(B0x2,B0x1,Bx2x1,x2groups,x1groups,x2,x1){
    px2x1<-dbinom(x1, size=(x1groups-1), prob=exp(B0x1+Bx2x1*x2)/(1+exp(B0x1+Bx2x1*x2)))*dbinom(x2, size=(x2groups-1), prob=exp(B0x2)/(1+exp(B0x2)))
    return(px2x1)
}


#------------------------------------------------------------------------------------------------------------
# Function to calculate P(x2)
#------------------------------------------------------------------------------------------------------------


Px2<-function(B0x2,x2groups,x2){
    dbinom(x2, size=(x2groups-1), prob=exp(B0x2)/(1+exp(B0x2)))
}

#------------------------------------------------------------------------------------------------------------
# Function to calculate P(x1)
#------------------------------------------------------------------------------------------------------------


Px1<-function(B0x2,B0x1,Bx2x1,x2groups,x1groups,x1){
    px1<-0
    for(j in 1:x2groups){
        px1<-dbinom(x1, size=(x1groups-1), prob=exp(B0x1+Bx2x1*(j-1))/(1+exp(B0x1+Bx2x1*(j-1))))*dbinom((j-1), size=(x2groups-1), prob=exp(B0x2)/(1+exp(B0x2)))+px1
    }
    return(px1)
}



#------------------------------------------------------------------------------------------------------------
# Function to calculate variances
#------------------------------------------------------------------------------------------------------------



variances<-function(N,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,x2,x1,y,w1,w2,w12,w12s,NstarX1,NstarX12,parnum){




    #-----------------------------------------------------------------------------------------------------
    # Calculating the fist term of the information as E[(f'/f)^2] in ascertainment sample
    #-----------------------------------------------------------------------------------------------------

    infopart1<-matrix(NA,parnum,parnum)
    for (i in 1:parnum){
        for (j in 1:parnum){
            infopart1[i,j]<-sum(psidata[,i]*psidata[,j]*w12)/N
        }
    }

    #-----------------------------------------------------------------------------------------------------
    # Calculating the second term of the information in ascertainment sample
    #-----------------------------------------------------------------------------------------------------



    #-----------------------------------------------------------------------------------------------------
    #Joint likelihood
    #-----------------------------------------------------------------------------------------------------

    infopart2joint<-matrix(NA,parnum,parnum)
    for (i in 1:parnum){
        for (j in 1:parnum){
            infopart2joint[i,j]<-sum(psidata[,i]*w12)/N*sum(psidata[,j]*w12)/N
        }
    }

    #-----------------------------------------------------------------------------------------------------
    #Prospective likelihood
    #-----------------------------------------------------------------------------------------------------

    infopart2prosp<-matrix(0,parnum,parnum)
    for (j in 1:x2groups){
        psidataj<-psidata[x2==j-1,]
        wj<-w12[x2==j-1]
        Nstarj<-sum(wj)
        for (i in 1:parnum){
            for (k in 1:parnum){
                infopart2prosp[i,k]<-(sum(psidataj[,i]*wj)/N*sum(psidataj[,k]*wj)/N)+infopart2prosp[i,k]
            }
        }
    }



    #-----------------------------------------------------------------------------------------------------
    # Retrospective likelihood X1
    #-----------------------------------------------------------------------------------------------------


    sumIyretroX1<-matrix(0,parnum,parnum)
    #------------------------------
    #for each sampled value of y
    #------------------------------
    for (j in 1:length(y)){
        Epsiy<-matrix(0,1,parnum)

        Epsiypsiy<-matrix(0,parnum,parnum)
        Epsiyx1etc<-matrix(0,parnum,parnum)

                    #------------------------------
                    #for each possible value of x2 and x1 (again) borde kunna flyttas ut ett steg
                    #------------------------------
                    for (q1 in 1:x1groups){
                        for (q2 in 1:x2groups){
                            psidatajs<-psi(1,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(q2-1),(q1-1),y[j])
                            Epsiy<-Pz(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(q2-1),(q1-1),y[j])/Py(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,y[j])*psidatajs+Epsiy
                         }
                    }

        #------------------------------
        #for each possible value of x1
        #------------------------------
        for (l in 1:x1groups){

            Epsiyx1<-matrix(0,1,parnum)
            Epsiyx1etc<-Pz(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(k-1),(l-1),y[j])/Pyx1(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(l-1),y[j])*t(Epsiyx1-Epsiy)%*%(Epsiyx1-Epsiy)+Epsiyx1etc

       }

        sumIyretroX1<-w1[j,]*Epsiyx1etc+sumIyretroX1
    }

    IretroX1<-sumIyretroX1/N






   #-----------------------------------------------------------------------------------------------------
   #Retrospective likelihood X2
   #-----------------------------------------------------------------------------------------------------



    sumIyretroX2<-matrix(0,parnum,parnum)
    #------------------------------
    #for each sampled value of y
    #------------------------------
    for (j in 1:length(y)){
        #------------------------------
        #for each possible value of x1
        #------------------------------
        for (l in 1:x1groups){
            Epsiyx1etc<-matrix(0,parnum,parnum)
            Epsiyx1<-matrix(0,1,parnum)
                    #------------------------------
                    #for each possible value of x2 and x1 (again) borde kunna flyttas ut ett steg
                    #------------------------------
                    for (q1 in 1:x1groups){
                        for (q2 in 1:x2groups){
                            psidatajs<-psi(1,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(q2-1),(q1-1),y[j])
                            Epsiyx1<-Pz(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(q2-1),(q1-1),y[j])/Py(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,y[j])*psidatajs+Epsiyx1
                         }
                    }
                  #------------------------------
                  #for each possible value of x2
                  #------------------------------
                for (k in 1:x2groups){
                    psidatajs<-psi(1,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(k-1),(l-1),y[j])
                    Epsiyx1etc<-Pz(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(k-1),(l-1),y[j])/Pyx1(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(l-1),y[j])*t(psidatajs-Epsiyx1)%*%(psidatajs-Epsiyx1)+Epsiyx1etc

                }
            sumIyretroX2<-w12s[j,l]*Pyx1(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(l-1),y[j])/Py(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,y[j])*Epsiyx1etc+sumIyretroX2
        }
    }








    #------------------------------------------------------------------------------------------------------------
    # Calculating information matrix from the parts above
    #------------------------------------------------------------------------------------------------------------


    IretroX2<-sumIyretroX2/N
    infojoint<-infopart1-infopart2joint
    infomulti<-(Iy+IretroX1+IretroX2)



    #------------------------------------------------------------------------------------------------------------
    # Inverting the information matrix to get variance
    #------------------------------------------------------------------------------------------------------------

    variancematrixjoint<-solve(infojoint)
    variancematrixmulti<-solve(infomulti)


    #------------------------------------------------------------------------------------------------------------
    # Output of function:variancematrix
    #------------------------------------------------------------------------------------------------------------

    return(list(variancematrixjoint=variancematrixjoint,variancematrixmulti=variancematrixmulti))

}



#============================================================================================================
# Running simulations
#============================================================================================================


#------------------------------------------------------------------------------------------------------------
# Simulate data
#------------------------------------------------------------------------------------------------------------


data<-simdata(N,B0x2,x2groups,B0x1,Bx2x1,x1groups,B0y,Bx2y,Bx1y,SDy)
x2<-data[,1]
x1<-data[,2]
y<-data[,3]


#------------------------------------------------------------------------------------------------------------
# Calculating psi using function "psi"
#------------------------------------------------------------------------------------------------------------


psidata<-psi(N,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,x2,x1,y)

#------------------------------------------------------------------------------------------------------------
# Calculating \hat{I}_{Y}(\theta)
#------------------------------------------------------------------------------------------------------------


sumIy<-matrix(0,parnum,parnum)
for (j in 1:length(y)){
    Epsiy<-matrix(0,1,parnum)
    for (k in 1:x2groups){
        for (l in 1:x1groups){
            psidatajs<-psi(1,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(k-1),(l-1),y[j])
            Epsiy<-Pz(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,(k-1),(l-1),y[j])/Py(B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,y[j])*psidatajs+Epsiy
         }
    }
    sumIy<-t(Epsiy)%*%Epsiy+ sumIy
}
Iy<-sumIy/N




#------------------------------------------------------------------------------------------------------------
# Running Monte Carlo simulations
#------------------------------------------------------------------------------------------------------------




varperpay0s<-matrix(NA,parnum,length(pAyhx1l)*length(PAyl))
varperpay0smulti<-matrix(NA,parnum,length(pAyhx1l)*length(PAyl))
pAhat1<-matrix(NA,1,length(PAyl))
pAhat12<-matrix(NA,1,length(pAyhx1l)*length(PAyl))

for (j in 1: length(PAyl)){
    w1<-weightsX1(y,cuty,PAyl[j],PAyh[j])
    pAhat1[,j]<-sum(w1)/N
    NstarX1<-N*pAhat1[,j]
    for (i in 1:length(pAyhx1l)){
        w2<-weightsX2(x1,y,cutx1,cuty,pAyhx1h[i]*PAyh[j],pAyhx1l[i]*PAyh[j],pAylx1h[i]*PAyl[j],pAylx1l[i]*PAyl[j])
        w12s<-matrix(c(w1,w1),N,2)*matrix(weightsX2(c(rep(0,N),rep(1,N)),c(y,y),cutx1,cuty,pAyhx1h[i]*PAyh[j],pAyhx1l[i]*PAyh[j],pAylx1h[i]*PAyl[j],pAylx1l[i]*PAyl[j]),N,2)
        w12<-w1*w2
        pAhat12[,((j-1)*length(pAyhx1l)+i)]<-sum(w12)/N
        # The sample size after ascertainment= sum of the ascertainment prob.
        NstarX12<-N*pAhat12[,((j-1)*length(pAyhx1l)+i)]
        varperpay0<-variances(N,B0x2,B0x1,Bx2x1,B0y,Bx2y,Bx1y,SDy,x2groups,x1groups,x2,x1,y,w1,w2,w12,w12s,NstarX1,NstarX12,parnum)
        varperpay0s[,((j-1)*length(pAyhx1l)+i)]<-diag(varperpay0$variancematrixjoint)
        varperpay0smulti[,((j-1)*length(pAyhx1l)+i)]<-diag(varperpay0$variancematrixmulti)

    }
}


#------------------------------------------------------------------------------------------------------------
# Saving the results to your working directory
#------------------------------------------------------------------------------------------------------------

results<-list(title="Efficiency of samling schemes for three-stage model", created=date(),N=N,B0x2=B0x2,x2groups=x2groups,B0x1=B0x1,Bx2x1=Bx2x1,B0y=B0y,Bx2y=Bx2y,Bx1y=Bx1y,SDy=SDy,cuty=cuty,cutx1=cutx1,PAyl=PAyl,PAyh=PAyh,pAyhx1h=pAyhx1h,pAyhx1l=pAyhx1l,pAylx1h=pAylx1h,pAylx1l=pAylx1l,varperpay0s=varperpay0s,varperpay0smulti=varperpay0smulti,pAhat1=pAhat1,pAhat12=pAhat12)
dput(results, file = "EfficiencyTreestage.txt")




#============================================================================================================
# Creating plots
#============================================================================================================

#------------------------------------------------------------------------------------------------------------
#Retrieving your results from file in your working directory
#------------------------------------------------------------------------------------------------------------

ResultsFromFile<-dget(file = "EfficiencyTreestage.txt")

#If you have trouble saving the file you can instead use the following line
#ResultsFromFile<-results



#------------------------------------------------------------------------------------------------------------
# Chosing your cost functions
#------------------------------------------------------------------------------------------------------------

C1<-0
C2<-0
C3<-1



#------------------------------------------------------------------------------------------------------------
#Preparation for plots
#------------------------------------------------------------------------------------------------------------


varBx2x1<-ResultsFromFile$varperpay0smulti[3,]
varBx2y<-ResultsFromFile$varperpay0smulti[5,]
varBx1y<-ResultsFromFile$varperpay0smulti[6,]


urval2<-ResultsFromFile$pAylx1l

# Function to calculate total avarage cost (TAC)
TAC<-function(C1,C2,C3){C1*(rep(1,times=length(ResultsFromFile$pAhat12))-rep(ResultsFromFile$pAhat1, each=length(ResultsFromFile$pAylx1l)))+C2*(rep(ResultsFromFile$pAhat1, each=length(ResultsFromFile$pAylx1l))-ResultsFromFile$pAhat12)+C3*ResultsFromFile$pAhat12}



#-----------------------------------------------------------------
# Threedimensional plots
#-----------------------------------------------------------------


x<-ResultsFromFile$PAyl
y<-ResultsFromFile$pAylx1l


cost<-TAC(C1,C2,C3)

costeffBx2x1<-(varBx2x1[length(varBx2x1)]/varBx2x1)/cost
costeffBx2y<-(varBx2y[length(varBx2y)]/varBx2y)/cost
costeffBx1y<-(varBx1y[length(varBx1y)]/varBx1y)/cost

effsummary<-((varBx2x1[length(varBx2x1)]/varBx2x1)^(-1)+(varBx2y[length(varBx2y)]/varBx2y)^(-1)+(varBx1y[length(varBx1y)]/varBx1y)^(-1))^(-1)
costeffsumamry<-effsummary/effsummary[length(varBx2y)]/cost
costeffs<-matrix(c(costeffBx2x1,costeffBx2y,costeffBx1y),length(cost),3)


z1<-(matrix(costeffs[,1],length(y),length(x)))
z2<-(matrix(costeffs[,2],length(y),length(x)))
z3<-(matrix(costeffs[,3],length(y),length(x)))
z4<-(matrix(costeffsumamry,length(y),length(x)))

#settings for plotting
op<-par(mfrow=c(2,4))
op <- par(bg = "white")
op<-par(cex=0.4, cex.main=2,cex.axis=1.2,cex.lab=1.8)


#hat(beta)[x2Y]
persp(y,x, z1, theta = 30, phi = 30, expand = 0.5, col = "darkolivegreen3",
      ltheta = 120, shade = 0.75, ticktype = "detailed",nticks=4,
      xlab="\n a        ",ylab="\n b      ",zlab="\n CE"
)
title(main=expression(paste('CE  ',hat(beta)[x2Y])),font.main = 4, cex.main = 3)


#hat(beta)[x2x1]
persp(y,x, z2, theta = 30, phi = 30, expand = 0.5, col = "darkolivegreen3",
      ltheta = 120, shade = 0.75, ticktype = "detailed",nticks=4,
      xlab="\n a        ",ylab="\n b      ",zlab="\n CE"
)
title(main=expression(paste('CE  ',hat(beta)[x2x1])),font.main = 4, cex.main = 3)

#hat(beta)[x1Y]
persp(y,x, z3, theta = 30, phi = 30, expand = 0.5, col = "darkolivegreen3",
      ltheta = 120, shade = 0.75, ticktype = "detailed",nticks=4,
      xlab="\n a        ",ylab="\n b      ",zlab="\n CE"
)
title(main=expression(paste('CE  ',hat(beta)[x1Y])),font.main = 4, cex.main = 3)

#summary measure of the three parameters
persp(y,x, z4, theta = 30, phi = 30, expand = 0.5, col = "orangered2",
      ltheta = 120, shade = 0.75, ticktype = "detailed",nticks=4,
      xlab="\n a        ",ylab="\n b      ",zlab="\n CE"
)
title(main=expression(paste('CE')),font.main = 4, cex.main = 3)


#----------------------------------------------------
# Two-dimensional plots with efficinecy illustrated by colour
#----------------------------------------------------

op<-par(cex=0.6, cex.main=1.6,cex.axis=1.2,cex.lab=1.2)

#hat(beta)[x2Y]
image(y, x, z1, col = terrain.colors(100), axes = TRUE, xlab="a",ylab="b")
contour(y, x, z1, levels = seq(0, 5, by = 0.1),
        add = TRUE, col = "peru")
box()

#hat(beta)[x2x1]
image(y, x, z2, col = terrain.colors(100), axes = TRUE, xlab="a",ylab="b")
contour(y, x, z2, levels = seq(0, 5, by = 0.1),
        add = TRUE, col = "peru")
box()

#hat(beta)[x1Y]
image(y, x, z3, col = terrain.colors(100), axes = TRUE, xlab="a",ylab="b")
contour(y, x, z3, levels = seq(0, 5, by = 0.1),
        add = TRUE, col = "peru")
box()

#summary measure of the three parameters
image(y, x, z4, col = heat.colors(100), axes = TRUE, xlab="a",ylab="b")
contour(y, x, z4, levels = seq(0, 5, by = 0.1),
        add = TRUE, col = "peru")
box()