install.packages("palmerpenguins")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
library(tidyverse)
library(palmerpenguins)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ── ✔ dplyr 1.2.1 ✔ readr 2.2.0 ✔ forcats 1.0.1 ✔ stringr 1.6.0 ✔ ggplot2 4.0.3 ✔ tibble 3.3.1 ✔ lubridate 1.9.5 ✔ tidyr 1.3.2 ✔ purrr 1.2.2 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag() ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors Attaching package: ‘palmerpenguins’ The following objects are masked from ‘package:datasets’: penguins, penguins_raw
head(penguins)| species | island | bill_length_mm | bill_depth_mm | flipper_length_mm | body_mass_g | sex | year |
|---|---|---|---|---|---|---|---|
| <fct> | <fct> | <dbl> | <dbl> | <int> | <int> | <fct> | <int> |
| Adelie | Torgersen | 39.1 | 18.7 | 181 | 3750 | male | 2007 |
| Adelie | Torgersen | 39.5 | 17.4 | 186 | 3800 | female | 2007 |
| Adelie | Torgersen | 40.3 | 18.0 | 195 | 3250 | female | 2007 |
| Adelie | Torgersen | NA | NA | NA | NA | NA | 2007 |
| Adelie | Torgersen | 36.7 | 19.3 | 193 | 3450 | female | 2007 |
| Adelie | Torgersen | 39.3 | 20.6 | 190 | 3650 | male | 2007 |
dim(penguins)- 344
- 8
install.packages("GGally")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘patchwork’, ‘ggstats’
library(GGally)
penguins_clean <- na.omit(penguins)
ggpairs(
penguins_clean,
columns = 3:6, # Select specific numeric columns (flipper length, bill length, etc.)
aes(color = species, alpha = 0.7)
) +
theme_classic(base_size = 12)
glimpse(penguins)
# Look at the raw dataset (17 columns with messy names)
glimpse(penguins_raw)Rows: 344 Columns: 8 $ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel… $ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse… $ bill_length_mm <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, … $ bill_depth_mm <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, … $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186… $ body_mass_g <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, … $ sex <fct> male, female, female, NA, female, male, female, male… $ year <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007… Rows: 344 Columns: 17 $ studyName <chr> "PAL0708", "PAL0708", "PAL0708", "PAL0708", "PAL… $ `Sample Number` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1… $ Species <chr> "Adelie Penguin (Pygoscelis adeliae)", "Adelie P… $ Region <chr> "Anvers", "Anvers", "Anvers", "Anvers", "Anvers"… $ Island <chr> "Torgersen", "Torgersen", "Torgersen", "Torgerse… $ Stage <chr> "Adult, 1 Egg Stage", "Adult, 1 Egg Stage", "Adu… $ `Individual ID` <chr> "N1A1", "N1A2", "N2A1", "N2A2", "N3A1", "N3A2", … $ `Clutch Completion` <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "No", … $ `Date Egg` <date> 2007-11-11, 2007-11-11, 2007-11-16, 2007-11-16,… $ `Culmen Length (mm)` <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34… $ `Culmen Depth (mm)` <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18… $ `Flipper Length (mm)` <dbl> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190,… $ `Body Mass (g)` <dbl> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 34… $ Sex <chr> "MALE", "FEMALE", "FEMALE", NA, "FEMALE", "MALE"… $ `Delta 15 N (o/oo)` <dbl> NA, 8.94956, 8.36821, NA, 8.76651, 8.66496, 9.18… $ `Delta 13 C (o/oo)` <dbl> NA, -24.69454, -25.33302, NA, -25.32426, -25.298… $ Comments <chr> "Not enough blood for isotopes.", NA, NA, "Adult…
tail(penguins_raw)| studyName | Sample Number | Species | Region | Island | Stage | Individual ID | Clutch Completion | Date Egg | Culmen Length (mm) | Culmen Depth (mm) | Flipper Length (mm) | Body Mass (g) | Sex | Delta 15 N (o/oo) | Delta 13 C (o/oo) | Comments |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <dbl> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <date> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <dbl> | <dbl> | <chr> |
| PAL0910 | 63 | Chinstrap penguin (Pygoscelis antarctica) | Anvers | Dream | Adult, 1 Egg Stage | N98A1 | Yes | 2009-11-19 | 45.7 | 17.0 | 195 | 3650 | FEMALE | 9.26715 | -24.31912 | NA |
| PAL0910 | 64 | Chinstrap penguin (Pygoscelis antarctica) | Anvers | Dream | Adult, 1 Egg Stage | N98A2 | Yes | 2009-11-19 | 55.8 | 19.8 | 207 | 4000 | MALE | 9.70465 | -24.53494 | NA |
| PAL0910 | 65 | Chinstrap penguin (Pygoscelis antarctica) | Anvers | Dream | Adult, 1 Egg Stage | N99A1 | No | 2009-11-21 | 43.5 | 18.1 | 202 | 3400 | FEMALE | 9.37608 | -24.40753 | Nest never observed with full clutch. |
| PAL0910 | 66 | Chinstrap penguin (Pygoscelis antarctica) | Anvers | Dream | Adult, 1 Egg Stage | N99A2 | No | 2009-11-21 | 49.6 | 18.2 | 193 | 3775 | MALE | 9.46180 | -24.70615 | Nest never observed with full clutch. |
| PAL0910 | 67 | Chinstrap penguin (Pygoscelis antarctica) | Anvers | Dream | Adult, 1 Egg Stage | N100A1 | Yes | 2009-11-21 | 50.8 | 19.0 | 210 | 4100 | MALE | 9.98044 | -24.68741 | NA |
| PAL0910 | 68 | Chinstrap penguin (Pygoscelis antarctica) | Anvers | Dream | Adult, 1 Egg Stage | N100A2 | Yes | 2009-11-21 | 50.2 | 18.7 | 198 | 3775 | FEMALE | 9.39305 | -24.25255 | NA |
install.packages("performance")
install.packages("see")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘bayestestR’, ‘insight’, ‘datawizard’
Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘correlation’, ‘effectsize’, ‘modelbased’, ‘patchwork’, ‘parameters’
library(performance)
library(see)model <- lm(body_mass_g ~ flipper_length_mm, data = penguins)model
Call:
lm(formula = body_mass_g ~ flipper_length_mm, data = penguins)
Coefficients:
(Intercept) flipper_length_mm
-5780.83 49.69 # 1. Install the missing Linux system library
system("apt-get update && apt-get install -y libfftw3-dev")
# 2. Now install qqplotr again
install.packages("qqplotr")
# 3. Load it up
library(qqplotr)Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependency ‘qqconf’
Attaching package: ‘qqplotr’
The following objects are masked from ‘package:ggplot2’:
stat_qq_line, StatQqLine
# Run this once to install the package
install.packages("qqplotr")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘bitops’, ‘iterators’, ‘DEoptimR’, ‘caTools’, ‘pracma’, ‘twosamples’, ‘doParallel’, ‘pbmcapply’, ‘foreach’, ‘robustbase’, ‘opdisDownsampling’, ‘qqconf’
Warning message in install.packages("qqplotr"):
“installation of package ‘qqconf’ had non-zero exit status”
Warning message in install.packages("qqplotr"):
“installation of package ‘qqplotr’ had non-zero exit status”options(repr.plot.width = 16, repr.plot.height = 10)
plot(check_model(model), n_columns = 2, theme = theme_classic(base_size = 16))Ignoring unknown labels: • size : ""
# Drop missing values
df <- penguins %>% drop_na(body_mass_g, flipper_length_mm)
# Fit model and add predictions/residuals to the dataframe
df_aug <- df %>% mutate(
predicted = predict(lm(body_mass_g ~ flipper_length_mm, data = .)),
residual = body_mass_g - predicted
)
# Plot the regression line AND the residual "stems"
ggplot(df_aug, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point(alpha = 0.5, color = "steelblue") +
geom_smooth(method = "lm", se = FALSE, color = "black") +
# This draws the vertical residual lines!
geom_segment(aes(xend = flipper_length_mm, yend = predicted), color = "red", linetype = "dashed") +
theme_classic(base_size = 14) +
labs(title = "Visualizing Residuals (Red Dashed Lines)",
x = "Flipper Length (mm)", y = "Body Mass (g)")`geom_smooth()` using formula = 'y ~ x'
# Fit a model and extract residuals
model <- lm(body_mass_g ~ flipper_length_mm, data = penguins)
res_data <- data.frame(residuals = residuals(model))
# Native ggplot2 Q-Q plot (no extra packages needed)
ggplot(res_data, aes(sample = residuals)) +
stat_qq(color = "darkblue", size = 2) +
stat_qq_line(color = "red", linetype = "dashed") +
theme_classic(base_size = 14) +
labs(
title = "Standard Q-Q Plot",
x = "Theoretical Quantiles",
y = "Sample Quantiles"
)
# Install and load
#install.packages("parameters")
library(parameters)
#library(palmerpenguins)
# Fit a model
model <- lm(body_mass_g ~ flipper_length_mm + species, data = penguins)
# Extract clean parameters
model_parameters(model)| Parameter | Coefficient | SE | CI | CI_low | CI_high | t | df_error | p |
|---|---|---|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | <dbl> |
| (Intercept) | -4031.4769 | 584.15134 | 0.95 | -5180.50685 | -2882.44693 | -6.901425 | 338 | 2.551453e-11 |
| flipper_length_mm | 40.7054 | 3.07102 | 0.95 | 34.66468 | 46.74612 | 13.254686 | 338 | 1.401600e-32 |
| speciesChinstrap | -206.5101 | 57.73065 | 0.95 | -320.06672 | -92.95352 | -3.577132 | 338 | 3.980986e-04 |
| speciesGentoo | 266.8096 | 95.26374 | 0.95 | 79.42513 | 454.19408 | 2.800747 | 338 | 5.391808e-03 |
library(performance)
# Get R-squared, AIC, RMSE, and other fit statistics
model_performance(model)| AIC | AICc | BIC | R2 | R2_adjusted | RMSE | Sigma | |
|---|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 5031.523 | 5031.702 | 5050.697 | 0.7826479 | 0.7807187 | 373.3325 | 375.5351 |
summary(model)
Call:
lm(formula = body_mass_g ~ flipper_length_mm + species, data = penguins)
Residuals:
Min 1Q Median 3Q Max
-927.70 -254.82 -23.92 241.16 1191.68
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4031.477 584.151 -6.901 2.55e-11 ***
flipper_length_mm 40.705 3.071 13.255 < 2e-16 ***
speciesChinstrap -206.510 57.731 -3.577 0.000398 ***
speciesGentoo 266.810 95.264 2.801 0.005392 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 375.5 on 338 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.7826, Adjusted R-squared: 0.7807
F-statistic: 405.7 on 3 and 338 DF, p-value: < 2.2e-16# Forces R to calculate a separate intercept for every species
model_no_intercept <- lm(body_mass_g ~ 0 + flipper_length_mm + species, data = penguins)
model_parameters(model_no_intercept)| Parameter | Coefficient | SE | CI | CI_low | CI_high | t | df_error | p |
|---|---|---|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <int> | <dbl> |
| flipper_length_mm | 40.7054 | 3.07102 | 0.95 | 34.66468 | 46.74612 | 13.254686 | 338 | 1.401600e-32 |
| speciesAdelie | -4031.4769 | 584.15134 | 0.95 | -5180.50685 | -2882.44693 | -6.901425 | 338 | 2.551453e-11 |
| speciesChinstrap | -4237.9870 | 603.09976 | 0.95 | -5424.28866 | -3051.68536 | -7.027008 | 338 | 1.169686e-11 |
| speciesGentoo | -3764.6673 | 667.84449 | 0.95 | -5078.32228 | -2451.01229 | -5.637042 | 338 | 3.652962e-08 |
install.packages("easystats")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘patchwork’, ‘bayestestR’, ‘correlation’, ‘datawizard’, ‘effectsize’, ‘insight’, ‘modelbased’, ‘parameters’, ‘performance’, ‘report’, ‘see’
library(tidyverse)
library(palmerpenguins)
library(performance)
library(see)
# Ensure clean data
df <- penguins %>% drop_na()
# 1. Define competing models
model1 <- lm(body_mass_g ~ flipper_length_mm, data = df)
model2 <- lm(body_mass_g ~ flipper_length_mm + species, data = df)
model3 <- lm(body_mass_g ~ flipper_length_mm + species + bill_length_mm, data = df)
# 2. Compare them all at once (and rank them automatically!)
comparison <- compare_performance(model1, model2, model3, rank = TRUE)
# 3. View as a clean table or markdown
comparison| Name | Model | R2 | R2_adjusted | RMSE | Sigma | AIC_wt | AICc_wt | BIC_wt | Performance_Score |
|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> |
| model3 | lm | 0.8243043 | 0.8221617 | 337.0077 | 339.5666 | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | 1.0000000 |
| model2 | lm | 0.7870333 | 0.7850913 | 371.0352 | 373.2839 | 3.335124e-14 | 3.461151e-14 | 2.238926e-13 | 0.2210171 |
| model1 | lm | 0.7620922 | 0.7613734 | 392.1603 | 393.3433 | 2.418503e-21 | 2.652517e-21 | 7.316942e-19 | 0.0000000 |
# Set plot dimensions for Colab
options(repr.plot.width = 9, repr.plot.height = 6)
# Generate a visual comparison plot
plot(comparison) + theme_classic(base_size = 14)
# Extract the exact coefficients for Model 3
model3 <- lm(body_mass_g ~ flipper_length_mm + species + bill_length_mm, data = df)
coef(model3)- (Intercept)
- -3864.07323411886
- flipper_length_mm
- 27.5442861236828
- speciesChinstrap
- -732.416673207492
- speciesGentoo
- 113.2541752724
- bill_length_mm
- 60.1173245626618
# Extract coefficients from your model
b <- coef(model3)
# Automatically print your equation text using base R
cat(sprintf(
"Body Mass = %.2f + (%.2f * Flipper) + (%.2f * Bill) + (%.2f * Chinstrap) + (%.2f * Gentoo)\n",
b["(Intercept)"],
b["flipper_length_mm"],
b["bill_length_mm"],
b["speciesChinstrap"],
b["speciesGentoo"]
))Body Mass = -3864.07 + (27.54 * Flipper) + (60.12 * Bill) + (-732.42 * Chinstrap) + (113.25 * Gentoo)# Install (run once)
install.packages("equatiomatic")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘listenv’, ‘parallelly’, ‘future’, ‘globals’, ‘coda’, ‘furrr’, ‘broom.mixed’
library(equatiomatic)
# Fit Model 3
model3 <- lm(body_mass_g ~ flipper_length_mm + species + bill_length_mm, data = df)
# Generate the equation automatically!
eq <- extract_eq(model3, wrap = TRUE, terms_per_line = 2)
cat(as.character(eq))\begin{aligned}
\operatorname{body\_mass\_g} &= \alpha + \beta_{1}(\operatorname{flipper\_length\_mm})\ + \\
&\quad \beta_{2}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{3}(\operatorname{species}_{\operatorname{Gentoo}})\ + \\
&\quad \beta_{4}(\operatorname{bill\_length\_mm}) + \epsilon
\end{aligned}\[ \begin{aligned} \operatorname{body\_mass\_g} &= \alpha + \beta_{1}(\operatorname{flipper\_length\_mm})\ + \\ &\quad \beta_{2}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{3}(\operatorname{species}_{\operatorname{Gentoo}})\ + \\ &\quad \beta_{4}(\operatorname{bill\_length\_mm}) + \epsilon \end{aligned} \]
df_gentoo <- penguins %>%
drop_na() %>%
mutate(species = relevel(factor(species), ref = "Gentoo"))
# 2. Verify the new order of levels
print(levels(df_gentoo$species)) # Should output: "Gentoo" "Adelie" "Chinstrap"
# 3. Fit the model again using the same formula
model_gentoo <- lm(body_mass_g ~ flipper_length_mm + species + bill_length_mm, data = df_gentoo)
# 4. View the new coefficients
coef(model_gentoo)[1] "Gentoo" "Adelie" "Chinstrap"- (Intercept)
- -3750.81905884646
- flipper_length_mm
- 27.5442861236828
- speciesAdelie
- -113.2541752724
- speciesChinstrap
- -845.670848479892
- bill_length_mm
- 60.1173245626618
# 1. Relevel species so Chinstrap is Level 1
df_chinstrap <- penguins %>%
drop_na() %>%
mutate(species = relevel(factor(species), ref = "Chinstrap"))
# 2. Verify the new order of levels
print(levels(df_chinstrap$species)) # Should output: "Chinstrap" "Adelie" "Gentoo"
# 3. Fit the model and view coefficients
model_chinstrap <- lm(body_mass_g ~ flipper_length_mm + species + bill_length_mm, data = df_chinstrap)
coef(model_chinstrap)[1] "Chinstrap" "Adelie" "Gentoo" - (Intercept)
- -4596.48990732635
- flipper_length_mm
- 27.5442861236828
- speciesAdelie
- 732.416673207492
- speciesGentoo
- 845.670848479891
- bill_length_mm
- 60.1173245626618
install.packages("MuMIn")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependency ‘insight’
if (!require("palmerpenguins")) install.packages("palmerpenguins")
if (!require("MuMIn")) install.packages("MuMIn")
library(tidyverse)
library(palmerpenguins)
library(MuMIn) # The gold standard package for AICc model selection
# 1. Prepare the data: Filter for Adelie and remove missing values in target/predictors
adelie_df <- penguins %>%
filter(species == "Adelie") %>%
drop_na(sex, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) %>%
mutate(sex = as.factor(sex)) # Ensure sex is binary factor
# 2. Fit the "Global Model" containing all candidate parameters
# Note: na.action = "na.fail" is strictly required for MuMIn to work
global_model <- glm(
sex ~ bill_length_mm + bill_depth_mm + flipper_length_mm + body_mass_g,
data = adelie_df,
family = binomial(link = "logit"),
na.action = "na.fail"
)
# 3. Generate Table 1: Dredge builds every possible combination of sub-models
# and ranks them automatically by AICc
model_table <- dredge(global_model, rank = "AICc")
# View the top models (equivalent to Table 1 in the paper)
print(model_table)
# 4. Generate Table 2: Model averaging and Relative Variable Importance (RWI)
# This calculates the cumulative Akaike weight for each predictor across all models
avg_model <- model.avg(model_table, subset = delta < 4) # usually delta < 2 or 4
summary(avg_model)Loading required package: palmerpenguins
Warning message in library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, :
“there is no package called ‘palmerpenguins’”
Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
Loading required package: MuMIn
Attaching package: ‘palmerpenguins’
The following objects are masked from ‘package:datasets’:
penguins, penguins_raw
Fixed term is "(Intercept)"
Global model call: glm(formula = sex ~ bill_length_mm + bill_depth_mm + flipper_length_mm +
body_mass_g, family = binomial(link = "logit"), data = adelie_df,
na.action = "na.fail")
---
Model selection table
(Int) bll_dpt_mm bll_lng_mm bdy_mss_g flp_lng_mm df logLik AICc
8 -8.509e+01 1.306 0.8404 0.007790 4 -27.012 62.3
16 -7.918e+01 1.336 0.8403 0.007861 -0.035710 5 -26.839 64.1
7 -5.873e+01 0.7132 0.008486 3 -33.665 73.5
15 -5.751e+01 0.7130 0.008523 -0.007109 4 -33.657 75.6
6 -4.598e+01 1.065 0.007184 3 -38.075 82.3
14 -4.345e+01 1.077 0.007316 -0.017070 4 -38.009 84.3
5 -2.875e+01 0.007814 2 -44.371 92.8
13 -2.857e+01 0.007821 -0.001104 3 -44.371 94.9
4 -6.461e+01 1.828 0.8007 3 -45.899 98.0
12 -7.044e+01 1.775 0.7767 0.040980 4 -45.462 99.2
10 -4.513e+01 1.576 0.085690 3 -65.056 136.3
11 -4.211e+01 0.6465 0.089920 3 -66.867 139.9
2 -2.982e+01 1.629 2 -68.454 141.0
3 -2.636e+01 0.6795 2 -70.156 144.4
9 -2.417e+01 0.127100 2 -91.279 186.6
1 7.400e-17 1 -101.199 204.4
delta weight
8 0.00 0.708
16 1.80 0.288
7 11.19 0.003
15 13.29 0.001
6 20.01 0.000
14 21.99 0.000
5 30.52 0.000
13 32.60 0.000
4 35.66 0.000
12 36.90 0.000
10 73.97 0.000
11 77.60 0.000
2 78.68 0.000
3 82.09 0.000
9 124.33 0.000
1 142.12 0.000
Models ranked by AICc(x)
Call:
model.avg(object = model_table, subset = delta < 4)
Component model call:
glm(formula = sex ~ <2 unique rhs>, family = binomial(link = "logit"),
data = adelie_df, na.action = na.fail)
Component models:
df logLik AICc delta weight
123 4 -27.01 62.31 0.0 0.71
1234 5 -26.84 64.11 1.8 0.29
Term codes:
bill_depth_mm bill_length_mm body_mass_g flipper_length_mm
1 2 3 4
Model-averaged coefficients:
(full average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) -83.380334 18.584118 18.740928 4.449 8.60e-06 ***
bill_depth_mm 1.314667 0.426794 0.430467 3.054 0.002258 **
bill_length_mm 0.840368 0.237291 0.239335 3.511 0.000446 ***
body_mass_g 0.007810 0.001964 0.001981 3.943 8.03e-05 ***
flipper_length_mm -0.010325 0.036570 0.036824 0.280 0.779176
(conditional average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) -83.380334 18.584118 18.740928 4.449 8.60e-06 ***
bill_depth_mm 1.314667 0.426794 0.430467 3.054 0.002258 **
bill_length_mm 0.840368 0.237291 0.239335 3.511 0.000446 ***
body_mass_g 0.007810 0.001964 0.001981 3.943 8.03e-05 ***
flipper_length_mm -0.035709 0.060981 0.061509 0.581 0.561539
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1install.packages("easystats")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘patchwork’, ‘bayestestR’, ‘correlation’, ‘datawizard’, ‘effectsize’, ‘modelbased’, ‘parameters’, ‘performance’, ‘report’, ‘see’
library(tidyverse)
library(palmerpenguins)
library(performance)
# 1. Prepare data for predicting sex
adelie_df <- penguins %>%
filter(species == "Adelie") %>%
drop_na(sex, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g)
# 2. Define your candidate models based on biological hypotheses
# Model A: Bill dimensions + Body Mass
model_a <- glm(sex ~ bill_length_mm + bill_depth_mm + body_mass_g,
data = adelie_df, family = binomial)
# Model B: Bill length + Body Mass (simpler)
model_b <- glm(sex ~ bill_length_mm + body_mass_g,
data = adelie_df, family = binomial)
# Model C: The "Global" model with everything
model_global <- glm(sex ~ bill_length_mm + bill_depth_mm + flipper_length_mm + body_mass_g,
data = adelie_df, family = binomial)
# 3. Compare them all at once using easystats
compare_performance(model_a, model_b, model_global, metrics = "AICc")| Name | Model | AICc | AICc_wt |
|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> |
| model_a | glm | 62.30747 | 0.708980296 |
| model_b | glm | 73.49996 | 0.002631568 |
| model_global | glm | 64.10651 | 0.288388136 |
install.packages("janitor")Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
also installing the dependency ‘snakecase’
library(performance)
# 1. Capture the performance comparison object
comp <- compare_performance(model_a, model_b, model_global, metrics = "AICc")
# 2. View it and manually calculate Delta AICc (AICc - min(AICc))
# Extract the AICc column and compute delta
comp_df <- as.data.frame(comp)
comp_df$Delta_AICc <- comp_df$AICc - min(comp_df$AICc)
# 3. Print the results with Delta AICc
print(comp_df[, c("Name", "AICc", "Delta_AICc", "AICc_wt")]) Name AICc Delta_AICc AICc_wt
1 model_a 62.30747 0.000000 0.708980296
2 model_b 73.49996 11.192496 0.002631568
3 model_global 64.10651 1.799041 0.288388136library(tidyverse)
library(palmerpenguins)
library(janitor)
library(MuMIn) # Required for automated multi-model comparison
# 1. Prepare and clean the raw Adelie data
adelie_raw <- penguins_raw %>%
clean_names() %>%
filter(str_detect(species, "Adelie")) %>%
select(sex, culmen_length_mm, culmen_depth_mm, flipper_length_mm, body_mass_g) %>%
drop_na() %>%
mutate(sex = as.factor(tolower(sex)))
# 2. Fit the global model (na.action = "na.fail" is mandatory for MuMIn)
global_model <- glm(
sex ~ culmen_length_mm + culmen_depth_mm + flipper_length_mm + body_mass_g,
data = adelie_raw,
family = binomial,
na.action = "na.fail"
)
# 3. Generate all model combinations ranked by AICc
# This automatically anchors the top model's Delta AICc at 0.00
table_results <- dredge(global_model, rank = "AICc")
# View the top models
print(table_results)Fixed term is "(Intercept)"
Global model call: glm(formula = sex ~ culmen_length_mm + culmen_depth_mm + flipper_length_mm +
body_mass_g, family = binomial, data = adelie_raw, na.action = "na.fail")
---
Model selection table
(Int) bdy_mss_g clm_dpt_mm clm_lng_mm flp_lng_mm df logLik AICc
8 -8.509e+01 0.007790 1.306 0.8404 4 -27.012 62.3
16 -7.918e+01 0.007861 1.336 0.8403 -0.035710 5 -26.839 64.1
6 -5.873e+01 0.008486 0.7132 3 -33.665 73.5
14 -5.751e+01 0.008523 0.7130 -0.007109 4 -33.657 75.6
4 -4.598e+01 0.007184 1.065 3 -38.075 82.3
12 -4.345e+01 0.007316 1.077 -0.017070 4 -38.009 84.3
2 -2.875e+01 0.007814 2 -44.371 92.8
10 -2.857e+01 0.007821 -0.001104 3 -44.371 94.9
7 -6.461e+01 1.828 0.8007 3 -45.899 98.0
15 -7.044e+01 1.775 0.7767 0.040980 4 -45.462 99.2
11 -4.513e+01 1.576 0.085690 3 -65.056 136.3
13 -4.211e+01 0.6465 0.089920 3 -66.867 139.9
3 -2.982e+01 1.629 2 -68.454 141.0
5 -2.636e+01 0.6795 2 -70.156 144.4
9 -2.417e+01 0.127100 2 -91.279 186.6
1 7.400e-17 1 -101.199 204.4
delta weight
8 0.00 0.708
16 1.80 0.288
6 11.19 0.003
14 13.29 0.001
4 20.01 0.000
12 21.99 0.000
2 30.52 0.000
10 32.60 0.000
7 35.66 0.000
15 36.90 0.000
11 73.97 0.000
13 77.60 0.000
3 78.68 0.000
5 82.09 0.000
9 124.33 0.000
1 142.12 0.000
Models ranked by AICc(x) library(tidyverse)
library(palmerpenguins)
library(janitor)
library(performance)
# 1. Filter penguins_raw for Adelie and the study's exact seasons
adelie_study <- penguins_raw %>%
clean_names() %>%
filter(str_detect(species, "Adelie"),
study_name %in% c("PAL0708", "PAL0809", "PAL0910")) %>%
select(sex, culmen_length_mm, culmen_depth_mm, flipper_length_mm, body_mass_g) %>%
drop_na() %>%
mutate(sex = as.factor(tolower(sex)))
# 2. Define the three specific candidate models
model_a <- glm(sex ~ culmen_length_mm + culmen_depth_mm + body_mass_g,
data = adelie_study, family = binomial)
model_b <- glm(sex ~ culmen_length_mm + body_mass_g,
data = adelie_study, family = binomial)
model_global <- glm(sex ~ culmen_length_mm + culmen_depth_mm + flipper_length_mm + body_mass_g,
data = adelie_study, family = binomial)
# 3. Compare performance using easystats and calculate Delta AICc
comp <- compare_performance(model_a, model_b, model_global, metrics = "AICc")
comp_df <- as.data.frame(comp)
comp_df$Delta_AICc <- comp_df$AICc - min(comp_df$AICc)
# Print the final clean comparison table
print(comp_df[, c("Name", "AICc", "Delta_AICc", "AICc_wt")]) Name AICc Delta_AICc AICc_wt
1 model_a 62.30747 0.000000 0.708980296
2 model_b 73.49996 11.192496 0.002631568
3 model_global 64.10651 1.799041 0.288388136library(tidyverse)
library(palmerpenguins)
library(performance)
# 1. Prepare data
adelie_df <- penguins %>%
filter(species == "Adelie") %>%
select(sex, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) %>%
drop_na() %>%
mutate(sex = as.factor(sex))
# 2. Fit models
model_a <- glm(sex ~ bill_length_mm + bill_depth_mm + body_mass_g,
data = adelie_df, family = binomial)
model_b <- glm(sex ~ bill_length_mm + body_mass_g,
data = adelie_df, family = binomial)
model_global <- glm(sex ~ bill_length_mm + bill_depth_mm + flipper_length_mm + body_mass_g,
data = adelie_df, family = binomial)
# 3. Get AICc and weights via easystats performance
comp <- compare_performance(model_a, model_b, model_global, metrics = "AICc")
summary_df <- as.data.frame(comp) %>%
mutate(Delta_AICc = AICc - min(AICc)) %>%
select(Name, AICc, Delta_AICc, AICc_wt)
# 4. Robust metric extractor using base R model properties
extract_metrics <- function(mod, data) {
# McFadden's R2 via deviance ratio (bulletproof for glm binomial)
r2_mf <- 1 - (mod$deviance / mod$null.deviance)
# Classification accuracy using a 0.5 probability threshold
preds <- predict(mod, type = "response")
pred_class <- ifelse(preds > 0.5, levels(data$sex)[2], levels(data$sex)[1])
accuracy <- mean(pred_class == data$sex) * 100
tibble(
r2_mf = r2_mf,
percent_correct = accuracy
)
}
# 5. Bind metrics cleanly by model name
models_list <- list(model_a = model_a, model_b = model_b, model_global = model_global)
extra_metrics_df <- map_dfr(models_list, extract_metrics, data = adelie_df, .id = "Model_Key")
# 6. Merge into final clean summary table
final_table <- summary_df %>%
mutate(Model_Key = c("model_a", "model_b", "model_global")) %>%
left_join(extra_metrics_df, by = "Model_Key") %>%
select(-Model_Key)
print(final_table) Name AICc Delta_AICc AICc_wt r2_mf percent_correct
1 model_a 62.30747 0.000000 0.708980296 0.7330827 91.78082
2 model_b 73.49996 11.192496 0.002631568 0.6673355 89.72603
3 model_global 64.10651 1.799041 0.288388136 0.7347915 93.15068