Build the emulator

We build the emulator for the output variable BALX. We use the built-in R function step to choose the emulator mean function by stepwise regression. step returns an object of class lm, whose list element terms contains the R formula for the mean function. That formula is passed to the DiceKriging function km, along with the design matrix and the response, to fit the Gaussian process emulator. See ?km for options to control the emulator; here we only specify to use a Matern covariance matrix (with nu=5/2) and let all parameters be estimated automatically:

# define response
response = tData$BALX

# stepwise regression to fit trend model
em_lm = step(lm(response ~ ., data=design), 
             direction='both', trace=0)

# fit GP using DiceKriging function `km`
em_km = km(formula  = em_lm$terms, 
           design   = design, 
           response = response, 
           covtype  = 'matern5_2',
           control  = list(trace=FALSE))

The object em_km contains the emulator. We can use predict to evaluate the emulator at new inputs, for example at the mean of all input points:

pred = predict(em_km, newdata=map_df(design, mean), type='UK')
pred[c('mean', 'sd', 'lower95', 'upper95')] %>% as.data.frame

##       mean       sd  lower95  upper95
## 1 6.230499 1.142164 3.935236 8.525763

Validate emulator

The emulator can be validated by leave-one-out cross validation, using the DiceKriging function leaveOneOut.km:

em_km_loo = leaveOneOut.km(em_km, type='UK') %>% as_data_frame
em_km_loo 

## # A tibble: 60 x 2
##          mean       sd
##         <dbl>    <dbl>
##  1 -1.0139507 1.731929
##  2  5.0540275 1.482308
##  3  0.4496219 1.619887
##  4  6.1947459 1.613111
##  5 -5.3986640 1.785493
##  6  6.6694612 1.610073
##  7 10.4919816 1.572071
##  8 -1.7387301 1.383712
##  9 -5.4089395 1.652246
## 10  7.7264413 1.680943
## # ... with 50 more rows

em_km_loo is a list with two vectors mean and sd of equal length, which is why it can be transformed to a data frame. To plot the cross validation results in a faceted panel plot with ggplot2, we have to set up the appropriate data frame with columns for the response, and lower and upper limit of the prediction interval. We also add a boolean variable inside and a corresponding colr variable to change the colour of the intervals if the response falls outside. The variable pit is the probability integral transform, which should be uniformly distributed if the emulator is reliable. At last we append the design matrix, and use gather to turn the data frame from wide into long form appropriate for facet plotting:

# dataframe for plotting
plot_df = em_km_loo %>% 
  mutate(response = response,
         index    = seq_along(response),
         lwr      = mean - 2 * sd, 
         upr      = mean + 2 * sd,
         pit      = pnorm(response, mean, sd),
         inside   = (response > lwr & response < upr),
         colr     = ifelse(inside, 'deepskyblue3', 'darkorange3')) %>%
  bind_cols(design) %>%
  gather('variable', 'value', names(design))

The following plots the emulator output and the actual value of the response over the value of the input in a separate panel for each input.

ggplot(plot_df, aes(x=value)) + facet_wrap( ~ variable, ncol=2) +
  geom_point(aes(y=mean), colour=plot_df$colr) + 
  geom_segment(aes(xend=value, y=lwr, yend=upr), colour=plot_df$colr) +
  geom_point(aes(y=response), cex=1.5) 

The histogram over the PIT values should be approximately flat:

ggplot(plot_df, aes(x=pit)) + 
  geom_histogram(breaks=seq(0,1,0.2), colour='black', fill='deepskyblue3')

Build emulator

response = tData$PCP
em_lm = step(lm(response ~ ., data=design), direction='both', trace=0)
em_km = km(formula=em_lm$terms, design=design, response=response, covtype='matern5_2', control=list(trace=FALSE))

Validate emulator

plot_df = leaveOneOut.km(em_km, type='UK') %>% as_data_frame %>% mutate(response = response, index = seq_along(response), lwr = mean - 2 * sd, upr = mean + 2 * sd, pit = pnorm(response, mean, sd), inside = (response > lwr & response < upr), colr = ifelse(inside, 'deepskyblue3', 'darkorange3')) %>% bind_cols(design) %>% gather('variable', 'value', names(design))

ggplot(plot_df, aes(x=value)) + facet_wrap(~variable, ncol=2) + geom_point(aes(y=mean), colour=plot_df$colr) + geom_segment(aes(xend=value, y=lwr, yend=upr), colour=plot_df$colr) + geom_point(aes(y=response), cex=1.5) 

ggplot(plot_df, aes(x=pit)) + geom_histogram(breaks=seq(0,1,0.2), colour='black', fill='deepskyblue3')