Fitting Hierarchical Generalized Linear Models

Description

This function fits hierarchical generalized linear models (HGLMs) using various approximation methods.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
HGLM(y = NULL, X = NULL, Z = NULL, X.disp = NULL, 
     family = gaussian(link = identity), 
     random.family = gaussian(link = identity), method = "EQL", 
     conv = 1e-04, maxit = 20, fixed = NULL, random = NULL, 
     disp = NULL, link.disp = "log", disp.random = NULL, 
     data = NULL, data.random = NULL, fix.disp = NULL, 
     Offset = NULL, Weights = NULL, disp.start = 0, binomial.N = NULL, 
     start.fixed = NULL, start.random = NULL, start.disp = NULL, 
     start.disp.random = NULL, info = TRUE, debug = FALSE, 
     contrasts = NULL)

Arguments

y

the dependent variable, only available when method = 'EQL'.

X

a design matrix for the fixed effects, only available when method = 'EQL'.

Z

an design matrix for the random effects, only available when method = 'EQL'.

X.disp

a design matrix for the fixed effects in the dispersion part of the model, only available when method = 'EQL'.

family

a description of the error distribution and link function to be used in the mean part of the model. (See family for details of family functions.)

random.family

a description of the error distribution and link function to be used in the variance part of the model.

method

estimation method, which can be 'EQL', 'HL01', or 'HL11', where 'EQL' can ONLY be used when ONLY ONE random effect term is specified. 'EQL' is the method of interconnected GLMs presented in Lee et al. (2006), and for 'HL01' and 'HL11', see Lee and Nelder (2001).

conv

convergence criterion, the default is 1e-4, for models with many random effects could be set less strict.

maxit

maximum number of iterations in the IWLS algorithm, only available when method = 'EQL'.

fixed

a formula specifying the fixed effects part of the model, and the format is Response ~ Fixed.Effect.1 + ... + Fixed.Effect.p.

random

a one-sided formula specifying the random effects part of the model, and the format is ~ (Random.Effect.1 | Subject.1) + ... + (Random.Effect.q | Subject.q).

disp

a one-sided formula specifying the fixed effects in the dispersion part of the model, and the format is ~ Effect.1 + ... + Effect.N.

link.disp

the link function for the dispersion part of the model, only available when method = 'EQL'.

disp.random

a list of one-sided formulae for the dispersion strucutre of each random effects, which has the format of list(one = ~ Effect.1.1 + ..., two = ~ Effect.2.1 + ..., three = ..., ...), only available when method = 'HL01' or 'HL11'.

data

the data frame to be used together with fixed and random.

data.random

a list of data.frames for disp.random, which has the format of list(one = data.Random.1, two = data.Random.2,...), only available when method = 'HL01' or 'HL11'.

Weights

prior weights to be specified in weighted regression, only available when method = 'EQL'.

fix.disp

a numeric value if the dispersion parameter of the mean model is known for example 1 for binomial and Poisson models.

Offset

an offset for the linear predictor of the mean model.

disp.start

(starting) values for the overdispersion structure - vector of length equal to the number of parameters in the overdispersion structure, only available when fix.disp = NULL and method = 'HL01' or 'HL11'.

binomial.N

the number of trials for each observation for binomial models.

start.fixed

optional starting values for fixed effects in the mean structure (one vector of numeric values).

start.random

optional starting values for random effects in the mean structure (one vector of numeric values).

start.disp

optional starting values for parameters of dispersion components of the residuals (one vector of numeric values).

start.disp.random

optional starting values for parameters of dispersion components of random effects (one vector of numeric values).

info

a request to display of iteration information if TRUE, only available when method = 'HL01' or 'HL11'.

debug

a request to display of iteration mechanism progress in detail if TRUE, only available when method = 'HL01' or 'HL11'.

contrasts

see lm, caution as it is currently not fully developed, only available when method = 'HL01' or 'HL11'.

Details

When method = 'EQL', all the model checking functions in the hglm-package are available on the object returned; Otherwise, all the model checking functions in the HGLMMM-package are available on the object returned.

Value

When method = 'EQL', an object of class hglm is returned, see hglm; Otherwise, an object of class HGLM is returned, see HGLMfit.

Note

The function provides a unified interface to the hglm-package developed by Moudud Alam, Lars Ronnegard and Xia Shen, and the HGLMMM-package developed by Marek Molas.

Author(s)

Xia Shen and Marek Molas

References

Lee, Y. and Nelder, J.A. (1996). Hierarchical generalized linear models (with discussion). Journal of the Royal Statistical Society. Series B (Methological) 58, 619-678.

Lee, Y. and 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). Generalized Linear Models with Random Effects. Boca Raton: Chapman & Hall/CRC.

Noh, M. and Lee, Y. (2007). REML estimation for binary data in GLMMs. Journal of Multivariate Analysis 98, 896-915.

Ronnegard, L., Shen, X. and Alam, M. (2010). hglm: A Package for Fitting Hierarchical Generalized Linear Models. The R Journal. (to appear)

Molas, M. and Lesaffre, E. (2010). Hierarchical Generalized Linear Models: the R Package HGLMMM. Submitted.

See Also

hglm-package, HGLMMM-package, hglm, HGLMfit.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
data(semiconductor)

# ----- use 'EQL'

h.gamma.normal <- HGLM(fixed = y ~ x1 + x3 + x5 + x6,
                       random = ~ 1|Device,
                       family = Gamma(link = log),
                       disp = ~ x2 + x3, data = semiconductor)
                       
summary(h.gamma.normal)

plot(h.gamma.normal, cex = .6, pch = 1,
     cex.axis = 1/.6, cex.lab = 1/.6,
     cex.main = 1/.6, mar = c(3, 4.5, 0, 1.5))
     
# ----- use 'HL(0,1)'

RSC <- data.frame(int = rep(1, 16))
h.gamma.normal <- HGLM(fixed = y ~ x1 + x3 + x5 + x6,
                       random = ~ 1|Device,
                       family = Gamma(link = log),
                       disp = ~ x2 + x3, data = semiconductor, 
                       method = 'HL01', disp.start = c(0, 0, 0),
                       disp.random = list(one = ~ 1), data.random = list(RSC))
                       
# ----- use 'HL(1,1)'

RSC <- data.frame(int = rep(1, 16))
h.gamma.normal <- HGLM(fixed = y ~ x1 + x3 + x5 + x6,
                       random = ~ 1|Device,
                       family = Gamma(link = log),
                       disp = ~ x2 + x3, data = semiconductor, 
                       method = 'HL11', disp.start = c(0, 0, 0),
                       disp.random = list(one = ~ 1), data.random = list(RSC))