The weight of Pokémon species - Regression analyses considering type
Lara Lusa, Sam Manski, Marianne Huebner
July 2026
Pokemon_clean$gender_cat <- factor(Pokemon_clean$gender_cat )
Pokemon_clean$gender_cat <-relevel(Pokemon_clean$gender_cat, ref = "Equal")
Pokemon_clean$happiness_cat <-relevel(factor(Pokemon_clean$happiness_cat), ref = "Normal")
mv_fit_step <- lm(formula = log2(weight_kg) ~ log2(height_m) + height_m + gender_cat +
I(hp/10) + I(defense/10) + I(sp_attack/10) + I(speed/10) +
special_status + happiness_cat, data = Pokemon_clean, x = T,
y = T)Introduction
The aim of this analysis project is to study the variability in the weight of Pokémon species. The statistical analysis plan (SAP) specifies that linear regression models will be used to estimate the association between height and weight, both univariable and multivariable. The document Regression_Report addresses the problem with univariable and multivariable analyses.
In this document we assess also the association of type with the outcome. We assume that all the models described in the document Regression_Report were fitted and analyzed.
Pokémon types
We start from the multivariable model selected with the backward AIC based method and we add the 18 Pokémon types as explanatory variables. Each type is defined as a binary variable, indicating if the species belongs or does not belong to the type. (See additional IDA report for more details about this variable).
mv_fit_type <- update(mv_fit_step, .~. + is_rock + is_dragon + is_bug + is_steel + is_ground+ is_water+is_fire+is_grass+is_ice+is_flying+is_fairy+is_ghost+is_poison+is_electric+is_normal+is_dark+is_psychic+is_fighting)
summary(mv_fit_type) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.31 | 3.81 – 4.81 | <0.001 |
| height m [log2] | 1.73 | 1.60 – 1.86 | <0.001 |
| height m | -0.31 | -0.40 – -0.21 | <0.001 |
|
gender cat [Mostly female] |
-0.24 | -0.56 – 0.08 | 0.140 |
| gender cat [Mostly male] | 0.17 | -0.05 – 0.38 | 0.129 |
| gender cat [Genderless] | 0.45 | 0.14 – 0.75 | 0.004 |
| hp/10 | 0.12 | 0.08 – 0.15 | <0.001 |
| defense/10 | 0.07 | 0.04 – 0.10 | <0.001 |
| sp attack/10 | -0.04 | -0.07 – -0.01 | 0.019 |
| speed/10 | -0.04 | -0.07 – -0.01 | 0.021 |
|
special status [sub-legendary] |
-0.44 | -0.86 – -0.02 | 0.042 |
|
special status [legendary] |
0.17 | -0.36 – 0.69 | 0.532 |
| special status [mythical] | -0.62 | -1.21 – -0.03 | 0.041 |
| happiness cat [Happy] | -0.38 | -0.84 – 0.08 | 0.103 |
| happiness cat [Unhappy] | 0.17 | -0.11 – 0.44 | 0.232 |
| is rock | 0.67 | 0.35 – 0.99 | <0.001 |
| is dragon | 0.13 | -0.19 – 0.46 | 0.430 |
| is bug | -0.04 | -0.32 – 0.25 | 0.799 |
| is steel | 0.81 | 0.47 – 1.14 | <0.001 |
| is ground | 0.40 | 0.10 – 0.71 | 0.010 |
| is water | -0.04 | -0.30 – 0.22 | 0.771 |
| is fire | 0.24 | -0.07 – 0.56 | 0.131 |
| is grass | -0.37 | -0.63 – -0.10 | 0.007 |
| is ice | 0.51 | 0.16 – 0.87 | 0.005 |
| is flying | -0.38 | -0.63 – -0.13 | 0.003 |
| is fairy | -0.62 | -0.96 – -0.28 | <0.001 |
| is ghost | -1.11 | -1.44 – -0.78 | <0.001 |
| is poison | -0.43 | -0.72 – -0.14 | 0.004 |
| is electric | 0.11 | -0.24 – 0.46 | 0.541 |
| is normal | 0.05 | -0.24 – 0.34 | 0.719 |
| is dark | 0.17 | -0.14 – 0.48 | 0.292 |
| is psychic | -0.04 | -0.33 – 0.26 | 0.806 |
| is fighting | 0.15 | -0.16 – 0.47 | 0.342 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.760 / 0.752 | ||
RESULTS: Including type increases the R^2 (from 0.71 to 0.75) and decreases the residual standard error. Many types have a statistically significant association with weight. Rock, ice, steel and ground are statistically significantly positively associated, grass, flying, fairy, ghost and poison are statistically significantly negatively associated. The estimated coefficients of other explanatory variables change compared to the model without type, and the associations are generally attenuated (as might be expected, as type is associated to the other variables).
We fit also a model where the interaction between height (linear + log2) and all types are considered. The joint distribution of type and height was examined during IDA (see additional IDA report); no type was found to be so sparse or so narrowly distributed in height as to make its interaction term inestimable, though estimates for less numerous types are expected to be less precise. To limit the number of terms, backward elimination in the interaction model is performed based on AIC. The backward elimination drops both interactions together for a given type (pairwise), or none. The procedure only allows dropping a main effect after both interactions are eliminated (hierarchical rule).
# first fit a model with all interactions between types and height
mv_fit_type_int <- lm(
log2(weight_kg) ~ log2(height_m) + height_m + gender_cat +
I(hp/10) + I(attack/10) + I(defense/10) +
I(sp_attack/10) + I(sp_defense/10) + I(speed/10) +
happiness_cat +
is_rock + is_ice + is_dragon + is_bug + is_steel + is_ground +
is_water + is_fire + is_grass + is_flying + is_fairy + is_ghost +
is_poison + is_electric + is_normal + is_dark + is_psychic + is_fighting +
# Interactions
is_rock*log2(height_m) + is_ice*log2(height_m) + is_dragon*log2(height_m) +
is_bug*log2(height_m) + is_steel*log2(height_m) + is_ground*log2(height_m) +
is_water*log2(height_m) + is_fire*log2(height_m) + is_grass*log2(height_m) +
is_flying*log2(height_m) + is_fairy*log2(height_m) + is_ghost*log2(height_m) +
is_poison*log2(height_m) + is_electric*log2(height_m) + is_normal*log2(height_m) +
is_dark*log2(height_m) + is_psychic*log2(height_m) + is_fighting*log2(height_m) +
is_rock*height_m + is_ice*height_m + is_dragon*height_m +
is_bug*height_m + is_steel*height_m + is_ground*height_m +
is_water*height_m + is_fire*height_m + is_grass*height_m +
is_flying*height_m + is_fairy*height_m + is_ghost*height_m +
is_poison*height_m + is_electric*height_m + is_normal*height_m +
is_dark*height_m + is_psychic*height_m + is_fighting*height_m,
data = Pokemon_clean
)## -------------------------------
## Blocked, hierarchical backward AIC for type–height interactions
## -------------------------------
# types you use as 0/1 indicators in the model:
type_names <- c(
"rock","ice","dragon","bug","steel","ground","water","fire","grass",
"flying","fairy","ghost","poison","electric","normal","dark","psychic","fighting"
)
# helper: given term labels from terms(model), detect presence of main & interactions (either order)
present_terms_for_type <- function(terms_in_model, type) {
main <- paste0("is_", type)
int_lin_a <- paste0("is_", type, ":height_m")
int_lin_b <- paste0("height_m:is_", type)
int_log_a <- paste0("is_", type, ":log2(height_m)")
int_log_b <- paste0("log2(height_m):is_", type)
list(
main_present = main %in% terms_in_model,
main_label = main,
lin_present = any(c(int_lin_a, int_lin_b) %in% terms_in_model),
lin_label = intersect(c(int_lin_a, int_lin_b), terms_in_model),
log2_present = any(c(int_log_a, int_log_b) %in% terms_in_model),
log2_label = intersect(c(int_log_a, int_log_b), terms_in_model)
)
}
# helper: update formula dropping a vector of exact term labels
drop_terms_from_formula <- function(fml, drops) {
if (length(drops) == 0) return(fml)
rhs <- paste(". ~ . -", paste(drops, collapse = " - "))
update(fml, rhs)
}
# main routine
step_blocked_types <- function(model, data, type_names,
strict_pairs = TRUE, # if TRUE, only drop both interactions together
verbose = TRUE) {
current_model <- model
current_aic <- AIC(current_model)
log_df <- data.frame(iter = integer(), type = character(),
dropped = character(), old_AIC = numeric(),
new_AIC = numeric(), dAIC = numeric(),
stringsAsFactors = FALSE)
improved <- TRUE
iter <- 0
while (improved) {
iter <- iter + 1
improved <- FALSE
best_model <- current_model
best_aic <- current_aic
best_drop <- NULL
best_type <- NULL
terms_in_model <- attr(terms(current_model), "term.labels")
for (tp in type_names) {
pres <- present_terms_for_type(terms_in_model, tp)
# Which drops are allowed under hierarchy?
# Case A: both interactions present -> propose dropping BOTH interactions together
if (pres$lin_present && pres$log2_present) {
drop_vec <- c(pres$lin_label, pres$log2_label)
drop_vec <- unique(as.character(drop_vec))
# Try dropping both interactions
cand_formula <- drop_terms_from_formula(formula(current_model), drop_vec)
cand_model <- lm(cand_formula, data = data)
cand_aic <- AIC(cand_model)
if (cand_aic < best_aic) {
best_model <- cand_model
best_aic <- cand_aic
best_drop <- drop_vec
best_type <- tp
improved <- TRUE
}
} else {
# Case B (optional): only one interaction present
if (!strict_pairs) {
# allow dropping whichever single interaction remains
single <- unique(as.character(c(pres$lin_label, pres$log2_label)))
single <- single[nchar(single) > 0]
if (length(single) == 1) {
cand_formula <- drop_terms_from_formula(formula(current_model), single)
cand_model <- lm(cand_formula, data = data)
cand_aic <- AIC(cand_model)
if (cand_aic < best_aic) {
best_model <- cand_model
best_aic <- cand_aic
best_drop <- single
best_type <- tp
improved <- TRUE
}
}
}
}
# Case C: interactions absent, main present -> main is eligible for drop
if (!pres$lin_present && !pres$log2_present && pres$main_present) {
drop_vec <- pres$main_label
cand_formula <- drop_terms_from_formula(formula(current_model), drop_vec)
cand_model <- lm(cand_formula, data = data)
cand_aic <- AIC(cand_model)
if (cand_aic < best_aic) {
best_model <- cand_model
best_aic <- cand_aic
best_drop <- drop_vec
best_type <- tp
improved <- TRUE
}
}
# Case D: interactions present but main absent -> OK, we already handle dropping interactions above.
# We never try to drop a main that isn't in the model.
} # end for tp
if (improved) {
if (verbose) {
cat(sprintf("Iter %d: drop (%s): %s ΔAIC = %.3f\n",
iter, best_type, paste(best_drop, collapse=", "),
best_aic - current_aic))
}
log_df <- rbind(log_df, data.frame(
iter = iter,
type = best_type,
dropped = paste(best_drop, collapse = ", "),
old_AIC = current_aic,
new_AIC = best_aic,
dAIC = best_aic - current_aic,
stringsAsFactors = FALSE
))
current_model <- best_model
current_aic <- best_aic
}
} # while
list(model = current_model, log = log_df)
}
# mv_fit_type_int (defined above) is the full model with all type×height interactions
# Example call (strict pair-drop on interactions):
res <- step_blocked_types(mv_fit_type_int, data = Pokemon_clean,
type_names = type_names,
strict_pairs = TRUE, # only drop both interactions together
verbose = FALSE)
final_model <- res$model
cat("\nFinal AIC:", AIC(final_model), "\n")##
## Final AIC: 3219.621summary(final_model)##
## Call:
## lm(formula = cand_formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4102 -0.5164 0.0581 0.6297 9.8383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.47045 0.24599 18.174 < 2e-16 ***
## log2(height_m) 1.83904 0.07838 23.462 < 2e-16 ***
## height_m -0.40111 0.06285 -6.382 2.68e-10 ***
## gender_catMostly female -0.26470 0.16188 -1.635 0.102335
## gender_catMostly male 0.13606 0.10412 1.307 0.191580
## gender_catGenderless 0.20753 0.13778 1.506 0.132302
## I(hp/10) 0.11007 0.01863 5.907 4.78e-09 ***
## I(attack/10) 0.01199 0.01755 0.683 0.494688
## I(defense/10) 0.04348 0.01866 2.330 0.020007 *
## I(sp_attack/10) -0.05899 0.01632 -3.614 0.000317 ***
## I(sp_defense/10) 0.04788 0.01908 2.509 0.012249 *
## I(speed/10) -0.04052 0.01597 -2.537 0.011337 *
## happiness_catHappy -0.62109 0.21441 -2.897 0.003853 **
## happiness_catUnhappy 0.21474 0.12564 1.709 0.087736 .
## is_rock 0.55701 0.31564 1.765 0.077926 .
## is_ice 0.45692 0.16773 2.724 0.006559 **
## is_steel 0.84240 0.16087 5.237 1.99e-07 ***
## is_ground 0.43987 0.14422 3.050 0.002350 **
## is_fire 0.22803 0.14014 1.627 0.104027
## is_grass -0.41087 0.11450 -3.588 0.000349 ***
## is_flying -0.40552 0.12161 -3.334 0.000886 ***
## is_fairy -0.68936 0.16127 -4.275 2.10e-05 ***
## is_ghost -2.25585 0.46530 -4.848 1.45e-06 ***
## is_poison -0.90020 0.19313 -4.661 3.57e-06 ***
## log2(height_m):is_rock -0.38153 0.26927 -1.417 0.156815
## log2(height_m):is_ghost -0.59590 0.26467 -2.251 0.024572 *
## log2(height_m):is_poison -0.34890 0.18125 -1.925 0.054525 .
## height_m:is_rock 0.06776 0.22237 0.305 0.760650
## height_m:is_ghost 0.84672 0.33430 2.533 0.011467 *
## height_m:is_poison 0.31471 0.10149 3.101 0.001985 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.146 on 995 degrees of freedom
## Multiple R-squared: 0.7635, Adjusted R-squared: 0.7566
## F-statistic: 110.8 on 29 and 995 DF, p-value: < 2.2e-16res$log # shows the elimination path## iter type dropped old_AIC
## 1 1 normal height_m:is_normal, log2(height_m):is_normal 3271.235
## 2 2 fighting height_m:is_fighting, log2(height_m):is_fighting 3267.327
## 3 3 ground height_m:is_ground, log2(height_m):is_ground 3263.581
## 4 4 fairy height_m:is_fairy, log2(height_m):is_fairy 3259.806
## 5 5 dark height_m:is_dark, log2(height_m):is_dark 3256.119
## 6 6 psychic height_m:is_psychic, log2(height_m):is_psychic 3252.718
## 7 7 dragon height_m:is_dragon, log2(height_m):is_dragon 3249.191
## 8 8 steel height_m:is_steel, log2(height_m):is_steel 3246.744
## 9 9 fire height_m:is_fire, log2(height_m):is_fire 3243.830
## 10 10 ice height_m:is_ice, log2(height_m):is_ice 3241.569
## 11 11 psychic is_psychic 3239.289
## 12 12 bug height_m:is_bug, log2(height_m):is_bug 3237.303
## 13 13 bug is_bug 3235.536
## 14 14 flying height_m:is_flying, log2(height_m):is_flying 3233.669
## 15 15 water height_m:is_water, log2(height_m):is_water 3232.105
## 16 16 water is_water 3229.614
## 17 17 normal is_normal 3227.642
## 18 18 grass height_m:is_grass, log2(height_m):is_grass 3226.007
## 19 19 fighting is_fighting 3224.724
## 20 20 dragon is_dragon 3223.596
## 21 21 dark is_dark 3222.600
## 22 22 electric height_m:is_electric, log2(height_m):is_electric 3221.519
## 23 23 electric is_electric 3220.426
## new_AIC dAIC
## 1 3267.327 -3.9078560
## 2 3263.581 -3.7463262
## 3 3259.806 -3.7747486
## 4 3256.119 -3.6873993
## 5 3252.718 -3.4008783
## 6 3249.191 -3.5272207
## 7 3246.744 -2.4467203
## 8 3243.830 -2.9145009
## 9 3241.569 -2.2606623
## 10 3239.289 -2.2798241
## 11 3237.303 -1.9858642
## 12 3235.536 -1.7675395
## 13 3233.669 -1.8669335
## 14 3232.105 -1.5641473
## 15 3229.614 -2.4908816
## 16 3227.642 -1.9722943
## 17 3226.007 -1.6348519
## 18 3224.724 -1.2826499
## 19 3223.596 -1.1280716
## 20 3222.600 -0.9961329
## 21 3221.519 -1.0813083
## 22 3220.426 -1.0920398
## 23 3219.621 -0.8053959summary(final_model) %>% tab_model()| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.47 | 3.99 – 4.95 | <0.001 |
| log2(height_m) | 1.84 | 1.69 – 1.99 | <0.001 |
| height_m | -0.40 | -0.52 – -0.28 | <0.001 |
| gender_catMostly female | -0.26 | -0.58 – 0.05 | 0.102 |
| gender_catMostly male | 0.14 | -0.07 – 0.34 | 0.192 |
| gender_catGenderless | 0.21 | -0.06 – 0.48 | 0.132 |
| I(hp/10) | 0.11 | 0.07 – 0.15 | <0.001 |
| I(attack/10) | 0.01 | -0.02 – 0.05 | 0.495 |
| I(defense/10) | 0.04 | 0.01 – 0.08 | 0.020 |
| I(sp_attack/10) | -0.06 | -0.09 – -0.03 | <0.001 |
| I(sp_defense/10) | 0.05 | 0.01 – 0.09 | 0.012 |
| I(speed/10) | -0.04 | -0.07 – -0.01 | 0.011 |
| happiness_catHappy | -0.62 | -1.04 – -0.20 | 0.004 |
| happiness_catUnhappy | 0.21 | -0.03 – 0.46 | 0.088 |
| is_rock | 0.56 | -0.06 – 1.18 | 0.078 |
| is_ice | 0.46 | 0.13 – 0.79 | 0.007 |
| is_steel | 0.84 | 0.53 – 1.16 | <0.001 |
| is_ground | 0.44 | 0.16 – 0.72 | 0.002 |
| is_fire | 0.23 | -0.05 – 0.50 | 0.104 |
| is_grass | -0.41 | -0.64 – -0.19 | <0.001 |
| is_flying | -0.41 | -0.64 – -0.17 | 0.001 |
| is_fairy | -0.69 | -1.01 – -0.37 | <0.001 |
| is_ghost | -2.26 | -3.17 – -1.34 | <0.001 |
| is_poison | -0.90 | -1.28 – -0.52 | <0.001 |
| log2(height_m):is_rock | -0.38 | -0.91 – 0.15 | 0.157 |
| log2(height_m):is_ghost | -0.60 | -1.12 – -0.08 | 0.025 |
| log2(height_m):is_poison | -0.35 | -0.70 – 0.01 | 0.055 |
| height_m:is_rock | 0.07 | -0.37 – 0.50 | 0.761 |
| height_m:is_ghost | 0.85 | 0.19 – 1.50 | 0.011 |
| height_m:is_poison | 0.31 | 0.12 – 0.51 | 0.002 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.764 / 0.757 | ||
RESULTS: The interactions retained in the model after variable selection are those for rock, ghost and poison. The main effects are those of rock, ice, steel, ground, fire, grass, flying, fairy, ghost and poison. The R2 is similar to that of the model without interactions.
The graphical presentation of the estimated association between height and weight among different types is shown below. The types that were omitted from variable selection are grouped in the Other class, the weight is estimated using the reference value for categorical variables and the median value for numerical variables.
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(scales)
library(ggrepel)
# --- 1. Define all types and those kept in final_model ---
all_types <- c("is_rock","is_ice","is_dragon","is_bug","is_steel","is_ground",
"is_water","is_fire","is_grass","is_flying","is_fairy","is_ghost",
"is_poison","is_electric","is_normal","is_dark","is_psychic","is_fighting")
kept_types <- intersect(all_types, attr(terms(final_model), "term.labels"))
# --- 2. Reference values for other predictors ---
ref_values <- Pokemon_clean %>%
summarise(across(where(is.numeric), median, na.rm = TRUE))
ref_values$gender_cat <- levels(Pokemon_clean$gender_cat)[1]
ref_values$happiness_cat <- levels(Pokemon_clean$happiness_cat)[1]
ref_values$PokemonPhase_5cat <- levels(Pokemon_clean$PokemonPhase_5cat)[1]
# --- 3. Build predictions per type within observed ranges ---
pred_list <- list()
for (t in all_types) {
obs_range <- Pokemon_clean %>%
filter(.data[[t]] == 1) %>%
summarise(min_h = min(height_m, na.rm = TRUE),
max_h = max(height_m, na.rm = TRUE))
if (nrow(obs_range) == 0 || is.na(obs_range$min_h)) next
height_seq <- seq(obs_range$min_h, obs_range$max_h, length.out = 100)
newdata <- data.frame(
height_m = height_seq,
log2_height_m = log2(height_seq)
)
for (col in setdiff(names(ref_values), c("height_m","log2_height_m"))) {
newdata[[col]] <- ref_values[[col]]
}
for (tt in all_types) {
newdata[[tt]] <- ifelse(tt == t, 1, 0)
}
newdata$type <- ifelse(t %in% kept_types, t, "Other")
newdata$pred <- predict(final_model, newdata)
newdata$pred_kg <- 2^(newdata$pred)
newdata$line_type <- ifelse(newdata$type == "Other", "dashed", "solid")
pred_list[[t]] <- newdata
}
pred_all <- bind_rows(pred_list)
# --- 4. Assign Dark2 colors ---
main_types <- unique(pred_all$type)
kept_colors <- brewer.pal(min(length(kept_types), 8), "Dark2") # max 8 colors in Dark2
# If more than 8 kept types, repeat colors
kept_colors <- rep(kept_colors, length.out = length(kept_types))
colors_vec <- c(kept_colors, "gray50") # last color for Other
names(colors_vec) <- c(kept_types, "Other")
# --- 5. Prepare labels at endpoints ---
end_points <- pred_all %>%
group_by(type) %>%
filter(height_m == max(height_m)) %>%
ungroup()
# --- 6. Plot ---
ggplot(pred_all, aes(x = height_m, y = pred_kg, color = type, linetype = line_type)) +
geom_line(size = 1.2) +
geom_text_repel(data = end_points,
aes(label = type),
nudge_x = 0.1,
segment.color = NA,
direction = "y",
hjust = 0,
size = 4,
box.padding = 0.2,
max.overlaps = 20) +
scale_x_continuous(trans = "log2",
breaks = c(0.25, 0.5, 1, 2, 4, 8),
labels = c("0.25","0.5","1","2","4","8")) +
scale_y_continuous(trans = "log2",
breaks = c(1, 2, 5, 10, 20, 50, 100, 200, 500),
labels = c("1","2","5","10","20","50","100","200","500")) +
scale_color_manual(values = colors_vec) +
scale_linetype_identity() +
labs(x = "Height (m)",
y = "Predicted Weight (kg)") +
theme_minimal(base_size = 14) +
theme(legend.position = "none")
RESULTS: The differences between types are substantial. Ghost type has a slightly different shape than the other types. The growth for some types is very steep (steel, ice, ground, fire, rock); other types as fairy, grass and flying increase weight at lower rate, as height increases. The distinctive shape of the poison type is determined by the extreme values of height and weight of Eternatus (poison and dragon type).