hglm2: Fitting Hierarchical Generalized Linear Models
In hglm: Hierarchical Generalized Linear Models
Description
Usage
Arguments
Details
Value
Author(s)
References
See Also
Examples
Description
hglm2 is used to fit hierarchical generalized linear models. hglm2 is used to fit hierarchical generalized linear models.
It extends the hglm function by allowing for several random effects, where the model is specified in lme4 convension,
and also by implementing sparse matrix techniques using the Matrix library.
Usage
1
2
3
4
5
6
7
hglm2(meanmodel, data = NULL, family = gaussian(link = identity),
rand.family = gaussian(link = identity), method = "EQL",
conv = 1e-6, maxit = 50, startval = NULL,
X.disp = NULL, disp = NULL, link.disp = "log",
weights = NULL, fix.disp = NULL, offset = NULL,
sparse = TRUE, vcovmat = FALSE, calc.like = FALSE,
RandC = NULL, bigRR = FALSE, verbose = FALSE, ...)
Arguments
meanmodel
formula. A two sided formula specifying the fixed and random terms in lme4 convention, e.g. y ~ x1 + (1|id) indicates y as response, x1 as the fixed effect and (1|id) represent a random intercept for each level of id.
data
data.frame. An optional data frame from where the variables in the meanmodel (and possibly disp) are to be obtained. It is expected that the data frame does not contain any missing value.
family
family. The 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.)
rand.family
family. The description of the distribution and link function to be used for the random effect.
method
character. Estimation method where EQL is the method of interconnected GLMs presented in Lee et al (2006). Apart from the default option EQL there is also an EQL1 option, which improves estimation for GLMMs (especially for Poisson models with a large number of levels in the random effects).
conv
numeric. The convergence criteria (change in linear predictor between iterations).
maxit
numeric. Maximum number of iterations in the hglm algorithm.
startval
numeric. A vector of starting values in the following order: fixed effects, random effect, variance of random effects, variance of residuals.
X.disp
matrix. The design matrix for the fixed effects in the dispersion part of the model.
disp
formula. A one-sided formula specifying the fixed effects in the dispersion part of the model.
link.disp
character. The link function for the dispersion part of the model.
weights
numeric. Prior weights to be specified in weighted regression.
fix.disp
numeric. A numeric value if the dispersion parameter of the mean model is known, e.g., 1 for binomial and Poisson model.
offset
An offset for the linear predictor of the mean model.
sparse
logical. If TRUE, the computation is to be carried out by using sparse matrix technique.
vcovmat
logical. If TRUE, the variance-covariance matrix is exported.
calc.like
logical. If TRUE, likelihoods will be computed at convergence and will be shown via the print or summary methods on the output object.
RandC
numeric. Necessary in old versions but can be neglected now. Integers (possibly a vector) specifying the number of column of Z to be used for each of the random-effect terms.
bigRR
logical. If TRUE, and only for the Gaussian model with one random effect term, a specific algorithm will be used for fast fitting high-dimensional (p >> n) problems. See Shen et al. (2013) for more details of the method.
verbose
logical. If TRUE, more information is printed during model fitting process.
...
not used.
Details
Models for hglm are either specified symbolically using formula
or by specifying the design matrices ( X, Z and X.disp). Currently, only the extended quasi likelihood (EQL)
method is available for the estimation of the model parameters. Only for the Gaussian-Gaussina linear mixed models, it
is REML. It should be noted that the EQL estimator can be biased and inconsistent in some special cases e.g. binary pair matched response. A higher order correction
might be useful to correct the bias of EQL (Lee et al. 2006). But, those currections are not implemented in the current version. By default, the dispersion
parameter is always estimated via EQL. If the dispersion parameter of the mean model is to be held constant, for example if it is
desired to be 1 for binomial and Poisson family, then fix.disp=value where, value=1 for the above example, should be used.
Value
It returns an object of class hglm consiting of the following values.
fixef
fixed effect estimates.
ranef
random effect estimates.
RandC
integers (possibly a vector) specified the number of column of Z to be used for each of the random-effect terms.
varFix
dispersion parameter of the mean model (residual variance for LMM).
varRanef
dispersion parameter of the random effects (variance of random effects for GLMM).
iter
number of iterations used.
Converge
specifies if the algorithm converged.
SeFe
standard errors of fixed effects.
SeRe
standard errors of random effects.
dfReFe
deviance degrees of freedom for the mean part of the model.
SummVC1
estimates and standard errors of the linear predictor in the dispersion model.
SummVC2
estimates and standard errors of the linear predictor for the dispersion parameter of the random effects.
dev
individual deviances for the mean part of the model.
hv
hatvalues for the mean part of the model.
resid
studentized residuals for the mean part of the model.
fv
fitted values for the mean part of the model.
disp.fv
fitted values for the dispersion part of the model.
disp.resid
standardized deviance residuals for the dispersion part of the model.
link.disp
link function for the dispersion part of the model.
vcov
the variance-covariance matrix.
likelihood
a list of log-likelihood values for model selection purposes, where $hlik is -2 times the log-h-likelihood, $pvh -2 times the adjusted profile log-likelihood profiled over random effects, $pbvh -2 times the adjusted profile log-likelihood profiled over fixed and random effects, and $cAIC the conditional AIC.
bad
the index of the influential observation.
Author(s)
Moudud Alam, Xia Shen, Lars Ronnegard
References
Lars Ronnegard, Xia Shen and Moudud Alam (2010). hglm: A Package for Fitting Hierarchical Generalized Linear Models. The R Journal, 2(2), 20-28.
Youngjo Lee, John A Nelder and Yudi Pawitan (2006) Generalized Linear Models with Random Effect: a unified analysis via h-likelihood. Chapman and Hall/CRC.
Xia Shen, Moudud Alam, Freddy Fikse and Lars Ronnegard (2013). A novel generalized ridge regression method for quantitative genetics. Genetics.
Moudud Alam, Lars Ronnegard, Xia Shen (2014). Fitting conditional and simultaneous autoregressive spatial models in hglm. Submitted.
See Also
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
# Find more examples and instructions in the package vignette:
# vignette('hglm')
require(hglm)
# --------------------- #
# semiconductor example #
# --------------------- #
data(semiconductor)
m11 <- hglm(fixed = y ~ x1 + x3 + x5 + x6,
random = ~ 1|Device,
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor)
summary(m11)
plot(m11, cex = .6, pch = 1,
cex.axis = 1/.6, cex.lab = 1/.6,
cex.main = 1/.6, mar = c(3, 4.5, 0, 1.5))
# ------------------- #
# redo it using hglm2 #
# ------------------- #
m12 <- hglm2(y ~ x1 + x3 + x5 + x6 + (1|Device),
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor)
summary(m12)
# -------------------------- #
# redo it using matrix input #
# -------------------------- #
attach(semiconductor)
m13 <- hglm(y = y, X = model.matrix(~ x1 + x3 + x5 + x6),
Z = kronecker(diag(16), rep(1, 4)),
X.disp = model.matrix(~ x2 + x3),
family = Gamma(link = log))
summary(m13)
# --------------------- #
# verbose & likelihoods #
# --------------------- #
m14 <- hglm(fixed = y ~ x1 + x3 + x5 + x6,
random = ~ 1|Device,
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor,
verbose = TRUE, calc.like = TRUE)
summary(m14)
# --------------------------------------------- #
# simulated example with 2 random effects terms #
# --------------------------------------------- #
## Not run:
set.seed(911)
x1 <- rnorm(100)
x2 <- rnorm(100)
x3 <- rnorm(100)
z1 <- factor(rep(LETTERS[1:10], rep(10, 10)))
z2 <- factor(rep(letters[1:5], rep(20, 5)))
Z1 <- model.matrix(~ 0 + z1)
Z2 <- model.matrix(~ 0 + z2)
u1 <- rnorm(10, 0, sqrt(2))
u2 <- rnorm(5, 0, sqrt(3))
y <- 1 + 2*x1 + 3*x2 + Z1%*%u1 + Z2%*%u2 + rnorm(100, 0, sqrt(exp(x3)))
dd <- data.frame(x1 = x1, x2 = x2, x3 = x3, z1 = z1, z2 = z2, y = y)
(m21 <- hglm(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1, Z2),
RandC = c(10, 5)))
summary(m21)
plot(m21)
(m22 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), data = dd, vcovmat = TRUE))
image(m22$vcov, main = 'Variance-covariance Matrix')
summary(m22)
plot(m22)
m31 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), disp = ~ x3, data = dd)
print (m31)
summary(m31)
plot(m31)
## End(Not run)
Example output
Loading required package: Matrix
Loading required package: MASS
Loading required package: hglm.data
hglm: Hierarchical Generalized Linear Models
Version 2.1-1 (2015-08-28) installed
Authors: Moudud Alam, Lars Ronnegard, Xia Shen
Maintainer: Xia Shen <xia.shen@ki.se>
Use citation("hglm") to know how to cite our work.
Discussion: https://r-forge.r-project.org/forum/?group_id=558
BugReports: https://r-forge.r-project.org/tracker/?group_id=558
VideoTutorials: http://www.youtube.com/playlist?list=PLn1OmZECD-n15vnYzvJDy5GxjNpVV5Jr8
Call:
hglm.formula(family = Gamma(link = log), fixed = y ~ x1 + x3 +
x5 + x6, random = ~1 | Device, disp = ~x2 + x3, data = semiconductor)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
as.factor(Device)1 0.2724 0.1787
as.factor(Device)2 0.0097 0.1787
as.factor(Device)3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
dev.new(): using pdf(file="Rplots1.pdf")
dev.new(): using pdf(file="Rplots2.pdf")
Call:
hglm2.formula(meanmodel = y ~ x1 + x3 + x5 + x6 + (1 | Device),
data = semiconductor, family = Gamma(link = log), disp = ~x2 +
x3)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| Device:1 0.2724 0.1787
(Intercept)| Device:2 0.0097 0.1787
(Intercept)| Device:3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
Call:
hglm.default(X = model.matrix(~x1 + x3 + x5 + x6), y = y, Z = kronecker(diag(16),
rep(1, 4)), family = Gamma(link = log), X.disp = model.matrix(~x2 +
x3))
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
Z.1 0.2724 0.1787
Z.2 0.0097 0.1787
Z.3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
-------------------
Iteration 1
-------------------
Sum of squared linear predictor: 1437.809
Convergence: 0.001198849
-------------------
Iteration 2
-------------------
Sum of squared linear predictor: 1441.496
Convergence: 9.225213e-05
-------------------
Iteration 3
-------------------
Sum of squared linear predictor: 1441.963
Convergence: 1.078477e-06
-------------------
Iteration 4
-------------------
Sum of squared linear predictor: 1442.064
Convergence: 8.388541e-08
Call:
hglm.formula(family = Gamma(link = log), fixed = y ~ x1 + x3 +
x5 + x6, random = ~1 | Device, disp = ~x2 + x3, data = semiconductor,
calc.like = TRUE, verbose = TRUE)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
as.factor(Device)1 0.2724 0.1787
as.factor(Device)2 0.0097 0.1787
as.factor(Device)3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
---------------
LOG-LIKELIHOODS
---------------
h-likelihood: 307.1976
Adjusted profile likelihood
Profiled over random effects: 286.9197
Profiled over fixed and random effects: 277.9467
Conditional AIC: -577.9375
EQL estimation converged in 4 iterations.
Call:
hglm.default(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1,
Z2), RandC = c(10, 5))
---------------------------
Estimates of the mean model
---------------------------
Fixed effects:
x1 x2
1.593223 2.163230 2.875127
Random effects:
z1A z1B z1C
0.5876582 0.7925926 -1.3893579
...
z1I z1J
-1.071900 1.615409
NOTE: to show all the random effects estimates, use print(hglm.object, print.ranef = TRUE).
Random effects:
z2a z2b z2c z2d z2e
1.0972865 -1.2008286 0.5139767 -0.8425195 0.4320849
Dispersion parameter for the mean model: 1.574819
Dispersion parameter for the random effects: 2.079812 1.597671
Estimation converged in 4 iterations
Call:
hglm.default(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1,
Z2), RandC = c(10, 5))
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
1.5932 0.7410 2.15 0.0343 *
x1 2.1632 0.1321 16.38 <2e-16 ***
x2 2.8751 0.1210 23.77 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 88 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
z1A 0.5877 0.8920
z1B 0.7926 0.8911
z1C -1.3894 0.8918
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
z2a 1.0973 0.9263
z2b -1.2008 0.9254
z2c 0.5140 0.9254
z2d -0.8425 0.9255
z2e 0.4321 0.9252
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Dispersion parameter for the mean model:
[1] 1.574819
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
0.4541 0.1503
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 2.080 1.598
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.7323 0.5707
.|Random2
Estimate Std. Error
0.4685 0.9162
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
dev.new(): using pdf(file="Rplots3.pdf")
dev.new(): using pdf(file="Rplots4.pdf")
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, vcovmat = TRUE)
---------------------------
Estimates of the mean model
---------------------------
Fixed effects:
(Intercept) x1 x2
1.593223 2.163230 2.875127
Random effects:
(Intercept)| z1:A (Intercept)| z1:B (Intercept)| z1:C
0.5876582 0.7925926 -1.3893579
...
(Intercept)| z1:I (Intercept)| z1:J
-1.071900 1.615409
NOTE: to show all the random effects estimates, use print(hglm.object, print.ranef = TRUE).
Random effects:
(Intercept)| z2:a (Intercept)| z2:b (Intercept)| z2:c (Intercept)| z2:d
1.0972865 -1.2008286 0.5139767 -0.8425195
(Intercept)| z2:e
0.4320849
Dispersion parameter for the mean model: 1.574819
Dispersion parameter for the random effects: 2.079812 1.597671
Estimation converged in 4 iterations
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, vcovmat = TRUE)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.5932 0.7410 2.15 0.0343 *
x1 2.1632 0.1321 16.38 <2e-16 ***
x2 2.8751 0.1210 23.77 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 88 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z1:A 0.5877 0.8920
(Intercept)| z1:B 0.7926 0.8911
(Intercept)| z1:C -1.3894 0.8918
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z2:a 1.0973 0.9263
(Intercept)| z2:b -1.2008 0.9254
(Intercept)| z2:c 0.5140 0.9254
(Intercept)| z2:d -0.8425 0.9255
(Intercept)| z2:e 0.4321 0.9252
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Dispersion parameter for the mean model:
[1] 1.574819
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
0.4541 0.1503
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 2.080 1.598
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.7323 0.5707
.|Random2
Estimate Std. Error
0.4685 0.9162
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
dev.new(): using pdf(file="Rplots5.pdf")
dev.new(): using pdf(file="Rplots6.pdf")
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, disp = ~x3)
---------------------------
Estimates of the mean model
---------------------------
Fixed effects:
(Intercept) x1 x2
1.676184 2.148076 2.989964
Random effects:
(Intercept)| z1:A (Intercept)| z1:B (Intercept)| z1:C
0.1914595 1.3496934 -1.2882400
...
(Intercept)| z1:I (Intercept)| z1:J
-1.327242 1.596106
NOTE: to show all the random effects estimates, use print(hglm.object, print.ranef = TRUE).
Random effects:
(Intercept)| z2:a (Intercept)| z2:b (Intercept)| z2:c (Intercept)| z2:d
0.8394165 -0.9463058 0.5146495 -0.5542017
(Intercept)| z2:e
0.1464415
------------------------------------------
Estimates of the residual dispersion model
------------------------------------------
Link = log
Effects:
(Intercept) x3
-0.1065797 0.9969769
Dispersion parameter for the random effects: 2.041197 1.067874
Estimation converged in 6 iterations
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, disp = ~x3)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.67618 0.65467 2.56 0.0122 *
x1 2.14808 0.08218 26.14 <2e-16 ***
x2 2.98996 0.07374 40.55 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 88 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z1:A 0.1915 0.8087
(Intercept)| z1:B 1.3497 0.8047
(Intercept)| z1:C -1.2882 0.8179
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z2:a 0.8394 0.8051
(Intercept)| z2:b -0.9463 0.8055
(Intercept)| z2:c 0.5146 0.8051
(Intercept)| z2:d -0.5542 0.8049
(Intercept)| z2:e 0.1464 0.8070
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -0.1066 0.1516
x3 0.9970 0.1485
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 2.041 1.068
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.7135 0.5449
.|Random2
Estimate Std. Error
0.0657 0.9867
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 6 iterations.
dev.new(): using pdf(file="Rplots7.pdf")
dev.new(): using pdf(file="Rplots8.pdf")
dev.new(): using pdf(file="Rplots9.pdf")
hglm documentation built on May 2, 2019, 7:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Author(s) References See Also Examples
Description
hglm2 is used to fit hierarchical generalized linear models. hglm2 is used to fit hierarchical generalized linear models.
It extends the hglm function by allowing for several random effects, where the model is specified in lme4 convension,
and also by implementing sparse matrix techniques using the Matrix library.
Usage
1 2 3 4 5 6 7 | hglm2(meanmodel, data = NULL, family = gaussian(link = identity),
rand.family = gaussian(link = identity), method = "EQL",
conv = 1e-6, maxit = 50, startval = NULL,
X.disp = NULL, disp = NULL, link.disp = "log",
weights = NULL, fix.disp = NULL, offset = NULL,
sparse = TRUE, vcovmat = FALSE, calc.like = FALSE,
RandC = NULL, bigRR = FALSE, verbose = FALSE, ...)
|
Arguments
meanmodel |
|
data |
|
family |
|
rand.family |
|
method |
|
conv |
|
maxit |
|
startval |
|
X.disp |
|
disp |
|
link.disp |
|
weights |
|
fix.disp |
|
offset |
An offset for the linear predictor of the mean model. |
sparse |
|
vcovmat |
logical. If |
calc.like |
logical. If |
RandC |
|
bigRR |
logical. If |
verbose |
logical. If |
... |
not used. |
Details
Models for hglm are either specified symbolically using formula
or by specifying the design matrices ( X, Z and X.disp). Currently, only the extended quasi likelihood (EQL)
method is available for the estimation of the model parameters. Only for the Gaussian-Gaussina linear mixed models, it
is REML. It should be noted that the EQL estimator can be biased and inconsistent in some special cases e.g. binary pair matched response. A higher order correction
might be useful to correct the bias of EQL (Lee et al. 2006). But, those currections are not implemented in the current version. By default, the dispersion
parameter is always estimated via EQL. If the dispersion parameter of the mean model is to be held constant, for example if it is
desired to be 1 for binomial and Poisson family, then fix.disp=value where, value=1 for the above example, should be used.
Value
It returns an object of class hglm consiting of the following values.
fixef |
fixed effect estimates. |
ranef |
random effect estimates. |
RandC |
integers (possibly a vector) specified the number of column of Z to be used for each of the random-effect terms. |
varFix |
dispersion parameter of the mean model (residual variance for LMM). |
varRanef |
dispersion parameter of the random effects (variance of random effects for GLMM). |
iter |
number of iterations used. |
Converge |
specifies if the algorithm converged. |
SeFe |
standard errors of fixed effects. |
SeRe |
standard errors of random effects. |
dfReFe |
deviance degrees of freedom for the mean part of the model. |
SummVC1 |
estimates and standard errors of the linear predictor in the dispersion model. |
SummVC2 |
estimates and standard errors of the linear predictor for the dispersion parameter of the random effects. |
dev |
individual deviances for the mean part of the model. |
hv |
hatvalues for the mean part of the model. |
resid |
studentized residuals for the mean part of the model. |
fv |
fitted values for the mean part of the model. |
disp.fv |
fitted values for the dispersion part of the model. |
disp.resid |
standardized deviance residuals for the dispersion part of the model. |
link.disp |
link function for the dispersion part of the model. |
vcov |
the variance-covariance matrix. |
likelihood |
a list of log-likelihood values for model selection purposes, where |
bad |
the index of the influential observation. |
Author(s)
Moudud Alam, Xia Shen, Lars Ronnegard
References
Lars Ronnegard, Xia Shen and Moudud Alam (2010). hglm: A Package for Fitting Hierarchical Generalized Linear Models. The R Journal, 2(2), 20-28.
Youngjo Lee, John A Nelder and Yudi Pawitan (2006) Generalized Linear Models with Random Effect: a unified analysis via h-likelihood. Chapman and Hall/CRC.
Xia Shen, Moudud Alam, Freddy Fikse and Lars Ronnegard (2013). A novel generalized ridge regression method for quantitative genetics. Genetics.
Moudud Alam, Lars Ronnegard, Xia Shen (2014). Fitting conditional and simultaneous autoregressive spatial models in hglm. Submitted.
See Also
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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 | # Find more examples and instructions in the package vignette:
# vignette('hglm')
require(hglm)
# --------------------- #
# semiconductor example #
# --------------------- #
data(semiconductor)
m11 <- hglm(fixed = y ~ x1 + x3 + x5 + x6,
random = ~ 1|Device,
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor)
summary(m11)
plot(m11, cex = .6, pch = 1,
cex.axis = 1/.6, cex.lab = 1/.6,
cex.main = 1/.6, mar = c(3, 4.5, 0, 1.5))
# ------------------- #
# redo it using hglm2 #
# ------------------- #
m12 <- hglm2(y ~ x1 + x3 + x5 + x6 + (1|Device),
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor)
summary(m12)
# -------------------------- #
# redo it using matrix input #
# -------------------------- #
attach(semiconductor)
m13 <- hglm(y = y, X = model.matrix(~ x1 + x3 + x5 + x6),
Z = kronecker(diag(16), rep(1, 4)),
X.disp = model.matrix(~ x2 + x3),
family = Gamma(link = log))
summary(m13)
# --------------------- #
# verbose & likelihoods #
# --------------------- #
m14 <- hglm(fixed = y ~ x1 + x3 + x5 + x6,
random = ~ 1|Device,
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor,
verbose = TRUE, calc.like = TRUE)
summary(m14)
# --------------------------------------------- #
# simulated example with 2 random effects terms #
# --------------------------------------------- #
## Not run:
set.seed(911)
x1 <- rnorm(100)
x2 <- rnorm(100)
x3 <- rnorm(100)
z1 <- factor(rep(LETTERS[1:10], rep(10, 10)))
z2 <- factor(rep(letters[1:5], rep(20, 5)))
Z1 <- model.matrix(~ 0 + z1)
Z2 <- model.matrix(~ 0 + z2)
u1 <- rnorm(10, 0, sqrt(2))
u2 <- rnorm(5, 0, sqrt(3))
y <- 1 + 2*x1 + 3*x2 + Z1%*%u1 + Z2%*%u2 + rnorm(100, 0, sqrt(exp(x3)))
dd <- data.frame(x1 = x1, x2 = x2, x3 = x3, z1 = z1, z2 = z2, y = y)
(m21 <- hglm(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1, Z2),
RandC = c(10, 5)))
summary(m21)
plot(m21)
(m22 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), data = dd, vcovmat = TRUE))
image(m22$vcov, main = 'Variance-covariance Matrix')
summary(m22)
plot(m22)
m31 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), disp = ~ x3, data = dd)
print (m31)
summary(m31)
plot(m31)
## End(Not run)
|
Example output
Loading required package: Matrix
Loading required package: MASS
Loading required package: hglm.data
hglm: Hierarchical Generalized Linear Models
Version 2.1-1 (2015-08-28) installed
Authors: Moudud Alam, Lars Ronnegard, Xia Shen
Maintainer: Xia Shen <xia.shen@ki.se>
Use citation("hglm") to know how to cite our work.
Discussion: https://r-forge.r-project.org/forum/?group_id=558
BugReports: https://r-forge.r-project.org/tracker/?group_id=558
VideoTutorials: http://www.youtube.com/playlist?list=PLn1OmZECD-n15vnYzvJDy5GxjNpVV5Jr8
Call:
hglm.formula(family = Gamma(link = log), fixed = y ~ x1 + x3 +
x5 + x6, random = ~1 | Device, disp = ~x2 + x3, data = semiconductor)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
as.factor(Device)1 0.2724 0.1787
as.factor(Device)2 0.0097 0.1787
as.factor(Device)3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
dev.new(): using pdf(file="Rplots1.pdf")
dev.new(): using pdf(file="Rplots2.pdf")
Call:
hglm2.formula(meanmodel = y ~ x1 + x3 + x5 + x6 + (1 | Device),
data = semiconductor, family = Gamma(link = log), disp = ~x2 +
x3)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| Device:1 0.2724 0.1787
(Intercept)| Device:2 0.0097 0.1787
(Intercept)| Device:3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
Call:
hglm.default(X = model.matrix(~x1 + x3 + x5 + x6), y = y, Z = kronecker(diag(16),
rep(1, 4)), family = Gamma(link = log), X.disp = model.matrix(~x2 +
x3))
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
Z.1 0.2724 0.1787
Z.2 0.0097 0.1787
Z.3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
-------------------
Iteration 1
-------------------
Sum of squared linear predictor: 1437.809
Convergence: 0.001198849
-------------------
Iteration 2
-------------------
Sum of squared linear predictor: 1441.496
Convergence: 9.225213e-05
-------------------
Iteration 3
-------------------
Sum of squared linear predictor: 1441.963
Convergence: 1.078477e-06
-------------------
Iteration 4
-------------------
Sum of squared linear predictor: 1442.064
Convergence: 8.388541e-08
Call:
hglm.formula(family = Gamma(link = log), fixed = y ~ x1 + x3 +
x5 + x6, random = ~1 | Device, disp = ~x2 + x3, data = semiconductor,
calc.like = TRUE, verbose = TRUE)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -4.71168 0.06696 -70.368 < 2e-16 ***
x1 0.20979 0.06638 3.160 0.00263 **
x3 0.32893 0.06696 4.913 9.34e-06 ***
x5 -0.17314 0.06638 -2.608 0.01185 *
x6 -0.35690 0.06633 -5.380 1.80e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 52 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
as.factor(Device)1 0.2724 0.1787
as.factor(Device)2 0.0097 0.1787
as.factor(Device)3 -0.2697 0.1584
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -2.5887 0.1972
x2 -0.6861 0.1971
x3 -0.5024 0.1971
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 0.0486
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
-3.0242 0.5172
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
---------------
LOG-LIKELIHOODS
---------------
h-likelihood: 307.1976
Adjusted profile likelihood
Profiled over random effects: 286.9197
Profiled over fixed and random effects: 277.9467
Conditional AIC: -577.9375
EQL estimation converged in 4 iterations.
Call:
hglm.default(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1,
Z2), RandC = c(10, 5))
---------------------------
Estimates of the mean model
---------------------------
Fixed effects:
x1 x2
1.593223 2.163230 2.875127
Random effects:
z1A z1B z1C
0.5876582 0.7925926 -1.3893579
...
z1I z1J
-1.071900 1.615409
NOTE: to show all the random effects estimates, use print(hglm.object, print.ranef = TRUE).
Random effects:
z2a z2b z2c z2d z2e
1.0972865 -1.2008286 0.5139767 -0.8425195 0.4320849
Dispersion parameter for the mean model: 1.574819
Dispersion parameter for the random effects: 2.079812 1.597671
Estimation converged in 4 iterations
Call:
hglm.default(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1,
Z2), RandC = c(10, 5))
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
1.5932 0.7410 2.15 0.0343 *
x1 2.1632 0.1321 16.38 <2e-16 ***
x2 2.8751 0.1210 23.77 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 88 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
z1A 0.5877 0.8920
z1B 0.7926 0.8911
z1C -1.3894 0.8918
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
z2a 1.0973 0.9263
z2b -1.2008 0.9254
z2c 0.5140 0.9254
z2d -0.8425 0.9255
z2e 0.4321 0.9252
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Dispersion parameter for the mean model:
[1] 1.574819
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
0.4541 0.1503
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 2.080 1.598
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.7323 0.5707
.|Random2
Estimate Std. Error
0.4685 0.9162
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
dev.new(): using pdf(file="Rplots3.pdf")
dev.new(): using pdf(file="Rplots4.pdf")
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, vcovmat = TRUE)
---------------------------
Estimates of the mean model
---------------------------
Fixed effects:
(Intercept) x1 x2
1.593223 2.163230 2.875127
Random effects:
(Intercept)| z1:A (Intercept)| z1:B (Intercept)| z1:C
0.5876582 0.7925926 -1.3893579
...
(Intercept)| z1:I (Intercept)| z1:J
-1.071900 1.615409
NOTE: to show all the random effects estimates, use print(hglm.object, print.ranef = TRUE).
Random effects:
(Intercept)| z2:a (Intercept)| z2:b (Intercept)| z2:c (Intercept)| z2:d
1.0972865 -1.2008286 0.5139767 -0.8425195
(Intercept)| z2:e
0.4320849
Dispersion parameter for the mean model: 1.574819
Dispersion parameter for the random effects: 2.079812 1.597671
Estimation converged in 4 iterations
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, vcovmat = TRUE)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.5932 0.7410 2.15 0.0343 *
x1 2.1632 0.1321 16.38 <2e-16 ***
x2 2.8751 0.1210 23.77 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 88 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z1:A 0.5877 0.8920
(Intercept)| z1:B 0.7926 0.8911
(Intercept)| z1:C -1.3894 0.8918
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z2:a 1.0973 0.9263
(Intercept)| z2:b -1.2008 0.9254
(Intercept)| z2:c 0.5140 0.9254
(Intercept)| z2:d -0.8425 0.9255
(Intercept)| z2:e 0.4321 0.9252
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Dispersion parameter for the mean model:
[1] 1.574819
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
0.4541 0.1503
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 2.080 1.598
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.7323 0.5707
.|Random2
Estimate Std. Error
0.4685 0.9162
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 4 iterations.
dev.new(): using pdf(file="Rplots5.pdf")
dev.new(): using pdf(file="Rplots6.pdf")
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, disp = ~x3)
---------------------------
Estimates of the mean model
---------------------------
Fixed effects:
(Intercept) x1 x2
1.676184 2.148076 2.989964
Random effects:
(Intercept)| z1:A (Intercept)| z1:B (Intercept)| z1:C
0.1914595 1.3496934 -1.2882400
...
(Intercept)| z1:I (Intercept)| z1:J
-1.327242 1.596106
NOTE: to show all the random effects estimates, use print(hglm.object, print.ranef = TRUE).
Random effects:
(Intercept)| z2:a (Intercept)| z2:b (Intercept)| z2:c (Intercept)| z2:d
0.8394165 -0.9463058 0.5146495 -0.5542017
(Intercept)| z2:e
0.1464415
------------------------------------------
Estimates of the residual dispersion model
------------------------------------------
Link = log
Effects:
(Intercept) x3
-0.1065797 0.9969769
Dispersion parameter for the random effects: 2.041197 1.067874
Estimation converged in 6 iterations
Call:
hglm2.formula(meanmodel = y ~ x1 + x2 + (1 | z1) + (1 | z2),
data = dd, disp = ~x3)
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.67618 0.65467 2.56 0.0122 *
x1 2.14808 0.08218 26.14 <2e-16 ***
x2 2.98996 0.07374 40.55 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 88 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z1:A 0.1915 0.8087
(Intercept)| z1:B 1.3497 0.8047
(Intercept)| z1:C -1.2882 0.8179
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
(Intercept)| z2:a 0.8394 0.8051
(Intercept)| z2:b -0.9463 0.8055
(Intercept)| z2:c 0.5146 0.8051
(Intercept)| z2:d -0.5542 0.8049
(Intercept)| z2:e 0.1464 0.8070
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
(Intercept) -0.1066 0.1516
x3 0.9970 0.1485
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 2.041 1.068
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.7135 0.5449
.|Random2
Estimate Std. Error
0.0657 0.9867
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 6 iterations.
dev.new(): using pdf(file="Rplots7.pdf")
dev.new(): using pdf(file="Rplots8.pdf")
dev.new(): using pdf(file="Rplots9.pdf")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.