############################################################################ ## Code by Vanja Dukic, U of Colorado at Boulder ## Bayesian Statistics and Computing 2019 ## ## LECTURE 12 - ## GLMM - frequentist approaches ############################################### install.packages("hglm", repos="http://R-Forge.R-project.org") require("hglm") # From the HGLM authors: 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. # # Generalized linear mixed models (GLMM) have previously been implemented in several R (R # Development Core Team 2009) function, such as the glmer() function in the lme4 library # and in the glmmPQL() function in the MASS library. In GLMM, the random effects are # assumed to be Gaussian whereas the hglm() function allow for other distributions for the # random effect. The hglm() function also extends the fitting algorithm of Gordon Smyth’s # dglm package by including random effects in the linear predictor for the mean. # Example: # Pump reliability data set from Gaver and O'Muircheartaigh (1987) # Robust Empirical Bayes Analyses of Event Rates, Technometrics 29(1), 1-15 # # The number of failures for 10 pumps in a neuclear plant: # System: The system number. # S: Number of pumps failures. # t: Time (in thousand hours) of operation. # Gr: Pump groups; two levels: 1 = operated continuously, # 0 = operated intermittently. # # data(pump) pump # System S t Gr #1 1 5 94.320 1 #2 2 1 15.720 0 #3 3 5 62.880 1 #4 4 14 125.760 1 #5 5 3 5.240 0 #6 6 19 31.440 1 #7 7 1 1.048 0 #8 8 1 1.048 0 #9 9 4 2.096 0 #10 10 22 10.480 0 #################### Normally distirbuted random effects pump.model.poinorm <- hglm(fixed = S ~ factor(Gr), random = ~1 | System, offset = log(t), fix.disp = 1, family = poisson(), rand.family = gaussian(link = identity), data = pump) print(pump.model.poinorm, print.ranef=TRUE) #Fixed effects: #(Intercept) factor(Gr)1 # -0.2431341 -1.7609144 # #Random effects: # as.factor(System)1 as.factor(System)2 as.factor(System)3 as.factor(System)4 # -0.8005086 -1.6210937 -0.4508527 -0.1799417 # as.factor(System)5 as.factor(System)6 as.factor(System)7 as.factor(System)8 # -0.2447570 1.4313030 0.1017428 0.1017428 # as.factor(System)9 as.factor(System)10 # 0.7165092 0.9458558 #################### Compare to GLMER -- also normally distirbuted random effects pump.model.poinorm.glmer = glmer(S ~ factor(Gr) + offset(log(t)) + (1 | System), family = poisson, data=pump) summary(pump.model.poinorm.glmer) #Generalized linear mixed model fit by maximum likelihood (Laplace # AIC BIC logLik deviance df.resid # 65.8 66.7 -29.9 59.8 7 # #Random effects: # Groups Name Variance Std.Dev. # System (Intercept) 0.9019 0.9497 # #Fixed effects: # Estimate Std. Error z value Pr(>|z|) #(Intercept) -0.3596 0.4808 -0.748 0.4545 #factor(Gr)1 -1.6877 0.6961 -2.425 0.0153 * # #Correlation of Fixed Effects: # (Intr) #factor(Gr)1 -0.686 # Compare random effect estimates cbind(pump.model.poinorm$ranef, random.effects(pump.model.poinorm.glmer)$System) # pump.model.poinorm$ranef (Intercept) # as.factor(System)1 -0.8005086 -0.7383270 # as.factor(System)2 -1.6210937 -1.4407563 # as.factor(System)3 -0.4508527 -0.3995894 # as.factor(System)4 -0.1799417 -0.1372133 # as.factor(System)5 -0.2447570 -0.1456537 # as.factor(System)6 1.4313030 1.4549246 # as.factor(System)7 0.1017428 0.1417471 # as.factor(System)8 0.1017428 0.1417471 # as.factor(System)9 0.7165092 0.7668877 # as.factor(System)10 0.9458558 1.0469829 # Compare fitted values cbind(pump.model.poinorm$fv, fitted(pump.model.poinorm.glmer)) # [,1] [,2] #1 5.7095270 5.818661 #2 2.4368493 2.597519 #3 5.3996111 5.443067 #4 14.1594905 14.152143 #5 3.2169392 3.161502 #6 17.7313717 17.386771 #7 0.9098207 0.842830 #8 0.9098207 0.842830 #9 3.3649256 3.149670 #10 21.1616453 20.839099 # RESIDUALS on the scale of the data cbind(pump$S,pump$S)-cbind(pump.model.poinorm$fv, fitted(pump.model.poinorm.glmer)) [,1] [,2] 1 -0.70952699 -0.8186614 2 -1.43684928 -1.5975192 3 -0.39961114 -0.4430672 4 -0.15949048 -0.1521429 5 -0.21693921 -0.1615017 6 1.26862834 1.6132291 7 0.09017927 0.1571700 8 0.09017927 0.1571700 9 0.63507437 0.8503296 10 0.83835466 1.1609009 # SUM OF SQUARED RESIDUALS colSums(cbind(pump$S,pump$S)-cbind(pump.model.poinorm$fv, fitted(pump.model.poinorm.glmer)))^2 # [1] 1.383842e-12 5.866138e-01 ############################################################################################### # Poisson model with Gamma distributed random effects # # From HGLM authors: # E(Yi | β, u) = exp(Xi β + Zi v) # where level j in the random effect v is given by vj = log( uj ) and uj are iid with gamma # distribution having mean and variance: E(uj ) = 1, var(uj ) = λ. pump.model.poigam <- hglm(fixed = S ~ factor(Gr), random = ~1 | System, offset = log(t), fix.disp = 1, family = poisson(), rand.family = Gamma(link = log), data = pump) print(pump.model.poigam, print.ranef=TRUE) ## Fixed effects: ## (Intercept) factor(Gr)1 ## 0.07479 -1.66527 ## ## Random effects: pump.model.poigam$ranef # # as.factor(System)1 as.factor(System)2 as.factor(System)3 as.factor(System)4 # 0.2950977 0.1092488 0.4324137 0.5624661 # as.factor(System)5 as.factor(System)6 as.factor(System)7 as.factor(System)8 # 0.5990728 2.7100225 0.9379283 0.9379283 # as.factor(System)9 as.factor(System)10 # 1.5417814 1.8740406 cbind(exp(pump.model.poinorm$ranef),pump.model.poigam$ranef) # [,1] [,2] #as.factor(System)1 0.4491005 0.2950977 #as.factor(System)2 0.1976824 0.1092488 #as.factor(System)3 0.6370847 0.4324137 #as.factor(System)4 0.8353189 0.5624661 #as.factor(System)5 0.7828947 0.5990728 #as.factor(System)6 4.1841476 2.7100225 #as.factor(System)7 1.1070987 0.9379283 #as.factor(System)8 1.1070987 0.9379283 #as.factor(System)9 2.0472742 1.5417814 #as.factor(System)10 2.5750162 1.8740406 > cbind(pump$S,pump.model.poinorm$fv,pump.model.poigam$fv) # [,1] [,2] [,3] # [1,] 5 5.7095270 5.673258 # [2,] 1 2.4368493 1.850764 # [3,] 5 5.3996111 5.542107 # [4,] 14 14.1594905 14.417893 # [5,] 3 3.2169392 3.382929 # [6,] 19 17.7313717 17.366742 # [7,] 1 0.9098207 1.059285 # [8,] 1 0.9098207 1.059285 # [9,] 4 3.3649256 3.482540 #[10,] 22 21.1616453 21.165196