offset+ X beta + Z v, where
X is the design matrix of fixed effects and Z is a design matrix of random effects.
The function also handles a linear predictor (with only fixed effects) for the residual variance.
HLfit(formula, data, family = gaussian(), rand.family = gaussian(), resid.model = ~1, resid.formula, REMLformula = NULL, verbose = c(warn = TRUE, trace = FALSE, summary = FALSE), HLmethod = "HL(1,1)", control.HLfit = list(), control.glm = list(), init.HLfit = list(), ranFix = list(), etaFix = list(), prior.weights = NULL, processed = NULL)
## see 'rand.family' argument for inverse.Gammafamily object describing the distribution of the response variable.
Possible values include the gaussian, poisson, binomial, Gamma, negbin and
COMPoisson families. Possible combinations of family and link are those allowed by each of these families
(see family for the first four, and specific documentation pages for the last two).
family object describing the distribution of the random effect, or a list of
family objects for different random effects (see Examples). Possible options are
gaussian(), Gamma(log), Gamma(identity) (see Details), Beta(logit), inverse.Gamma(-1/mu), and inverse.Gamma(log).
For discussion of these alternatives see Lee and Nelder 2001 or Lee et al. 2006, p. 178-.
Here the family gives the distribution of a random effect $u$
and the link gives v as function of $u$ (see Details).
If there are several random effects and only one family is given, this family holds for all random effects.
phi of the residual error, or a list
with such a formula as element formula, and the optional element family for the model family of the residual model.
Currently the formula can only contain fixed effects, including an offset. A log link is assumed by default, but an identity link can be requested by setting resid.model$family to Gamma(identity).
resid.model.
formula that allows the estimation of dispersion parameters, and
computation of restricted likelihood (p_bv) under a model different from the predictor formula. For example, if only random effects are included in REMLformula, an ML fit is performed and p_bv equals
the marginal likelihood (or its approximation), p_v. This ML fit can be performed more simply by setting
HLmethod="ML" and leaving REMLformula at its default NULL value.
trace controls various diagnostic (possibly messy) messages about the iterations.
summary controls whether a summary of the fit is called by HLfit.
warn is for programming purposes and best ignored.
"REML", "ML", "EQL-" and "EQL+" for all models;
"PQL" (="REPQL") and "PQL/L" for GLMMs only;
and (only for those curious to experiment) expressions of the
form "HL(<...>)", "ML(<...>)" and "RE(<...>)". HL and RE are equivalent.
The default behaviour is RE(1,1), which by default performs REML (standard REML for LMMs,
an extended definition for other models). But it can also perform non-standard forms of REML
(indeed including ML), depending on the REMLformula given. EQL stands for the EQL method of Lee and Nelder (2001). The '+' version includes the d v/ d tau correction
described p. 997 of that paper, and the '-' version ignores it.
PQL can be seen as the version of EQL- for GLMMs.
PQL/L is PQL without the leverage corrections that define REML estimation of random-effect parameters.
For GLMs, the default is still REML. For binomial and Poisson GLMs, this cannot be distinguished from an exact ML fit.
Note that the glm function performs an EQL analysis, which will differ from the REML (and ML) one for Gamma GLMs.
See Details for the more general syntax.
resid.family allows one to change the link for modeling of residual variance $\phi$, which is "log" by default. The family is always Gamma, so the non-default possible values of resid.family are Gamma(identity) or Gamma(inverse). Only the default value ensures that the fitted $\phi$ is positive.
Controls for the fitting algorithms should be ignored in routine use. They are
conv.threshold,
a convergence threshold for the iterative algorithm, which controls whether
linear predictor terms (fixed effects and inferred random effects), and dispersion parameter estimates have converged. Defaults to 1e-05;
break_conv_logL, a boolean specifying whether the iterative algorithm should terminate when log-likelihood appears to have convergence (roughly, when its relative variation over on iteration is lower than 1e-8). Default is FALSE.
iter.mean.dispFix, the number of iterations of the iterative algorithm for coefficients of the linear predictor,
if no dispersion parameters are estimated by the iterative algorithm. Defaults to 200;
iter.mean.dispVar, the number of iterations of the iterative algorithm for coefficients of the linear predictor,
if some dispersion parameter(s) is estimated by the iterative algorithm. Defaults to 50;
max.iter, the number of iterations of the iterative algorithm for joint estimation of dispersion parameters and
of coefficients of the linear predictor. Defaults to 200. This is typically much more than necessary,
unless there is little information to separately estimate $\lambda$ and $\phi$ parameters.
glm.control; e.g. control.glm=list(maxit=100). See glm.control for further details.
fixef for fixed effect estimates (beta),
v_h for random effects vector v in the linear predictor,
lambda for the parameter determining the variance of random effects $u$ as drawn from the rand.family distribution
phi for the residual variance.
However, this argument can be ignored in routine use.
lambda, and also phi for gaussian and Gamma HGLMs.
Inhibits the estimation of these parameters.
beta. etaFix$beta should be a vector with names matching (a subset of) coefficient names of a fit without fixed values. It provides a convenient interface for fixing (some of) the fixed-effects coefficients ($\beta$). In contrast to an offset specification, it affects the REML correction for estimation of dispersion parameters, which depends only on which $\beta$ coefficients are estimated. However, for non-standard use, REML can still be performed as if all $\beta$ coefficients were estimated, by adding attribute keepInREML=TRUE to etaFix$beta (see Examples). These different behaviours will be overridden whenever a non-null REMLformula is provided.
glm. This fits the data to a model with residual variance phi/prior.weights, so that increasing the weights by a constant factor f will yield (Intercept) estimates of phi also increased by f (this effect cannot be generally achieved if a non-trivial resid.formula with log link is used). This is not necessarily the way prior weights are interpreted in widely used packages, but this is consistent with what glm.
corrHLfit).
HLfit, actually a list with many elements, several of which represent input arguments.
Some elements may be undocumented.A few extractor functions are available (see extractors),
and should be used as far as possible as they should be backward-compatible from version 1.4 onwards, while the structure of the return object may still evolve. The following information will be useful for extracting further elements of the object.Elements describing the fit include descriptors of the fit:Information about the input is contained in output elements named as HLfit or corrHLfit arguments (data,family,resid.family,ranFix,prior.weights), with the following notable exceptions or modifications:Further miscellaneous diagnostics and descriptors of model structure:Finally, the object includes programming tools: call, spaMM.version, get_w_h_coeffs, get_beta_w_cov, get_invColdoldList, get_logdispObject.Fitting methods: Many approximations for likelihood have been defined to fit mixed models (e.g. Noh and Lee (2007) for some overview), and this function only considers a subset of them, but it adds a new complication in terms of REML methods. For example, PQL as originally defined by Breslow and Clayton uses REML to estimate dispersion parameters, but this function allows one to use ML instead. Moreover, it allows some non-standard specification of the model formula that determines the conditional distribution used in REML.
In the more general syntax for HLmethod, used as e.g. HLmethod="RE(1,1)"
the first '1' means that a first order Laplace approximation to the likelihood is used to estimate fixed effects
(a '0' would instead mean that the h likelihood is used as the objective function).
The second '1' means that a first order Laplace approximation to the likelihood or restricted likelihood
is used to estimate dispersion parameters, including the dv/d tau term specifically discussed by Lee & Nelder 2001, p. 997
(a '0' would instead mean that these terms are ignored).
It is possible to enforce the EQL approximation for estimation of dispersion parameter by adding a third index with value 0.
"HL(0,1,0)" is Lee & Nelder's (2001) method, i.e. "EQL+".
For a Gamma GLM with log link, ML and EQL results will differ in their phi estimates, and the EQL estimate will match that from the glm function.
Random effects are constructed in several steps. first, a vector u of independent and identically distributed (iid) random effects is drawn from some distribution;
second, a transformation v=f(u) is applied to each element (this defines v which elements are still iid); third, correlated random effects are obtained as Lv
where L is the “square root” of a correlation matrix (this may be meaningful only for Gaussian random effects).
Finally, a matrix Z (or sometimes ZA, see Predictor) allows to specify how the correlated random effects
affect the response values. In particular, Z is the identity matrix if there is a single observation (response) for each location, but otherwise
its elements $z_{ji}$ are 1 for the $j$th observation in the $i$th location.
The design matrix for v is then of the form ZL.
The specification of the random effects u and v handles the following cases:
The standard errors reported may sometimes be misleading. For each set of parameters among $\beta$, $\lambda$, and $\phi$ parameters these are computed assuming that the other parameters are known without error. This is why they are labelled Cond. SE (conditional standard error). This is most uninformative in the unusual case where $\lambda$ and $\phi$ are not separately estimable parameters. Further, the SEs for $\lambda$ and $\phi$ are rough approximations as discussed in particular by Smyth et al. (2001; $V_1$ method).
The extractor function get_any_IC can compute information criteria (“AIC”) and effective degrees of freedom from the HLfit results. See the get_any_IC documentation for more details.
Cox, D. R. and Donnelly C. A. (2011) Principles of Applied Statistics. Cambridge Univ. Press.
Ha, I. D., Lee, Y. and MacKenzie, G. (2007) Model selection for multi-component frailty models. Statistics in Medicine 26: 4790-4807.
Lee, Y., Nelder, J. A. (2001) Hierarchical generalised linear models: A synthesis of generalised linear models, random-effect models and structured dispersions. Biometrika 88, 987-1006.
Lee, Y., Nelder, J. A. and Pawitan, Y. (2006). Generalised linear models with random effects: unified analysis via h-likelihood. Chapman & Hall: London.
Noh, M., and Lee, Y. (2007). REML estimation for binary data in GLMMs, J. Multivariate Anal. 98, 896-915.
Smyth GK, Huele AF, Verbyla AP (2001). Exact and approximate REML for heteroscedastic regression. Statistical Modelling 1, 161-175.
Vaida, F., and Blanchard, S. 2005. Conditional Akaike information for mixed-effects models. Biometrika 92, 351-370.
HLCor for estimation with given spatial correlation parameters;
corrHLfit for joint estimation with spatial correlation parameters.
data(wafers)
## Gamma GLMM with log link
HLfit(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
resid.model = ~ X3+I(X3^2) ,data=wafers)
## Gamma - inverseGamma HGLM with log link
HLfit(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
HLmethod="HL(1,1)",rand.family=inverse.Gamma(log),
resid.model = ~ X3+I(X3^2) ,data=wafers)
Run the code above in your browser using DataLab