Fit second-order response-surface model to model objective function, see
rsm; and evaluate on an independent sample. Uses the
rsm package.
runRSM(modx, ..., objective = hydromad.getOption("objective"))
evalRSM(
modx,
rsm.object,
n = 100,
objective = hydromad.getOption("objective"),
method = "latin.hypercube"
)
# S3 method for summary.rsm.hydromad
print(x, digits = max(3L, getOption("digits") - 3L), ...)a model specification created by hydromad. It
should not be fully specified, i.e one or more parameters should be defined
by ranges of values rather than exact values.
The response surface for this model and objective within these bounds should be able to be approximated by a quadratic, which may only be the case close to the optimum.
findUnivariateBounds can be used as a heuristic method to
narrow down bounds
Passed to ccd to select an experimental design.
objective function to use to generate response surface,
given as a function(Q, X, ...). See objFunVal.
rsm object produced by runRSM
Number of independent samples to draw for evaluation of fit
Sampling scheme for evalRSM. See
parameterSets
Placeholder
Placeholder
For runRSM, an rsm object.
For evalRSM, a list with elements:
objective values from rsm object
objective values from hydromad runs
coded.data describing new sampled points
Degrees of freedom, see summary.lm
the
square root of the estimated variance of the random error, see
summary.lm
R^2, the 'fraction of variance
explained by the model', see summary.lm
the above R^2 statistic 'adjusted', penalizing for higher p
Wrapper around ccd, evalPars and
rsm to respectively generate a response-surface design, run a
hydromad model and fit a second-order response-surface model.
Parameters for ccd need to be determined before hand by the
user, e.g. using ccd.pick.
The fit of the resulting response surface model should evaluated, e.g. using
summary.lm. Fit might be improved with better design parameters
ccd or by reducing the ranges of parameters in modx to a
region that can be better approximated by a quadratic.
evalRSM randomly samples an independent set of parameters, evaluates
the objective using the hydromad model and response surface at each point
and calculates measures of fit similarly to summary.lm,
displaying them using the print method.
Shin, Mun-Ju, Joseph H.A. Guillaume, Barry F.W. Croke, and Anthony J. Jakeman. 2015. "A Review of Foundational Methods for Checking the Structural Identifiability of Models: Results for Rainfall-Runoff." Journal of Hydrology 520(November): 1-16. doi:10.1016/j.jhydrol.2014.11.040
Lenth RV (2009) "Response-Surface Methods in R, Using rsm", Journal of Statistical Software, 32(7), 1-17. http://www.jstatsoft.org/v32/i07/
G.E.P. Box, N.R. Draper. 1987. Empirical Model-building and Response Surfaces. John Wiley and Sons, New York
A. Abusam, K.J. Keesman, G. van Straten, H. Spanjers, K. Meinema. 2001. "Sensitivity analysis in oxidation ditch modelling: the effect of variations in stoichiometric, kinetic and operating parameters on the performance indices" J. Chem. Technol. Biot., 76(4):430-438. doi: http://dx.doi.org/10.1002/jctb.39810.1002/jctb.398
the rsm package;
findUnivariateBounds to help reduce bounds to a quadratic
region, ccd.pick to help determine parameters of design
(passed as ...)
## This method depends on the rsm package
library(rsm)
################################################################################
## Load data and define model
data(Cotter)
x <- Cotter[1:1000]
## IHACRES CWI model with exponential unit hydrograph
## an unfitted model, with ranges of possible parameter values
modx <- hydromad(x,
sma = "cwi", routing = "expuh",
tau_s = c(2, 100), v_s = c(0, 1)
)
################################################################################
## Find bounds around best guess
## Best fit, used as initial feasible solution
## fitx=fitBySCE(modx,control=list(ncomplex=20))
## Stored solution
fitx <- hydromad(x,
sma = "cwi", routing = "expuh",
tw = 99.9999996914052, f = 4.82636111146549, scale = 0.00129007217632751,
l = 0, p = 1, t_ref = 20, tau_s = 25.2327996372133, v_s = 0.925195113458179
)
## Identify bounds
thres <- objFunVal(fitx) - 0.05 ## Look for models within 0.05 of best fit
bounds <- findUnivariateBounds(modx, fitx, thres)
#> Error in uniroot(function(val) { :
#> f() values at end points not of opposite sign
#> Error in uniroot(function(val) { :
#> f() values at end points not of opposite sign
#> Error in uniroot(function(val) { :
#> f() values at end points not of opposite sign
## If fit of RSM is poor, probably need to narrow bounds further to concentrate on quadratic region
modx <- update(modx, newpars = bounds)
################################################################################
## Evaluate designs with 4 variables, because there are 4 free parameters
getFreeParsRanges(modx)
#> $tw
#> [1] 54.19842 100.00000
#>
#> $f
#> [1] 2.559639 8.000000
#>
#> $tau_s
#> [1] 12.72902 46.54968
#>
#> $v_s
#> [1] 0.760626 1.000000
#>
ccd.pick(4)
#> n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s N alpha.rot alpha.orth
#> 1 16 2 1 8 1 1 1 1 27 2 2.000000
#> 2 16 4 1 8 2 1 1 1 30 2 2.000000
#> 3 16 6 1 8 3 1 1 1 33 2 2.000000
#> 4 16 8 1 8 4 1 1 1 36 2 2.000000
#> 5 16 10 1 8 5 1 1 1 39 2 2.000000
#> 6 16 9 1 8 5 1 1 1 38 2 2.039608
#> 7 16 9 1 8 4 1 1 1 37 2 1.959592
#> 8 16 7 1 8 4 1 1 1 35 2 2.043016
#> 9 16 7 1 8 3 1 1 1 34 2 1.956039
#> 10 16 5 1 8 3 1 1 1 32 2 2.047065
## Design 5 with n0.c=10 and n0.s=5 is both rotatable and orthogonal,
## as alpha.rot and alpha.orth are identical
## and it uses N=39 model evaluations
## This seems suitable
################################################################################
## Run RSM with n0.c=10 and n0.s=5
set.seed(10) # make replicable
# evals.rsm <- runRSM(modx, n0 = c(10, 5))
#
# ## evaluate fit
# summary.lm(evals.rsm)
# ## Residual standard error of 0.008575
# ## and R2=0.9235 could be ok,
# ## but has room for improvement.
#
# ## Standard results from ?rsm
# summary(evals.rsm)
#
# ## Evaluate on an independent sample
# evalRSM(modx, evals.rsm, n = 100)
#
# ## Plot eigen plots
# eigen.plot(evals.rsm, fixed.axis = TRUE)
# ## Wide/tall ellipses indicate relatively lower identifiability
# ## Rotated ellipses indicate correlation between pairs of parameters
# dev.new()
# eigen.plot(evals.rsm, fixed.axis = FALSE)
# ## Rotation of ellipses is more visible with fixed axes