Computer experiments
We start by generating a maximin Latin hypercube, following the run of thumb of having 10 times as many runs as input variables. The function maximinSA_LHS requires a starting design for the optimisation; we will use a random Latin hypercube of the correct size. We can visualise our final design by looking at one-dimensional projections (which should be uniform) and two-dimensional projections (the design criterion doesn’t consider these, but they hopefully will look OK).

library(DiceDesign)
start.design <- lhsDesign(80, 8)$design
maximin.d <- maximinSA_LHS(start.design)
colnames(maximin.d$design) <- c("rw", "r", "Tu", "Hu", "Tl", "Hl", "L", "Kw")
hist(maximin.d$design[, 1])

pairs(maximin.d$design)

To run the experiment, we need to scale our design to the appropriate ranges (DiceDesign finds designs on the scale \([0,1]\)) and pass each row of the design to the borehole function.

library(apts.doe)
l <- c(.05, 100, 63070, 990, 63.1, 700, 1120, 9855)
u <- c(.15, 50000, 115600, 1110, 116, 820, 1680, 12045)
center <- -l / (u - l)
scalev <- 1 / (u - l)
bhdesign <- scale(maximin.d$design, center = center, scale = scalev)
bhdesign
##               rw          r        Tu        Hu        Tl       Hl        L
##  [1,] 0.11537239 29219.0282  75382.12 1068.3173 100.50402 709.3176 1215.375
##  [2,] 0.14546237 33514.5571  85695.92  993.0483 104.21359 801.9650 1439.911
##  [3,] 0.05535548 10359.4135  65444.77 1049.6694  71.83889 813.4733 1474.682
##  [4,] 0.12662909 46364.5476  67959.53 1018.0191  83.26376 744.0609 1226.499
##  [5,] 0.14826760 17918.0545 102395.51 1043.3945  88.60545 755.9380 1466.979
##  [6,] 0.07673155 37013.2550 103611.68 1078.5740  94.62532 703.0235 1149.465
##  [7,] 0.13140425  5431.1244  71252.47 1060.4337 104.08928 748.3192 1377.999
##  [8,] 0.09267319 45991.7186  92676.63 1048.3820  86.82722 753.9493 1556.455
##  [9,] 0.06828096  7684.0230  82597.72 1001.8662  77.86824 793.1268 1199.162
## [10,] 0.10024538 17454.0793 106630.97 1029.5972  74.62280 791.0753 1610.302
## [11,] 0.10745694 23907.8701  74227.83 1093.4344  64.56479 784.1258 1645.538
## [12,] 0.14046010 25035.2237  94832.93 1000.0361  73.79977 798.8818 1320.860
## [13,] 0.12378778  6696.8456 105347.17 1070.7515 112.68672 720.7161 1536.508
## [14,] 0.08256554 35750.2627  90150.01 1046.2768 109.84955 806.0497 1335.705
## [15,] 0.10843667 43652.1238  80056.60 1041.6872  63.89793 819.7385 1409.281
## [16,] 0.11915562   216.4491  90991.28 1007.6579  85.54829 736.6169 1444.136
## [17,] 0.08726074 32474.8338  96346.25 1087.1908  83.61669 789.0849 1171.394
## [18,] 0.05024158 41304.5375  63979.75 1060.7550 110.53604 817.4315 1597.155
## [19,] 0.10336697  4582.1316 109748.44 1051.8512  87.94302 723.3986 1551.503
## [20,] 0.11625584 13937.5755  69929.21 1081.5701 101.84764 802.1231 1155.981
## [21,] 0.09530599 38821.8122 113310.53 1015.1067  67.68323 773.3843 1249.737
## [22,] 0.11026803 48228.7805  89150.53 1078.3435 106.01184 727.9497 1384.578
## [23,] 0.12584459 34047.9232  91695.07 1098.9901  65.62162 715.6370 1271.196
## [24,] 0.05690062 12267.1576  76991.21 1004.0658  84.85474 732.5157 1414.009
## [25,] 0.05323357 18909.4550 104705.24 1025.9688  80.74662 816.8315 1431.135
## [26,] 0.12325839 47955.7351 108796.01 1106.1315 106.93599 781.4195 1206.240
## [27,] 0.13068331 22908.0898 111607.97  994.8169  81.17870 724.4976 1490.263
## [28,] 0.06568001 19807.5055  98743.60  991.1691  70.30063 712.5901 1134.984
## [29,] 0.14462486 39755.8378  78584.60 1030.8242  96.43516 788.3104 1310.067
## [30,] 0.08135119  9676.3349 114436.25 1026.7778  85.88069 758.4766 1306.977
## [31,] 0.14883191 11695.1538  85029.02 1085.3877  79.69802 779.8023 1164.253
## [32,] 0.13807450 31158.4296  67064.50 1108.3249 113.38561 747.2191 1391.777
## [33,] 0.07036292 32882.4044  87082.03 1033.3715  63.58912 795.9810 1584.049
## [34,] 0.14740286 15756.0164  80528.32 1067.0936  90.90622 800.2501 1502.210
## [35,] 0.08556261 49572.4357 113804.02 1034.4461  93.15134 778.3113 1520.349
## [36,] 0.09918574 26619.5901  63528.95  998.1053  68.73569 711.3410 1327.345
## [37,] 0.11429947 36424.4928  95442.48 1083.0503  87.48217 777.3974 1669.518
## [38,] 0.10468956 45326.4390  79146.04 1062.2473  69.07220 707.6860 1513.953
## [39,] 0.05499937 28357.9052 101749.60 1091.2152  89.46995 716.5142 1652.387
## [40,] 0.07286393 31450.7610  74278.94 1075.9705  79.04955 719.2500 1276.114
## [41,] 0.09825578 20667.2211  89889.52 1039.4895  65.88083 737.7054 1219.017
## [42,] 0.09208811 40974.9448  69006.64 1024.3730 106.27045 729.3686 1511.314
## [43,] 0.09656920 12720.2762 100994.63 1005.1530  77.32677 811.6852 1405.582
## [44,] 0.06262735 27124.6050  99465.08 1037.0555 110.73302 746.4228 1143.578
## [45,] 0.08110154  8460.7936  72227.38 1095.0119 114.06277 721.2008 1345.387
## [46,] 0.14297965 35396.1110  77732.18 1053.4431  70.57050 765.2836 1295.570
## [47,] 0.08953049 40625.6981  83868.31 1009.8812 115.01395 700.5589 1179.720
## [48,] 0.11275790  1422.6345  88657.36 1028.2201  73.47946 754.8292 1617.818
## [49,] 0.12888255 41903.5233  87806.97 1019.3981  75.74385 706.8625 1125.641
## [50,] 0.09035706 44222.8789  64845.37 1081.4846  89.78920 770.4595 1364.803
## [51,] 0.05919128 15646.4149 102482.22 1103.1121  92.70216 713.8839 1262.364
## [52,] 0.06175129 16678.5523  94089.25 1016.3615  90.49860 733.3580 1572.307
## [53,] 0.07542878 13607.2687  83226.52 1088.1101  94.15689 759.0658 1609.237
## [54,] 0.10884000 25414.8398  66545.16 1099.7812 103.06002 702.5487 1592.742
## [55,] 0.12197877 38658.5840 104181.05 1035.3087 102.51564 743.3940 1342.045
## [56,] 0.13619583 21272.8531 115399.87 1097.5579  75.64264 796.6270 1254.332
## [57,] 0.05185551 27572.1276  81283.59 1105.1241 105.29440 749.6732 1626.149
## [58,] 0.11786273 23287.1938 105979.37 1109.8261 113.34630 791.8956 1526.425
## [59,] 0.13357524 18555.2707 107552.39 1093.7044  71.14706 751.2677 1395.319
## [60,] 0.06107568 11012.2957  86576.20 1008.5761 111.43871 806.5267 1282.055
## [61,] 0.05844912 42783.6494  98336.54 1057.8017  78.33010 730.8198 1366.445
## [62,] 0.12122857 21859.8515 109093.34 1012.3711 108.34387 771.4973 1635.609
## [63,] 0.07837248 14477.9223  70794.39 1056.3205  82.05365 768.8164 1566.171
## [64,] 0.09381159 34880.8278 107959.41 1051.4048  98.96874 775.3314 1288.599
## [65,] 0.13904817  9375.8857  93330.19 1054.5696  95.06026 739.8963 1239.122
## [66,] 0.07926744 26056.8854  81759.95 1040.9956  98.75985 815.1655 1453.896
## [67,] 0.10155085  2274.8016  84508.31 1066.1041 108.93848 766.4887 1661.411
## [68,] 0.10530698  1150.7812  97625.51 1071.7642  98.10596 803.6990 1427.224
## [69,] 0.06930397  7492.2642  97099.05 1090.2116 101.35303 785.7663 1185.073
## [70,] 0.14191578 49182.9315 100272.85 1064.5311  96.88066 810.4880 1478.020
## [71,] 0.08811211  6277.1382  76088.29 1101.0572  66.98848 774.3849 1190.959
## [72,] 0.06674825 29765.0298  68383.32 1044.1148 115.68107 741.8344 1545.613
## [73,] 0.13473304 22164.6973  65854.30 1073.8928  68.32198 782.6980 1233.671
## [74,] 0.11179465 38069.6631  76504.51  996.9301  72.45594 760.2474 1494.353
## [75,] 0.07217883 30415.7865 112836.86  992.6525  82.30438 763.2768 1641.383
## [76,] 0.08425932  3941.1619  73513.91 1020.9802  95.76458 705.5891 1129.240
## [77,] 0.07486668  3695.3500 112285.99 1074.1623  76.39219 761.5849 1357.239
## [78,] 0.06421231 46931.9994  92319.58 1013.7209 107.91794 735.9758 1575.166
## [79,] 0.12818286 44915.9181  72855.72 1021.5587  99.57134 726.1911 1678.686
## [80,] 0.13744039  3093.0436 110982.75 1002.3575  92.08545 808.8613 1460.183
##              Kw
##  [1,] 12018.159
##  [2,] 10659.358
##  [3,] 10371.356
##  [4,] 11097.064
##  [5,] 10159.971
##  [6,] 10768.605
##  [7,] 10238.328
##  [8,]  9975.980
##  [9,] 11793.153
## [10,] 11848.767
## [11,] 10377.094
## [12,] 10969.724
## [13,] 11927.574
## [14,] 11592.299
## [15,] 11478.995
## [16,] 12007.658
## [17,] 11886.046
## [18,] 10451.922
## [19,] 10874.200
## [20,] 11440.417
## [21,] 11055.051
## [22,] 10721.943
## [23,] 11712.842
## [24,] 10423.417
## [25,] 11028.423
## [26,] 10612.965
## [27,] 11189.162
## [28,] 11547.664
## [29,] 11936.648
## [30,] 11407.783
## [31,] 10472.416
## [32,] 11314.454
## [33,] 10581.374
## [34,] 11062.217
## [35,] 11164.030
## [36,] 10993.274
## [37,] 11990.246
## [38,] 11234.917
## [39,] 10691.952
## [40,] 10303.823
## [41,] 10548.950
## [42,] 11869.202
## [43,]  9867.160
## [44,] 11466.005
## [45,] 11751.526
## [46,] 10087.467
## [47,] 10931.112
## [48,] 10043.607
## [49,] 10140.466
## [50,] 11273.097
## [51,] 11673.577
## [52,] 11739.460
## [53,] 10115.595
## [54,] 10634.730
## [55,] 11634.032
## [56,]  9884.088
## [57,] 11572.273
## [58,] 10270.492
## [59,] 10900.585
## [60,] 10866.570
## [61,] 11809.990
## [62,] 11215.913
## [63,] 11648.332
## [64,] 10007.112
## [65,] 11377.770
## [66,]  9945.304
## [67,] 11347.897
## [68,]  9919.965
## [69,] 10798.987
## [70,] 11283.628
## [71,] 10061.167
## [72,] 10198.403
## [73,] 11509.160
## [74,] 10332.810
## [75,] 10487.524
## [76,] 10817.402
## [77,] 10229.048
## [78,] 10745.249
## [79,] 10521.850
## [80,] 11128.587
## attr(,"scaled:center")
## [1] -0.500000000 -0.002004008 -1.200647249 -8.250000000 -1.192816635
## [6] -5.833333333 -2.000000000 -4.500000000
## attr(,"scaled:scale")
## [1] 1.000000e+01 2.004008e-05 1.903674e-05 8.333333e-03 1.890359e-02
## [6] 8.333333e-03 1.785714e-03 4.566210e-04
y <- apply(bhdesign, 1, borehole)
y
##  [1] 147.24706  93.16504  15.96210 123.48550 136.12401  64.82985 124.97501
##  [8]  50.70603  29.96537  55.07848  70.27019 102.17133 130.05917  44.49503
## [15]  66.14203 100.00442  71.93234  12.61441  76.92288 116.53137  60.55221
## [22] 103.05302 173.32664  20.32551  14.32303 135.26979 107.86412  38.26139
## [29] 143.42115  48.49855 188.62181 174.35139  24.58111 132.77665  43.09265
## [36]  72.86306  89.48699  89.90764  23.00756  47.87238  78.59548  61.47283
## [43]  39.59209  35.82351  67.27342 142.10412  71.87257  67.45110 145.10412
## [50]  65.57051  39.53133  25.26330  36.88517  98.24134 117.35111 136.60791
## [57]  21.34307  92.89031 148.23353  20.03475  30.25394  75.83410  41.12257
## [64]  59.01433 173.69731  30.40927  66.08662  64.67514  41.75676 121.48120
## [71]  66.94430  27.87099 152.65194  63.77259  23.92452  67.12967  41.35870
## [78]  24.49958  94.93653  86.83445

We use DiceKriging and the km function to fit a GP model to this data. We use our original design on the interval \([0,1]\) so that the correlation parameters are comparable (they are all on the same scale).

library(DiceKriging)
gp <- km(formula = ~., design = maximin.d$design, response = y, control = list(trace = F))
gp
## 
## Call:
## km(formula = ~., design = maximin.d$design, response = y, control = list(trace = F))
## 
## Trend  coeff.:
##                Estimate
##  (Intercept)    21.8785
##           rw   144.1632
##            r     0.8479
##           Tu    -1.8313
##           Hu    35.0752
##           Tl     1.6642
##           Hl   -30.1165
##            L   -32.6835
##           Kw    14.6344
## 
## Covar. type  : matern5_2 
## Covar. coeff.:
##                Estimate
##    theta(rw)     0.4887
##     theta(r)     1.9782
##    theta(Tu)     1.9749
##    theta(Hu)     1.4027
##    theta(Tl)     1.9694
##    theta(Hl)     1.4490
##     theta(L)     1.1735
##    theta(Kw)     1.9644
## 
## Variance estimate: 163.1625

The dominant variable is \(r_w\), with a large linear coefficient and a small \(\theta\) (correlation parameter).

We can graphically explore a bit using leave-one-out predictions. If we plot these against the response, we would hope to see a line close to the diagonal, with narrow intervals mostly crossing the line (so the response is usually within the interval). We can also work out an approximate \(R^2\) value as the square of the correlation between the response and the mean prediction. The model looks pretty good on these measures.

library(ggplot2)
loo <- leaveOneOut.km(gp, type = "UK", trend.reestim = T)
summary_df <- data.frame(y = y, pred = loo$mean, ymin = loo$mean - 2 * loo$sd, ymax = loo$mean + 2 * loo$sd)
ggplot(summary_df, aes(x = y, y = pred, color = 100 * maximin.d$design[, "rw"])) + geom_point(size = 2) +
  geom_abline(slope = 1, intercept = 0) +
  geom_errorbar(aes(ymin = ymin, ymax = ymax)) + labs(color = "r_w")

with(summary_df, cor(y, pred)) ^ 2
## [1] 0.9979813