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.formula = ~1, REMLformula = NULL,
verbose=c(warn=TRUE,trace=FALSE,summary=FALSE),
HLmethod = "HL(1,1)", control.HLfit = list(), init.HLfit = list(),
ranFix = list(), etaFix = list(), prior.weights = rep(1, nobs), processed = NULL)
## see 'rand.family' argument for inverse.Gamma
## see 'control.HLfit' argument for AICfamily object describing the distribution of the response variable.
Possible values include the gaussian, poisson, binomial and Gamma families. Possible links are identity, log,
inverse, logit, probit, and cloglog. Possible combinafamily 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), Beta(lphi of the residual error. Currently can only contain fixed effects, including an offset. A log link is assumed by default, but an identity link can be used (see control.HLfit).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 artrace 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"REML", "ML", "EQL-" and "EQL+" for all models;
"PQL" (="REPQL") and "PQL/L" for GLMMs only;
and (only for those curious to experiment) eAIC=TRUE provides several information criteria and a related FALSfixef for fixed effect estimates (beta),
v_h for random effects vector v in the linear predictor,
lambdalambda, 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 provglm. 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 (corrHLfit code).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 (to further deal with this, the return object includes a version tag as element spaMM.version). The following information will be useful for extracting further elements of the object.
Elements describing the fit include:fixef function)ranef function)fitted function)p_v and the (approximate) restricted likelihood p_bv (the latter two available
through the logLik function).family object corresponding to the family inputfamily object corresponding to the rand.family inputranFix inputHLmethod, 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 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:
[object Object],[object Object],[object Object],[object Object]
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).
Effective degrees of freedom for the random effects: this number (approximately) characterizes the expectation of a goodness of fit statistic discussed by Lee and Nelder (2001), which gave a general formula for it in HGLMs. It is computed if AIC=TRUE and returned as the object's $APHLs$GoFdf.
AIC: The literature considers several information criteria for mixed models. The conditional AIC (Vaida and Blanchard 2005) is notable in involving the conditional likelihood and the effective degrees of freedom. Lee et al. (2006) and Ha et al (2007) similarly define a corrected AIC [i.e., AIC(D*) in their eq. 7]. The conditional AIC returned by HLfit includes both this effective df, df for estimated fixed effects, and df for estimated parameters of the variance of random effects. HLfit will also return the classical AIC (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.formula = ~ 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.formula = ~ X3+I(X3^2) ,data=wafers)Run the code above in your browser using DataLab