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 + special_status +
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: 3220.826summary(final_model)##
## Call:
## lm(formula = cand_formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5434 -0.4987 0.0593 0.6100 9.7873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.54322 0.24120 18.836 < 2e-16 ***
## log2(height_m) 1.86438 0.07675 24.292 < 2e-16 ***
## height_m -0.42692 0.06363 -6.709 3.28e-11 ***
## gender_catMostly female -0.25145 0.16052 -1.566 0.117560
## gender_catMostly male 0.15309 0.10425 1.468 0.142305
## gender_catGenderless 0.37678 0.15529 2.426 0.015428 *
## I(hp/10) 0.11848 0.01805 6.565 8.39e-11 ***
## I(defense/10) 0.06845 0.01624 4.214 2.74e-05 ***
## I(sp_attack/10) -0.04549 0.01573 -2.891 0.003919 **
## I(speed/10) -0.03462 0.01531 -2.261 0.023966 *
## happiness_catHappy -0.37912 0.22732 -1.668 0.095674 .
## happiness_catUnhappy 0.25039 0.13138 1.906 0.056960 .
## special_statussub-legendary -0.37872 0.21211 -1.785 0.074494 .
## special_statuslegendary 0.09786 0.26297 0.372 0.709888
## special_statusmythical -0.59544 0.29641 -2.009 0.044825 *
## is_rock 0.54828 0.31453 1.743 0.081611 .
## is_ice 0.46672 0.16784 2.781 0.005527 **
## is_steel 0.81677 0.16128 5.064 4.88e-07 ***
## is_ground 0.40066 0.14365 2.789 0.005387 **
## is_fire 0.22605 0.14011 1.613 0.106972
## is_grass -0.39055 0.11434 -3.416 0.000662 ***
## is_flying -0.39868 0.12177 -3.274 0.001097 **
## is_fairy -0.62814 0.15963 -3.935 8.90e-05 ***
## is_ghost -2.12191 0.46886 -4.526 6.75e-06 ***
## is_poison -0.88010 0.19400 -4.537 6.41e-06 ***
## log2(height_m):is_rock -0.37506 0.26858 -1.396 0.162882
## log2(height_m):is_ghost -0.52935 0.26593 -1.991 0.046802 *
## log2(height_m):is_poison -0.34475 0.18137 -1.901 0.057611 .
## height_m:is_rock 0.05965 0.22157 0.269 0.787800
## height_m:is_ghost 0.75234 0.33709 2.232 0.025846 *
## height_m:is_poison 0.31191 0.10161 3.070 0.002202 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.146 on 994 degrees of freedom
## Multiple R-squared: 0.7637, Adjusted R-squared: 0.7566
## F-statistic: 107.1 on 30 and 994 DF, p-value: < 2.2e-16res$log # shows the elimination path## iter type dropped old_AIC
## 1 1 psychic height_m:is_psychic, log2(height_m):is_psychic 3269.676
## 2 2 dark height_m:is_dark, log2(height_m):is_dark 3265.832
## 3 3 dragon height_m:is_dragon, log2(height_m):is_dragon 3262.110
## 4 4 normal height_m:is_normal, log2(height_m):is_normal 3258.506
## 5 5 ground height_m:is_ground, log2(height_m):is_ground 3254.810
## 6 6 fighting height_m:is_fighting, log2(height_m):is_fighting 3251.738
## 7 7 fairy height_m:is_fairy, log2(height_m):is_fairy 3248.676
## 8 8 fire height_m:is_fire, log2(height_m):is_fire 3245.802
## 9 9 steel height_m:is_steel, log2(height_m):is_steel 3243.768
## 10 10 ice height_m:is_ice, log2(height_m):is_ice 3241.832
## 11 11 psychic is_psychic 3239.830
## 12 12 bug height_m:is_bug, log2(height_m):is_bug 3237.899
## 13 13 bug is_bug 3236.010
## 14 14 normal is_normal 3234.055
## 15 15 dragon is_dragon 3232.461
## 16 16 fighting is_fighting 3231.400
## 17 17 dark is_dark 3230.207
## 18 18 electric height_m:is_electric, log2(height_m):is_electric 3228.957
## 19 19 electric is_electric 3227.641
## 20 20 flying height_m:is_flying, log2(height_m):is_flying 3226.147
## 21 21 water height_m:is_water, log2(height_m):is_water 3224.981
## 22 22 grass height_m:is_grass, log2(height_m):is_grass 3223.398
## 23 23 water is_water 3222.066
## new_AIC dAIC
## 1 3265.832 -3.844584
## 2 3262.110 -3.721907
## 3 3258.506 -3.604121
## 4 3254.810 -3.695140
## 5 3251.738 -3.072149
## 6 3248.676 -3.061840
## 7 3245.802 -2.874496
## 8 3243.768 -2.033720
## 9 3241.832 -1.936009
## 10 3239.830 -2.002363
## 11 3237.899 -1.930950
## 12 3236.010 -1.888966
## 13 3234.055 -1.954656
## 14 3232.461 -1.594553
## 15 3231.400 -1.060380
## 16 3230.207 -1.192989
## 17 3228.957 -1.249958
## 18 3227.641 -1.316759
## 19 3226.147 -1.494029
## 20 3224.981 -1.165555
## 21 3223.398 -1.583144
## 22 3222.066 -1.331807
## 23 3220.826 -1.240028summary(final_model) %>% tab_model()| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.54 | 4.07 – 5.02 | <0.001 |
| log2(height_m) | 1.86 | 1.71 – 2.01 | <0.001 |
| height_m | -0.43 | -0.55 – -0.30 | <0.001 |
| gender_catMostly female | -0.25 | -0.57 – 0.06 | 0.118 |
| gender_catMostly male | 0.15 | -0.05 – 0.36 | 0.142 |
| gender_catGenderless | 0.38 | 0.07 – 0.68 | 0.015 |
| I(hp/10) | 0.12 | 0.08 – 0.15 | <0.001 |
| I(defense/10) | 0.07 | 0.04 – 0.10 | <0.001 |
| I(sp_attack/10) | -0.05 | -0.08 – -0.01 | 0.004 |
| I(speed/10) | -0.03 | -0.06 – -0.00 | 0.024 |
| happiness_catHappy | -0.38 | -0.83 – 0.07 | 0.096 |
| happiness_catUnhappy | 0.25 | -0.01 – 0.51 | 0.057 |
| special_statussub-legendary | -0.38 | -0.79 – 0.04 | 0.074 |
| special_statuslegendary | 0.10 | -0.42 – 0.61 | 0.710 |
| special_statusmythical | -0.60 | -1.18 – -0.01 | 0.045 |
| is_rock | 0.55 | -0.07 – 1.17 | 0.082 |
| is_ice | 0.47 | 0.14 – 0.80 | 0.006 |
| is_steel | 0.82 | 0.50 – 1.13 | <0.001 |
| is_ground | 0.40 | 0.12 – 0.68 | 0.005 |
| is_fire | 0.23 | -0.05 – 0.50 | 0.107 |
| is_grass | -0.39 | -0.61 – -0.17 | 0.001 |
| is_flying | -0.40 | -0.64 – -0.16 | 0.001 |
| is_fairy | -0.63 | -0.94 – -0.31 | <0.001 |
| is_ghost | -2.12 | -3.04 – -1.20 | <0.001 |
| is_poison | -0.88 | -1.26 – -0.50 | <0.001 |
| log2(height_m):is_rock | -0.38 | -0.90 – 0.15 | 0.163 |
| log2(height_m):is_ghost | -0.53 | -1.05 – -0.01 | 0.047 |
| log2(height_m):is_poison | -0.34 | -0.70 – 0.01 | 0.058 |
| height_m:is_rock | 0.06 | -0.38 – 0.49 | 0.788 |
| height_m:is_ghost | 0.75 | 0.09 – 1.41 | 0.026 |
| height_m:is_poison | 0.31 | 0.11 – 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]
ref_values$special_status <- levels(Pokemon_clean$special_status)[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).