plot.hglm: Plot Hierarchical Generalized Linear Model Objects

Description Usage Arguments Details Author(s) Examples

Description

Plots residuals for the mean and dispersion models, individual deviances and hatvalues for hglm objects

Usage

1
2
3
## S3 method for class 'hglm'
plot(x, pch = "+", pcol = 'slateblue', lcol = 2, 
                    device = NULL, name = NULL, ...)

Arguments

x

the hglm object to be plotted

pch

symbol used in the plots

pcol

color of points

lcol

color of lines

device

if NULL, plot on screen devices, if 'pdf', plot to PDF files in the current working directory.

name

a string gives the main name of the PDF file when device = 'pdf'.

...

graphical parameters

Details

A S3 generic plot method for hglm objects. It produces a set of diagnostic plots for a hierarchical model.

Author(s)

Xia Shen

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
# --------------------- #
# semiconductor example #
# --------------------- #

data(semiconductor)

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))

# ------------------- #
# redo it using hglm2 #
# ------------------- #

m1 <- hglm2(y ~ x1 + x3 + x5 + x6 + (1|Device),
            family = Gamma(link = log),
            disp = ~ x2 + x3, data = semiconductor)
summary(m1)
plot(m1, cex = .6, pch = 1,
     cex.axis = 1/.6, cex.lab = 1/.6,
     cex.main = 1/.6, mar = c(3, 4.5, 0, 1.5))

# --------------------------------------------- #  
# 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)

(m2.1 <- hglm(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1, Z2), 
              RandC = c(10, 5)))
summary(m2.1)
plot(m2.1)

(m2.2 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), data = dd, vcovmat = TRUE))
image(m2.2$vcov)
summary(m2.2)
plot(m2.2)

m3 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), disp = ~ x3, data = dd)
print (m3)
summary(m3)
plot(m3)

## End(Not run)

Example output

Loading required package: Matrix
Loading required package: MASS
Loading required package: hglm.data
Loading required package: sp

hglm: Hierarchical Generalized Linear Models
Version 2.2-1 (2019-04-04) 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.
dev.new(): using pdf(file="Rplots3.pdf")
dev.new(): using pdf(file="Rplots4.pdf")
dev.new(): using pdf(file="Rplots5.pdf")
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="Rplots6.pdf")
dev.new(): using pdf(file="Rplots7.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="Rplots8.pdf")
dev.new(): using pdf(file="Rplots9.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="Rplots10.pdf")
dev.new(): using pdf(file="Rplots11.pdf")
dev.new(): using pdf(file="Rplots12.pdf")

hglm documentation built on May 2, 2019, 6:11 p.m.