Fits a generalized linear model with empirical likelihood.
el_glm(
formula,
family = gaussian,
data,
weights = NULL,
na.action,
start = NULL,
etastart = NULL,
mustart = NULL,
offset,
control = el_control(),
...
)An object of class formula (or one that can be coerced to
that class): a symbolic description of the model to be fitted.
A description of the error distribution and link function to be used in the model. Only the result of a call to a family function is supported. See ‘Details’.
An optional data frame, list or environment (or object coercible
by as.data.frame() to a data frame) containing the variables in the
formula. If not found in data, the variables are taken from
environment(formula).
An optional numeric vector of weights to be used in the
fitting process. Defaults to NULL, corresponding to identical weights. If
non-NULL, weighted empirical likelihood is computed.
A function which indicates what should happen when the data
contain NAs. The default is set by the na.action setting of
options, and is na.fail if that is unset.
Starting values for the parameters in the linear predictor.
Defaults to NULL and is passed to glm.fit().
Starting values for the linear predictor. Defaults to NULL
and is passed to glm.fit().
Starting values for the vector of means. Defaults to NULL
and is passed to glm.fit().
An optional expression for specifying an a priori known
component to be included in the linear predictor during fitting. This
should be NULL or a numeric vector or matrix of extents matching those of
the response. One or more offset terms can be included in the formula
instead or as well, and if more than one are specified their sum is used.
An object of class ControlEL constructed by
el_control().
Additional arguments to be passed to glm.control().
An object of class of GLM.
Suppose that we observe \(n\) independent random variables \({Z_i} \equiv {(X_i, Y_i)}\) from a common distribution, where \(X_i\) is the \(p\)-dimensional covariate (including the intercept if any) and \(Y_i\) is the response. A generalized linear model specifies that \({\textrm{E}(Y_i | X_i)} = {\mu_i}\), \({G(\mu_i)} = {X_i^\top \theta}\), and \({\textrm{Var}(Y_i | X_i)} = {\phi V(\mu_i)}\), where \(\theta = (\theta_0, \dots, \theta_{p-1})\) is an unknown \(p\)-dimensional parameter, \(\phi\) is an optional dispersion parameter, \(G\) is a known smooth link function, and \(V\) is a known variance function.
With \(H\) denoting the inverse link function, define the quasi-score
$${g_1(Z_i, \theta)} =
\left\{
H^\prime(X_i^\top \theta) \left(Y_i - H(X_i^\top \theta)\right) /
\left(\phi V\left(H(X_i^\top \theta)\right)\right)
\right\}
X_i.$$
Then we have the estimating equations
\(\sum_{i = 1}^n g_1(Z_i, \theta) = 0\).
When \(\phi\) is known, the (profile) empirical likelihood ratio for a
given \(\theta\) is defined by
$$R_1(\theta) =
\max_{p_i}\left\{\prod_{i = 1}^n np_i :
\sum_{i = 1}^n p_i g_1(Z_i, \theta) = 0,\
p_i \geq 0,\
\sum_{i = 1}^n p_i = 1
\right\}.$$
With unknown \(\phi\), we introduce another estimating function based on
the squared residuals. Let \({\eta} = {(\theta, \phi)}\) and
$${g_2(Z_i, \eta)} =
\left(Y_i - H(X_i^\top \theta)\right)^2 /
\left(\phi^2 V\left(H(X_i^\top \theta)\right)\right) - 1 / \phi.$$
Now the empirical likelihood ratio is defined by
$$R_2(\eta) =
\max_{p_i}\left\{\prod_{i = 1}^n np_i :
\sum_{i = 1}^n p_i g_1(Z_i, \eta) = 0,\
\sum_{i = 1}^n p_i g_2(Z_i, \eta) = 0,\
p_i \geq 0,\
\sum_{i = 1}^n p_i = 1
\right\}.$$
el_glm() first computes the parameter estimates by calling glm.fit()
(with ... if any) with the model.frame and model.matrix obtained from
the formula. Note that the maximum empirical likelihood estimator is the
same as the the quasi-maximum likelihood estimator in our model. Next, it
tests hypotheses based on asymptotic chi-square distributions of the
empirical likelihood ratio statistics. Included in the tests are overall
test with
$$H_0: \theta_1 = \theta_2 = \cdots = \theta_{p-1} = 0,$$
and significance tests for each parameter with
$$H_{0j}: \theta_j = 0,\ j = 0, \dots, p-1.$$
The available families and link functions are as follows:
gaussian: "identity", "log", and "inverse".
binomial: "logit", "probit", and "log".
poisson: "log", "identity", and "sqrt".
quasipoisson: "log", "identity", and "sqrt".
Chen SX, Cui H (2003). ``An Extended Empirical Likelihood for Generalized Linear Models.'' Statistica Sinica, 13(1), 69--81.
Kolaczyk ED (1994). ``Empirical Likelihood for Generalized Linear Models.'' Statistica Sinica, 4(1), 199--218.
EL, GLM, el_lm(), elt(),
el_control()
data("warpbreaks")
fit <- el_glm(wool ~ .,
family = binomial, data = warpbreaks, weights = NULL, na.action = na.omit,
start = NULL, etastart = NULL, mustart = NULL, offset = NULL
)
summary(fit)
#>
#> Empirical Likelihood
#>
#> Model: glm (binomial family with logit link)
#>
#> Call:
#> el_glm(formula = wool ~ ., family = binomial, data = warpbreaks,
#> weights = NULL, na.action = na.omit, start = NULL, etastart = NULL,
#> mustart = NULL, offset = NULL)
#>
#> Number of observations: 54
#> Number of parameters: 4
#>
#> Parameter values under the null hypothesis:
#> (Intercept) breaks tensionM tensionH
#> 6.337e-08 0.000e+00 0.000e+00 0.000e+00
#>
#> Lagrange multipliers:
#> [1] 1.81250 -0.05261 -0.41072 -0.71437
#>
#> Maximum EL estimates:
#> (Intercept) breaks tensionM tensionH
#> 1.67965 -0.04677 -0.44662 -0.67062
#>
#> logL: -217.4 , logLR: -1.97
#> Chisq: 3.94, df: 3, Pr(>Chisq): 0.268
#> Constrained EL: converged
#>
#> Coefficients:
#> Estimate Chisq Pr(>Chisq)
#> (Intercept) 1.67965 6.226 0.0126 *
#> breaks -0.04677 3.941 0.0471 *
#> tensionM -0.44662 0.568 0.4511
#> tensionH -0.67062 1.628 0.2020
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Dispersion for binomial family: 1
#>