estimate.default: Influence function based inference

Description

Primary tool for obtaining parameter estimates with robust (sandwich) standard errors, applying the delta method, and testing linear hypotheses. The function returns an object of class estimate which serves as a general container for parameter estimates and their influence functions (IFs). Three calling conventions are supported:

Usage

# S3 method for default
estimate(
  x = NULL,
  f = NULL,
  ...,
  data,
  id,
  coef,
  IC = TRUE,
  vcov,
  stack = TRUE,
  average = FALSE,
  subset,
  keep,
  use,
  regex = FALSE,
  ignore.case = FALSE,
  print = NULL,
  labels,
  label.width,
  contrast,
  null,
  level = NULL,
  type = NULL,
  var.adj = NULL,
  df = NULL,
  back.transform = NULL
)

Value

Object of class estimate with the following elements:

coef: Named vector of parameter estimates.
vcov: Variance-covariance matrix.
IC: Influence function matrix (observations x parameters).
coefmat: Formatted coefficient table (estimate, std.err, confidence limits, p-value).
id: Cluster/id variable used.
ncluster: Number of clusters.
n: Number of observations.
compare: (When null or contrasts are specified) Wald test result.

Arguments

x

model object (glm, lvmfit, ...) or an existing estimate object. When two model objects are supplied (e.g., estimate(g, g0)) a likelihood-ratio test is performed.

f

transformation of model parameters. Accepts several input types: - A function f(p) or f(p, data): applies the delta method. When f returns a named list the names are used as parameter labels. - A matrix: used as a contrast (linear combination) matrix. - A numeric vector of parameter indices: converted to a contrast that selects and differences those parameters. - A list of indices: each element selects one parameter. - Character expressions: supports wildcards ("?", "*") and arithmetic on parameter names (e.g., "z" - "x", 2 * "z" - 3 * "x").

...

additional arguments to lower level functions

data

data.frame used by f when the transformation depends on covariates (see average). Defaults to model.frame(x).

id

(optional) cluster identifier. Can be a vector of cluster IDs, a one-sided formula (evaluated in data), a single character column name, or a logical scalar (TRUE for one-to-one matching, FALSE for independence). When supplied, the IF is aggregated within clusters to produce cluster-robust standard errors.

coef

(optional) named parameter vector. Used instead of coef(x) when constructing an estimate object without a model.

IC

if TRUE (default) the influence function matrix is estimated and stored in the returned object (extract with the IC method). Can also be a user-supplied IF matrix (one row per observation, one column per parameter), which is used directly instead of estimating it from x.

vcov

(optional) covariance matrix of parameter estimates, or a logical. If TRUE, stats::vcov is used to obtain the (model-based) covariance matrix from x, yielding non-robust standard errors. If a matrix is supplied it is used directly. When omitted or FALSE, robust standard errors are computed from the influence function.

stack

if TRUE (default) the influence function contributions are summed within each cluster defined by id. Set to FALSE to keep the un-stacked (per-observation) decomposition.

average

if TRUE the function computes the standardized (marginalized) estimate $\hat\Psi = P_n f(X; \hat\theta)$, i.e., the empirical mean of f(p, data) over all rows of data. The influence function accounts for both the empirical averaging and the parameter estimation uncertainty (see Details).

subset

(optional) logical vector, expression evaluated in data, or column name. When used together with average = TRUE, the average is conditioned on the subpopulation where subset is TRUE, yielding a conditional marginalized estimate.

keep

(optional) index of parameters to keep from final result. Accepts integer indices, character names, or (with regex = TRUE) perl-compatible regular expressions.

use

(optional) index of parameters to use in calculations. The selected parameters are first extracted (via keep) and then the remaining arguments (f, contrast, etc.) are applied to this subset.

regex

if TRUE use perl-compatible regular expressions for keep and use arguments

ignore.case

ignore case in regular expressions

print

(optional) custom print function for the resulting estimate object

labels

(optional) character vector of coefficient names

label.width

(optional) max display width of labels

contrast

(deprecated, use summary method) contrast matrix for a final Wald test. When supplied together with null, tests $H_0: B\theta = b_0$.

null

(deprecated, use summary method) null hypothesis vector $b_0$ to test against (default 0)

level

(deprecated, use summary method) level of confidence limits (default 0.95)

type

(deprecated, use summary method) type of small-sample correction for cluster-robust variance. One of: - "robust" (default): no correction.

"df": applies $n/(n-p)$ correction (Mancl & DeRouen, 2001). - "mbn": Morel-Bokossa-Neerchal (2003) correction. - "hc3": leverage-adjusted HC3-type correction (blended with var.adj). - "hc4": Cribari-Neto (2004) leverage-adjusted correction.

var.adj

(deprecated, use summary method) blending parameter for the HC3 leverage adjustment (default 0.25). Controls the weight between observation-level empirical leverage and the average leverage $p/n$.

df

(deprecated, use summary method) degrees of freedom for t-based inference (default: NULL for Gaussian approximation; when set, confidence intervals and p-values use the t-distribution with df degrees of freedom).

back.transform

(deprecated, use summary method) function applied to the point estimates and confidence interval bounds after inference is performed on the original scale. Useful for variance-stabilizing transformations, e.g., compute CIs on the atanh (Fisher z) scale and back-transform with tanh.

Influence functions and robust standard errors

An estimator $\widehat{\theta}$ is regular and asymptotically linear (RAL) when it admits the iid decomposition $$\sqrt{n}(\widehat{\theta}-\theta) = \frac{1}{\sqrt{n}}\sum_{i=1}^n \mathrm{IC}(Z_i; P) + o_p(1)$$ where $\mathrm{IC}$ is the unique influence function satisfying $E\{\mathrm{IC}(Z; P)\} = 0$. By the central limit theorem $$\sqrt{n}(\widehat{\theta}-\theta) \overset{d}{\longrightarrow} N(0,\; \mathrm{Var}\{\mathrm{IC}(Z; P)\})$$ and the asymptotic variance is consistently estimated by the empirical variance of the plugin IF estimate, yielding robust (sandwich) standard errors. The estimated IF can be extracted with the IC method.

Parameter transformations (delta method)

When f is a function $\phi: R^p \to R^m$, the delta method is applied: $$\sqrt{n}\{\phi(\widehat{\theta}) - \phi(\theta)\} = \frac{1}{\sqrt{n}}\sum_{i=1}^n \nabla\phi(\theta)\,\mathrm{IC}(Z_i; P) + o_p(1)$$ Derivatives are computed numerically via numDeriv::jacobian unless the function returns an attribute "grad" with the analytic Jacobian.

Alternatively, estimate objects support direct arithmetic operations (e.g., a * b, exp(a), a^b) which apply the delta method with exact (analytical) derivatives computed automatically. This influence function calculus allows building complex transformations from simple building blocks without numerical differentiation. See the last example section ("influence function calculus") and vignette("influencefunction", package = "lava") for details.

Averaging and marginalization

When average = TRUE and f(p, data) depends on covariates, the target parameter is the standardized (marginalized) estimate $\Psi = E\{f(X;\theta)\}$. The IF for the averaged estimate accounts for both the empirical averaging and parameter estimation uncertainty: $$\mathrm{IC}_\Psi(Z; P) = f(X;\theta) - \Psi + [E\nabla_\theta f(X;\theta)]\,\phi(Z; P)$$ When subset is also specified, the average is conditioned on the subpopulation, yielding a conditional marginalized estimate.

Cluster-robust standard errors

When id is supplied, the per-observation IF contributions are summed within clusters (when stack = TRUE), producing the cluster-level IF $\widetilde{\mathrm{IC}}(Z_i; P) = \sum_{k=1}^{N_i} \frac{n}{N}\mathrm{IC}(Z_{ik}; P)$. The resulting variance estimate is equivalent to the GEE working independence sandwich estimator.

For full theoretical background and worked examples see vignette("influencefunction", package = "lava").

Details

estimate(x, ...) -- extract estimates from a model object
estimate(coef=, IC=, ...) -- construct from coefficients and IF matrix
estimate(coef=, vcov=, ...) -- construct from coefficients and covariance matrix

Examples

Run this code


## Simulation from logistic regression model
m <- lvm(y~x+z);
distribution(m,y~x) <- dist_bernoulli("logit")
d <- sim(m,1000)
g <- glm(y~z+x,data=d,family=binomial())
g0 <- glm(y~1,data=d,family=binomial())

## LRT
estimate(g, g0)


estimate(g)

## Testing contrasts
summary(estimate(g), null=0)
estimate(g, rbind(c(1,1,0), c(1,0,2)))
summary(estimate(g, rbind(c(1,1,0), c(1,0,2))), null=c(1,2))
estimate(g, 2:3) ## same as cbind(0,1,-1)
estimate(g, as.list(2:3)) ## same as rbind(c(0,1,0),c(0,0,1))
## Alternative syntax
estimate(g, "z", "z"-"x", 2*"z"-3*"x")
estimate(g, "?")  ## Wildcards
estimate(g, "*Int*", "z")
summary(estimate(g, "1", "2"-"3"), null = c(0,1))
estimate(g, 2, 3)

## Usual (non-robust) confidence intervals
estimate(g, vcov=TRUE)
estimate(g, vcov=vcov(g))

## Transformations
estimate(g, function(p) p[1]+p[2])

## Multiple parameters
e <- estimate(g, function(p) c(p[1]+p[2], p[1]*p[2]))
e
vcov(e)

## Label new parameters
estimate(g, function(p) list("a1"=p[1]+p[2], "b1"=p[1]*p[2]))
#'
## Multiple group
m <- lvm(y~x)
m <- baptize(m)
d2 <- d1 <- sim(m,50,seed=1)
e <- estimate(list(m,m),list(d1,d2))
estimate(e) ## Wrong
ee <- estimate(e, id=rep(seq(nrow(d1)), 2)) ## Clustered
ee
estimate(lm(y~x,d1))

## Marginalize / standardization
f <- function(p,data)
  list(p0=expit(p["(Intercept)"] + p["z"]*data[,"z"]),
       p1=expit(p["(Intercept)"] + p["x"] + p["z"]*data[,"z"]))
e <- estimate(g, f, average=TRUE)
e
estimate(e,diff)
estimate(e,cbind(1,1))

## Clusters and subset (conditional marginal effects)
d$id <- rep(seq(nrow(d)/4),each=4)
estimate(g,function(p,data)
         list(p0=expit(p[1] + p["z"]*data[,"z"])),
         subset=d$z>0, id=d$id, average=TRUE)

## More examples with clusters:
m <- lvm(c(y1,y2,y3)~u+x)
d <- sim(m,10)
l1 <- glm(y1~x,data=d)
l2 <- glm(y2~x,data=d)
l3 <- glm(y3~x,data=d)

## Some random id-numbers
id1 <- c(1,1,4,1,3,1,2,3,4,5)
id2 <- c(1,2,3,4,5,6,7,8,1,1)
id3 <- seq(10)

## Un-stacked and stacked i.i.d. decomposition
IC(estimate(l1,id=id1,stack=FALSE))
IC(estimate(l1,id=id1))

## Combined i.i.d. decomposition
e1 <- estimate(l1,id=id1)
e2 <- estimate(l2,id=id2)
e3 <- estimate(l3,id=id3)
(a2 <- merge(e1,e2,e3))

## If all models were estimated on the same data we could use the
## syntax:
## Reduce(merge,estimate(list(l1,l2,l3)))

## Same:
IC(a1 <- merge(l1,l2,l3,id=list(id1,id2,id3)))

IC(merge(l1,l2,l3,id=TRUE)) # one-to-one (same clusters)
IC(merge(l1,l2,l3,id=FALSE)) # independence


# ------ influence function calculus -------
a <- estimate(coef = c("a" = 0.5), IC = scale(rnorm(10), scale=FALSE), id = 1:10)
b <- estimate(coef = c("b" = 0.8), IC = scale(rnorm(10), scale=FALSE), id = 1:10)

e <- c(a, b) # merge
merge(a, b)
c(e1=a, b) # naming of par
labels(e, c("p1", "p2")) # renaming parameters
e["a"] # subset
subset(e, "a")

# pipes
# c(a, b) |>
#  transform(function(x) x^2) |>
#  subset("a") |>
#  labels("sq")

# Parameter transformation with automatic calculation of derivatives
a * b
(3 * cos(a) / sqrt(b) + 1) / a
expit(c(a,b))
c(sum=sum(e), sum2=a+b,
  prod=prod(e), prod2=a*b)
e %*% e # inner prod.
c(1, 2) %*% e
c(pow = a^b)
a^c(0.5, 2)
c(b=e["a"] * e["b"] / a, also.b=e["b"])

B <- rbind(c(1,-1), c(1,0), c(0,1))
B %*% e
e == 1 # wald-test, null-hypothesis H0: b=1
e == c(1,2)
B %*% e == 1

Run the code above in your browser using DataLab