Chapter A03: Methods available in glmbayes

Kjell Nygren

2026-04-30

The glmb function and related method functions that handle the output are designed to be Bayesian versions of the glm function and many of its method functions. This vignette shows how the basic setup/calling of the functions compare and then walks through how the method functions for glmb can be called to generate similar outputs to those from the glm functions. GLM background is standard (McCullagh and Nelder 1989; Venables and Ripley 2002); the glmb interface parallels glm as described in (Hastie and Pregibon 1992).

Dobson Randomized Control Data

To understand how the outputs of the glmb function mirrors those for the glm function, it is useful to take a look at the first portion of the example that is provided with the glm function. The data are randomized controlled trial counts from (Dobson 1990). Here is a view of the data:

## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
#>   treatment outcome counts
#> 1         1       1     18
#> 2         1       2     17
#> 3         1       3     15
#> 4         2       1     20
#> 5         2       2     10
#> 6         2       3     20
#> 7         3       1     25
#> 8         3       2     13
#> 9         3       3     12

Calling the two functions

The example code for the glm function specifies a Poisson regression model for this data. The below code chunks shows how the glmb function can be used in much the same fashion but with some extra requirements. In particular, we need to provide a prior distribution (in this case a Multivariate Normal prior) and related constants (in this case, a prior Mean and and a prior Variance-Covariance matrix) for the coefficients of interest. Optionally, we can also tell the function how many draws to make from the posterior distribution (the default used here generates 1000 iid draws).

Here is the call to the classical glm function:

## Call to glm
glm.D93 <- glm(counts ~ outcome + treatment, 
              family = poisson())

To use the glmb function, we first use a call to a function Prior_Setup to initialize the prior and to get the variable names needed. We use the default prior from the Prior_Setup function (discussed further in a separate vignette).

## Using glmb
## Step 1: Set up Default Prior 
ps=Prior_Setup(counts ~ outcome + treatment,family=poisson())
#> Using default pwt = 0.01 (low-d default).
mu=ps$mu
V=ps$Sigma

We now call the glmb function and include the default prior distribution in addition to the two required arguments for the glm function.

# Step 3: Call the glmb function
glmb.D93<-glmb(counts ~ outcome + treatment, family=poisson(), pfamily=dNormal(mu=mu,Sigma=V))

In the above, it is worth noting the additional step that we went through when specifying the prior. We used a call to the function Prior_Setup to get the correct dimensions for the mean and variance-covariance matrices and to initialize the constants. The Prior_Setup function also provided information on the Variable names in the design matrix (which also corresponds to the names of the coefficients eventually estimated) so we can make informed changes to the prior if so desired.

The next couple of vignette’s will discuss the prior specification in more details. Here we instead focus on reviewing the output from the function and how it can be used.

Printing the output

Taking a look at the basic printed output, we can see that the two closely mirror each other with the glmb posterior means replacing the glm maximum likelihood estimates.

## Printed view of the output from the glm function 
print(glm.D93)
#> 
#> Call:  glm(formula = counts ~ outcome + treatment, family = poisson())
#> 
#> Coefficients:
#> (Intercept)     outcome2     outcome3   treatment2   treatment3  
#>   3.045e+00   -4.543e-01   -2.930e-01    1.218e-15    8.438e-16  
#> 
#> Degrees of Freedom: 8 Total (i.e. Null);  4 Residual
#> Null Deviance:       10.58 
#> Residual Deviance: 5.129     AIC: 56.76
## Printed view of the output from the glmb function 
print(glmb.D93)
#> 
#> Call:  glmb(formula = counts ~ outcome + treatment, family = poisson(), 
#>     pfamily = dNormal(mu = mu, Sigma = V))
#> 
#> Posterior Mean Coefficients:
#> (Intercept)     outcome2     outcome3   treatment2   treatment3  
#>    3.045375    -0.469607    -0.303245    -0.019815    -0.007108  
#> 
#> Effective Number of Parameters: 4.918062 
#> Expected Residual Deviance: 10.09891 
#> DIC: 56.64915

In addition to the posterior means, the glmb printed output also returned three pieces of information that are similar to (but not quite the same) as the classical output. The “Effective Number of Parameters” should in general be close to (but not exactly the same) as the number of parameters estimated while the Expected Residual Deviance in general will be higher than the corresponding maximum likelihood estimate for the Residual Deviance (since the latter is designed to minimize it). Finally, the DIC is essentially a Bayesian version of the AIC. These measures will be discussed in greater detail in one of our later vignettes.

Methods Available

In addition to the basic print function output, the glmb function returns an object with an assigned class “glmb” for which a number of generic functions (or methods) are available. The class “glmb” inherits from “glm” and “lm” and as such many functions for those classes work directly for “glmb”. For some of the instances where the inherited methods fail and/or could produce incorrect results, we have implemented methods specifically for the glmb class. The methods for classes lm, glm, and glmb are listed below. We will use many of these later in this vignette.

## Methods for class "lm"
methods(class="lm")
#>  [1] add1           addterm        alias          anova          boxcox        
#>  [6] case.names     coerce         confint        cooks.distance deviance      
#> [11] dfbeta         dfbetas        dffits         drop1          dropterm      
#> [16] dummy.coef     effects        extractAIC     family         formula       
#> [21] hatvalues      influence      initialize     kappa          labels        
#> [26] logLik         logtrans       model.frame    model.matrix   nobs          
#> [31] plot           predict        print          proj           qr            
#> [36] residuals      rstandard      rstudent       show           simulate      
#> [41] slotsFromS3    summary        variable.names vcov          
#> see '?methods' for accessing help and source code
## Methods for class "glm"
methods(class="glm")
#>  [1] add1           addterm        anova          coerce         confint       
#>  [6] cooks.distance deviance       dfbetas        dffits         drop1         
#> [11] dropterm       effects        extractAIC     family         formula       
#> [16] gamma.shape    influence      initialize     logLik         model.frame   
#> [21] nobs           predict        print          profile        residuals     
#> [26] rstandard      rstudent       show           sigma          slotsFromS3   
#> [31] summary        vcov           weights       
#> see '?methods' for accessing help and source code
## Methods for class "glmb"
methods(class="glmb")
#>  [1] anova          case.names     confint        cooks.distance dfbetas       
#>  [6] dummy.coef     extractAIC     influence      logLik         plot          
#> [11] predict        print          residuals      rstandard      rstudent      
#> [16] simulate       summary        variable.names vcov          
#> see '?methods' for accessing help and source code

The summary functions

Let’s take a closer look at the outputs of the summary functions.

glm summary output

In turn, we see the Call, a list of Deviance residuals, and the estimated coefficients. The coefficients are then followed by some additional model related information.

## summary output for the "glm" class
summary(glm.D93)
#> 
#> Call:
#> glm(formula = counts ~ outcome + treatment, family = poisson())
#> 
#> Coefficients:
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  3.045e+00  1.709e-01  17.815   <2e-16 ***
#> outcome2    -4.543e-01  2.022e-01  -2.247   0.0246 *  
#> outcome3    -2.930e-01  1.927e-01  -1.520   0.1285    
#> treatment2   1.217e-15  2.000e-01   0.000   1.0000    
#> treatment3   8.438e-16  2.000e-01   0.000   1.0000    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 10.5814  on 8  degrees of freedom
#> Residual deviance:  5.1291  on 4  degrees of freedom
#> AIC: 56.761
#> 
#> Number of Fisher Scoring iterations: 4

glmb summary output

The summary for the glmb function follow a similar structure but adds a table containing information related to the prior and the maximum likelihood above the table with the means for the estimated Bayesian coefficients. The output below the main table with coefficients is also modified to contain similar (but slightly different) pieces of information (the details of which are discussed elsewhere).

## summary output for the "glm" class
summary(glmb.D93)
#> Call
#> glmb(formula = counts ~ outcome + treatment, family = poisson(), 
#>     pfamily = dNormal(mu = mu, Sigma = V))
#> 
#> Expected Residuals:
#>        Min         1Q     Median         3Q        Max 
#> -0.9276586 -0.6953010 -0.1542997  0.8542582  1.1469176 
#> 
#> Prior and Maximum Likelihood Estimates with Standard Deviations
#> 
#>              Null Mode Prior Mean   Prior.sd  Max Like. Like.sd
#> (Intercept)  2.813e+00  2.813e+00  1.700e+00  3.045e+00   0.171
#> outcome2     0.000e+00  0.000e+00  2.012e+00 -4.543e-01   0.202
#> outcome3     0.000e+00  0.000e+00  1.918e+00 -2.930e-01   0.193
#> treatment2   0.000e+00  0.000e+00  1.990e+00  1.217e-15   0.200
#> treatment3   0.000e+00  0.000e+00  1.990e+00  8.438e-16   0.200
#> 
#> Bayesian Estimates Based on 1000 iid draws
#> 
#>              Post.Mode  Post.Mean    Post.Sd   MC Error Pr(Null_tail) SE(tail)
#> (Intercept)  3.042e+00  3.045e+00  1.629e-01  5.151e-03     8.300e-02    0.009
#> outcome2    -4.497e-01 -4.696e-01  2.000e-01  6.325e-03     9.000e-03    0.003
#> outcome3    -2.901e-01 -3.032e-01  1.857e-01  5.871e-03     4.900e-02    0.007
#> treatment2  -1.395e-05 -1.982e-02  2.048e-01  6.476e-03     4.390e-01    0.016
#> treatment3   1.209e-06 -7.108e-03  2.013e-01  6.365e-03     4.990e-01    0.016
#>             Pr(Prior_tail)   
#> (Intercept)          0.083 . 
#> outcome2             0.009 **
#> outcome3             0.049 * 
#> treatment2           0.439   
#> treatment3           0.499   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>   Directional Tail Summaries:
#> 
#>                Metric vs Null vs Prior
#>  Mahalanobis Distance  2.4364   2.4364
#>      Tail Probability  0.0100   0.0100
#>   [Tail probabilities are P(delta^T * Z <= 0) in whitened space]
#> 
#> 
#> Distribution Percentiles
#> 
#>                  1.0%      2.5%      5.0%    Median     95.0%     97.5%  99.0%
#> (Intercept)  2.662145  2.721820  2.775757  3.045359  3.308822  3.338802  3.399
#> outcome2    -0.910834 -0.852322 -0.793186 -0.474863 -0.146508 -0.074720 -0.027
#> outcome3    -0.761660 -0.666781 -0.602717 -0.295056 -0.006935  0.057228  0.139
#> treatment2  -0.488061 -0.415362 -0.343864 -0.029525  0.327688  0.390714  0.457
#> treatment3  -0.484902 -0.414396 -0.352729  0.001169  0.324980  0.376591  0.440
#> 
#> Effective Number of Parameters: 4.918062 
#> Expected Residual Deviance: 10.09891 
#> DIC: 56.64915 
#> 
#> Expected Mean dispersion: 1 
#> Sq.root of Expected Mean dispersion: 1 
#> 
#> Mean Likelihood Subgradient Candidates Per iid sample: 1.827

Model Fit, Predictions, Deviance Residuals, Covariance Matrices, and Confidence/Credible Intervals

Let’s next take a look at the outputs from these additional methods to see how they compare. Note that the Bayesian version of these contain random draws tied to the underlying distributions so the column means are mostly used in these comparisons.

Fitted values

## fitted outputs for the glm function
fitted(glm.D93)
#>        1        2        3        4        5        6        7        8 
#> 21.00000 13.33333 15.66667 21.00000 13.33333 15.66667 21.00000 13.33333 
#>        9 
#> 15.66667
## mean of fitted outputs for the glm function
## works without a "glmb" class specific generic function
colMeans(fitted(glmb.D93))
#>        1        2        3        4        5        6        7        8 
#> 21.29688 13.38418 15.77647 20.90986 13.14965 15.48002 21.17498 13.31044 
#>        9 
#> 15.68789

Predictions(linear predictors)

## predictions for the glm function
predict(glm.D93)
#>        1        2        3        4        5        6        7        8 
#> 3.044522 2.590267 2.751535 3.044522 2.590267 2.751535 3.044522 2.590267 
#>        9 
#> 2.751535
## predictions for the glmb function
colMeans(glmb.D93$linear.predictors) # no current predict function
#>        1        2        3        4        5        6        7        8 
#> 3.045375 2.575768 2.742129 3.025559 2.555953 2.722314 3.038267 2.568660 
#>        9 
#> 2.735021
colMeans(predict(glmb.D93)) 
#>        1        2        3        4        5        6        7        8 
#> 3.045375 2.575768 2.742129 3.025559 2.555953 2.722314 3.038267 2.568660 
#>        9 
#> 2.735021

Residuals

## residuals for the glm function
residuals(glm.D93)
#>           1           2           3           4           5           6 
#> -0.67124923  0.96272360 -0.16964662 -0.21998507 -0.95552353  1.04938637 
#>           7           8           9 
#>  0.84715368 -0.09167147 -0.96656372
## residuals for the glmb function
colMeans(residuals(glmb.D93))
#>           1           2           3           4           5           6 
#> -0.69530096  0.99567823 -0.15429968 -0.15633497 -0.86121772  1.14691762 
#>           7           8           9 
#>  0.85425816 -0.03755228 -0.92765863

vcov

## vcov for the glm function
vcov(glm.D93)
#>             (Intercept)      outcome2      outcome3    treatment2    treatment3
#> (Intercept)  0.02920635 -1.587302e-02 -1.587302e-02 -2.000000e-02 -2.000000e-02
#> outcome2    -0.01587302  4.087301e-02  1.587302e-02 -7.889764e-18 -7.435533e-18
#> outcome3    -0.01587302  1.587302e-02  3.714961e-02 -5.991636e-18 -6.847584e-18
#> treatment2  -0.02000000 -7.889764e-18 -5.991636e-18  4.000000e-02  2.000000e-02
#> treatment3  -0.02000000 -7.435533e-18 -6.847584e-18  2.000000e-02  4.000000e-02
## vcov for the glmb function
vcov(glmb.D93)
#>             (Intercept)      outcome2      outcome3    treatment2    treatment3
#> (Intercept)  0.02653307 -0.0150206964 -0.0141197005 -0.0194907847 -0.0186977072
#> outcome2    -0.01502070  0.0400022842  0.0143093362  0.0009934275  0.0003620475
#> outcome3    -0.01411970  0.0143093362  0.0344740844 -0.0004603668  0.0001177139
#> treatment2  -0.01949078  0.0009934275 -0.0004603668  0.0419343750  0.0200409058
#> treatment3  -0.01869771  0.0003620475  0.0001177139  0.0200409058  0.0405073918

confint

Confidence intervals

## confint for the glm function
confint(glm.D93)
#> Waiting for profiling to be done...
#>                  2.5 %      97.5 %
#> (Intercept)  2.6958215  3.36655581
#> outcome2    -0.8577018 -0.06255840
#> outcome3    -0.6753696  0.08244089
#> treatment2  -0.3932548  0.39325483
#> treatment3  -0.3932548  0.39325483
## confint for the glmb function
confint(glmb.D93)
#>                   2.5%       97.5%
#> (Intercept)  2.7218196  3.33880223
#> outcome2    -0.8523221 -0.07471989
#> outcome3    -0.6667814  0.05722808
#> treatment2  -0.4153618  0.39071396
#> treatment3  -0.4143964  0.37659097

AIC/DIC, Deviance, and the Log-Likelihood

These model statistics are useful when comparing different model specifications. The Bayesian versions of these will be discussed in greater detail in a separate Vignette.

AIC/DIC

## AIC for the glm function (equivalent degrees of freedom and the AIC)
extractAIC(glm.D93)
#> [1]  5.00000 56.76132
## DIC for the glmb function (estimated effective number of parameters and the DIC)
extractAIC(glmb.D93)
#>        pD       DIC 
#>  4.918062 56.649153

Deviance

## Deviance for the glm function
deviance(glm.D93)
#> [1] 5.129141
## Deviance for the glmb function
## works without a "glmb" class specific generic function
mean(deviance(glmb.D93))
#> [1] 10.09891

Log-Likelihoods

## Deviance for the glm function
logLik(glm.D93)
#> 'log Lik.' -23.38066 (df=5)
## Deviance for the glmb function
mean(logLik(glmb.D93))
#> [1] -25.86555

Model Frame, formula, family, nobs, and show

These are mostly useful for understanding various aspects of the model.

Model Frames

## Model Frame for the glm function
model.frame(glm.D93)
#>   counts outcome treatment
#> 1     18       1         1
#> 2     17       2         1
#> 3     15       3         1
#> 4     20       1         2
#> 5     10       2         2
#> 6     20       3         2
#> 7     25       1         3
#> 8     13       2         3
#> 9     12       3         3
## Model Frame for the glmb function
model.frame(glmb.D93$glm)
#>   counts outcome treatment
#> 1     18       1         1
#> 2     17       2         1
#> 3     15       3         1
#> 4     20       1         2
#> 5     10       2         2
#> 6     20       3         2
#> 7     25       1         3
#> 8     13       2         3
#> 9     12       3         3

formula

## formula for the glm function
formula(glm.D93)
#> counts ~ outcome + treatment
## formula for the glmb function
formula(glmb.D93)
#> counts ~ outcome + treatment

family

## family for the glm function
family(glm.D93)
#> 
#> Family: poisson 
#> Link function: log
## family for the glmb function
family(glmb.D93$glm)
#> 
#> Family: poisson 
#> Link function: log

nobs

## nobs for the glm function
nobs(glm.D93)
#> [1] 9
## nobs for the glmb function
nobs(glmb.D93)
#> [1] 9

show

The show method returns the same output as the print function

References

Dobson, A. J. 1990. An Introduction to Generalized Linear Models. Chapman; Hall.
Hastie, T. J., and D. Pregibon. 1992. “Generalized Linear Models.” Chap. 6 in Statistical Models in S, edited by J. M. Chambers and T. J. Hastie. Wadsworth & Brooks/Cole.
McCullagh, P., and J. A. Nelder. 1989. Generalized Linear Models. Chapman; Hall.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. Springer.