Description Usage Arguments Details Author(s) Examples
Plots residuals for the mean and dispersion models, individual deviances and hatvalues for hglm objects
| 1 2 3 | 
| x | the  | 
| pch | symbol used in the plots | 
| pcol | color of points | 
| lcol | color of lines | 
| device | if  | 
| name | a string gives the main name of the PDF file when  | 
| ... | graphical parameters | 
A S3 generic plot method for hglm objects. It produces a set of diagnostic plots for a hierarchical model.
Xia Shen
| 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)
 | 
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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.