1 Introduction

This analysis uses the penguins_bayes data to study the relationship between penguin flipper length in millimeters (flipper_length_mm) and body mass in grams (body_mass_g).

Let:

\[ Y_i=\text{flipper length of penguin }i \]

and:

\[ X_i=\text{body mass of penguin }i. \]

The Bayesian simple Normal regression model is:

\[ Y_i\mid\beta_0,\beta_1,\sigma \stackrel{ind}{\sim} N(\beta_0+\beta_1X_i,\sigma^2). \]

data(penguins_bayes)

penguins_mass <- penguins_bayes %>%
  drop_na(
    flipper_length_mm,
    body_mass_g
  )

nrow(penguins_bayes)
## [1] 344
nrow(penguins_mass)
## [1] 342

The original dataset contains 344 penguins. Observations missing flipper length or body mass are omitted from the fitted model.

2 1. Interpretation of the Slope Prior

The researchers specify:

\[ \beta_1\sim N(0.01,0.002^2). \]

The prior mean of 0.01 indicates that researchers expect typical flipper length to increase by 0.01 mm for every additional gram of body mass.

For an additional 1,000 grams, the expected increase is:

\[ 1000(0.01)=10\text{ mm}. \]

An approximate 95% prior interval is:

\[ 0.01\pm1.96(0.002) = (0.00608,0.01392). \]

Expressed per 1,000 grams, researchers believe that typical flipper length likely increases by approximately 6.08 to 13.92 mm.

Therefore, researchers are quite certain that body mass and flipper length have a positive relationship.

prior_slopes <- seq(
  0.002,
  0.018,
  length.out = 1000
)

plot(
  prior_slopes,
  dnorm(
    prior_slopes,
    mean = 0.01,
    sd = 0.002
  ),
  type = "l",
  lwd = 3,
  col = "steelblue",
  main = "Prior Model for the Body-Mass Coefficient",
  xlab = expression(
    beta[1]~"(mm per gram)"
  ),
  ylab = "Density"
)

abline(
  v = 0.01,
  col = "darkred",
  lty = 2,
  lwd = 2
)

3 2. Observed Relationship

ggplot(
  penguins_mass,
  aes(
    x = body_mass_g,
    y = flipper_length_mm
  )
) +
  geom_point(
    alpha = 0.65,
    color = "steelblue"
  ) +
  geom_smooth(
    method = "lm",
    se = TRUE,
    color = "darkred"
  ) +
  labs(
    title = "Penguin Flipper Length and Body Mass",
    x = "Body mass (grams)",
    y = "Flipper length (mm)"
  ) +
  theme_minimal(base_size = 12)

observed_correlation <- cor(
  penguins_mass$body_mass_g,
  penguins_mass$flipper_length_mm
)

classical_mass_model <- lm(
  flipper_length_mm ~ body_mass_g,
  data = penguins_mass
)

observed_slope <- coef(
  classical_mass_model
)[["body_mass_g"]]

observed_correlation
## [1] 0.8712018
observed_slope
## [1] 0.01527592
summary(classical_mass_model)
## 
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g, data = penguins_mass)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.7626  -4.9138   0.9891   5.1166  16.6392 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.367e+02  1.997e+00   68.47   <2e-16 ***
## body_mass_g 1.528e-02  4.668e-04   32.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.913 on 340 degrees of freedom
## Multiple R-squared:  0.759,  Adjusted R-squared:  0.7583 
## F-statistic:  1071 on 1 and 340 DF,  p-value: < 2.2e-16

The scatterplot displays a strong, positive, approximately linear association. Heavier penguins generally have longer flippers.

The observed correlation is 0.871. The ordinary least-squares slope is approximately 0.01528 mm per gram, which is equivalent to approximately 15.28 mm per 1,000 grams.

The visible grouping may partly reflect differences among penguin species. Nevertheless, simple Normal regression provides a useful initial model of the overall relationship.

4 3. Comparing Residual Variability

The parameter \(\sigma\) measures the typical amount by which observed flipper lengths differ from the regression mean among penguins with the same predictor value.

A stronger predictor should leave less unexplained variation and therefore produce a smaller residual standard deviation.

Both preliminary models below use the same complete observations.

penguins_comparison <- penguins_bayes %>%
  drop_na(
    flipper_length_mm,
    body_mass_g,
    bill_length_mm
  )

bill_lm <- lm(
  flipper_length_mm ~ bill_length_mm,
  data = penguins_comparison
)

mass_lm <- lm(
  flipper_length_mm ~ body_mass_g,
  data = penguins_comparison
)

sigma_comparison <- tibble(
  Predictor = c(
    "Bill length",
    "Body mass"
  ),
  Correlation = c(
    cor(
      penguins_comparison$bill_length_mm,
      penguins_comparison$flipper_length_mm
    ),
    cor(
      penguins_comparison$body_mass_g,
      penguins_comparison$flipper_length_mm
    )
  ),
  Residual_SD = c(
    sigma(bill_lm),
    sigma(mass_lm)
  )
)

knitr::kable(
  sigma_comparison,
  digits = 3,
  caption = "Bill Length and Body Mass Comparison"
)
Bill Length and Body Mass Comparison
Predictor Correlation Residual_SD
Bill length 0.656 10.627
Body mass 0.871 6.913

Body mass has the stronger absolute correlation with flipper length and produces the smaller residual standard deviation.

Therefore, \(\sigma\) should be larger when bill length is the predictor and smaller when body mass is the predictor. Body mass explains more of the differences in flipper length, leaving less unexplained variation around the regression line.

5 4. Bayesian Posterior Simulation

The researchers’ informative slope prior is:

\[ \beta_1\sim N(0.01,0.002^2). \]

Weakly informative priors are used for the centered intercept and residual standard deviation.

The slope prior uses autoscale = FALSE so that it remains exactly \(N(0.01,0.002^2)\). The other priors use autoscale = TRUE so that rstanarm adjusts them to the observed variable scales.

mass_model <- stan_glm(
  flipper_length_mm ~ body_mass_g,
  data = penguins_mass,
  family = gaussian,
  prior_intercept = normal(
    location = 200,
    scale = 2.5,
    autoscale = TRUE
  ),
  prior = normal(
    location = 0.01,
    scale = 0.002,
    autoscale = FALSE
  ),
  prior_aux = exponential(
    rate = 1,
    autoscale = TRUE
  ),
  chains = 4,
  iter = 10000,
  seed = 84735,
  refresh = 0
)

prior_summary(mass_model)
## Priors for model 'mass_model' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 200, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 200, scale = 35)
## 
## Coefficients
##  ~ normal(location = 0.01, scale = 0.002)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.071)
## ------
## See help('prior_summary.stanreg') for more details

Each of the four MCMC chains contains 10,000 iterations.

5.1 MCMC Diagnostics

print(
  summary(mass_model),
  digits = 3
)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      flipper_length_mm ~ body_mass_g
##  algorithm:    sampling
##  sample:       20000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 342
##  predictors:   2
## 
## Estimates:
##               mean    sd      10%     50%     90%  
## (Intercept) 137.874   1.954 135.365 137.866 140.397
## body_mass_g   0.015   0.000   0.014   0.015   0.016
## sigma         6.935   0.269   6.598   6.925   7.288
## 
## Fit Diagnostics:
##            mean    sd      10%     50%     90%  
## mean_PPD 200.912   0.532 200.236 200.913 201.588
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse  Rhat  n_eff
## (Intercept)   0.014 1.000 20525
## body_mass_g   0.000 1.000 20583
## sigma         0.002 1.000 19425
## mean_PPD      0.004 1.000 19900
## log-posterior 0.013 1.001  8582
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
mcmc_trace(
  as.array(mass_model),
  pars = c(
    "(Intercept)",
    "body_mass_g",
    "sigma"
  )
)

The trace plots should show stable, overlapping chains without long-term trends or separation. Effective sample sizes should be sufficiently large, and split-\(\widehat R\) values should be close to 1.00.

These results indicate whether the MCMC simulation mixed adequately and can be used for posterior inference.

5.2 Posterior Summary

posterior_summary <- tidy(
  mass_model,
  conf.int = TRUE,
  conf.level = 0.90
)

knitr::kable(
  posterior_summary,
  digits = 4,
  caption = "Posterior Summary with 90% Credible Intervals"
)
Posterior Summary with 90% Credible Intervals
term estimate std.error conf.low conf.high
(Intercept) 137.866 1.9522 134.6795 141.0915
body_mass_g 0.015 0.0005 0.0142 0.0158 
slope_summary <- posterior_summary %>%
  filter(term == "body_mass_g")

6 5. Posterior Body-Mass Coefficient

posterior_draws <- as.data.frame(
  mass_model
)

posterior_slopes <- posterior_draws$body_mass_g

posterior_slope_median <- median(
  posterior_slopes
)

posterior_slope_interval <- quantile(
  posterior_slopes,
  probs = c(0.05, 0.95)
)

posterior_probability_positive <- mean(
  posterior_slopes > 0
)

posterior_slope_median
## [1] 0.01500204
posterior_slope_interval
##         5%        95% 
## 0.01424679 0.01575491
posterior_probability_positive
## [1] 1
ggplot(
  data.frame(
    body_mass_coefficient = posterior_slopes
  ),
  aes(x = body_mass_coefficient)
) +
  geom_density(
    fill = "orange",
    color = "darkorange4",
    alpha = 0.40,
    linewidth = 1
  ) +
  geom_vline(
    xintercept = posterior_slope_median,
    color = "darkred",
    linetype = "dashed",
    linewidth = 1
  ) +
  labs(
    title = "Posterior Model of the Body-Mass Coefficient",
    x = expression(
      beta[1]~"(mm per gram)"
    ),
    y = "Density"
  ) +
  theme_minimal(base_size = 12)

The posterior median body-mass coefficient is 0.015 mm per gram.

The central 90% posterior credible interval is:

\[ \left( 0.01425, \; 0.01575 \right) \]

millimeters per gram.

The posterior median implies that an additional 1,000 grams of body mass is associated with an increase of approximately 15 mm in typical flipper length.

There is a 90% posterior probability that the corresponding increase per 1,000 grams lies between approximately 14.25 and 15.75 mm, conditional on the model and prior assumptions.

The posterior probability that \(\beta_1>0\) is 1. Therefore, the posterior provides ample evidence of a positive association between body mass and flipper length.

6.1 Prior and Posterior Comparison

set.seed(84735)

prior_draws <- rnorm(
  40000,
  mean = 0.01,
  sd = 0.002
)

comparison_data <- bind_rows(
  tibble(
    slope = prior_draws,
    Model = "Prior"
  ),
  tibble(
    slope = posterior_slopes,
    Model = "Posterior"
  )
)

ggplot(
  comparison_data,
  aes(
    x = slope,
    color = Model,
    fill = Model
  )
) +
  geom_density(
    alpha = 0.20,
    linewidth = 1
  ) +
  labs(
    title = "Prior and Posterior Body-Mass Coefficients",
    x = expression(
      beta[1]~"(mm per gram)"
    ),
    y = "Density"
  ) +
  theme_minimal(base_size = 12)

The researchers began with a strong prior belief that the slope was centered at 0.01 mm per gram. The observed data also demonstrate a positive association.

A change in the center from 0.01 indicates that the data modified the researchers’ estimate of the magnitude of the relationship. A narrower posterior distribution indicates that the observed data increased their certainty.

Overall, the posterior preserves the researchers’ belief in a positive relationship while updating the size and uncertainty of the body-mass coefficient using the observed penguin data.

7 Conclusion

The observed data show a strong positive association between penguin body mass and flipper length. Heavier penguins tend to have longer flippers.

Body mass produces less residual variation than bill length in a comparable simple regression, indicating that body mass is the stronger of these two individual predictors of flipper length.

The Bayesian Normal regression combines the positive informative prior with the observed evidence. The posterior slope distribution quantifies both the likely size of the relationship and the remaining uncertainty.

This analysis identifies an association and does not prove that increased body mass causes longer flippers.

8 AI-Use Disclosure

I used OpenAI’s ChatGPT/Codex to assist with this assignment. The tool helped me consult Exercises 9.18–9.20 in Chapter 9 of Bayes Rules!, interpret the questions, construct the required stan_glm() model, develop the R code, organize the R Markdown document, and formulate the statistical explanations.

I reviewed the proposed calculations, model specification, code, plots, and interpretations. I remain responsible for verifying the knitted results and the accuracy of the submitted work.

8.1 Interaction Transcript

User: Complete Exercise 9.20 given at the end of Chapter 9 of the Bayes Rules! book. Consult Exercises 9.19 and 9.18 for additional information. Submit an HTML file knitted from an Rmd file and the original Rmd containing all code and explanations. Do not submit a zipped archive. If an AI tool is used, provide the complete transcript of the interactions and explain how the solution was reached using the tool’s help.

AI assistant: I consulted the Chapter 9 exercise descriptions. Exercise 9.18 provides the earlier posterior analysis of flipper length by bill length. Exercise 9.19 extends that model to prediction. Exercise 9.20 changes the predictor to body mass. I structured the solution to interpret the informative slope prior, plot and assess the observed relationship, compare residual variation using bill length and body mass, fit the requested four-chain Bayesian Normal regression, check the MCMC simulation, and examine the posterior body-mass coefficient.

User: Exercise 9.18 asks for a posterior analysis of flipper length by bill length.

AI assistant: I explained that Exercise 9.18 provides information about the bill-length model and its residual standard deviation, which is needed for the comparison in Exercise 9.20.

User: Exercise 9.19 asks for posterior prediction of Pablo’s flipper length when his bill is 51 mm long.

AI assistant: I explained the distinction between the posterior model for a typical flipper length and the posterior predictive model for one individual penguin.

User: Exercise 9.20 requests a regression of flipper length by body mass using the informative slope prior and weakly informative priors for the centered intercept and residual standard deviation.

AI assistant: I created the complete R Markdown analysis shown above. The researchers’ exact \(N(0.01,0.002^2)\) prior is retained by setting autoscale = FALSE. Weakly informative autoscaled priors are used for the centered intercept and residual standard deviation. The model uses four chains with 10,000 iterations each.

Any subsequent AI interactions used to troubleshoot or revise this document should be added to the transcript before submission.

9 Reference

Johnson, A. A., Ott, M. Q., and Dogucu, M. Bayes Rules! An Introduction to Applied Bayesian Modeling. Chapter 9, Exercise 9.20.
https://www.bayesrulesbook.com/chapter-9#exercises-8