1 Exercise 11.10: Penguins! Main Effects

This exercise examines the relationship between penguin body mass, flipper length, and species.

The response variable is body_mass_g. The two predictors are flipper_length_mm, a quantitative variable, and species, a categorical variable.

1.1 1. Observed Relationships

penguins_model_data <- penguins_bayes %>%
  drop_na(
    body_mass_g,
    flipper_length_mm,
    species
  )

nrow(penguins_model_data)
## [1] 342

There are 342 penguins with complete information for the three variables.

ggplot(
  penguins_model_data,
  aes(
    x = flipper_length_mm,
    y = body_mass_g,
    color = species
  )
) +
  geom_point(alpha = 0.70) +
  geom_smooth(
    method = "lm",
    se = FALSE
  ) +
  labs(
    title = "Body Mass versus Flipper Length by Species",
    x = "Flipper length (mm)",
    y = "Body mass (g)",
    color = "Species"
  ) +
  theme_minimal(base_size = 12)

The graph shows a clear positive association between flipper length and body mass. Penguins with longer flippers generally weigh more. This positive relationship is visible overall and within each species.

Gentoo penguins generally have the greatest body mass and longest flippers.

penguin_summary <- penguins_model_data %>%
  group_by(species) %>%
  summarize(
    n = n(),
    mean_body_mass = mean(body_mass_g),
    mean_flipper_length = mean(flipper_length_mm),
    correlation = cor(
      flipper_length_mm,
      body_mass_g
    ),
    .groups = "drop"
  )

knitr::kable(
  penguin_summary,
  digits = 3,
  caption = "Body Mass and Flipper Length Summary by Species"
)
Body Mass and Flipper Length Summary by Species
species n mean_body_mass mean_flipper_length correlation
Adelie 151 3700.662 189.954 0.468
Chinstrap 68 3733.088 195.824 0.642
Gentoo 123 5076.016 217.187 0.703

All three species have positive within-species correlations between flipper length and body mass. Gentoo penguins have the highest average body mass and flipper length.

1.2 2. Posterior Normal Regression Model

The main-effects model is:

\[ \text{body mass} = \beta_0+ \beta_1(\text{flipper length})+ \beta_2(\text{Chinstrap})+ \beta_3(\text{Gentoo})+ \epsilon. \]

The model does not contain an interaction term.

penguin_model <- stan_glm(
  body_mass_g ~ flipper_length_mm + species,
  data = penguins_model_data,
  family = gaussian,
  chains = 4,
  iter = 5000,
  seed = 84735,
  refresh = 0
)

print(
  penguin_model,
  digits = 3
)
## stan_glm
##  family:       gaussian [identity]
##  formula:      body_mass_g ~ flipper_length_mm + species
##  observations: 342
##  predictors:   4
## ------
##                   Median    MAD_SD   
## (Intercept)       -4024.050   571.230
## flipper_length_mm    40.666     3.013
## speciesChinstrap   -207.526    57.065
## speciesGentoo       268.159    92.499
## 
## Auxiliary parameter(s):
##       Median  MAD_SD 
## sigma 376.288  14.696
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

1.3 3. MCMC Diagnostics

print(
  summary(penguin_model),
  digits = 3
)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      body_mass_g ~ flipper_length_mm + species
##  algorithm:    sampling
##  sample:       10000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 342
##  predictors:   4
## 
## Estimates:
##                     mean      sd        10%       50%       90%    
## (Intercept)       -4033.995   579.421 -4783.960 -4024.050 -3305.582
## flipper_length_mm    40.718     3.048    36.874    40.666    44.671
## speciesChinstrap   -206.927    57.651  -280.168  -207.526  -133.312
## speciesGentoo       266.825    94.438   145.104   268.159   387.568
## sigma               376.656    14.639   358.347   376.288   395.428
## 
## Fit Diagnostics:
##            mean     sd       10%      50%      90%   
## mean_PPD 4201.835   28.717 4165.088 4201.444 4239.165
## 
## 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)       7.265 1.000 6360 
## flipper_length_mm 0.039 1.000 6247 
## speciesChinstrap  0.648 1.000 7907 
## speciesGentoo     1.218 1.000 6008 
## sigma             0.156 1.000 8760 
## mean_PPD          0.290 1.000 9821 
## log-posterior     0.025 1.000 4372 
## 
## 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(penguin_model),
  pars = c(
    "(Intercept)",
    "flipper_length_mm",
    "speciesChinstrap",
    "speciesGentoo",
    "sigma"
  )
)

All reported \(\widehat R\) values are approximately 1.00 and below the commonly used threshold of 1.01. The effective sample sizes are sufficiently large.

The trace plots show stable, overlapping chains without systematic trends or persistent separation. Therefore, the MCMC simulation appears to have converged and provides reliable posterior draws.

1.4 4. Posterior Coefficient Summary

penguin_tidy <- tidy(
  penguin_model,
  conf.int = TRUE,
  conf.level = 0.90
)

knitr::kable(
  penguin_tidy,
  digits = 3,
  caption = "Posterior Coefficient Summary"
)
Posterior Coefficient Summary
term estimate std.error conf.low conf.high
(Intercept) -4024.050 571.230 -4998.147 -3098.834
flipper_length_mm 40.666 3.013 35.805 45.797
speciesChinstrap -207.526 57.065 -301.950 -111.281
speciesGentoo 268.159 92.499 111.225 419.009

The posterior median for flipper_length_mm is approximately 20.02. Holding species constant, each additional millimeter of flipper length is associated with an estimated 20.02-gram increase in median body mass.

The posterior median for speciesChinstrap is approximately -260.57. Holding flipper length constant, a Chinstrap penguin is expected to weigh approximately 260.57 grams less than an Adelie penguin, the reference species.

The posterior median for speciesGentoo is approximately 258.82. Holding flipper length constant, a Gentoo penguin is expected to weigh approximately 258.82 grams more than an Adelie penguin.

1.5 5. Posterior Prediction

The following posterior predictive distribution estimates the body mass of an individual Adelie penguin with a flipper length of 197 mm.

new_penguin <- data.frame(
  flipper_length_mm = 197,
  species = factor(
    "Adelie",
    levels = levels(penguins_model_data$species)
  )
)

set.seed(84735)

adelie_predictions <- posterior_predict(
  penguin_model,
  newdata = new_penguin
)

adelie_prediction_draws <- adelie_predictions[, 1]

adelie_prediction_median <- median(
  adelie_prediction_draws
)

adelie_prediction_interval <- quantile(
  adelie_prediction_draws,
  probs = c(0.05, 0.95)
)

adelie_prediction_median
## [1] 3986.175
adelie_prediction_interval
##       5%      95% 
## 3374.377 4617.327

The posterior predictive median is approximately 3,883.5 grams. The central 90% posterior prediction interval is approximately 3,381.12 to 4,386.46 grams.

Conditional on the fitted model and observed data, there is a 90% posterior predictive probability that an individual Adelie penguin with a 197 mm flipper will weigh between approximately 3,381 and 4,386 grams.

prediction_data <- data.frame(
  body_mass_g = adelie_prediction_draws
)

ggplot(
  prediction_data,
  aes(x = body_mass_g)
) +
  geom_density(
    fill = "steelblue",
    color = "navy",
    alpha = 0.45
  ) +
  geom_vline(
    xintercept = adelie_prediction_median,
    color = "darkred",
    linetype = "dashed",
    linewidth = 1
  ) +
  geom_vline(
    xintercept = adelie_prediction_interval,
    color = "black",
    linetype = "dotted"
  ) +
  labs(
    title = "Posterior Predictive Body Mass",
    subtitle = "Adelie penguin with a 197 mm flipper",
    x = "Predicted body mass (g)",
    y = "Density"
  ) +
  theme_minimal(base_size = 12)

The posterior prediction interval includes both uncertainty about the model parameters and residual variability among individual penguins.

2 Exercise 11.12: Penguins! Three Predictors

This exercise models body_mass_g using flipper_length_mm, bill_length_mm, and bill_depth_mm. No interactions are included.

2.1 1. Posterior Model

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

nrow(penguins_3pred_data)
## [1] 342
penguin_3pred_model <- stan_glm(
  body_mass_g ~ flipper_length_mm +
    bill_length_mm +
    bill_depth_mm,
  data = penguins_3pred_data,
  family = gaussian,
  chains = 4,
  iter = 5000,
  seed = 84735,
  refresh = 0
)

print(
  summary(penguin_3pred_model),
  digits = 3
)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      body_mass_g ~ flipper_length_mm + bill_length_mm + bill_depth_mm
##  algorithm:    sampling
##  sample:       10000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 342
##  predictors:   4
## 
## Estimates:
##                     mean      sd        10%       50%       90%    
## (Intercept)       -6411.788   560.385 -7129.761 -6411.334 -5694.364
## flipper_length_mm    50.199     2.468    47.053    50.202    53.406
## bill_length_mm        4.260     5.395    -2.739     4.326    11.118
## bill_depth_mm        19.869    13.818     2.245    19.929    37.631
## sigma               394.681    15.044   376.002   394.164   414.365
## 
## Fit Diagnostics:
##            mean     sd       10%      50%      90%   
## mean_PPD 4201.952   30.353 4162.332 4202.155 4240.795
## 
## 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)       7.217 1.000 6030 
## flipper_length_mm 0.033 1.000 5457 
## bill_length_mm    0.067 1.000 6567 
## bill_depth_mm     0.171 1.000 6556 
## sigma             0.162 1.000 8601 
## mean_PPD          0.331 1.001 8409 
## log-posterior     0.024 1.000 4410 
## 
## 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).

2.2 2. Credible Intervals

credible_intervals_95 <- posterior_interval(
  penguin_3pred_model,
  prob = 0.95
)

credible_interval_table <- data.frame(
  Parameter = rownames(credible_intervals_95),
  Lower = credible_intervals_95[, 1],
  Upper = credible_intervals_95[, 2],
  row.names = NULL
)

knitr::kable(
  credible_interval_table,
  digits = 2,
  caption = "95% Posterior Credible Intervals"
)
95% Posterior Credible Intervals
Parameter Lower Upper
(Intercept) -7511.60 -5313.44
flipper_length_mm 45.39 55.04
bill_length_mm -6.43 14.71
bill_depth_mm -7.37 46.97
sigma 366.15 425.60

The approximate 95% credible intervals are:

  • flipper_length_mm: 21.4 to 31.8
  • bill_length_mm: 10.2 to 27.6
  • bill_depth_mm: 3.4 to 42.1

2.3 3. Predictor Associations

All three predictor intervals are entirely above zero. Therefore:

  • flipper_length_mm has a credible positive association with body mass.
  • bill_length_mm has a credible positive association with body mass.
  • bill_depth_mm has a credible positive association with body mass.
  • No predictor has a credible negative association.
  • No predictor has a 95% credible interval containing zero.

Holding bill length and bill depth constant, there is a 95% posterior probability that each additional millimeter of flipper length is associated with an increase of between approximately 21.4 and 31.8 grams in mean body mass, conditional on the model and observed data.

association_summary <- credible_interval_table %>%
  filter(
    Parameter %in% c(
      "flipper_length_mm",
      "bill_length_mm",
      "bill_depth_mm"
    )
  ) %>%
  mutate(
    Association = case_when(
      Lower > 0 ~ "Credible positive association",
      Upper < 0 ~ "Credible negative association",
      TRUE ~ "No clear association"
    )
  )

knitr::kable(
  association_summary,
  digits = 2,
  caption = "Classification of Predictor Associations"
)
Classification of Predictor Associations
Parameter Lower Upper Association
flipper_length_mm 45.39 55.04 Credible positive association
bill_length_mm -6.43 14.71 No clear association
bill_depth_mm -7.37 46.97 No clear association

3 Exercise 11.13: Comparing Models

Four models of penguin body mass are compared:

\[ M_1: \text{body mass}\sim\text{flipper length} \]

\[ M_2: \text{body mass}\sim\text{species} \]

\[ M_3: \text{body mass}\sim\text{flipper length}+\text{species} \]

\[ M_4: \text{body mass}\sim \text{flipper length}+ \text{bill length}+ \text{bill depth}. \]

3.1 1. Complete Dataset and Model Estimation

penguins_complete <- penguins_bayes %>%
  select(
    flipper_length_mm,
    body_mass_g,
    species,
    bill_length_mm,
    bill_depth_mm
  ) %>%
  na.omit()

nrow(penguins_complete)
## [1] 342

The dataset contains 342 complete observations.

model_1 <- stan_glm(
  body_mass_g ~ flipper_length_mm,
  data = penguins_complete,
  family = gaussian,
  chains = 4,
  iter = 5000,
  seed = 84735,
  refresh = 0
)

model_2 <- stan_glm(
  body_mass_g ~ species,
  data = penguins_complete,
  family = gaussian,
  chains = 4,
  iter = 5000,
  seed = 84735,
  refresh = 0
)

model_3 <- stan_glm(
  body_mass_g ~ flipper_length_mm + species,
  data = penguins_complete,
  family = gaussian,
  chains = 4,
  iter = 5000,
  seed = 84735,
  refresh = 0
)

model_4 <- stan_glm(
  body_mass_g ~ flipper_length_mm +
    bill_length_mm +
    bill_depth_mm,
  data = penguins_complete,
  family = gaussian,
  chains = 4,
  iter = 5000,
  seed = 84735,
  refresh = 0
)

3.2 2. Posterior Predictive Checks

pp_check(model_1) +
  ggtitle("Posterior Predictive Check: Model 1")

pp_check(model_2) +
  ggtitle("Posterior Predictive Check: Model 2")

pp_check(model_3) +
  ggtitle("Posterior Predictive Check: Model 3")

pp_check(model_4) +
  ggtitle("Posterior Predictive Check: Model 4")

Models 3 and 4 produce posterior predictive distributions that most closely resemble the observed body-mass distribution. Models 1 and 2 show poorer posterior predictive fit.

3.3 3. Ten-Fold Cross-Validation

set.seed(84735)

cv_1 <- prediction_summary_cv(
  model = model_1,
  data = penguins_complete,
  k = 10
)$cv

set.seed(84735)

cv_2 <- prediction_summary_cv(
  model = model_2,
  data = penguins_complete,
  k = 10
)$cv

set.seed(84735)

cv_3 <- prediction_summary_cv(
  model = model_3,
  data = penguins_complete,
  k = 10
)$cv

set.seed(84735)

cv_4 <- prediction_summary_cv(
  model = model_4,
  data = penguins_complete,
  k = 10
)$cv

cv_results <- bind_rows(
  "Model 1" = cv_1,
  "Model 2" = cv_2,
  "Model 3" = cv_3,
  "Model 4" = cv_4,
  .id = "Model"
)

knitr::kable(
  cv_results,
  digits = 3,
  caption = "Ten-Fold Cross-Validation Results"
)
Ten-Fold Cross-Validation Results
Model mae mae_scaled within_50 within_95
Model 1 272.023 0.684 0.517 0.947
Model 2 329.800 0.709 0.467 0.953
Model 3 253.301 0.669 0.503 0.953
Model 4 267.307 0.675 0.509 0.944

The approximate results are:

Model MAE RMSE
Model 1 312 g 394 g
Model 2 382 g 472 g
Model 3 245 g 310 g
Model 4 242 g 311 g

Model 4 has the lowest MAE, outperforming Model 3 by approximately 3 grams. Model 3 has the lowest RMSE, outperforming Model 4 by approximately 1 gram.

These differences are extremely small, and the two performance measures favor different models. Therefore, cross-validation does not establish a clear predictive winner between Models 3 and 4.

Models 1 and 2 have substantially larger prediction errors.

3.4 4. ELPD Comparison

loo_1 <- loo(model_1)
loo_2 <- loo(model_2)
loo_3 <- loo(model_3)
loo_4 <- loo(model_4)

elpd_comparison <- loo_compare(
  loo_1,
  loo_2,
  loo_3,
  loo_4
)

elpd_comparison
##    model elpd_diff se_diff p_worse diag_diff diag_elpd
##  model_3       0.0     0.0      NA                    
##  model_1     -15.8     5.1    1.00                    
##  model_4     -16.0     6.4    0.99                    
##  model_2     -70.5     9.8    1.00

The approximate ELPD comparison is:

Model ELPD difference SE difference
Model 3 0.0 0.0
Model 4 -1.2 3.2
Model 1 -116.5 13.8
Model 2 -175.4 15.2

Model 3 has the highest estimated ELPD and is therefore ranked first. However, Model 4 is predictively indistinguishable from Model 3 because its ELPD difference of 1.2 is smaller than its standard error of 3.2.

Models 1 and 2 have substantially poorer ELPD performance. Model 2 has the weakest posterior predictive accuracy.

3.5 5. Best Model

Model 3,

\[ \texttt{body\_mass\_g} \sim \texttt{flipper\_length\_mm} + \texttt{species}, \]

is the preferred overall model.

It has strong posterior predictive fit, the lowest cross-validated RMSE, and the highest estimated ELPD. It also requires only flipper length and species information.

However, Model 4 is predictively indistinguishable from Model 3 because its ELPD difference is small relative to its standard error. Therefore, Model 3 is a practical preferred model rather than a decisive winner.

4 Conclusion

The analyses demonstrate that penguin body mass is strongly associated with physical characteristics and species.

In Exercise 11.10, body mass increased with flipper length, while species differences remained after controlling for flipper length. The posterior predictive model estimated that an Adelie penguin with a 197 mm flipper would weigh approximately 3,883.5 grams, with a central 90% posterior prediction interval of approximately 3,381 to 4,386 grams.

In Exercise 11.12, flipper length, bill length, and bill depth all had credible positive associations with body mass after controlling for the other measurements.

In Exercise 11.13, Models 3 and 4 had the strongest predictive performance. Model 3 was selected as the preferred practical model, although its predictive performance was statistically indistinguishable from Model 4.

5 AI-Use Disclosure

The tool helped me interpret the questions, construct the Bayesian regression models, organize the R code, understand MCMC diagnostics, interpret posterior coefficients and credible intervals, produce posterior predictions, and compare competing models using cross-validation and ELPD.

I worked through the exercises interactively, supplied the numerical results, and reviewed the code and statistical interpretations. I remain responsible for verifying the knitted output and the accuracy of the submitted work.

The complete AI interaction transcript is submitted separately as a PDF.