Inverse estimation is a classical and well-known problem in regression. In simple terms, it involves the use of an observed value of the response to make inference on the corresponding unknown value of the explanatory variable. To our knowledge, however, statistical software is somewhat lacking the capabilities for analyzing these types of problems. In this paper, we introduce investr (which stands for inverse estimation in R), a package for solving inverse estimation problems in both linear and nonlinear regression models.
Consider the regression model
There are three main functions available in the investr package (Greenwell 2014):
calibrate
;
invest
;
plotFit
.
calibrate
operates on objects of class "lm"
and can only be used
when the expectation function has the form
invest
, which calculates the point
estimate and confidence intervals for uniroot
from the stats package. The function plotFit
produces a
scatterplot of the data with fitted regression curve and the option to
add confidence/prediction bands for the response (pointwise or
adjusted). It can be used with single-predictor objects of class "lm"
or "nls"
; however, for objects of class "nls"
, confidence/prediction
bands are based on the linear approximation and can be misleading
(Bates and Watts 1988 65). The development version of investr can
be found on GitHub at https://github.com/w108bmg/investr. To report
bugs or issues, contact the main author or submit them to
https://github.com/w108bmg/investr/issues.
Consider the most common calibration model,
where
It can be shown (see, for example, Draper and Smith 1998 47–51) that
an exact
where plotFit
function in the investr package can be used for drawing
scatterplots with the fitted model and confidence/prediction bands. To
calculate the inversion interval for the linear calibration problem, we
use the calibrate
function and specify the option
interval = "inversion"
(the default) as in the following example.
The data frame arsenic
contains the true amounts of arsenic present in
32 water samples (Graybill and Iyer 1994). Also present, is the
amount of arsenic measured by some field test, which is subject to
error. A new water sample is obtained and subjected to the field test
producing a reading of 3.0 calibrate
function to compute the ML estimate and corresponding 90%
calibration interval based on Equation (2).
library(investr)
mod <- lm(measured ~ actual, data = arsenic)
(res <- calibrate(mod, y0 = 3, interval = "inversion", level = 0.9))
## Figure 1
plotFit(mod, interval = "prediction", level = 0.9, shade = TRUE, col.pred = "skyblue")
abline(h = 3, v = c(res$lower, res$estimate, res$upper), lty = 2)
Running the above block of code produces Figure 1 and the following output to the R console:
estimate lower upper
2.9314 2.6035 3.2587
where estimate
is the ML estimate (1), and
lower
/upper
correspond to the lower/upper bounds of the 90%
inversion interval for calibrate
with
a vector of these values as in
calibrate(mod, y0 = c(3.17, 3.09, 3.16), interval = "inversion", level = 0.9)
If interval = "inversion"
, and the slope of the model is not
ignificant at the specified badfit
is
an "lm"
object for which the slope is not significant at the
calibrate(badfit, y0 = 10, level = 0.9)
will either produce two semi-infinite intervals, e.g.,
Error: The calibration line is not well determined. The resulting
confidence region is the union of two semi-infinite intervals:
( -Inf , -282.0006 ) U ( 393.1267 , Inf)
or the entire real line, e.g.,
estimate lower upper
-97.5987 -Inf Inf
Warning message:
The calibration line is not well determined.
Another common approach to computing calibration intervals is to use the delta method (Dorfman 1938). It is easy to show that an approximate standard error for the ML estimator (1), based on a first-order Taylor series approximation, is given by
Assuming large sample normality for
This is equivalent to taking interval = "Wald"
in the call
to calibrate
. For instance,
calibrate(mod, y0 = 3, interval = "Wald", level = 0.9)
produces the output
estimate lower upper se
2.9314 2.6040 3.2589 0.1929
The point estimate remains unchanged (as expected) but the calibration
interval is slightly different. The major benefit of using the delta
method to compute calibration intervals is that it always exists and
provides us with an asymptotic estimate of the standard error, which in
this case gives
Instead of inferring the value of
which is the same as Equation (1) but with
where mean.response = TRUE
.
To illustrate, consider the crystal growth data taken from
Graybill and Iyer (1994 119). These data are from an experiment in
which the weight in grams of 14 crystals were recorded after letting the
crystals grow for different (predetermined) amounts of time in hours.
The weight of each crystal is plotted against time in
Figure 3. Suppose we want to estimate the growth
time in hours that corresponds to an average weight of 8 grams; that is,
we want to estimate
mod <- lm(weight ~ time, data = crystal)
(res <- calibrate(mod, y0 = 8, mean.response = TRUE))
## Figure 3
plotFit(mod, interval = "confidence", pch = 19, shade = TRUE, col.conf = "plum",
extend.range = TRUE)
abline(h = 8, v = c(res$lower, res$estimate, res$upper), lty = 2)
The output for this code chunk should be
estimate lower upper
15.8882 14.6590 17.1596
Thus, in order to produce crystals with an average weight of 8 grams,
they should be grown for an estimated 15.8882 hours. A 95% confidence
interval for the growth time is (14.6590, 17.1596). Obviously, if
mean.response = TRUE
, then y0
can only take a single value;
otherwise, an error will be displayed as in the following:
calibrate(mod, y0 = c(8, 9), mean.response = TRUE)
which displays the message
Error in calibrate.default(cbind(x, y), ...) :
Only one mean response value allowed.
This type of calibration problem is similar to computing the median
effective dose (glm
function from the stats package. The function dose.p
from the
MASS package
(Venables and Ripley 2002) can then be used to compute the ML estimate of
The calibration intervals discussed so far are one-at-a-time intervals.
If adjust
option
which can take any of the following arguments: "none"
(the default),
"bonferroni"
, or "scheffe"
. The value
In application, many relationships are nonlinear (e.g., dose-response
curves). The classical estimator along with the inversion and Wald
intervals can easily be extended to such nonlinear calibration curves.
However, classical inference in these models (such as prediction
intervals) are based on large samples and linear approximations (see Bates and Watts 1988 2). Thus, for nonlinear calibration curves,
the inversion interval provides only an approximate
The invest
function can be used for inverse estimation with any
univariate regression model in R that inherits from class "lm"
or
"nls"
(with the exception of weighted fits). For instance, consider
the regression model
invest
function computes predict
and uniroot
to solve
Equation (2) gives a closed-form expression for the inversion interval for the case of simple linear regression. In more complicated cases, such an expression is not available and the interval must be found numerically. It can be shown (see, for example, Seber and Wild 2003 245–246) that
where
Essentially, invest
finds the lower and upper limits for this interval
by solving the equations
The Wald interval for
Since the point estimate and confidence intervals for uniroot
, invest
has the additional options
lower
, upper
, tol
, and maxiter
. See the reference manual for
details.
To get an idea of how invest
works, consider the data in
Figure 4 which were generated from the model
invest
function.
set.seed(101) # for reproducibility
x <- rep(seq(from = 0, to = 25, by = 2), each = 2)
d <- data.frame(x, y = 5 + x - sin(x) + rnorm(length(x), sd = 1.5))
mod <- lm(y ~ I(x - sin(x)), data = d)
res <- invest(mod, y0 = 22, interval = "Wald")
## Figure 4
plotFit(mod, interval = "prediction", shade = TRUE, col.pred = "seagreen1",
extend.range = TRUE)
abline(h = 22, v = res$estimate, lty = 2)
The point estimate is invest
for obtaining the
point estimate is essentially
uniroot(function(x) predict(mod, newdata = list("x" = x)) - 22,
lower = min(d$x), upper = max(d$x))$root
The code for obtaining the standard error used in the Wald interval is slightly more involved:
## Write x0.hat as function of parameters (theta1.hat, theta2.hat, Y0)
x0.fun <- function(params, object = mod) {
object$coefficients <- params[1:2]
uniroot(function(x) predict(object, list("x" = x)) - params[3],
lower = 0, upper = 25, tol = 1e-10)$root
}
## Variance-covariance matrix of (theta1.hat, theta2.hat, Y0)'
covmat <- diag(3)
covmat[1:2, 1:2] <- vcov(mod)
covmat[3, 3] <- summary(mod)$sigma^2
## Numerically evaluate gradient
params <- c(coef(mod), y0 = 22)
grad <- attr(numericDeriv(quote(x0.fun(params)), "params"), "gradient")
## Calculate standard error
(se <- sqrt(grad %*% covmat %*% t(grad)))
The following example uses the nasturtium data from the drc package.
These data were analyzed in Racine-Poon (1988) using an approximate
Bayesian approach. Bioassays were performed on a type of garden cress
called nasturtium. Six replicates of the response (plant weight in mg
after the third week of growth) were measured at seven preselected
concentrations of an agrochemical. The weights corresponding to three
new soil samples, all sharing the same (unknown) concentration
## Load package containing nasturtium data
library(drc)
## Fit log-logistic model
mod <- nls(weight ~ theta1/(1 + exp(theta2 + theta3 * log(conc))),
start = list(theta1 = 1000, theta2 = -1, theta3 = 1),
data = nasturtium)
plotFit(mod, interval = "prediction") # figure not shown
## Compute approximate 95% inversion interval
invest(mod, y0 = c(309, 296, 419), interval = "inversion")
The interval obtained is (1.7722, 2.9694) with a point estimate of
2.2639. We can check the point estimate manually. Some algebra gives
A Wald interval can also be obtained by running the following line of code:
invest(mod, y0 = c(309, 296, 419), interval = "Wald")
the output for which is
estimate lower upper se
2.2639 1.6889 2.8388 0.2847
In certain cases (i.e., when deltaMethod
function from
the car package
(Fox and Weisberg 2011) to compute the approximate standard error used in the Wald
interval. For the nasturtium example, the minimal code to obtain
## Using the deltaMethod function in the car package
library(car)
covmat <- diag(4)
covmat[1:3, 1:3] <- vcov(mod)
covmat[4, 4] <- summary(mod)$sigma^2 / 3 # since length(y0) = 3
(se <- deltaMethod(c(coef(mod), y0.bar = mean(c(309, 296, 419))),
g = "exp((log(theta1 / y0.bar - 1) - theta2) / theta3)",
vcov. = covmat)$SE)
which produces an estimated standard error of 0.2847019, the same as
that obtained automatically by invest
.
As indicated by the bootstrap distribution obtained in the next section,
the symmetric Wald interval obtained here seems unrealistic for this
problem. In the next section, we show how to obtain a bias-corrected and
accelerated (
The bootstrap (Efron 1979) provides an alternative means for
computing calibration intervals. This is useful in nonlinear calibration
problems where inference traditionally relies on large samples,
approximate normality, and linear approximations. A practical guide to
the bootstrap is provided by Davison and Hinkley (1997) and the companion R
package boot, and bootstrap resampling for controlled calibration is
discussed in Jones and Rocke (1999). Although it is likely for a
"bootstrap"
option to appear in a future release of investr, it is
quite simple to set up using the recommended boot package. First, we
discuss a naive approach to calculating bootstrap calibration intervals.
A simple, but ultimately wrong approach to resampling in controlled calibration is demonstrated in the following example for the nasturtium data:
library(boot)
## Function to compute estimate of x0
x0.fun <- function(object, y) {
theta <- unname(coef(object))
exp((log(theta[1] / mean(y) - 1) - theta[2]) / theta[3])
}
## Bootstrap setup
y0 <- c(309, 296, 419)
res <- resid(mod) - mean(resid(mod)) # center the residuals
n <- length(res)
boot.data <- data.frame(nasturtium, res = res, fit = fitted(mod))
boot.fun <- function(data, i) {
boot.mod <- nls(fit + res[i] ~ theta1 / (1 + exp(theta2 + theta3 * log(conc))),
start = list(theta1 = 1000, theta2 = -1, theta3 = 1), data = data)
## Make sure the original estimate also gets returned
if (all(i == 1:n)) x0.fun(mod, y0) else x0.fun(boot.mod, y0)
}
## Run bootstrap simulation (takes about 50s on a standard laptop)
set.seed(123) # for reproducibility
res <- boot(boot.data, boot.fun, R = 9999) # collect 9,999 bootstrap samples
boot.ci(res, type = "bca") # obtain BCa confidence interval for x0
This produces an approximate 95%
invest(mod, y0 = mean(c(309, 296, 419)), mean.response = TRUE)
We can get a hold of this extra source of variability by again sampling
with replacement from the centered residuals. In particular, we resample
the responses boot.fun
(the changes are colored magenta):
boot.fun <- function(data, i) {
boot.mod <- nls(fit + res[i] ~ theta1 / (1 + exp(theta2 + theta3 * log(conc))),
start = list(theta1 = 1000, theta2 = -1, theta3 = 1), data = data)
## Simulate the correct variance
Y0 <- y0 + sample(data$res, size = length(y0), replace = TRUE)
## Make sure the original estimate also gets returned
if (all(i == 1:n)) x0.fun(mod, y0) else x0.fun(boot.mod, Y0)
}
Rerunning the previous example with the modified boot.fun
function
produces a more realistic 95% confidence interval for
We introduced the investr package for computing point estimates along
with inversion and Wald-based confidence intervals for linear and
nonlinear calibration problems with constant variance. We also showed
how the boot package can be used for constructing approximate
"lme"
from the recommended nlme package;
Pinheiro et al. (2014)]), and even multivariate calibration problems (e.g.,
objects of class "mlm"
).
The authors would like to thank two anonymous reviewers and the Editor for their helpful comments and suggestions.
Agriculture, ChemPhys, Distributions, Econometrics, Environmetrics, Finance, MixedModels, NumericalMathematics, Optimization, Psychometrics, Robust, Survival, TeachingStatistics, TimeSeries
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Greenwell & Kabban, "investr: An R Package for Inverse Estimation", The R Journal, 2014
BibTeX citation
@article{RJ-2014-009, author = {Greenwell, Brandon M. and Kabban, Christine M. Schubert}, title = {investr: An R Package for Inverse Estimation}, journal = {The R Journal}, year = {2014}, note = {https://rjournal.github.io/}, volume = {6}, issue = {1}, issn = {2073-4859}, pages = {90-100} }