### R function for computing the partial expected value of perfect information when inputs are correlated

### As described in Strong, M. and Oakley, J.E. 
### An efficient method for computing partial expected value of perfect information for correlated inputs

### September 2011
### m.strong@sheffield.ac.uk

### Assume P input parameters, D decicion options, and sample size of S

### Function has the following arguments:
### 1. Sample from the inputs as matrix `inputs', sample sets row-wise. S rows, P columns.
### 2. Input of interest. Must be a single input.
### 3. Sample from costs output of model. Costs for each decision column-wise. S rows, D columns.
### 4. Sample from effects output of model. Costs for each decision column-wise. S rows, D columns.
### 5. Willingness to pay for one unit health effect.
### 6. Sample set is partitioned into K blocks of J samples each. Total sample size must be J*K.


partial.evpi.function<-function(inputs,input.of.interest,costs,effects,lambda,J,K)
{
  S <- nrow(inputs) # number of samples
  if(J*K!=S) stop("The number of samples does not equal J times K")
  D <- ncol(costs) # number of decision options

  nb <- lambda*effects-costs
  baseline <- max(colMeans(nb))
  perfect.info <- mean(apply(nb,1,max))
  evpi <- perfect.info-baseline

  sort.order <- order(inputs[,input.of.interest])
  sort.nb <- nb[sort.order,]

  nb.array <- array(sort.nb,dim=c(J,K,D))
  mean.k <- apply(nb.array,c(2,3),mean)
  partial.info <- mean(apply(mean.k,1,max))
  partial.evpi <- partial.info-baseline
  partial.evpi.index <- partial.evpi/evpi

  return(list(
    baseline = baseline,
    perfect.info = perfect.info,
    evpi = evpi,
    partial.info = partial.info,
    partial.evpi = partial.evpi,
    partial.evpi.index = partial.evpi.index
  ))
}