################################################################################################ ## Code by Vanja Dukic, U of Colorado at Boulder ## Bayesian Statistics and Computing 2019 ## ## Week 10 - Examples, Gibbs sampler ## ## NUCLEAR REACTOR PUMPS EXAMPLE # # Pump reliability data set from Gaver and O'Muircheartaigh (1987) # # The "pump" data set presents the failures of pumps in several # systems of the water reactor neuclear plant Farley 1. # # The "pump" data set contains 4 columns and 10 rows: # # 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. # # Gaver, D P. and O'Muircheartaigh, I. G. 1987. Robust Empirical # Bayes Analyses of Event Rates, Technometrics 29(1),1-15 #################################################### # Full model # Si ~ Poisson( ri ti ) # log(ri) = log(lambdai) + d Gi # lambdai ~ Gamma( a , b ), i = 1,...,10 # # # ti = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5) # Si = c( 5, 1, 5, 14, 3, 19, 1, 1, 4, 22) require("hglm") 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 ############################################################################################### ##################################################### ## Bayesian model # Si ~ Poisson( ri * ti ) # log(ri) = log(lambdai) + d Gi # lambdai ~ Gamma( a , b ) # a ~ Exponential(1.0) # b ~ Gamma(1, 1.0) # d ~ Norm(0,1) # ######################################################3 ### Gibbs sampler ## ## Simpler model, no covariate ## #S ~ Poi(deltais*time) #deltais ~ G(alpha,beta) #beta ~ Gamma(gamma1,gamma2)) alpha=1.8 gamma1 = 0.01 gamma2 = 1 Niter = 10000 n=dim(pump)[1] set.seed(1234567) deltais=matrix(0.1,Niter,n) deltais[1,]=c(rep(0.1,n)) #gamma1s=matrix(0.1,Niter,1) #gamma2s=matrix(0.1,Niter,1) #alphas=matrix(1,Niter,1) betas=matrix(1,Niter,1) #par.chain = cbind(deltais,gamma1s,gamma2s,alphas,betas) par.chain = cbind(deltais,betas) for(i in 2:Niter){ for(j in 1:n) deltais[i,j]=rgamma(1,sh=pump$S[j]+alpha,ra=pump$t[j]+betas[i-1]) #alpha[i]=rgamma(1,sh=gamma+n*alpha,ra=delta+sum(lambda[i,])) betas[i]=rgamma(1,sh=gamma1+n*alpha,ra=gamma2+sum(deltais[i,])) } par.chain = cbind(deltais,betas) par.thin=10 par.thinned = par.chain[seq(from=1, to=Niter,by=par.thin),] acf(par.thinned) plot(as.mcmc.list(mcmc(par.thinned))) summary(as.mcmc.list(mcmc(par.thinned))) plot(as.mcmc(par.chain)) summary.chain = summary(as.mcmc(par.chain)) 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE [1,] 0.06911 0.02563 0.0008104 0.0008104 [2,] 0.15605 0.09409 0.0029753 0.0029753 [3,] 0.10476 0.04024 0.0012726 0.0011500 [4,] 0.12377 0.03119 0.0009864 0.0009864 [5,] 0.62134 0.28180 0.0089114 0.0089114 [6,] 0.60887 0.13411 0.0042409 0.0039974 [7,] 0.83049 0.51678 0.0163420 0.0163420 [8,] 0.84224 0.53600 0.0169498 0.0169498 [9,] 1.29646 0.63051 0.0199385 0.0199385 [10,] 1.84468 0.37727 0.0119303 0.0119303 [11,] 2.45862 0.71081 0.0224777 0.0209236 ######################################################################################3 ## Compare with the HGLM output pump.model.simple.poigam <- hglm(fixed = S ~ 1, random = ~1 | System, offset = log(t), fix.disp = 1, family = poisson(), rand.family = Gamma(link = log), data = pump) print(pump.model.simple.poigam, print.ranef=TRUE) Fixed effects: (Intercept) -0.399823 Random effects: as.factor(System)1 as.factor(System)2 as.factor(System)3 as.factor(System)4 0.08809069 0.14559857 0.13149197 0.17218741 as.factor(System)5 as.factor(System)6 as.factor(System)7 as.factor(System)8 0.87602389 0.90423080 1.22389527 1.22389527 as.factor(System)9 as.factor(System)10 2.27768400 2.95690214 # compare fitted values: > cbind(pump.model.simple.poigam$fv,pump$t*(summary.chain$statistics[1:10,1] )) [,1] [,2] [1,] 5.570483 6.6236025 [2,] 1.534506 2.4194366 [3,] 5.543331 6.5404675 [4,] 14.517873 15.5049636 [5,] 3.077558 3.3086325 [6,] 19.059912 19.2615212 [7,] 0.859933 0.8592895 [8,] 0.859933 0.8749703 [9,] 3.200691 2.7361967 [10,] 20.775778 19.3202894 #rough approximation for st dev from hglm: sqrt((pump.model.simple.poigam$varRanef)) [1] 1.264272 # better intervals, based on poisson pmf: cbind(qpois(0.025,pump.model.simple.poigam$fv), qpois(0.975,pump.model.simple.poigam$fv)) # complete uncertainty estimate cbind(pump.model.simple.poigam$fv, qpois(0.025,pump.model.simple.poigam$fv), qpois(0.975,pump.model.simple.poigam$fv), pump$t*(summary.chain$statistics[1:10,1]), pump$t*(summary.chain$quantiles[1:10,c(1,5)] )) 2.5% 97.5% var1 5.570483 1 11 6.6236025 2.6435443 12.416257 var2 1.534506 0 4 2.4194366 0.4896116 6.055911 var3 5.543331 1 11 6.5404675 2.5656267 12.407775 var4 14.517873 8 22 15.5049636 8.9100527 24.143862 var5 3.077558 0 7 3.3086325 0.9936318 7.082288 var6 19.059912 11 28 19.2615212 11.8381543 28.389860 var7 0.859933 0 3 0.8592895 0.1526564 2.216199 var8 0.859933 0 3 0.8749703 0.1556191 2.294269 var9 3.200691 0 7 2.7361967 0.9273745 5.658913 var10 20.775778 12 30 19.3202894 12.0460900 28.063262 ############################################################################ # With covariate # # # From the HGLM manual: # Poisson model with Gamma distributed random effects # # For dependent count data it is common to model a Poisson distributed response with a gamma # distributed random effect (Lee et al. 2006): # E(yi | β, u) = exp(Xi β + Zi v) # where a 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 ) = λ. # This model is also possible to fit with the hglm package and extends other GLMM functions # (e.g. glmer()) to allow for non-normal distributions for the random effect. 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 pump.model.poigam$fv # [1] 5.673258 1.850764 5.542107 14.417893 3.382929 17.366742 1.059285 # [8] 1.059285 3.482540 21.165196 ##############################################################################################################3333