## R code to demonstrate algorithm for estimating bias and CI width in partial EVPI calculations
## Accompanies the paper "Simulation Sample Sizes for Monte Carlo Partial EVPI Calculations" 
## by Jeremy E. Oakley, Alan Brennan, Paul Tappenden and Jim Chilcott
## 09/11/09

library(MASS) # required for multivariate normal sampling

## CASE STUDY 1


################################################################################################
################################################################################################
## Net benefit functions for the two treatment options
## X is a matrix, with each column corresponding to one input and each row
## corresponding to  a separate run of the economic model
NB1 <- function(X, lambda){
    lambda * (X[,5] * X[,6] * X[,7] + X[,8] * X[,9] * X[,10]) - (X[,1] + X[,2] *X[,3] * X[,4])
}
    
NB2 <- function(X, lambda){
    lambda * (X[,14] * X[,15] * X[,16] + X[,17] * X[,18] * X[,19]) - (X[,11] +X[,12] * X[,13] * X[,4])
}

################################################################################################
################################################################################################




################################################################################################
################################################################################################
## Input distributions
## Means and standard deviations of the 19 uncertain inputs
meanvector<-c(1000.00, 0.10, 5.20, 400.00, 0.70, 0.30, 3.00, 0.25, -0.10, 0.50, 1500.00, 0.08, 6.10, 0.80, 0.30, 3.00, 0.20, -0.10, 0.50)
stdvector<-c(1, .02, 1, 200, .1, .1, .5,0.1, .02, .2, 1, .02, 1, .1,.05,  1, .05, .02, .2)
n.inputs<-19
################################################################################################
################################################################################################



################################################################################################
################################################################################################
### STAGE 1
### Function to sample values of X_i and estimate mu(X_i) and V(X_i) for each sampled X_i.

Estimate.Mu.and.V<-function(inputnumber,meanvector,stdvector,n.treatments,n.inputs,outerlength,innerlength){
    # Sample X_i
    
    # Option 1: Simple random sampling
    X.i<-rnorm(outerlength,meanvector[inputnumber],stdvector[inputnumber])
    
    # Option 2: Stratified sampling
    #X.i.quantiles<-seq(from=0,to=1-1/outerlength,length=outerlength)+runif(outerlength,0,1/outerlength)
    #X.i<-qnorm(X.i.quantiles,meanvector[inputnumber],stdvector[inputnumber])
   

    # Placeholders for the economic model outputs
    NB1matrix<-matrix(0,innerlength,outerlength)
    NB2matrix<-matrix(0,innerlength,outerlength)

    for(k in 1:outerlength){
        # Sample X_{-i} from conditional distribution of X_{-i} | X_i=X.i[k]
        # In case study 1, this is just distribution of X_{-i}, as inputs are independent
        inputs<-matrix(rnorm(n.inputs*innerlength,meanvector,stdvector),nrow=innerlength,byrow=T) 

        # Fix X_i at its k-th sampled value
        inputs[,inputnumber]<-rep(X.i[k],innerlength) 

        # Evaluate the economic model outputs
        NB1matrix[,k]<-NB1(inputs,lambda=10000)
        NB2matrix[,k]<-NB2(inputs,lambda=10000)
    }

    # Estimate mu(X_i)
    mu.tilde<-matrix(0,n.treatments,outerlength)
    mu.tilde[1,]<-apply(NB1matrix,2,mean) # t=1
    mu.tilde[2,]<-apply(NB2matrix,2,mean) # t=2

    # Estimate V(X_i)
    V.tilde<-array(0,c(2,2,outerlength))
    for(k in 1:outerlength){
        V.tilde[,,k]<-cov(cbind(NB1matrix[,k],NB2matrix[,k]))
    }   
    # Function outputs
    list(mu.tilde=mu.tilde,V.tilde=V.tilde)
}
################################################################################################
################################################################################################



################################################################################################
################################################################################################
# STAGES 2 AND 3
# Function to estimate bias and CI width for given J and K
bias.and.CI<-function(mu.tilde,V.tilde,J,K,N=10000,outerlength){
    mean.bias.vector<-rep(0,outerlength)
    var.bias.vector<-rep(0,outerlength)

    for(k in 1:outerlength){
        # Approximate distribution of \hat{\mu}(X_i) by multivariate normal
        # and sample random mu vectors
        random.mu<-mvrnorm(N,mu.tilde[,k],V.tilde[,,k]/J)
        bias<-apply(random.mu,1,max)-max(mu.tilde[,k])
        mean.bias.vector[k]<-mean(bias)
        var.bias.vector[k]<-var(bias)
    }

    # Estimate overall bias
    overall.bias<-mean(mean.bias.vector)
    
    # Estimate width of 95% confidence interval
    m.X<-apply(mu.tilde,2,max)
    CI.width<-2*1.96*((mean(var.bias.vector)+var(m.X))/K)^.5

    # Function outputs
    list(overall.bias=overall.bias,CI.width=CI.width)
}
################################################################################################
################################################################################################



################################################################################################
################################################################################################
# EXAMPLE: bias and CI width for EVPI(X_5)

# Estimate mu(X_5) and V(X_5) for 21 values of X_5.
X5.mu.and.V<-Estimate.Mu.and.V(inputnumber=5,meanvector=meanvector,stdvector=stdvector,n.treatments=2,n.inputs=19,outerlength=21,innerlength=30)

# Obtain bias and CI width for J=10 and K=10
raw.bias.and.CI<-bias.and.CI(X5.mu.and.V$mu.tilde,X5.mu.and.V$V.tilde,J=10,K=10,N=10000,outerlength=21)

print(raw.bias.and.CI$overall.bias/1319*100) # index to overall EVPI=1319
print(raw.bias.and.CI$CI.width/1319*100) # index to overall EVPI=1319
################################################################################################
################################################################################################