# Chapter 08: Estimating Bayesian generalized linear models

#### Kjell Nygren

#### 2026-07-26

# General Discussion

As we saw in the previous chapter, using the Bayesian version glmb of the classical function glm requires the specification of an additional prior distribution (a pfamily). Envelope-based sampling for these models builds on (Nygren and Nygren 2006). For Poisson and Binomial families, the most appropriate of the currently implemented priors would be the dNormal pfamily which requires specifying a prior mean and a prior variance-covariance matrix for the regression coefficients. Some models (covered more in a later Vignette) also requires providing either a constant or a prior for the dispersion (variance) parameter associated with the model. Here we focus on the case where only the regression coefficients require a prior and walk through an example where we specify a multivariate prior for a model with two coeffcients related to the Age of menarche (i.e., the age at which girls experience their first period). As part of the setup, we discuss the approach taken and give some general advice in regards to setting of prior distributions.

# The Age of Menarche Data

The data we will be using here is the menarche data covering the Age of Menarche for girls in Warsaw. This data appears in (Venables and Ripley 2002) and in their supporting MASS package. It is provided (and renamed) here so as to allow examples to run during build without loading of the MASS package.

In the data, there are observations related to 25 cohorts of girls at ages ranging from an average age of 9.21 to an average age of 17.58. The average age for each cohort is captured in a column named Age. For each of the cohorts, there is also a column (Total) giving the total number of girls in the cohort and a column Menarche giving the number of girls out of the cohort who had experienced Menarche at that age.

The first 5 observations are printed below. It is important to note that in Bayesian analyses (at least those that seek to incorporate prior/external information), the priors should ideally be based on information outside of the data being studied, so we will hold off on a more complete view of the data until we have mapped out our prior. The analysis we will be exploring here is a model of the relationship between age and the percentage of girls who have had their first period using various formulations of binomial regression.

    ## Load menarche data
    data(menarche,package="MASS")
    head(menarche, 5)
    #>     Age Total Menarche
    #> 1  9.21   376        0
    #> 2 10.21   200        0
    #> 3 10.58    93        0
    #> 4 10.83   120        2
    #> 5 11.08    90        2

# Variable Transformation and Preliminary Setup

Because of the need for a prior distribution, it is important that coefficients in Bayesian models can be interpreted so that analysts can set reasonable priors. This is true regardless of whether analysts are looking to give a model a weak prior (that will limit any move of the estimates away from those driven purely by the data for that specific model) or if the analyst is looking to give the model a strong prior that incorporates meaningful external information and/or prior knowledge. If the analysts can’t interpret the coefficients, then is is difficult to know whether a prior distribution is weak or strong and whether it is likely to move estimates meaningfully away from classical estimates.

A useful approach to making coefficients more interpretable involves centering variables closer to the center (or at least within the support of) the data. To facilitate the interpretation of the coefficients in our model, we will therefore pick a reference age ref_age1 within the support of the data (we pick 13 but you can feel free to pick a different reference age or to use the mean across the observations) and define a new variable Age2=Age-ref_age1. This will allow our prior for the intercept to represent a point estimate (on the chosen link function scale) for the percentage of girls who have had their first period by age 13. It also happens to be an age close to the average age of menarche in other sources (Hamilton-Fairley 2004, 29) which reported the UK average age of menarche in the UK as 12.9.

In addition to the intercept, our model will contain a slope coefficient for the modified age variable that will capture the change in the percentage of girls who have had their first period (on the relevant link scale). To be able to formulate a prior for this, we will set up a second reference age ref_age2 (our choice will be 15) at which we also will formulate a point estimate for the percentage who have had their first period. A point estimate for the slope on the link scale will then be derived as the difference between the two estimates (on the link scale) divided by the Age gap between the two reference ages. Below is code doing the variable transformation and a preliminary setup to prepare for the prior specification.

    ## Number of variables in model
    Age=menarche$Age
    nvars=2
    set.seed(333)

    ## Reference Ages for setting of priors and Age_Difference
    ref_age1=13  # user can modify this
    ref_age2=15  ## user can modify this

    ## Define variables used later in analysis
    Age2=Age-ref_age1
    Age_Diff=ref_age2-ref_age1

# Setting prior point estimates, associated credible intervals, and correlation

We are now ready to specify our two point estimates with associated prior credible intervals and an assumed correlation (on the link scale). For the first reference age (13), we will provide a point estimate of 0.5 (50%) for the percentage of girls who have had their first period by age 13. To account for the possibility that this percentage was a bit lower in Warsaw at the time data was collected, we will assume that the lower end of a 90% credible interval is at 0.3 (30%).

For the second reference age (15), we have not accessed as much meaningful data from other sources to inform the prior estimate, but we have a general sense that “most” girls likely have had their period by that age. Based on this, we select 0.9 (90%) as a prior point estimate and select 0.7 (70%) as the lower end of our associated credible interval.

To complete the information needed to map this information to both prior means and a fully populated Variance-Covariance matrix, we account for the likely correlation between the percentage at the two ages. If one of our two point estimates turn out to be either too high (or low) for the Warsaw data, it is likely that the other estimate also could be too high (or low) for the data as well. To that end, we assume a prior 0.4 correlation between the two estimates on the relevant link scale.

Here is code setting these assumptions up.

    ## Point estimates at reference ages
    m1=0.5
    m2=0.9

    ## Lower bound of prior credible intervals for point estimates
    m1_lower=0.3
    m2_lower=0.7

    ## Assumed correlation between the two (on link scale)
    m_corr=0.4

# Mapping the prior to the link scale

To be able to use the prior information in our model, the information has to be mapped to the link scale and the point estimates (with associated credible intervals) have to be mapped to a point estimate and uncertainty for the slope. In order to do this, we first get the required link function and set up the matrices into which we will store the prior information (mu1 and V1).

    ## Set up link function and initialize prior mean and Variance-Covariance matrices
    bi_logit <- binomial(link="logit")
    mu1<-matrix(0,nrow=nvars,ncol=1)
    rownames(mu1)=c("Intercept","Age2")
    colnames(mu1)=c("Prior Mean")
    V1<-1*diag(nvars)
    rownames(V1)=c("Intercept","Age2")
    colnames(V1)=c("Intercept","Age2")

Setting the prior mean

This piece of code now maps the prior assumptions to the logit scale and sets the prior mean vector. The mean for the intercept is just the logit transformation of the prior point estimate while the mean for the slope uses the difference between the logit transformed point estimates divided by the Age_Difference.

    ## Prior mean for intercept is set to point estimate
    ## at reference age1 (on logit scale)
    mu1[1,1]=bi_logit$linkfun(m1)

    ## Prior mean for slope is set to difference in point estimates
    ## on logit scale divided by Age_Diff

    mu1[2,1]=(bi_logit$linkfun(m2) -bi_logit$linkfun(m1))/Age_Diff
    print(mu1)
    #>           Prior Mean
    #> Intercept   0.000000
    #> Age2        1.098612

Setting the prior Variance-Covariance Matrix

The piece of code setting the prior variance-covariance matrix for this prior specification is a bit more complex. It proceeds as follows:

- Calculate the prior standard deviation for each of the point estimates by leveraging the lower bounds
- Calculate the implied standard deviation for the slope using the formula for the variance of differences
- Calculate the implied covariance between the intercept and the slope
- Populate the Prior Variance-Covariance matrix using the derived information

    ## Implied standard deviations for point estimates on logit scale

    sd_m1= (bi_logit$linkfun(m1) -bi_logit$linkfun(m1_lower))/1.96
    sd_m2= (bi_logit$linkfun(m2) -bi_logit$linkfun(m2_lower))/1.96

    ## Also compute implied estimate for upper bound of confidence intervals

    m1_upper=bi_logit$linkinv(bi_logit$linkfun(m1)+sd_m1*1.96)
    m2_upper=bi_logit$linkinv(bi_logit$linkfun(m2)+sd_m2*1.96)
    print("m1_upper is:")
    #> [1] "m1_upper is:"
    m1_upper
    #> [1] 0.7
    print("m2_upper is:")
    #> [1] "m2_upper is:"
    m2_upper
    #> [1] 0.972

    ## Implied Standard deviation for slope (using variance formula for difference between two variables)
    a=(1/Age_Diff)
    sd_slope=sqrt((a*sd_m1)^2+(a*sd_m2)^2-2*a*a*(sd_m1*sd_m2*m_corr))

    #Cov(m1,slope)=cov(m1, a*(m2-m1)) =a*E[(m1-E[m1])((m2-m1)-E[m2-m1])]
    #   =a*E[(m1-E[m1])(m2-E[m2])]- a* E[(m1-E[m1])(m1-E[m1])]
    ##   =a*Cov[m1,m2] - a*Var[m1]
    ##  =a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1
    cov_V1=a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1

    # Set covariance matrix
    V1[1,1]=sd_m1^2
    V1[2,2]=sd_slope^2
    V1[1,2]=cov_V1
    V1[2,1]=V1[1,2]
    print("V1 is:")
    #> [1] "V1 is:"
    print(V1)
    #>             Intercept        Age2
    #> Intercept  0.18687882 -0.03389182
    #> Age2      -0.03389182  0.10576217

# Running the Model

We are now ready to run a Logit model using the above specification. As a best practice, we recommend creating a data frame holding all the variables used in the model formula and to use it as the “data” for the model. We will show how this data can be used an an input to the prediction function in the next Chapter.

The below code performs the step creating the Model data frame and then submits the model.

    Menarche_Model_Data=data.frame(Age=menarche$Age,Total=as.integer(menarche$Total),
    Menarche=as.integer(menarche$Menarche),Age2)

    glmb.out1<-glmb(n=1000,cbind(Menarche, as.integer(Total-Menarche)) ~Age2,family=binomial(logit),
    pfamily=dNormal(mu=mu1,Sigma=V1),data=Menarche_Model_Data)

Reviewing the printed output

To get a first read of the results, we now print the results of the model and re-print the prior mean specification so that we can do an initial comparison of the prior to the estimated posterior means.

    # Print model output
    print(glmb.out1)
    #>
    #> Call:  glmb(formula = cbind(Menarche, as.integer(Total - Menarche)) ~
    #>     Age2, family = binomial(logit), pfamily = dNormal(mu = mu1,
    #>     Sigma = V1), n = 1000, data = Menarche_Model_Data)
    #>
    #> Posterior Mean Coefficients:
    #> (Intercept)         Age2
    #>    -0.01547      1.61809
    #>
    #> Effective Number of Parameters: 2.037394
    #> Expected Residual Deviance: 28.80323
    #> DIC: 114.8924

    # Print prior mean as comparison
    print(t(mu1))
    #>            Intercept     Age2
    #> Prior Mean         0 1.098612

In reviewing the coefficients, we see that the intercepts look quite close, so our prior estimate at age 13 does not seem to have been a bad estimate for what the Warsaw data suggests. Our prior estimate for the slope parameter seems to be a bit lower than the posterior mean estimate but without information on the standard error or credible intervals connected to the posterior mean, we can’t tell how off it might have been. We ignore the other printed output for now (it will be discussed in the next chapter) and instead turn to a more extensive model summary.

Summarizing the model

To get a more detailed view of the model output, we simply use the summary function and its print method much like you would to generate a summary for output from the glm function. This results in a printed view with several sections of outputs

- Expected Deviance Residuals (Discussed more in the next Chapter)
- Prior and Maximum Likelihood Estimates with Standard Deviations
- Bayesian Estimates Based on iid draws from the posterior (in this case 1000 draws)
- Distribution Percentiles (Credible Interval type information)
- DIC related information (Discussed more in the next chapter)
- Information related to the number of candidates per iid draw (discussed more in one of our Appendix chapters)

We discuss several of these after the print view of the output from the call to the summary function.

    summary(glmb.out1)
    #> Call
    #> glmb(formula = cbind(Menarche, as.integer(Total - Menarche)) ~ Age2, family = binomial(logit), pfamily = dNormal(mu = mu1, Sigma = V1), n = 1000, data = Menarche_Model_Data)
    #>
    #> Expected Residuals:
    #>        Min         1Q     Median         3Q        Max
    #> -2.0786671 -0.9942019 -0.4341926  0.7324099  1.3458398
    #>
    #> Prior and Maximum Likelihood Estimates with Standard Deviations
    #>
    #>             Null Mode Prior Mean Prior.sd Max Like. Like.sd
    #> (Intercept)   0.36015    0.00000  0.43229  -0.01081   0.063
    #> Age2          0.00000    1.09861  0.32521   1.63197   0.059
    #>
    #> Bayesian Estimates Based on 1000 iid draws
    #>
    #>             Post.Mode Post.Mean   Post.Sd  MC Error Pr(Null_tail) SE(tail)
    #> (Intercept) -0.013688 -0.015468  0.064810  0.002049      0.000000    0.016
    #> Age2         1.614377  1.618086  0.057518  0.001819      0.000000    0.000
    #>             Pr(Prior_tail)
    #> (Intercept)           0.42
    #> Age2                <2e-16 ***
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #>   Directional Tail Summaries:
    #>
    #>                Metric vs Null vs Prior
    #>  Mahalanobis Distance 28.7720   9.0922
    #>      Tail Probability  0.0000   0.0000
    #>   [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) -0.15923 -0.13802 -0.12209 -0.01379  0.08861  0.10770 0.134
    #> Age2         1.48741  1.51361  1.52336  1.61863  1.71254  1.73036 1.759
    #>
    #> Effective Number of Parameters: 2.037394
    #> Expected Residual Deviance: 28.80323
    #> DIC: 114.8924
    #>
    #> Expected Mean dispersion: 1
    #> Sq.root of Expected Mean dispersion: 1
    #>
    #> Mean Likelihood Subgradient Candidates Per iid sample: 1.272

Interpreting the “Prior and Maximum Likelihood Estimates with Standard Deviations” table

This table gives a side-by-side comparison of the prior and the data driven maximum likelihood estimates (that should match the ouput from the glm function for the same model). In addition to revealing any differences in the estimates themselves, the relative magnitude of the standard deviations give insight into how “strong” the prior specification is relative to the information contained in the data. The unstandardized weights on the prior and data respectively are in most cases likely to be close to (1/sd_prior)^2 and (1/sd_data)^2 respectively. The fact that the prior standard deviations are much larger in this case tells us that the posterior mean estimates are likely to be much closer to the maximum likelihood estimates than to the prior estimates.

Interpreting the “Bayesian Estimates Based on 1000 iid draws table

This table gives more information related to the estimates seen in the print output. We now have a view of both the posterior mode and posterior mean. If these were substantially different, it would indicate that the posterior density had substantial asymmetry. The posterior standard deviation will generally be smaller than both the prior and the likelihood standard deviations. In this case, the posterior standard deviation for the intercept appears to be nearly the same as the likelihood standard deviation while the posterior standard deviation for Age2 is substantially smaller (suggesting some weight on the prior).

The tail probability gives the probability (under the posterior density) of the true coefficient being at least as far away from the center as the prior point estimate. In this case, there is ~44% probability of the intercept being greater than 0 but only a ~0.01% chance that the Age2 coefficient is smaller than the prior point estimate of 1.098. The latter is flagged as highly significant suggesting that the true slope coefficient is substantially larger than the prior point estimate.

The MC error column gives a sense of any error in the estimates driven by the fact that we are using monte carlo simulation to generate our estimates. In this case, we used 1000 draws which appears to have been enough to make this error small. The user can experiment with this to see how the errors increase if we base the estimates on a smaller number of draws. In general, it is not recommended to interpret model output based on small number of draws. On the other hand, the fact that estimates are iid draws as opposed to Markov Chain Monte Carlo outputs eliminates any need for convergence diagnostics and burn-in phases.

Interpreting the “Distribution Percentile” table

This table providence credible interval information at several significance levels for the model. It gives the user insight into the range of most of the support for the posterior density and how it compares to both the prior and maximum likelihood point estimates. In this case, we can see that the prior point estimate for the intercept is close to the center of the density while the prior point estimate is far away from it.

Brief discussion of Findings from model

As we noted above, the posterior estimates for the intercept appear broadly consistent with the prior specification, so in a sense our prior beliefs are “validated” for the 13 year old girls. As the slope is substantially larger than in our prior specification, however, it seems like we underestimated the % of girls at age 15 who had already had their first period. The percentage is apparently substantially higher than the 90% we specificed using our prior and if we are Bayesians we should hence update our beliefs substantially based on the data. To understand exactly how much requires us to look at the model predictions(in particular predictions on the inverse link scale). We will do so in a later chapter.

# Appendix A. Bayes Rules! companion — logistic priors on `weather_perth`

Book: (Johnson et al. 2022), Chapter 13 — Logistic regression (model 13.5)
Data: `bayesrules::weather_perth`
Informative book priors (used here): \(\beta_0 \sim N(-1.4, 0.7^2)\), \(\beta_1 \sim N(0.07, 0.035^2)\) on the logit scale.
glmbayes: `glmb()` with the same prior means and SDs in `dNormal(mu, Sigma)`.

> Requires suggested package `bayesrules`. Book comparison uses 80% posterior intervals from Ch. 13 (`posterior_interval(rain_model_1)`); midpoints are shown as a rough book “mean.”

    library(bayesrules)
    weather <- bayesrules::weather_perth
    weather$raintomorrow <- as.integer(weather$raintomorrow == "Yes")

    mu_w <- matrix(c(-1.4, 0.07), nrow = 1)
    colnames(mu_w) <- c("(Intercept)", "humidity9am")
    Sigma_w <- diag(c(0.7^2, 0.035^2))
    dimnames(Sigma_w) <- list(colnames(mu_w), colnames(mu_w))

    ## Bayes Rules! Ch. 13: 80% posterior intervals on log-odds scale (rain_model_1)
    book_br08 <- data.frame(
      parameter = c("(Intercept)", "humidity9am"),
      book_lo   = c(-5.08785, 0.04147),
      book_hi   = c(-4.13450, 0.05487),
      book_mid  = c(-4.611175, 0.04817),
      check.names = FALSE
    )

    set.seed(2026)
    glmb_w <- glmb(
      raintomorrow ~ humidity9am,
      family  = binomial(),
      pfamily = dNormal(mu = mu_w, Sigma = Sigma_w),
      data    = weather,
      n       = 2000
    )
    #> [glmb_Standardize_Model][NOTE] Posterior Hessian is moderately ill-conditioned.
    #>   kappa(H) = 67316.9
    print(glmb_w)
    #>
    #> Call:  glmb(formula = raintomorrow ~ humidity9am, family = binomial(),
    #>     pfamily = dNormal(mu = mu_w, Sigma = Sigma_w), n = 2000,
    #>     data = weather)
    #>
    #> Posterior Mean Coefficients:
    #> (Intercept)  humidity9am
    #>    -3.94489      0.03889
    #>
    #> Effective Number of Parameters: 1.821911
    #> Expected Residual Deviance: 874.6306
    #> DIC: 876.4525

    br08_compare <- data.frame(
      parameter = book_br08$parameter,
      `Book 80% lo` = book_br08$book_lo,
      `Book 80% hi` = book_br08$book_hi,
      `Book midpoint` = book_br08$book_mid,
      `glmb Post.Mean` = as.numeric(glmb_w$coef.means[book_br08$parameter]),
      `glmb Post.Sd`   = sapply(book_br08$parameter, function(p)
        sd(glmb_w$coefficients[, p, drop = TRUE])),
      check.names = FALSE
    )
    knitr::kable(br08_compare, digits = 4,
      caption = "Bayes Rules! Ch. 13 (80% intervals) vs. glmb() with informative priors")

Bayes Rules! Ch. 13 (80% intervals) vs. glmb() with informative priors
| parameter | Book 80% lo | Book 80% hi | Book midpoint | glmb Post.Mean | glmb Post.Sd
--- | --- | --- | --- | --- | --- | ---
(Intercept) | (Intercept) | -5.0879 | -4.1345 | -4.6112 | -3.9449 | 0.3159
humidity9am | humidity9am | 0.0415 | 0.0549 | 0.0482 | 0.0389 | 0.0046

The prior intercept \(-1.4\) on the logit scale corresponds to \(\mathrm{logit}^{-1}(-1.4) \approx 0.20\), matching the book’s “about 20% chance of rain on a typical day.” `glmb Post.Mean` for `humidity9am` should fall between the book’s 80% bounds (0.0415, 0.0549) and below the prior mean 0.07.

# Concluding Discussion

In this note, we walked through an example illustrating a process for setting of the prior in the case of the Age of Menarche model based on the data in (Venables and Ripley 2002). There are several approaches taken here the we believe analysts can find useful in setting reasonable priors. These include:

- Centering variables so that the intercept can be given a meaningful prior and so that meaningful priors for the other variables can be formulated based on sensitivity in the support of data
- Incorporating external information if available (like our information regarding the average age of Menarche in the UK)
- Making reasonable assumptions based on other knowledge/insight (e.g., our prior point estimate at age 15)

In the next section, we will take a take a look at how the analysts can check whether the data seems consistent with the specified prior.

# References

Hamilton-Fairley, Diana. 2004. Lecture Notes: Obstetrics and Gynaecology. Blackwell Publishing.
Johnson, Alicia A., Miles Q. Ott, and Mine Dogucu. 2022. Bayes Rules! An Introduction to Applied Bayesian Modeling. CRC Press. https://www.bayesrulesbook.com.
Nygren, K. N., and L. M. Nygren. 2006. “Likelihood Subgradient Densities.” Journal of the American Statistical Association 101 (475): 1144–56. https://doi.org/10.1198/016214506000000357.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. Springer.