Package 'deli'

Description

Data from AIDS Clinical Trials Group Study 175, used to demonstrate confidence bands for treatment effects.

Usage

actg175

Format

A data frame.

References

Hammer SM, et al. (1996). A trial comparing nucleoside monotherapy with combination therapy in HIV-infected adults with CD4 cell counts from 200 to 500 per cubic millimeter. New England Journal of Medicine, 335(15), 1081-1090.

Description

Constructs the expanded design matrix for generalized additive models by appending spline basis terms according to per-column specifications. Each column in X keeps its linear term; columns with a non-NULL specification get additional spline basis columns appended.

This function mirrors additive_design_matrix() in Python delicatessen, so code translated from Python can keep its shape. There is no base R equivalent of this construction, so this is the interface for it in deli as well.

Usage

additive_design_matrix(X, specifications, return_penalty = FALSE)

Arguments

Value

If return_penalty = FALSE, a numeric matrix. If TRUE, a list with elements X (the design matrix) and penalty (numeric vector).

Examples

set.seed(42)
X <- cbind(rnorm(50), rnorm(50))

# The first column stays linear; the second also gets penalized splines
specs <- list(NULL, list(knots = c(-1, 0, 1), penalty = 5))
out <- additive_design_matrix(X, specs, return_penalty = TRUE)

# Linear terms are unpenalized, spline terms carry the requested penalty
out$penalty

dim(out$X)

Description

Computes predicted survival analysis measures and point-wise confidence intervals from an accelerated failure time model for a single covariate pattern across a set of time points. The point estimates mirror aft_predictions_individual(); the variance is obtained with the delta method (see delta_method()) and Wald-type intervals are formed on the resulting standard errors. Meant to be used after fitting ee_aft() with MEstimator(), typically to draw a measure and its confidence band over time.

Usage

aft_predictions_function(
  X,
  times,
  theta,
  covariance,
  distribution,
  measure = "survival",
  alpha = 0.05,
  deriv_method = "capprox",
  dx = 1e-09
)

Arguments

Details

The survival and hazard for a covariate pattern $X$ are

$S(t) = S_{\epsilon}\left( \frac{\log(t) - X \beta^T}{\sigma} \right)$

$h(t) = (\sigma t)^{-1} h_{\epsilon}\left( \frac{\log(t) - X \beta^T}{\sigma} \right)$

where $S_{\epsilon}$ and $h_{\epsilon}$ are the error survival and hazard functions for the chosen distribution. The requested measure is derived from these through convert_survival_measures().

Value

A data frame with one row per time point and columns: time, predicted, variance, lower, upper.

Length-one times

A single time point is supported and returns a one-row data frame equal to the corresponding row of a multi-time call. This is a deliberate improvement over Python Delicatessen, whose aft_predictions_function raises on a scalar or length-one times because it takes the diagonal of a scalar delta-method covariance.

Examples

# Weibull AFT fit, then a survival curve for one covariate pattern
set.seed(1)
n <- 200
x <- rbinom(n, 1, 0.5)
Xd <- cbind(1, x)
eps <- log(-log(runif(n)))
t_event <- exp(2 + 0.5 * x + 0.8 * eps)
t_censor <- rexp(n, rate = 0.02)
t_obs <- pmin(t_event, t_censor)
delta <- as.numeric(t_event <= t_censor)

psi <- function(theta) {
  ee_aft(theta, X = Xd, time = t_obs, event = delta, distribution = "weibull")
}
m <- m_estimate(
  stacked_equations = psi,
  init = c(mean(log(t_obs)), 0, 0),
  solver = "nleqslv"
)

aft_predictions_function(
  X = matrix(c(1, 1), nrow = 1), times = c(5, 10, 20),
  theta = coef(m), covariance = vcov(m),
  distribution = "weibull", measure = "risk"
)

Description

Computes individual-level predicted survival measures from an accelerated failure time model at specified time points. Meant to be used after fitting ee_aft() with MEstimator().

Usage

aft_predictions_individual(X, times, theta, distribution, measure = "survival")

Arguments

Value

A data frame with n rows and one column per time point.

Examples

# Weibull AFT fit, then individual-level survival for the first four people
set.seed(1)
n <- 200
x <- rbinom(n, 1, 0.5)
Xd <- cbind(1, x)
eps <- log(-log(runif(n)))
t_event <- exp(2 + 0.5 * x + 0.8 * eps)
t_censor <- rexp(n, rate = 0.02)
t_obs <- pmin(t_event, t_censor)
delta <- as.numeric(t_event <= t_censor)

psi <- function(theta) {
  ee_aft(theta, X = Xd, time = t_obs, event = delta, distribution = "weibull")
}
m <- m_estimate(
  stacked_equations = psi,
  init = c(mean(log(t_obs)), 0, 0),
  solver = "nleqslv"
)

# Rows are individuals and columns are the requested times. The first two
# people share a covariate pattern, as do the third and fourth, so their
# predictions agree.
aft_predictions_individual(
  X = Xd[1:4, ], times = c(5, 10, 20),
  theta = coef(m), distribution = "weibull", measure = "survival"
)

Description

Collapses unit-level estimating function contributions into group-level contributions for clustered or grouped data. Uses an independent working correlation structure (summing within groups).

This function mirrors aggregate_efuncs() in Python delicatessen, so code translated from Python can keep its shape. There is no base R equivalent that operates on estimating function contributions, so this is the interface for them in deli as well.

Usage

aggregate_efuncs(est_funcs, group)

Arguments

Details

This function should be called inside the psi function after computing unit-level estimating equations but before returning them to MEstimator(). This changes the effective sample size used by the empirical sandwich variance estimator.

Value

A p-by-m matrix, where m is the number of unique groups. Row names are those of est_funcs, since the rows are the same parameters. Columns are ordered by the sorted unique values of group and are labeled with those values, as character. A factor group is coerced with as.vector() to its character labels before sorting, so its columns sort lexically by label rather than by factor-level order and carry those labels; a level with no observations contributes no column and so no label.

Examples

# Fifty clusters of four observations, sharing a cluster-level shift in y
set.seed(42)
n <- 200
group <- rep(1:50, each = 4)
cluster_effect <- rnorm(50, sd = 2)
y <- cluster_effect[group] + rnorm(n)

psi <- function(theta) aggregate_efuncs(ee_mean(theta, y = y), group = group)

m <- m_estimate(stacked_equations = psi, init = mean(y))

# Cluster-robust standard error, larger than the naive independence version
sqrt(diag(vcov(m)))

Description

Adverse events case study data from Bonate (2011), used to demonstrate generalized linear models (logistic, probit, complementary log-log).

Usage

bonate_adverse

Format

A data frame.

References

Bonate PL. (2011). Pharmacokinetic-Pharmacodynamic Modeling and Simulation. 2nd Ed. Springer.

Description

Survival data from Collett (2015) on the survival times of 45 women with breast cancer in Middlesex Hospital July 1987 (Table 1.2).

Usage

breast_cancer

Format

A data frame with 45 rows and 3 columns:

delta

Event indicator (1 = event, 0 = censored)

times

Observation time, in months

stain

Tumor HPA stain indicator (1 = positive, 0 = negative)

References

Collett, D. (2015). Survival analysis. In Modelling survival data in medical research. 3rd Ed. Chapman and Hall/CRC. pg 6-7

Description

Bladder cancer recurrence data from Collett (2015), used to demonstrate pooled logistic regression for survival analysis.

Usage

collett_bladder

Format

A data frame with 86 rows and 6 columns:

patient

Patient ID

time

Follow-up time, in months

delta

Event indicator

treat

Treatment group

init

Number of initial tumors

size

Size of largest initial tumor

References

Collett, D. (2015). Modelling survival data in medical research. 3rd Ed. Chapman and Hall/CRC.

Description

Compute confidence bands from theta and covariance

Usage

compute_confidence_bands(
  theta,
  covariance,
  alpha = 0.05,
  method = "supt",
  n_draws = 100000L,
  seed = NULL
)

Arguments

Value

A p-by-2 matrix with columns "lower" and "upper". Rows take their names from theta, when it has any.

Examples

fit <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                  model = "linear")

# Bands from the estimates and covariance alone, without the fitted object
compute_confidence_bands(coef(fit), covariance = vcov(fit),
                         method = "supt", seed = 1)

Description

Computes the empirical sandwich variance estimator directly from a set of estimating equations and a vector of parameter estimates. Unlike MEstimator(), this function does not solve for the parameters; it assumes that theta is already the root of the estimating equations and only assembles the sandwich covariance at that point.

Usage

compute_sandwich(
  stacked_equations,
  theta,
  deriv_method = "capprox",
  dx = 1e-09,
  allow_pinv = TRUE,
  finite_correction = NULL,
  summed_equations = NULL,
  check_summed_equations = TRUE
)

Arguments

Details

The sandwich is built from a bread matrix and a meat matrix. The bread is the negative Jacobian of the summed estimating equations, and the meat is the cross-product of the equation evaluations. Each is scaled by 1/n internally, so the returned matrix is on the asymptotic scale (see Value). The bread Jacobian is obtained either by finite differences or by forward-mode automatic differentiation.

Value

A p-by-p covariance matrix on the asymptotic scale. The bread and meat are each divided by n internally, so the returned matrix is the variance that corresponds to the standard deviation. Dividing it by the number of observations gives the standard-error-scale variance, whose square-rooted diagonal is the vector of standard errors.

A covariance matrix is the only thing this function returns. Where the bread has no inverse, and so no sandwich can be assembled from it, the call raises an error carrying the class deli_bread_not_invertible rather than returning something that has to be tested for. That covers a bread holding NA, and, under allow_pinv = FALSE, a rectangular bread and one the solve could not invert. estimate() makes the other choice from the same matrices: a fit whose bread holds NA warns and comes back with no variance, since it still carries the estimates.

References

Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference (pp. 297-337). Springer, New York, NY.

See Also

MEstimator() and GMMEstimator(), which solve for theta and report this variance internally, and delta_method() for the variance of a transformation of the parameters.

Examples

# A generic data set for estimating a mean and variance
y <- c(1, 2, 4, 1, 2, 3, 1, 5, 2)

# The mean and variance are the roots of ee_mean_variance, so they can be
# computed directly rather than solved for
theta <- c(mean(y), stats::var(y) * (length(y) - 1) / length(y))

# Wrap the built-in estimating equation as a function of theta alone
psi <- function(theta) ee_mean_variance(theta, y = y)

# compute_sandwich() returns the asymptotic-scale variance, so dividing by
# n puts it on the standard-error scale
sandwich <- compute_sandwich(psi, theta = theta) / length(y)
sandwich

# The diagonal square roots are the standard errors
sqrt(diag(sandwich))

# The bread only ever needs the sums of the estimating equations, so an
# equation whose sums have a closed form can supply them directly
summed <- function(theta) {
  c(
    sum(y) - length(y) * theta[1],
    sum((y - theta[1])^2) - length(y) * theta[2]
  )
}
compute_sandwich(psi, theta = theta, summed_equations = summed) / length(y)

Description

Computes simultaneous confidence bands that provide coverage for the entire parameter vector, adjusting for multiple comparisons. The formula is:

$\hat{\theta} \pm \hat{c}_{\alpha/2} \times \widehat{SE}(\hat{\theta})$

where $\hat{c}$ is the adjusted critical value.

Usage

confidence_bands(
  object,
  alpha = 0.05,
  method = "supt",
  n_draws = 100000L,
  seed = NULL,
  subset = NULL,
  covariance = NULL,
  ...
)

Arguments

Value

A p-by-2 matrix with columns "lower" and "upper". For a fitted estimator the rows are named for the parameters, as in confint(). For a numeric object they take their names from it, when it has any.

Examples

# Two independent samples, each contributing one mean parameter
set.seed(42)
n <- 200
y1 <- rnorm(n, 2)
y2 <- rnorm(n, 3)

psi <- function(theta) {
  rbind(y1 - theta[1], y2 - theta[2])
}

m <- m_estimate(stacked_equations = psi, init = c(0, 0))

# Simultaneous bands cover both means at once, so they are wider than the
# pointwise intervals from confint(m)
confidence_bands(m, method = "supt", seed = 1)

Description

Computes two-sided Wald-type $(1 - \alpha) \times 100\%$ confidence intervals using the point estimates and sandwich variance: $\hat{\theta} \pm c_{\alpha/2} \times \widehat{SE}(\hat{\theta})$. The critical value $c_{\alpha/2}$ comes from the standard normal distribution by default. When a finite_correction is set on the fit, it comes instead from the t-distribution with $n - p$ degrees of freedom, matching the finite-sample adjustment of the variance.

This function mirrors m.confidence_intervals() in Python delicatessen, so code translated from Python can keep its shape.

Usage

confidence_intervals(object, alpha = 0.05, ...)

Arguments

Value

A p-by-2 matrix with columns "lower" and "upper".

See Also

confint(), the standard R accessor for the same intervals. It returns identical values but is parameterized by level = 0.95 where this function takes alpha = 0.05.

Examples

psi <- function(theta) {
  y <- c(1, 2, 3, 4, 5)
  matrix(y - theta[1], nrow = 1)
}
m <- m_estimate(stacked_equations = psi, init = 0)
confidence_intervals(m)

Description

Converts a survival probability (and optionally a hazard) to other survival analysis metrics.

Usage

convert_survival_measures(survival, hazard = NULL, measure)

Arguments

Value

Numeric value or vector of the requested measure.

Examples

# Risk is the complement of survival
convert_survival_measures(0.8, measure = "risk")

# The density is the product of the hazard and the survival, so the
# "density" and "hazard" measures need the hazard as well.
convert_survival_measures(
  c(0.9, 0.8, 0.6),
  hazard = c(0.05, 0.06, 0.08),
  measure = "density"
)

Description

Data on violent crime rates by US state from the 2005 Statistical Abstract of the United States, used to demonstrate robust regression.

Usage

crime

Format

A data frame with 51 rows (50 states + DC).

References

Agresti A & Finlay B. (2009). Statistical Methods for the Social Sciences, 4th Edition. Pearson.

Description

Pharmacodynamic data from Cutler et al. (1995) on acetylcholinesterase inhibition, used to demonstrate E-max dose-response models.

Usage

cutler1995

Format

A data frame.

References

Cutler NR et al. (1995). Dose-dependent CSF acetylcholinesterase inhibition by SDZ ENA 713 in Alzheimer's disease. Acta Neurologica Scandinavica, 91(5), 347-351.

Description

Evaluates the digamma function, the first derivative of the log gamma function.

This is the equivalent of base::digamma(), and unlike most of the deli utilities it is not required anywhere: digamma() belongs to the Math group generic and has a tangent rule, so it differentiates under deriv_method = "exact" just as deli_digamma() does. The two agree value for value, including at the poles and on missing input. The only difference is that deli_digamma() returns NaN at the non-positive integers, where the digamma function has poles, without base R's ⁠NaNs produced⁠ warning. deli_digamma() exists so that code translated from Python delicatessen, where the function is called digamma(), can keep its shape.

Usage

deli_digamma(z)

Arguments

Value

Numeric digamma values.

See Also

base::digamma(), which is usable in the same positions, and deli_polygamma() for the higher-order derivatives. inverse_logit() has the full list of deli utilities and their base R counterparts, including the five that base R cannot stand in for under exact differentiation.

Examples

deli_digamma(1)
deli_digamma(c(0.5, 1, 2))

# The same values as the base R counterpart
all.equal(deli_digamma(c(0.5, 1, 2)), digamma(c(0.5, 1, 2)))

# NaN at the poles, where base R warns, and missing values pass through
deli_digamma(c(-1, 0, 2.5, NA))

Description

Evaluates the nth derivative of the digamma function.

This is the equivalent of base::psigamma() and returns identical values for numeric input, except that it carries derivatives. Exact differentiation (deriv_method = "exact") propagates a tangent alongside each value, and deli_polygamma() recognizes a tangent-carrying argument and applies the analytic rule itself, the polygamma function of order n + 1. psigamma() hands its argument to compiled code without dispatching, and errors on such an argument. Use deli_polygamma() inside estimating equations and inside transforms passed to delta_method(), and psigamma() for ordinary numeric work.

The two take their arguments in opposite orders. deli_polygamma(n, x) takes the derivative order first, matching Python's scipy.special.polygamma(n, x); psigamma(x, deriv = n) takes it second. deli_polygamma(1, 3) is psigamma(3, deriv = 1), not psigamma(1, deriv = 3). A positional substitution between the two computes a different quantity and raises no error.

Usage

deli_polygamma(n, x)

Arguments

Value

Numeric polygamma values.

Exact differentiation

deriv_method = "exact" is forward-mode automatic differentiation: it replaces each value with an object carrying both the value and its derivative. deli supports those objects through S3 methods: the Ops, Math, and Summary group generics, plus non-group methods such as [, %*%, t(), c(), and mean(). standard_normal_cdf(), standard_normal_pdf(), deli_polygamma(), and deli_digamma() recognize a tangent-carrying argument themselves and apply their own analytic rule. Support within the group generics is partial; see vignette("getting-started") for the operations deli differentiates.

plogis(), qlogis(), pnorm(), dnorm(), and psigamma() take none of these paths. Each hands its argument straight to compiled code through .Call() or .Internal() without dispatching, so the tangent reaches C code that requires a plain number. deli catches the resulting failure and raises its own error, naming the function that stopped the computation and the deli function to write in its place. The same applies to every other distribution function in stats, such as qnorm(), which is named in the error even though deli exports no counterpart for it.

Each deli utility below returns the same values as its base R counterpart for numeric input. What separates them is whether the counterpart survives exact mode:

For the first five rows, use the deli function inside estimating equations and inside transforms passed to delta_method(), and the base R function everywhere else: simulating data, post-fit display, plain numeric work. For the last two rows the base R function is usable everywhere, because digamma() is a Math group member with a tangent rule and identity() passes its argument through untouched.

The polygamma row is the one place where the arguments do not line up. deli_polygamma(n, x) takes the derivative order first and psigamma(x, deriv = n) takes it second, so a positional substitution between the two computes a different quantity and raises no error.

See Also

base::psigamma(), the base R equivalent for ordinary numeric work, which takes its arguments in the other order, and deli_digamma() for the order-zero case.

Examples

deli_polygamma(0, 1)
deli_polygamma(1, c(1, 2, 5))

# The same values as the base R counterpart, with the arguments swapped
all.equal(deli_polygamma(1, c(1, 2, 5)), psigamma(c(1, 2, 5), deriv = 1))

# A gamma GLM, whose last parameter is the log of the shape
m <- m_estimate(
  mpg ~ wt,
  data = mtcars,
  .ee = ee_glm,
  distribution = "gamma",
  link = "log",
  init = c(3, 0, 0)
)

# Variance of the trigamma function at the estimated shape. Writing
# `psigamma(exp(theta[3]), deriv = 1)` here instead would error, because
# exact differentiation hands the transform a tangent-carrying argument.
delta_method(
  m,
  transform = function(theta) deli_polygamma(1, exp(theta[3])),
  deriv_method = "exact"
)

Description

Generates polynomial spline terms for a numeric vector at pre-specified knot locations. Default is restricted (natural) cubic splines, but unrestricted splines with different polynomial terms can also be generated.

This function mirrors spline() in Python delicatessen, so code translated from Python can keep its shape. Despite the name it has nothing to do with stats::spline(), which interpolates a curve through data points and returns the interpolated values. deli_spline() builds a truncated power basis, a matrix of design columns to be used as covariates in a regression. The deli_ prefix is there so that the two names do not collide.

No base R function returns these truncated power columns, but the splines package, which installs with every copy of R, builds the same spline spaces in a B-spline parameterization. splines::bs() is the counterpart to restricted = FALSE: for knots k, cbind(1, x, x^2, x^3, deli_spline(x, k, restricted = FALSE)) and cbind(1, splines::bs(x, knots = k)) span the same space and give identical fitted values. splines::ns() is the nearest counterpart to restricted = TRUE but not an exact one, since it also constrains the fit to be linear beyond the boundary knots and so spans a subspace of the columns here. deli_spline() is what deli offers for parity with Python delicatessen, for the truncated power form itself, and because additive_design_matrix() builds on it.

Usage

deli_spline(x, knots, power = 3, restricted = TRUE, normalized = FALSE)

Arguments

Details

Unrestricted splines for knot $k$ are:

$s_k(X) = I(X > k) (X - k)^a$

Restricted (natural) splines subtract the last spline term:

$r_k(X) = s_k(X) - s_K(X)$

where $K$ is the largest knot. Restricted splines return one fewer column than the number of knots.

Value

A matrix with length(x) rows. Number of columns is length(knots) for unrestricted or length(knots) - 1 for restricted.

Examples

# Restricted quadratic splines at four knots, so three basis columns
s <- deli_spline(1:59, knots = c(10, 20, 30, 40), power = 2)
dim(s)

# Each term stays at zero below its knot and grows above it
s[c(5, 15, 25, 35, 45, 55), ]

Description

An augment() method for MEstimator and GMMEstimator objects fitted through the formula interface. It returns the model frame the fit was built from, or newdata when supplied, with the fitted values, their standard errors, a Wald confidence interval, and the residuals as columns beside it.

Usage

## S3 method for class ''deli::deli_estimator''
augment(
  x,
  newdata = NULL,
  type.predict = c("link", "response"),
  conf.level = 0.95,
  ...
)

Arguments

Details

The columns added are .fitted, .se.fit, .lower, .upper, and, when newdata is not supplied, .resid. They are exactly what predict() and residuals() return for the same fit, so .fitted is predict(), .lower and .upper are predict(interval = "confidence"), and .resid is residuals().

.fitted is on the link scale by default, matching both predict() and broom::augment() on a glm, and type.predict = "response" puts it and its interval on the scale of the response. .resid is the response residual either way, since a residual measured against a linear predictor would be a response minus a quantity the response is not measured in.

Rows the fit dropped for missing data are not reported, so the result has one row per nobs() and its row names are those of the retained rows. A newdata row with a missing value is kept, with NA in the added columns, so that the result lines up with the rows handed in.

augment() covers the estimating equations whose linear predictor predict() forms, and refuses the same fits with the same reasons; see deli-predict. It has no counterpart to that method's times argument: a survival measure is one value per row of the data and time rather than one per row, so the predictions do not go beside the data as columns.

Value

A data frame: the model frame, or newdata, followed by the columns .fitted, .se.fit, .lower, .upper, and .resid. The last is absent when newdata is supplied.

See Also

deli-predict for the predictions themselves and the equations they are available for, deli-tidiers for the parameter-level and model-level summaries, and reexports for the generic itself, which deli re-exports so that augment() resolves with deli alone attached.

Examples

fit <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                  model = "linear")

head(augment(fit))

# New covariate patterns come back with no residual column, since they carry
# no response to residualize against.
augment(fit, newdata = data.frame(wt = c(2, 3, 4), hp = 110))

Description

Where a caller might reasonably want to answer an error or a warning rather than let it stop the work, deli attaches a condition class to it, so that base::tryCatch(), base::withCallingHandlers(), or rlang::try_fetch() can match on the class rather than on the wording of a message. This page lists those classes and says what each of them reports on.

Details

Every class below belongs to an error unless the entry says otherwise. A condition carrying two of them leads with the narrower one, so a handler for either the family or the particular fault catches it.

A condition a user's own estimating function raises reaches the caller with whatever class that function gave it. Nothing here renames one.

What an estimating function returned

• deli_psi_return_error: the estimating function returned something no sandwich can be built from, either at the starting values or at the point handed to compute_sandwich(). Every refusal of a return carries it, so recognizing one is a match on this rather than on four messages.

• deli_psi_shape_error: the number of estimating equations returned cannot be solved against the number of parameters. Leads deli_psi_return_error.

• deli_psi_list_unsupported: the estimating function returned a list of one element per equation. compute_bread() and compute_sandwich() take that shape, and no fit can solve it, so estimate() refuses it rather than failing on the summed equations. Leads deli_psi_return_error.

The estimating equations deli supplies

• deli_mlogit_binary_outcome: the indicator matrix handed to ee_mlogit() holds two columns. Two categories leave a single non-reference residual whose reference residual is its exact negative, so the equations cancel to zero at every value of the parameters and a fit of them reports the values it started from with no variance. A two-level outcome is logistic regression, which ee_glm() fits with distribution = "binomial".

The formula interface

• deli_formula_ee_lookup_error: .ee is neither a function nor the name of one. Either it arrived as some other type entirely, or it arrived as a character vector naming no one function, because no function of that name was found or because the vector is not of length one. A handler that means to read the name it failed to resolve should confirm it has a character vector first.

• deli_formula_ee_signature_error: the arguments of .ee leave something the interface fills nowhere to go: theta, the model matrix it passes as X, the response it passes by position, or the offset an offset() term in the formula supplies. It also covers an equation that writes formals after its ..., which can be filled by name and by nothing else, so the positional response would fall into the dots and leave them unfilled.

• deli_formula_ee_argument_error: something in the ... the caller wrote cannot be forwarded to .ee. On its own it is a name matching no argument of .ee exactly, which is what keeps an abbreviation from reaching an argument the caller did not mean.

• deli_formula_ee_unnamed_argument: an argument was forwarded with no name, and so by position, into whichever argument the response did not take. Leads deli_formula_ee_argument_error.

• deli_formula_ee_reserved_argument: a name the caller supplied is one the interface fills itself, either theta, X, or the argument the response is passed to. Leads deli_formula_ee_argument_error.

• deli_formula_auto_init_error: the init the formula interface generated itself, one zero per model-matrix column, does not fit the estimating equation.

Solving

• deli_solver_not_converged, a warning: the solver stopped without solving the estimating equations, either because it reported a failure of its own or because the point it returned does not solve them.

• deli_nested_solver_error: a rootSolve::multiroot() solve was asked for while another one was running, which it cannot be. An error rather than a warning, because the solve cannot be attempted at all.

• deli_gmm_moments_rejected, a warning: the Hansen J-statistic of an over-identified GMM fit stands so far out against its reference distribution as to all but rule the fit out, which usually means the moment conditions cannot all hold at one value of the parameters.

• deli_gmm_moments_dependent, a warning: the moment conditions are linearly dependent at the estimated values, so the two-step weight matrix is a pseudo-inverse and what comes back is the fit of the independent conditions alone. It also carries the rarer covariance that has no inverse for want of conditioning rather than of rank, which the message tells apart.

• deli_meat_not_invertible: the covariance of the moment conditions cannot be inverted, so the two-step GMM weight update has no weight matrix to take, and allow_pinv = FALSE refused the pseudo-inverse. The counterpart of deli_bread_not_invertible for the other half of the sandwich.

The sandwich variance

• deli_bread_na, a warning: the bread matrix holds NA, so no inverse of it exists. A fit records no variance and carries on; compute_sandwich() has nothing but a matrix to return and converts it into the error below.

• deli_bread_not_invertible: the bread cannot be inverted and allow_pinv = FALSE refused the pseudo-inverse, or it is rectangular and has no inverse at any rank, or it holds NA.

• deli_summed_equations_error: the summed_equations handed to compute_sandwich() or compute_bread(), or carried as a property of MEstimator() or GMMEstimator(), is not one the bread can be differentiated from. It is neither NULL nor a function, or its return at the point in hand is not numeric, or that return does not hold one value per estimating equation. Every refusal of the argument carries it, as does a check_summed_equations that is not a single TRUE or FALSE.

• deli_summed_equations_disagree: the reduction is a function of the right shape and does not sum the stacked_equations it is paired with, either by differing from their row sums or by returning a value that is not a number where they sum to one. The bread is differentiated from the one and the meat is built from the other, so a sandwich assembled from two systems would carry the shape of a covariance matrix without the meaning, and a fit would report the estimates of the other system as well. Leads deli_summed_equations_error.

Differentiation

• deli_exact_unsupported_function: under deriv_method = "exact", a function reached with a tangent-carrying argument has no rule, either because it hands its argument to compiled code without dispatching or because deli declines to differentiate it. auto_differentiation() names the replacements to use instead.

• deli_exact_unsupported_shape: the estimating equations arrived in a container the summing step has no rule for, either because the tangents survived in a shape it cannot reduce or because a per-equation list does not hold one element per parameter, each of one length. Nothing was lost; the shape is the problem.

• deli_exact_tangent_lost: derivative information is gone or would be, either because a tangent-carrying value was asked to become a plain double or because a result arrived with no tangent and evidence that one had been dropped on the way.

• deli_finite_difference_lost, a warning: a finite-difference step changed the function by less than the floating-point spacing of its values, so an entry of the Jacobian carries none of the digits of the values it came from.

Parameter names

• deli_param_name_collision: filling the unnamed entries of a partly named parameter vector with positional theta_1, theta_2, ... labels would repeat a label another parameter in the same vector already carries.

Description

Methods for base::print() and base::summary() so that a fitted estimator shows its coefficients at the console and reports its parameter-level inference as a table.

Usage

## S3 method for class ''deli::deli_estimator''
print(x, ..., subset = NULL)

## S3 method for class ''deli::deli_estimator''
summary(object, alpha = 0.05, subset = NULL, ...)

Arguments

Details

print() reports the parameter and observation counts followed by the estimates, rounded to four decimal places. An estimator that has not been through estimate() has no estimates to show, so it reports its parameter count and says so.

summary() collects the estimates, their standard errors, Z-scores, confidence limits, P-values, and S-values into one object, which prints as a table with one row per parameter. The values are those confint(), z_scores(), p_values(), and s_values() return individually, so alpha here means what it means there: 0.05 gives 95% limits.

An over-identified GMMEstimator fit reports Hansen's J-statistic above the table, with its degrees of freedom and its P-value, because it judges the fit as a whole rather than any one parameter. See GMMEstimator() for what it means and where its reference distribution holds. No other fit has one, so no other output carries the line.

The S column reads Inf when a P-value underflows to exactly zero, for the reason s_values() gives. The P column reports that same underflow as ⁠<2e-16⁠: base::format.pval() stops printing digits below the eps it is given, and the table gives it 2.2e-16. tidy() returns the literal 0 instead.

subset restricts which parameters are displayed and nothing else. The reported parameter count and every reported value are computed from the whole fit, so displaying a subset of a stacked estimator is a way to read the parameters of interest without the nuisance parameters, not a way to refit without them. Row labels keep the names the full fit gave them.

Value

• print(): its input, invisibly.

• summary(): an object holding the estimates, standard errors, Z-scores, confidence limits, P-values, and S-values, which prints as a table.

See Also

deli-generics for the accessors that return these quantities as plain vectors and matrices, and deli-tidiers for the same results as a data frame.

Examples

fit <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                  model = "linear")

fit

# `subset` is an argument of the method rather than of the fit, so showing it
# here means calling the generic by name.
print(fit, subset = 2:3)

summary(fit)

# `summary()` takes `subset` as well, alongside `alpha` for the width of the
# reported limits.
summary(fit, subset = 2:3, alpha = 0.1)

Description

Methods for base R generics stats::coef(), stats::vcov(), stats::confint(), stats::nobs(), stats::df.residual(), stats::fitted(), stats::residuals(), stats::weights(), stats::model.frame(), stats::model.matrix(), stats::formula(), stats::terms(), stats::sigma(), stats::logLik(), and stats::deviance() so that deli estimator objects interoperate with the broader R modeling ecosystem.

Usage

## S3 method for class ''deli::deli_estimator''
coef(object, ...)

## S3 method for class ''deli::deli_estimator''
vcov(object, ...)

## S3 method for class ''deli::deli_estimator''
confint(object, parm, level = 0.95, ...)

## S3 method for class ''deli::deli_estimator''
nobs(object, ...)

## S3 method for class ''deli::deli_estimator''
df.residual(object, ...)

## S3 method for class ''deli::deli_estimator''
fitted(object, ...)

## S3 method for class ''deli::deli_estimator''
residuals(object, type = "response", ...)

## S3 method for class ''deli::deli_estimator''
weights(object, ...)

## S3 method for class ''deli::deli_estimator''
model.frame(
  formula,
  data = NULL,
  subset = NULL,
  na.action,
  drop.unused.levels = FALSE,
  xlev = NULL,
  ...
)

## S3 method for class ''deli::deli_estimator''
model.matrix(object, data = NULL, ...)

## S3 method for class ''deli::deli_estimator''
formula(x, ...)

## S3 method for class ''deli::deli_estimator''
terms(x, ...)

## S3 method for class ''deli::deli_estimator''
sigma(object, ...)

## S3 method for class ''deli::deli_estimator''
logLik(object, ...)

## S3 method for class ''deli::deli_estimator''
deviance(object, ...)

Arguments

Details

fitted(), residuals(), weights(), model.frame(), model.matrix(), formula(), and terms() report on the model a fit was specified as, which only the formula interface records, so each of them is an error for a fit built from a stacked_equations function. fitted() and residuals() are also an error for a formula fit of an estimating equation predict() does not support, since a fitted value is a prediction; see deli-predict for the equations it covers.

fitted() is the conditional mean of the response, so it is on the response scale rather than the link scale, matching stats::fitted() on a glm object. residuals() is the response minus that mean, the residual a GLM calls a response residual, and the only type it takes. Deviance and Pearson residuals are not offered: the first needs a likelihood and the second a variance function, and an M-estimator need not state either, while the response residual is defined the same way for every equation predict() supports. Under a non-identity link it is heteroscedastic by construction.

model.frame() returns the frame the fit was solved on, or builds the frame of data when one is supplied, through the terms and factor levels the fit recorded. A data-dependent term such as poly(x, 2) is therefore evaluated at its fitted basis and a factor gets its fitted levels, as they are under predict(). data has to carry the response, since a model frame holds every variable the formula names; covariate values that carry no response are what predict() and deli-augment take a newdata for.

model.matrix() returns the design the fit was solved on, or the design of data when one is supplied, coded with the contrasts and factor levels the fit recorded rather than with whatever getOption("contrasts") says when the call is made. A design matrix is a property of the fit, so a fit made under one setting of that option answers with the coding it was solved on under any other, which is what predict() does as well. A design names no response, so data needs only the predictors.

df.residual() is the number of observations less the number of parameters. Both counts belong to the fit rather than to a specification, so it answers for a fit built from a stacked_equations function as readily as for a formula one. weights() returns the observation weights the fit was solved with, which reach a fit only through the formula interface. A formula fit specified without weights returns NULL, as stats::weights() does for an unweighted stats::lm() fit: a vector of ones would be indistinguishable from a fit weighted by ones.

sigma(), logLik(), and deviance() are errors for every fit, and say why. An M-estimator is defined by its estimating equations, which need come from no likelihood at all, so there is no log-likelihood to return and no deviance to take from one; stats::AIC() and stats::BIC() reach a fit through logLik() and are refused with it. A residual standard deviation belongs to the model an equation states rather than to the solve: a linear equation has one, a logistic equation has none, and nothing in the estimates says which was solved. Reporting an error is the only answer that does not invent a quantity the fit never had.

Value

• coef(): Named numeric vector of parameter estimates.

• vcov(): Named variance-covariance matrix.

• confint(): Matrix with columns "lower" and "upper".

• nobs(): Integer number of observations.

• df.residual(): Integer residual degrees of freedom, the number of observations less the number of parameters.

• fitted(): Named numeric vector of fitted values on the response scale, one per observation the fit was solved on.

• residuals(): Named numeric vector of response residuals.

• weights(): The observation weights the fit was solved with, as they were recorded, or NULL for a formula fit specified without any. Weights that vary over time are reported as the matrix they were supplied as rather than as a vector: ee_plogit() takes an n-by-K matrix carrying one column per time interval, and an observation weighted differently in each interval has no one weight for a vector to hold.

• model.frame(): The model frame the fit was built from, with the rows dropped for missing data already removed, or the model frame of data when one is supplied.

• model.matrix(): The design matrix the fit was solved on, or the design matrix of data when one is supplied, coded with the contrasts and factor levels the fit recorded.

• formula(): The model formula.

• terms(): The terms object of the model frame, carrying the response index, any offset, and the predvars of a data-dependent term.

• sigma(), logLik(), deviance(): Nothing. Each raises an error saying that an M-estimator states no likelihood and records no residual scale.

See Also

deli-predict for predictions at new covariate values, and deli-augment for the fitted values, intervals, and residuals as columns beside the data.

Examples

fit <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                  model = "linear")

coef(fit)

vcov(fit)

confint(fit)

nobs(fit)

df.residual(fit)

head(fitted(fit))

head(residuals(fit))

formula(fit)

# Weights reach a fit through the formula interface, and `weights()` reports
# the vector the fit was solved with.
weighted <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                       model = "linear", weights = 1 / hp)

head(weights(weighted))

Description

A method for stats::predict() that predicts from a fit made through the formula interface, with sandwich standard errors and Wald confidence intervals. It forms the linear predictor of a regression fit on either the link or the response scale, and evaluates a survival measure at a set of times for a fit of ee_aft() or ee_plogit().

Usage

## S3 method for class ''deli::deli_estimator''
predict(
  object,
  newdata = NULL,
  type = c("link", "response"),
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  times = NULL,
  measure = "survival",
  deriv_method = "capprox",
  dx = 1e-09,
  ...
)

Arguments

Value

Without times, a named numeric vector of predictions, one per row of newdata or of the fitted design, and with interval = "confidence" a matrix with columns "fit", "lwr", and "upr".

With times, a data frame with one row per row of the design and time, every time for the first row before any time for the second, and columns .row, the row label of the design, time, and fit. With interval = "confidence" it also has "lwr" and "upr".

With se.fit = TRUE, either shape becomes a list whose fit element is whichever of the two the other arguments call for and whose se.fit element is a numeric vector of standard errors.

Two surfaces

times chooses what is predicted. Without it, predict() returns the linear predictor on the scale type names. With it, predict() returns the survival measure measure names at each of the times, one prediction per row of the design at each time.

The two cover disjoint sets of fits. An equation reaches the first surface by having a linear predictor that is a conditional mean of the response, and the two equations on the second surface are there because theirs is not: ee_aft() puts its linear predictor on the log-time scale, and ee_plogit() has one linear predictor per person and time interval rather than one per person. Supplying type beside times, or measure, deriv_method, or dx without times, is an error rather than an argument silently ignored, since no fit takes both sets.

The linear predictor

The design for newdata is rebuilt through the terms, factor levels, and contrasts the fit recorded, so a factor whose new values cover only some of the fitted levels still produces the fitted set of columns, and a data-dependent term such as poly(x, 2) or scale(x) is evaluated with the coefficients it was fitted with rather than refitted to newdata. A factor level the fit never saw, or a predictor newdata does not carry, is an error.

The standard error is the delta-method standard error of the linear predictor, $\sqrt{\mathrm{diag}(X \hat{V} X^{T})}$, formed from the coefficient block of the sandwich variance. On the response scale it is scaled by the derivative of the inverse link, which is exact because the inverse link is applied elementwise. The intervals are Wald intervals on the scale asked for, so they are symmetric about fit on that scale rather than transformed from the link scale.

predict() needs to know which parameters are coefficients on the design and what takes the linear predictor to the mean of the response. It supports ee_regression(), ee_glm(), ee_robust_regression(), ee_beta_regression(), and the five penalized regressions ee_bridge_regression(), ee_ridge_regression(), ee_lasso_regression(), ee_dlasso_regression(), and ee_elasticnet_regression(), whose parameters are one coefficient per design column followed by at most one parameter of the outcome distribution. Any other estimating equation, and any fit built from a stacked_equations function, is an error naming the reason; regression_predictions() takes a design, estimates, and a covariance matrix directly and can be used wherever this method declines.

Survival measures at a set of times

times is supported for a fit of ee_aft() or of ee_plogit(), and measure names any of the measures convert_survival_measures() defines. The point estimates are those of aft_predictions_individual() and plogit_predict() for the same fit.

ee_survival_model() has no surface here because the formula interface cannot drive it: it takes no design matrix, and the interface always passes the one it built. Predict from such a fit with survival_predictions().

Each covariate pattern has its own variance, so each row gets its own interval. The variance of the measure for a row at a time is the delta-method variance $G \hat{V} G^{T}$ of that one prediction, where $G$ is its row of the Jacobian of the whole grid with respect to the parameters, built by deriv_method. A pooled logistic fit is predicted through plogit_predict(), whose matrix products cannot carry tangents, so deriv_method = "exact" is available for an AFT fit alone. Asking for it on a pooled logistic fit is an error whether or not a standard error is wanted, since the Jacobian that could not be built is only built when one is.

See Also

deli-augment, which returns the linear-predictor predictions as columns beside the data; regression_predictions(), which computes the same quantities from a design matrix, a vector of estimates, and a covariance matrix, for fits this method does not cover; and aft_predictions_individual(), aft_predictions_function(), plogit_predict(), and survival_predictions(), which compute the survival measures from estimates directly.

Examples

fit <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                  model = "linear")

head(predict(fit))

# New covariate patterns, with sandwich standard errors.
at <- data.frame(wt = c(2, 3, 4), hp = 110)

predict(fit, newdata = at, se.fit = TRUE)

predict(fit, newdata = at, interval = "confidence")

# A logistic fit predicts log-odds by default and probabilities on the
# response scale.
vs_fit <- m_estimate(vs ~ wt, data = mtcars, .ee = ee_regression,
                     model = "logistic")

predict(vs_fit, newdata = data.frame(wt = c(2, 3)), type = "response")

# A pooled logistic fit has no linear predictor to put on a scale, and
# predicts a survival measure at a set of times instead. Every row of the
# design gets its own interval.
bladder <- collett_bladder
bladder$novel <- bladder$treat - 1
k <- length(unique(bladder$time[bladder$delta == 1]))

plogit_fit <- m_estimate(
  time ~ novel + init + size - 1, data = bladder, .ee = ee_plogit,
  event = delta, init = c(rep(0, 3), -4, rep(0, k - 1))
)

head(predict(plogit_fit, times = c(12, 24), interval = "confidence"))

Description

tidy() and glance() methods for MEstimator and GMMEstimator objects. These allow deli results to flow into tidyverse pipelines.

Usage

## S3 method for class ''deli::deli_estimator''
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class ''deli::deli_estimator''
glance(x, ...)

Arguments

Value

• tidy(): A data.frame with columns term, estimate, std.error, statistic, p.value, s.value. If conf.int = TRUE, also includes conf.low and conf.high. A p.value that underflows to exactly zero is reported as 0 alongside an infinite s.value; see s_values().

• glance(): A single-row data.frame with model-level summaries: nobs, npar, estimator, finite_correction, and the Hansen J-statistic of an over-identified GMM fit in j_statistic, j_df and j_p_value. The three J columns are present on every fit and hold the typed missing value of their own type where there is no such statistic, which is every M-estimation fit and every just-identified or subset GMM fit; see GMMEstimator() for what the statistic reads and where it is left unset.

See Also

deli-augment, the third broom generic, which returns the observation-level fitted values, intervals, and residuals, and reexports for the generics themselves, which deli re-exports so that tidy() and glance() resolve with deli alone attached.

Examples

fit <- m_estimate(mpg ~ wt + hp, data = mtcars, .ee = ee_regression,
                  model = "linear")

tidy(fit, conf.int = TRUE)

glance(fit)

Description

Computes the variance-covariance matrix for a transformation of parameters using the Delta Method:

$Var[g(\theta)] \approx G \Sigma G^T$

where $G$ is the Jacobian of $g$ and $\Sigma$ is the covariance matrix of $\theta$.

Usage

delta_method(
  object,
  transform,
  covariance = NULL,
  deriv_method = "capprox",
  dx = 1e-09,
  ...
)

Arguments

Value

A covariance matrix for transform(theta).

Examples

fit <- m_estimate(vs ~ mpg, data = mtcars, .ee = ee_regression,
                  model = "logistic")

# Variance of the odds ratio for mpg, exponentiating the log-odds coefficient
delta_method(fit, transform = function(theta) exp(theta[2]))

# The same variance from the estimates and covariance alone
delta_method(coef(fit), transform = function(theta) exp(theta[2]),
             covariance = vcov(fit))

Description

Estimates the causal effect using two-stage least squares for instrumental variable analysis. The first stage regresses the treatment on the instruments (and exogenous variables), and the second stage regresses the outcome on predicted treatment (and exogenous variables).

Usage

ee_2sls(theta, y, A, Z, W = NULL, weights = NULL)

Arguments

Value

A ⁠(1+b+2c)⁠-by-n matrix of estimating equation contributions. The second-stage rows are named stage2_A for the fitted treatment and stage2_W_1 through stage2_W_c; the first-stage rows are named stage1_Z_1 through stage1_Z_b and stage1_W_1 through stage1_W_c. The columns of W appear in both stages, so each name carries the stage it belongs to.

Examples

# A continuous treatment confounded by an unmeasured U, a strong instrument
# Z, and one measured exogenous covariate. The true causal effect is 2.
set.seed(42)
n <- 500
W1 <- rnorm(n)
Z <- cbind(rbinom(n, 1, 0.5))
U <- rnorm(n)
A <- 1.5 * Z[, 1] + 0.3 * W1 + U + rnorm(n, sd = 0.5)
Y <- 2 * A - U + 0.5 * W1 + rnorm(n)
W <- cbind(1, W1)

# At a first stage of all zeros the fitted treatment is identically zero, so
# the second stage degenerates: its leading design column vanishes and the
# coefficient on it, the causal effect, has nothing left to be estimated
# from. Seed the starting values with the ordinary least squares fit of each
# stage instead: the first stage regresses A on the instrument and the
# exogenous covariates, the second regresses Y on the fitted A and the same
# covariates.
alpha_init <- as.numeric(coef(lm(A ~ cbind(Z, W) - 1)))
a_hat <- as.numeric(cbind(Z, W) %*% alpha_init)
beta_init <- as.numeric(coef(lm(Y ~ cbind(a_hat, W) - 1)))

psi <- function(theta) ee_2sls(theta, y = Y, A = A, Z = Z, W = W)

# theta holds the three second-stage coefficients followed by the three
# first-stage coefficients.
m <- m_estimate(
  stacked_equations = psi,
  init = c(beta_init, alpha_init)
)

# theta_1 is the coefficient on the fitted treatment, the causal effect.
coef(m)

Description

Generalized Additive Model via L2-penalized splines. Internally expands X using additive_design_matrix() and delegates to ee_bridge_regression() with gamma = 2 (ridge penalty). The penalty only applies to the spline basis terms, not to the original linear terms.

Usage

ee_additive_regression(
  theta,
  X,
  y,
  specifications,
  model,
  weights = NULL,
  offset = NULL
)

Arguments

Value

A p-by-n matrix, where p is the number of columns in the expanded additive design matrix.

Examples

set.seed(42)
n <- 200
x <- runif(n, -3, 3)
y <- sin(x) + rnorm(n, sd = 0.3)
X <- cbind(1, x)

# No spline on the intercept column, a penalized spline on x.
specs <- list(NULL, list(knots = c(-2, -1, 0, 1, 2), penalty = 5))

psi <- function(theta) {
  ee_additive_regression(
    theta,
    X = X,
    y = y,
    specifications = specs,
    model = "linear"
  )
}

# One parameter per column of the expanded design matrix.
m <- m_estimate(
  stacked_equations = psi,
  init = rep(0, ncol(additive_design_matrix(X, specs)))
)
coef(m)

Description

Returns a p-by-n matrix of estimating equation contributions for accelerated failure time (AFT) models. Supports Weibull, exponential, log-logistic, and log-normal distributions.

Usage

ee_aft(theta, X, time, event, distribution, weights = NULL, offset = NULL)

Arguments

Details

The AFT model uses the parameterization where $Z_i = (\log(t_i) - X_i \beta) / \sigma$ and the last element of theta is $\log(1/\sigma)$.

Value

A p-by-n matrix where p is the number of parameters. For the exponential distribution, whose parameters are all regression coefficients, the rows are unnamed. For every other distribution the regression rows are named X_1 through X_b for the columns of X and the final row is named log_inv_scale.

Examples

# A Weibull AFT model for times generated from one binary covariate, with
# some of the times right censored by an independent exponential.
set.seed(1)
n <- 200
x <- rbinom(n, 1, 0.5)
Xd <- cbind(1, x)
eps <- log(-log(runif(n)))
t_event <- exp(2 + 0.5 * x + 0.8 * eps)
t_censor <- rexp(n, rate = 0.02)
t_obs <- pmin(t_event, t_censor)
delta <- as.numeric(t_event <= t_censor)

psi <- function(theta) {
  ee_aft(theta, X = Xd, time = t_obs, event = delta, distribution = "weibull")
}

# The default rootSolve solver does not converge here, so nleqslv is used.
# The last parameter is log(1/sigma), not a regression coefficient.
m <- m_estimate(
  stacked_equations = psi,
  init = c(mean(log(t_obs)), 0, 0),
  solver = "nleqslv"
)
coef(m)

Description

Estimates the average causal effect using AIPW, which combines a propensity score model and an outcome model for doubly-robust estimation.

Usage

ee_aipw(theta, y, A, W, X, X1, X0, truncate = NULL, force_continuous = FALSE)

Arguments

Value

A (3+b+c)-by-n matrix of estimating equation contributions, with the first three rows named ACE, E[Y^1], and E[Y^0], the propensity score rows named W_1 through W_b, and the outcome model rows named X_1 through X_c.

Examples

# A binary treatment, two confounders, and a continuous outcome whose true
# average causal effect is 1.5.
set.seed(42)
n <- 1000
W1 <- rnorm(n)
W2 <- rbinom(n, 1, 0.4)
A <- rbinom(n, 1, inverse_logit(-0.5 + 0.5 * W1 + 0.3 * W2))
Y <- 2 + 1.5 * A + W1 - 0.5 * W2 + rnorm(n)

W_ps <- cbind(1, W1, W2) # Propensity score design matrix
X <- cbind(1, A, W1, W2) # Outcome model, observed design matrix
X1 <- cbind(1, 1, W1, W2) # Outcome model, everyone treated
X0 <- cbind(1, 0, W1, W2) # Outcome model, everyone untreated

psi <- function(theta) {
  ee_aipw(theta, y = Y, A = A, W = W_ps, X = X, X1 = X1, X0 = X0)
}

# theta holds the average causal effect, the mean under treatment, and the
# mean under no treatment, followed by the three propensity score
# coefficients and the four outcome model coefficients.
m <- m_estimate(stacked_equations = psi, init = rep(0, 10))
summary(m, subset = 1:3)

Description

Beta regression for outcomes in (0, 1) using mean-precision parameterization. The last element of theta is log(phi), the log precision parameter.

Usage

ee_beta_regression(theta, X, y, weights = NULL, offset = NULL)

Arguments

Value

A (b+1)-by-n matrix. The rows are named X_1 through X_b for the columns of X, and the final row is named log_phi.

Examples

set.seed(42)
n <- 50
W <- rnorm(n)
mu <- 1 / (1 + exp(-(0.5 + 0.3 * W)))
y <- rbeta(n, shape1 = mu * 10, shape2 = (1 - mu) * 10)
X <- cbind(1, W)

psi <- function(theta) ee_beta_regression(theta, X = X, y = y)

# The last parameter is log(phi), started here at a precision of 10. The
# default solver diverges on this equation, so use nleqslv.
m <- m_estimate(
  stacked_equations = psi,
  init = c(0, 0, log(10)),
  solver = "nleqslv"
)
coef(m)

Description

Returns a p-by-n matrix for bridge penalized regression. Bridge is the general case: ridge is gamma = 2, LASSO approximation is gamma = 1 + epsilon.

Usage

ee_bridge_regression(
  theta,
  X,
  y,
  model,
  penalty,
  gamma,
  weights = NULL,
  center = 0,
  offset = NULL
)

Arguments

Value

A p-by-n matrix.

Examples

# A penalty vector gives one value per column of the design matrix. A scalar
# penalty would shrink the intercept along with the slopes. The estimating
# equation carries the penalty's derivative, which varies like
# |theta|^(gamma - 1) near the penalty center. That derivative is itself
# differentiable only once gamma reaches 2, so gamma = 2.3 issues no
# warning while a value below 2 would.
fit <- m_estimate(
  mpg ~ wt + hp,
  data = mtcars,
  .ee = ee_bridge_regression,
  model = "linear",
  penalty = c(0, 5, 5),
  gamma = 2.3
)
coef(fit)

Description

Uses a smooth approximation to the L1 penalty based on the standard normal CDF and PDF.

Usage

ee_dlasso_regression(
  theta,
  X,
  y,
  model,
  penalty,
  s = 1e-06,
  weights = NULL,
  center = 0,
  offset = NULL
)

Arguments

Value

A p-by-n matrix.

Examples

# A penalty vector gives one value per column of the design matrix. A scalar
# penalty would shrink the intercept along with the slopes. The estimating
# equation carries the penalty's derivative. Here that derivative is smooth,
# so the estimating equation is differentiable and no warning is issued,
# unlike ee_lasso_regression() whose derivative has unbounded slope at the
# penalty center.
fit <- m_estimate(
  mpg ~ wt + hp,
  data = mtcars,
  .ee = ee_dlasso_regression,
  model = "linear",
  penalty = c(0, 5, 5)
)
coef(fit)

Description

Combines L1 (approximate LASSO via bridge) and L2 (ridge) penalties at a given ratio. When ratio = 1, this is LASSO; when ratio = 0, ridge.

Usage

ee_elasticnet_regression(
  theta,
  X,
  y,
  model,
  penalty,
  ratio,
  epsilon = 0.003,
  weights = NULL,
  center = 0,
  offset = NULL
)

Arguments

Value

A p-by-n matrix.

Examples

# A penalty vector gives one value per column of the design matrix. A scalar
# penalty would shrink the intercept along with the slopes.
#
# The L1 half of the penalty enters the estimating equation as its own
# derivative, and that derivative has unbounded slope at the penalty center.
# The estimating equation is therefore not differentiable there, so the
# bread matrix is undefined and the fit warns once that the sandwich
# variance should not be trusted here.
fit <- m_estimate(
  mpg ~ wt + hp,
  data = mtcars,
  .ee = ee_elasticnet_regression,
  model = "linear",
  penalty = c(0, 5, 5),
  ratio = 0.5
)
coef(fit)

Description

Implements the hyperbolic E-max (Hill) model:

$R_i = \theta_0 + \frac{\theta_{max} D_i}{\theta_{50} + D_i}$

Usage

ee_emax(theta, dose, response, loss = NULL, k = NULL)

Arguments

Value

A 3-by-n matrix, with rows named e0, emax, and ed50.

See Also

ee_emax_ed() for the effective dose at a given level, which is stacked with this equation to give it a sandwich standard error.

Examples

# Dose-response of a herbicide on ryegrass root length. The response falls
# with dose, so the maximum change in response is initialized negative.
psi <- function(theta) {
  ee_emax(theta, dose = inderjit$dose, response = inderjit$response)
}

m <- m_estimate(stacked_equations = psi, init = c(8, -8, 2))

# Zero-dose response, maximum change in response, and ED50
coef(m)

Description

Computes the effective dose at level delta from the E-max model. Should be stacked with ee_emax().

Usage

ee_emax_ed(theta, dose, delta, ed50)

Arguments

Value

A 1-by-n matrix.

See Also

ee_emax(), the dose-response equation this one is stacked with.

Examples

# This equation carries no information on its own, so stack it with the
# E-max equation that supplies the ED50. Here delta = 0.9 requests the ED90.
psi <- function(theta) {
  emax <- ee_emax(
    theta[1:3],
    dose = inderjit$dose,
    response = inderjit$response
  )
  ed90 <- ee_emax_ed(
    theta[4],
    dose = inderjit$dose,
    delta = 0.9,
    ed50 = theta[3]
  )
  rbind(emax, ed90)
}

m <- m_estimate(stacked_equations = psi, init = c(8, -8, 2, 10))

# theta_4 is the ED90. Stacking is what gives it a standard error of its own,
# so summary() rather than coef() is what shows the gain.
summary(m)

Description

Estimates the parameters of a structural nested mean model via g-estimation. Supports both inefficient (X = NULL) and efficient (X provided) g-estimators, and both linear and log-linear (Poisson) structural mean models.

Usage

ee_gestimation_snmm(
  theta,
  y,
  A,
  W,
  V,
  X = NULL,
  model = "linear",
  weights = NULL
)

Arguments

Value

A matrix of estimating equation contributions. The structural mean model rows are named ⁠SNM phi_0⁠ through ⁠SNM phi_(b-1)⁠, matching the zero-based subscripts the literature gives those parameters. The propensity score rows are named W_1 through W_c, and, for the efficient g-estimator, the outcome model rows are named X_1 through X_d.

Examples

# A confounded binary treatment whose true effect on the outcome is -2.
set.seed(42)
n <- 500
W <- rbinom(n, 1, 0.5)
A <- rbinom(n, 1, 0.25 + 0.5 * W)
Y <- 5 + 2 * W - 2 * A + rnorm(n)

W_ps <- cbind(1, W) # Propensity score design matrix

# An intercept-only structural mean model gives a single causal contrast.
# Build it with rep() so the column has one entry per observation.
V <- cbind(rep(1, n))

psi <- function(theta) {
  ee_gestimation_snmm(theta, y = Y, A = A, W = W_ps, V = V, model = "linear")
}

# theta holds the structural mean model coefficient, which is the causal
# effect, followed by the two propensity score coefficients.
m <- m_estimate(stacked_equations = psi, init = rep(0, 3))
coef(m)

Description

Returns a stacked set of estimating equations for the g-formula. When X0 = NULL, estimates a single causal mean under the plan encoded by X1. When X0 is provided, estimates the average causal effect (difference between two plans).

Usage

ee_gformula(theta, y, X, X1, X0 = NULL, force_continuous = FALSE)

Arguments

Value

A matrix of estimating equation contributions. When X0 = NULL the first row is named causal_mean; when X0 is provided the first three rows are named ACE, E[Y^1], and E[Y^0], where 1 and 0 index the two plans. The outcome model rows are named X_1 through X_p for the columns of X.

Examples

# A binary treatment, two confounders, and a continuous outcome whose true
# average causal effect is 1.5.
set.seed(42)
n <- 1000
W1 <- rnorm(n)
W2 <- rbinom(n, 1, 0.4)
A <- rbinom(n, 1, inverse_logit(-0.5 + 0.5 * W1 + 0.3 * W2))
Y <- 2 + 1.5 * A + W1 - 0.5 * W2 + rnorm(n)

X <- cbind(1, A, W1, W2) # Observed design matrix
X1 <- cbind(1, 1, W1, W2) # Everyone treated
X0 <- cbind(1, 0, W1, W2) # Everyone untreated

psi <- function(theta) ee_gformula(theta, y = Y, X = X, X1 = X1, X0 = X0)

# theta holds the average causal effect, the mean under treatment, and the
# mean under no treatment, followed by the four outcome model coefficients.
m <- m_estimate(stacked_equations = psi, init = rep(0, 7))
coef(m)[1:3]

Description

Returns a p-by-n matrix of estimating equation contributions for a GLM with specified distribution and link function:

$\psi_i(\theta) = \{Y_i - g^{-1}(X_i^T \theta)\} \frac{D(\theta)}{v(\theta)} X_i$

Usage

ee_glm(
  theta,
  X,
  y,
  distribution,
  link,
  hyperparameter = NULL,
  weights = NULL,
  offset = NULL
)

Arguments

Value

A p-by-n matrix. For the gamma and negative binomial distributions the matrix has ncol(X) + 1 rows, the final row being the nuisance equation for that distribution's extra parameter (the gamma shape or the negative binomial dispersion). Those two distributions name their rows X_1 through X_p for the columns of X, followed by log_shape or log_dispersion. Every other distribution estimates one coefficient per design column and leaves its rows unnamed.

Gamma

The gamma distribution estimates one additional shape parameter beyond the regression coefficients, so theta has length ncol(X) + 1. The final element is the log of the shape $\alpha$ (the reciprocal of the GLM dispersion $\phi = 1 / \alpha$), and the variance function is $\mu^2$. An extra nuisance estimating equation for $\log(\alpha)$ is appended:

$\left(1 - \frac{Y_i}{\mu_i}\right) + \log\left(\frac{\alpha Y_i}{\mu_i}\right) - \psi(\alpha)$

where $\psi$ is the digamma function. The returned matrix therefore has ncol(X) + 1 rows.

Negative binomial

The negative binomial distribution estimates one additional dispersion parameter beyond the regression coefficients, so theta has length ncol(X) + 1. The final element is the log of the dispersion $\alpha$, and the variance function is $\mu + \alpha \mu^2$. An extra nuisance estimating equation for $\log(\alpha)$ is appended so that the uncertainty in the dispersion is carried honestly through the sandwich variance:

$-\alpha^{-2}\psi\left(Y_i + \alpha^{-1}\right) + \alpha^{-2}\psi\left(\alpha^{-1}\right) + \frac{Y_i}{\alpha^2 \mu_i + \alpha} - \frac{\frac{\alpha \mu_i}{\alpha \mu_i + 1} + \log\left(\frac{1}{\alpha \mu_i + 1}\right)}{\alpha^2}$

where $\psi$ is the digamma function. The returned matrix therefore has ncol(X) + 1 rows.

Tweedie

The tweedie distribution uses the variance function $v(\mu) = \mu^p$, where the power $p$ is the fixed hyperparameter rather than an estimated parameter. No nuisance estimating equation is appended, so theta has length ncol(X) and the returned matrix has ncol(X) rows. The power recovers familiar special cases: $p = 1$ gives the Poisson variance $\mu$ and $p = 2$ gives the gamma variance $\mu^2$, so those choices share the corresponding beta score equations. The hyperparameter must be non-negative.

Examples

# Negative binomial GLM with a log link estimates the regression
# coefficients plus one log-dispersion parameter, so init has ncol(X) + 1
# entries.
set.seed(1)
n <- 200
X <- cbind(1, rnorm(n), rnorm(n))
mu <- exp(0.5 + 0.5 * X[, 2] - 0.3 * X[, 3])
y <- rnbinom(n, size = 2, mu = mu)

psi <- function(theta) {
  ee_glm(theta, X = X, y = y, distribution = "negative_binomial",
         link = "log")
}

m <- m_estimate(stacked_equations = psi, init = c(0, 0, 0, 0))
coef(m)

# Tweedie GLM with a log link. The power p is a fixed hyperparameter, not an
# estimated parameter, so init has ncol(X) entries. Here p = 1.5 sits in the
# compound Poisson-gamma regime, appropriate for non-negative data with a
# point mass at zero.
set.seed(2)
mu <- exp(0.5 + 0.5 * X[, 2] - 0.3 * X[, 3])
y_tw <- rpois(n, lambda = mu)

psi_tw <- function(theta) {
  ee_glm(theta, X = X, y = y_tw, distribution = "tweedie", link = "log",
         hyperparameter = 1.5)
}

m_tw <- m_estimate(stacked_equations = psi_tw, init = c(0, 0, 0))
coef(m_tw)

Description

Estimates the average causal effect using IPW with a logistic propensity score model.

Usage

ee_ipw(theta, y, A, W, truncate = NULL, weights = NULL)

Arguments

Value

A (3+b)-by-n matrix of estimating equation contributions, with the first three rows named ACE, E[Y^1], and E[Y^0] and the propensity score rows named W_1 through W_b for the columns of W.

Examples

# A binary treatment, two confounders, and a continuous outcome whose true
# average causal effect is 1.5.
set.seed(42)
n <- 1000
W1 <- rnorm(n)
W2 <- rbinom(n, 1, 0.4)
A <- rbinom(n, 1, inverse_logit(-0.5 + 0.5 * W1 + 0.3 * W2))
Y <- 2 + 1.5 * A + W1 - 0.5 * W2 + rnorm(n)

W_ps <- cbind(1, W1, W2) # Propensity score design matrix

psi <- function(theta) ee_ipw(theta, y = Y, A = A, W = W_ps)

# theta holds the average causal effect, the mean under treatment, and the
# mean under no treatment, followed by the three propensity score
# coefficients.
m <- m_estimate(stacked_equations = psi, init = rep(0, 6))
coef(m)[1:3]

Description

Estimates the parameters of a marginal structural model using inverse probability weighting with a logistic propensity score model.

Usage

ee_ipw_msm(
  theta,
  y,
  A,
  W,
  V,
  distribution,
  link,
  hyperparameter = NULL,
  truncate = NULL,
  weights = NULL
)

Arguments

Value

A (c+b)-by-n matrix of estimating equation contributions. The marginal structural model rows are named ⁠MSM alpha_0⁠ through ⁠MSM alpha_(c-1)⁠, matching the zero-based subscripts the literature gives those parameters. The propensity score rows are named W_1 through W_b for the columns of W.

Examples

# A confounded binary treatment whose true effect on the outcome is -2.
set.seed(42)
n <- 500
W <- rbinom(n, 1, 0.5)
A <- rbinom(n, 1, 0.25 + 0.5 * W)
Y <- 5 + 2 * W - 2 * A + rnorm(n)

W_ps <- cbind(1, W) # Propensity score design matrix
V <- cbind(1, A) # Marginal structural model design matrix

psi <- function(theta) {
  ee_ipw_msm(
    theta,
    y = Y,
    A = A,
    W = W_ps,
    V = V,
    distribution = "normal",
    link = "identity"
  )
}

# theta holds the two marginal structural model coefficients, whose slope is
# the causal effect, followed by the two propensity score coefficients.
m <- m_estimate(stacked_equations = psi, init = rep(0, 4))
coef(m)

Description

Estimates the causal effect using the usual IV / Wald estimator. The parameter of interest is the additive effect of treatment A on outcome Y leveraging instrument Z.

Usage

ee_iv_causal(theta, y, A, Z, weights = NULL)

Arguments

Value

A 2-by-n matrix of estimating equation contributions, with rows named causal_effect and mean_Z.

Examples

# An unmeasured confounder U biases the association between A and Y, but the
# instrument Z affects Y only through A. The true causal effect is 3.
set.seed(123)
n <- 500
Z <- rbinom(n, 1, 0.5)
U <- rnorm(n)
A <- rbinom(n, 1, inverse_logit(-1 + 3 * Z + U))
Y <- 3 * A - U + rnorm(n, sd = 0.5)

psi <- function(theta) ee_iv_causal(theta, y = Y, A = A, Z = Z)

# theta holds the causal effect followed by the mean of the instrument.
m <- m_estimate(stacked_equations = psi, init = c(0, 0.5))
coef(m)

Description

Uses the bridge penalty with gamma = 1 + epsilon to approximate LASSO (L1 penalty). The true LASSO is not differentiable at zero, so an approximation is used.

Usage

ee_lasso_regression(
  theta,
  X,
  y,
  model,
  penalty,
  epsilon = 0.003,
  weights = NULL,
  center = 0,
  offset = NULL
)

Arguments

Value

A p-by-n matrix.

Examples

# A penalty vector gives one value per column of the design matrix. A scalar
# penalty would shrink the intercept along with the slopes.
#
# The approximate L1 penalty enters the estimating equation as its own
# derivative, and that derivative has unbounded slope at the penalty center.
# The estimating equation is therefore not differentiable there, so the
# bread matrix is undefined and the fit warns once that the sandwich
# variance should not be trusted here.
fit <- m_estimate(
  mpg ~ wt + hp,
  data = mtcars,
  .ee = ee_lasso_regression,
  model = "linear",
  penalty = c(0, 5, 5)
)
coef(fit)

Description

Implements the 4-parameter log-logistic model:

$R_i = \theta_0 + \frac{\theta_m - \theta_0}{1 + \exp[\theta_s (\log D_i - \log \theta_{50})]}$

Usage

ee_loglogistic(theta, dose, response, loss = NULL, k = NULL)

Arguments

Value

A 4-by-n matrix, with rows named lower, upper, ed50, and steepness.

See Also

ee_loglogistic_ed() for the effective dose at a given level, which is stacked with this equation to give it a sandwich standard error.

Examples

# Four-parameter log-logistic dose-response for the ryegrass data.
psi <- function(theta) {
  ee_loglogistic(theta, dose = inderjit$dose, response = inderjit$response)
}

# The default rootSolve solver does not converge here, so nleqslv is used.
m <- m_estimate(
  stacked_equations = psi,
  init = c(0.2, 8, 2, 1),
  solver = "nleqslv"
)

# Lower limit, upper limit, ED50, and steepness
coef(m)

Description

Computes the effective dose at level delta from the log-logistic model. Should be stacked with ee_loglogistic().

Usage

ee_loglogistic_ed(theta, dose, delta, lower, upper, ed50, steepness)

Arguments

Value

A 1-by-n matrix.

See Also

ee_loglogistic(), the dose-response equation this one is stacked with.

Examples

# The lower limit is held at zero instead of being estimated: root length
# cannot fall below zero, and the five-parameter stack diverges on these
# data. The lower-limit row of the log-logistic equation is therefore
# dropped, leaving the upper limit, ED50, and steepness to be estimated
# alongside the ED90.
psi <- function(theta) {
  loglogistic <- ee_loglogistic(
    c(0, theta[1:3]),
    dose = inderjit$dose,
    response = inderjit$response
  )
  ed90 <- ee_loglogistic_ed(
    theta[4],
    dose = inderjit$dose,
    delta = 0.9,
    lower = 0,
    upper = theta[1],
    ed50 = theta[2],
    steepness = theta[3]
  )
  rbind(loglogistic[-1, , drop = FALSE], ed90)
}

# This stacked system is sensitive to its starting values, so a reasonable
# init matters more here than it does for most equations.
m <- m_estimate(stacked_equations = psi, init = c(8, 2, 1, 5))

# Upper limit, ED50, steepness, and the ED90. Stacking is what gives the ED90
# a standard error of its own, so summary() rather than coef() shows the
# gain.
summary(m)

Description

Returns a 1-by-n matrix of estimating equation contributions for the mean: $\psi_i(\theta) = w_i (Y_i - \theta)$.

Usage

ee_mean(theta, y, weights = NULL)

Arguments

Value

A 1-by-n matrix.

Examples

y <- c(1, 2, 3, 1, 4, 5, 3, 2, 6, 7)
psi <- function(theta) ee_mean(theta, y = y)
m <- m_estimate(stacked_equations = psi, init = 0)
coef(m)

Description

Returns a 1-by-n matrix of estimating equation contributions for the geometric mean. When log_theta = TRUE (default), solves $\log(Y_i) - \log(\theta)$. When log_theta = FALSE, solves $\log(Y_i) - \theta$ where $\theta$ is the log of the geometric mean.

Usage

ee_mean_geometric(theta, y, weights = NULL, log_theta = TRUE)

Arguments

Value

A 1-by-n matrix.

Examples

y <- c(1, 2, 3, 1, 4, 5, 3, 2, 6, 7)
psi <- function(theta) ee_mean_geometric(theta, y = y)
m <- m_estimate(stacked_equations = psi, init = 1)
coef(m)

Description

Returns a 1-by-n matrix of estimating equation contributions for the robust mean using the specified loss function: $\psi_i(\theta) = f_k(Y_i - \theta)$.

Usage

ee_mean_robust(theta, y, k, loss = "huber")

Arguments

Value

A 1-by-n matrix.

Examples

# Forty-nine standard normal observations plus one gross outlier.
set.seed(1)
y <- c(rnorm(49), 50)
psi <- function(theta) ee_mean_robust(theta, y = y, k = 1.345, loss = "huber")
m <- m_estimate(stacked_equations = psi, init = 0)

# The Huber estimate stays near the uncontaminated mean, unlike mean(y).
coef(m)
mean(y)

Description

Estimates the mean of an outcome with missing data using a weighted sensitivity analysis approach. Handles MCAR, MAR, and MNAR mechanisms by specifying a user-defined sensitivity function q_eval and a monotone increasing distribution function H_function.

Usage

ee_mean_sensitivity_analysis(theta, y, delta, X, q_eval, H_function)

Arguments

Value

A (1+b)-by-n matrix of estimating equation contributions, with the first row named corrected_mean and the missingness model rows named X_1 through X_b for the columns of X.

Examples

# An outcome observed for only part of the sample, with missingness driven by
# the measured covariate W.
set.seed(42)
n <- 500
W <- rbinom(n, 1, 0.5)
Y_full <- 200 - 35 * W + rnorm(n, sd = 5)
delta <- rbinom(n, 1, inverse_logit(2 + W))

# Missing outcomes never enter the estimating equation, so any placeholder
# value works; a zero keeps the arithmetic finite.
Y <- ifelse(delta == 1, Y_full, 0)
X <- cbind(1, W) # Missingness model design matrix

# A sensitivity function of zero everywhere assumes the outcome is missing at
# random given W. Nonzero values encode departures from that assumption.
psi <- function(theta) {
  ee_mean_sensitivity_analysis(
    theta,
    y = Y,
    delta = delta,
    X = X,
    q_eval = rep(0, n),
    H_function = inverse_logit
  )
}

# theta holds the corrected mean, started near the observed outcomes,
# followed by the two missingness model coefficients.
m <- m_estimate(stacked_equations = psi, init = c(180, 0, 0))
coef(m)