Routed cost curve of a design: the travel constant of Hansen, Hurwitz and Madow
Source:R/costcurve.R
routed_cost_curve.RdSelects samples of each size in n_grid, routes each one from the depot
(or bases), prices it with the cost model, and fits the expected cost as
G(n) = c0 + a n + b sqrt(n). This is the cost function of Hansen, Hurwitz
and Madow (1953, Vol. I, Chapter 6, Section 12), C0 sqrt(m) + C1 m, with b in
the role of their travel constant C0 and a of the cost per primary unit
C1. They derived C0 from points on a regular grid and called it a rough
approximation; here it is measured on the frame, with its roads (when the
matrix comes from OSRM), its depots and its route limits. The travel of a
tour through n scattered points grows like sqrt(n) (Beardwood, Halton
and Hammersley 1959), and fixed costs per visit and per route add the
linear term.
Usage
routed_cost_curve(
frame,
depot,
cost_model,
n_grid,
size = NULL,
strata = NULL,
selection = c("lpm", "systematic", "srs"),
replicates = 1,
matrix = NULL,
method = c("haversine", "euclidean", "osrm"),
max_length = Inf,
max_stops = Inf,
n_rep = 10,
iterations = 100,
seed = 1
)
# S3 method for class 'fieldopt_cost_curve'
predict(object, newdata, type = c("total", "average", "marginal"), ...)
# S3 method for class 'fieldopt_cost_curve'
tidy(x, ...)
# S3 method for class 'fieldopt_cost_curve'
glance(x, ...)
# S3 method for class 'fieldopt_cost_curve'
autoplot(object, ...)Arguments
- frame
Data frame of units as in
select_units(), including the depot row named indepot.- depot
Name of the depot unit in
frame$unit.- cost_model
- n_grid
Sample sizes to route (at least three distinct values for the three coefficients).
- size
Optional size column for probability-proportional-to-size selection.
- strata
Optional stratum column, passed to
select_units().- selection
Selection method of
select_units():"lpm","systematic"or"srs".- replicates
Replicates of a
"systematic"selection.- matrix
Optional
travel_matrix()of the frame (computed from the coordinates otherwise).- method
Distance method of
travel_matrix()whenmatrixisNULL("osrm"queries the road network once for the whole frame).- max_length, max_stops
Route limits passed to
route_fieldwork().- n_rep
Routed samples per sample size.
- iterations
GRASP iterations per routing.
- seed
Seed.
- newdata
Sample sizes (a numeric vector) for
predict().- type
"total"forG(n),"average"forG(n)/n,"marginal"forG'(n).- ...
Unused.
- x, object
A
fieldopt_cost_curve.
Value
An object of class fieldopt_cost_curve: a list with points (a
tibble with n, cost_mean, cost_sd, travel_mean, routes_mean per
routed size), coef (c0, a, b), r_squared, currency and the
settings. predict() gives the total, average or marginal cost at any
n; tidy(), glance() and autoplot() are available.
Details
The fit is constrained to be increasing and concave (a >= 0, b >= 0),
so the marginal cost G'(n) = a + b / (2 sqrt(n)) is positive and falls
with n. The marginal cost is the one the classical allocation formulas
need (see two_stage_design()).
References
Hansen, M. H., Hurwitz, W. N. and Madow, W. G. (1953). Sample Survey Methods and Theory, Vol. I, Chapter 6, Sections 11-15, and Vol. II, Chapter 6, Section 11. Wiley. Beardwood, J., Halton, J. H. and Hammersley, J. M. (1959). The shortest path through many points. Proc. Cambridge Phil. Soc., 55, 299-327.
Examples
set.seed(2)
cells <- expand.grid(x = 1:15, y = 1:15)
cells$unit <- paste0("c", seq_len(nrow(cells)))
frame <- rbind(data.frame(unit = "depot", x = 8, y = 8), cells)
model <- field_cost_model(per_travel = 3, per_unit = 40)
curve <- routed_cost_curve(frame, "depot", model,
n_grid = c(10, 20, 40, 80), method = "euclidean", n_rep = 2,
iterations = 20
)
curve
#>
#> ── Routed cost curve ───────────────────────────────────────────────────────────
#> G(n) = 3.9 + 39.17 n + 46.46 sqrt(n) (cost; 4 sample sizes from 10 to 80, 2
#> routed samples each, lpm selection; R-squared 1).
#> In Hansen, Hurwitz and Madow's notation: C0 = 46.46 (travel), C1 = 39.17 (per
#> primary unit).
#> At n = 30: average cost 47.79, marginal cost 43.41 per unit.
predict(curve, c(30, 60), type = "marginal")
#> [1] 43.41366 42.17136