expectation: Mean of (estimated) function via g-computation

Description

Estimates the mean parameter $\Psi = E{f(X; \theta)}$ by the empirical average of f(p, data) over the rows of data, where the parameter $\theta$ is estimated by object. The influence function of the estimate accounts for both the empirical averaging and the uncertainty of $\widehat\theta$. Standardization and g-computation are special applications of this function.

Usage

expectation(
  f,
  object,
  data,
  ...,
  id,
  subset,
  labels = NULL,
  details = FALSE,
  args = list()
)

Value

estimate object

Arguments

f: function with named argument p providing the parameter vector and the optional named arguments data and object (the model object, e.g. an estimate object). Should return a vector, matrix or a (named) list of vectors with one value per observation (row of data unless data is missing). If the result has an attribute "grad" it is used as the Jacobian of f with respect to p
object: model object (glm, lvmfit, ...) or estimate object providing the parameter estimates and their influence function. Alternatively, an evaluated expression of a function with values in data can be provided to calculate the simple mean (see examples).
data: data.frame over which f is averaged. Defaults to model.frame(object) (which is not available for estimate objects).
...: additional arguments passed to f
id: (optional) ids of the values returned by f (the g-computation part): a vector with one value per value of f, or a column name or one-sided formula evaluated in data (requiring one value of f per row of data). Defaults to the names of the values of f, and otherwise rownames(data). id = NULL uses the default ids but removes the id (index) from the returned object.
subset: (optional) subpopulation on which the average is conditioned (conditional marginal mean estimate). Either a logical vector (or numeric 0/1 vector) with one value per value of f, a column name of data, or a one-sided formula evaluated in data (e.g., ~ a == 1 & z > 0). Note that subset uses standard evaluation, i.e., subset = a == 1 is evaluated in the calling environment and not in data. As for model formulas (see stats::model.frame), the variables of a formula subset are looked up in the columns of data and then in the environment where the formula was created (not the environment from which expectation is called).
labels: (optional) character vector of parameter names
details: (optional) if TRUE extra information from the estimation procedure will be stored in the output as attributes: ic.theta contains the term $E[Df(theta) s]IC(theta)$ where $Df$ is gradient of f evaluted in the data and parameter estimate of object and s is the subset indicator. $IC(theta)$ is the influence function of object. The subset attribute holds the subset indicator if provided.
args: (list) arguments to f

Details

The influence function of the estimate is $$\mathrm{IC}_\Psi(Z; P) = f(X;\theta) - \Psi + [E\nabla_\theta f(X;\theta)]\,\phi(Z; P)$$ where $\phi$ is the influence function of $\widehat\theta$.

The two terms are identified by different ids: the first term by id (the rows of data), and the second term by the ids of the model, i.e., index(object) for estimate objects, or the row names of the influence function (for model objects such as glm the row names of the model frame). The model may therefore be estimated on another (e.g., smaller or partly overlapping) dataset than data. The two terms are aligned on the union of ids: each term is zero for ids outside its own support and rescaled by the inverse proportion of observed ids (as in merge.estimate). If there are no common ids, the model estimate and data are treated as independent. Use estimate(object, id=...) to assign ids (or clusters) to the model.

Examples

Run this code

m <- lvm(y ~ x + z)
distribution(m, y ~ x) <- dist_bernoulli("logit")
d <- sim(m, 1000, seed = 1)
fit <- glm(y ~ z + x, data = d, family = binomial())

## Standardization (g-computation)
f <- function(p, data)
  list(p0 = expit(p["(Intercept)"] + p["z"] * data[, "z"]),
       p1 = expit(p["(Intercept)"] + p["x"] + p["z"] * data[, "z"]))
e <- expectation(f, fit)
e
estimate(e, diff)

## Model estimated on a subset, standardized over the full data
d$id <- paste0("i", seq_len(nrow(d)))
d$w <- rbinom(nrow(d), 1, 0.5)
d1 <- subset(d, w == 1)
g1 <- glm(y ~ x + z, data = d1, family = binomial)
e1 <- estimate(g1, id = d1$id)
expectation(f, e1, data = d, id = "id")

## Conditional marginal mean effects (subset) and clusters
d$cl <- rep(seq(nrow(d) / 4), each = 4)
expectation(function(p, data) expit(p[1] + p["z"] * data[, "z"]),
            estimate(fit, id = d$cl),
            data = d, subset = ~ z > 0, id = "cl")

## Simple mean of known function
expectation(~ 1 / (1+z^2), data=d, subset = ~x, labels="p")