The weight of Pokémon species - Regression analyses
Lara Lusa, Sam Manski, Marianne Huebner
July 2026
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.
Before conducting the regression analysis, an initial data analysis report was created. Initial Data Analysis (IDA) is a systematic process to provide reliable knowledge about the data to determine the suitability of the data for the main data analysis. Information from the IDA findings was used to updated the statistical analysis plan (SAP v1 to SAP v2). Analyses shown in this document are from SAP v2 unless specifically indicated for comparison purposes.
COMMENT:
In this document we first summarize the dataset with a table that would be suitable for inclusion in a scientific rerport. Then we show results from a naïve linear regression model (SAP version 1) explaining the variation in weight with variation in height, not taking into account results from IDA. We compare this to an IDA-informed linear regression model using log-transformations for weight and height.
Then we conduct multivariable linear regression analyses as outlined in the statistical analysis plan that was updated after IDA.
Data can be read also online - use the code below.
# read data online
pokemon <- read.delim("https://stratosida.github.io/assets/Pokemon_clean.tsv",
strip.white = TRUE, na.string =c(NA, ""), blank.lines.skip = TRUE)
pokemon$PokemonPhase_5cat <- factor(pokemon$PokemonPhase_5cat, levels=c("Baby", "Base", "Phase1", "Phase2", "NotEvolving"))
pokemon$special_status <- factor(pokemon$special_status, levels=c("normal", "sub-legendary", "legendary", "mythical"))
pokemon <- dplyr::arrange(pokemon, pokedex_number)
pokemon$gender_cat<-factor(pokemon$gender_cat, levels= c("Equal", "Mostly females", "Mostly males", "Genderless"), labels=c("Equal", "Mostly female", "Mostly male","Genderless"))A0 - Description of the dataset
We summarize the main characteristics of the analysis sample. This table is more concise than the descriptive statistics presented in the IDA report and is intended to accompany the regression results, as would typically appear in a scientific paper.
### Table 1. Descriptive statistics of the analysis sample
s <- Hmisc::summaryM(
weight_kg + height_m + special_status + gender_cat + PokemonPhase_5cat +
hp + attack + defense + sp_attack + sp_defense + speed + happiness_cat ~ 1,
data = Pokemon_clean,
overall = FALSE,
test = FALSE
)
Hmisc::html(
s,
caption = 'Table 1. Descriptive statistics of variables included in the regression models',
exclude1 = TRUE,
npct = 'both',
digits = 2,
prmsd = TRUE,
brmsd = TRUE
)| Table 1. Descriptive statistics of variables included in the regression models. | |
| N=1025 |
|
|---|---|
| weight_kg | 8.5 28.0 70.0 67.0 ± 121.2 |
| height_m | 0.5 1.0 1.5 1.2 ± 1.2 |
| special_status : normal | 0.90 920⁄1025 |
| sub-legendary | 0.05 51⁄1025 |
| legendary | 0.03 31⁄1025 |
| mythical | 0.02 23⁄1025 |
| gender_cat : Equal | 0.62 631⁄1025 |
| Mostly female | 0.06 64⁄1025 |
| Mostly male | 0.17 176⁄1025 |
| Genderless | 0.15 154⁄1025 |
| PokemonPhase_5cat : Baby | 0.02 17⁄1025 |
| Base | 0.33 342⁄1025 |
| Phase1 | 0.35 354⁄1025 |
| Phase2 | 0.11 111⁄1025 |
| NotEvolving | 0.20 201⁄1025 |
| hp | 50 68 85 70 ± 27 |
| attack | 55 75 100 78 ± 30 |
| defense | 50 70 90 73 ± 29 |
| sp_attack | 47 65 90 70 ± 30 |
| sp_defense | 50 67 86 70 ± 27 |
| speed | 45 65 88 67 ± 29 |
| happiness_cat : Happy | 0.03 34⁄1025 |
| Normal | 0.81 835⁄1025 |
| Unhappy | 0.15 156⁄1025 |
| a b c represent the lower quartile a, the median b, and the upper quartile c for continuous variables. x ± s represents X ± 1 SD. | |
A1 - Univariable linear regression models
Here weight is the outcome and height is the explanatory variable.
SAP v1: Naïve linear regression model
Here we ignore that height and weight have skewed distributions.
naive_first_fit <- lm(weight_kg ~ height_m, data = Pokemon_clean, x=T, y=T)The results are summarized below.
summary(naive_first_fit) %>% tab_model()| weight kg | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -7.49 | -15.52 – 0.55 | 0.068 |
| height m | 61.38 | 56.76 – 66.00 | <0.001 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.399 / 0.399 | ||
#par(mfrow=c(2,2))
#plot(naive_first_fit)
autoplot(naive_first_fit, which = 1:6, ncol = 2) #ggfortify package
RESULTS:
If we compare two Pokémon that differ by 1 m, the one that is higher is on average 61.4 kg heavier.
The naïve fit violates some assumptions of linear regression. For example, the Residuals vs Fitted and Scale-Location plots show a violation of homogeneous residual variance, with residual variance increasing as fitted values increase. Further, a QQ-plot of residuals shows a violation of normally distributed residuals.
The Pokémon with largest residuals are those with very large height, whose weight is strongly underestimated by the model.
SAP v2: Linear regression model with log-transformed variables
We fit a second linear regression model, where weight and height are log2-transformed, and summarize its findings.
IDA_first_fit <- lm(log2(weight_kg) ~ log2(height_m), data = Pokemon_clean, x=T, y=T)summary(IDA_first_fit) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.78 | 4.70 – 4.86 | <0.001 |
| height m [log2] | 1.75 | 1.68 – 1.83 | <0.001 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.657 / 0.657 | ||
#par(mfrow=c(2,2))
#plot(IDA_first_fit)
autoplot(IDA_first_fit, which = 1:6, ncol = 2) #ggfortify package
RESULTS:
The model estimates that if we compare two Pokémon, where one is twice the height of the other, the weight of the taller Pokémon is, on average, 3.4 times greater.
When examining the IDA-informed fit, we see that the issues of heteroskedasticity and non-normal residual variance are largely resolved. Further, the regression output shows that the R^2 for the IDA-informed fit is much higher than that of the naïve fit (0.66 vs 0.40).
We notice that the Pokémon with largest residuals are Cosmoem (ID=790), Gastly and Haunter (ID=92 and 93), which have unusual weight/height combinations (or equivalently, BMI values).
In the graph shown below (in log2 scale on both axes) we display the estimated values of weight using original data (red), the log2-transformed weight and height (gray), and the estimated association using a smoothing function (in blue).
RESULTS: We observe that the model based on the non-transformed data substantially overestimates weights in the the lower range of heights. The fit of the model with log-transformed data is poor for large heights (after about 2m), where the weight is substantially overestimated, indicating that the relationship between height and weight cannot be captured well with assuming a linear relationship on the log2 scale with a univariable model.
To address the problem of overestimation of weight for large heights, we fit a model that includes also the height term (shown in yellow in the graph).
IDA_quad_fit <- lm(log2(weight_kg) ~ log2(height_m) + height_m, data = Pokemon_clean, x=T, y=T)
anova(IDA_first_fit, IDA_quad_fit)## Analysis of Variance Table
##
## Model 1: log2(weight_kg) ~ log2(height_m)
## Model 2: log2(weight_kg) ~ log2(height_m) + height_m
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1023 1892.4
## 2 1022 1849.5 1 42.895 23.703 1.302e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1x.pred <- data.frame(height_m = seq(0.1, 20, by=0.1))
data.plot <- data.frame(x=x.pred[,1],y=2^predict(IDA_quad_fit, x.pred))
# species to label inthe plot
rows_to_label <- c(790, 92, 93, 890, 977)
my.plot <- ggplot(Pokemon_clean, aes(height_m, weight_kg)) + geom_point(pos="jitter", alpha=.7) + scale_x_continuous(trans="log2")+ scale_y_continuous(trans="log2") + labs(x="Height (m), log 2 scale", y="Weight (kg) - log 2 scale") +
geom_smooth(aes(col="GAM smoother")) +
geom_line(aes(height_m, naive_first_fit$fitted.values, col="Weight~Height"), cex=2) + geom_line(aes(height_m, 2^IDA_first_fit$fitted.values, col="log2(Weight)~log2(Height)"), cex=2)+
geom_line(data=data.plot, aes(x,y,col="log2(Weight)~log2(Height)+Height"), cex=2) +
scale_colour_manual("Regression model",
#breaks = c("A", "B", "log Weight ~ log Height + Height ","Loess smoother"),
breaks = c("Weight~Height", "log2(Weight)~log2(Height)", "log2(Weight)~log2(Height)+Height","GAM smoother"),
values=c("Weight~Height"="red", "log2(Weight)~log2(Height)"="yellow", "log2(Weight)~log2(Height)+Height"="purple", "GAM smoother"="blue")
) + guides(colour = guide_legend(nrow = 2, byrow = TRUE)) + guides(colour = guide_legend(nrow = 2, byrow = TRUE)) +
theme_bw() + theme(
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9)
)
my.plot + geom_label_repel(data=Pokemon_clean[rows_to_label,], aes(label=name)); 
# fig UnivariableModels for the paper or supplement
######## paper figure
invisible(tiff("pokemon_plot_regression.tiff", width = 6, height = 6, units = "in", res = 600, compression = "lzw"))
# Build the plot
my.plot <- ggplot(Pokemon_clean, aes(x = height_m, y = weight_kg)) +
geom_point(position = position_jitter(width = 0.1, height = 0.1), alpha = 0.7, size = 1) +
scale_x_continuous(trans = "log2") +
scale_y_continuous(trans = "log2") +
labs(
x = "Height (m), log2 scale",
y = "Weight (kg), log2 scale"
) +
geom_smooth(aes(col = "GAM smoother"), method = "gam", se = FALSE, size = 0.8) +
geom_line(aes(y = naive_first_fit$fitted.values, col = "Weight~Height"), size = 0.8) +
geom_line(aes(y = 2^IDA_first_fit$fitted.values, col = "log2(Weight)~log2(Height)"), size = 0.8) +
geom_line(data = data.plot, aes(x = x, y = y, col = "log2(Weight)~log2(Height)+Height"), size = 0.8) +
scale_colour_manual(
name = "Regression model",
breaks = c(
"Weight~Height",
"log2(Weight)~log2(Height)",
"log2(Weight)~log2(Height)+Height",
"GAM smoother"
),
values = c(
"Weight~Height" = "red",
"log2(Weight)~log2(Height)" = "orange",
"log2(Weight)~log2(Height)+Height" = "purple",
"GAM smoother" = "blue"
)
) + guides(colour = guide_legend(nrow = 2, byrow = TRUE)) +
theme_bw(base_size = 11) +
theme(
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9)
) +
geom_label_repel(
data = Pokemon_clean[rows_to_label,],
aes(label = name),
size = 2.5,
max.overlaps = 10,
box.padding = 0.3,
point.padding = 0.2,
segment.size = 0.3
)
# Render the plot
print(my.plot);
# Close TIFF device
invisible(dev.off())Diagnostic plots
We present the diagnostic plots from the model.
#regression diagnostics
autoplot(IDA_quad_fit, which = 1:6, ncol = 2)
RESULTS: The diagnostic plots suggest that a few species (Cosmoem, Gastly, Haunter, Dondozo, Eternatus) may be influential in this univariable model, either due to unusual weight-for-height (BMI) combinations or extreme height values. We do not formally exclude them at this stage; their influence is assessed systematically in the multivariable model below (SA3, SA4).
Now we limit the estimation to the Pokémon with height within 5 meters (removing 16 Pokémon). As a preview of SA2 (formally conducted on the multivariable model below), we limit the estimation to the Pokémon with height within 5 meters (removing 16 Pokémon).
IDA_quad_fit_sens2 <- lm(log2(weight_kg) ~ log2(height_m) + height_m, data = Pokemon_clean[Pokemon_clean$height_m<5,], x=T, y=T)
#x.pred <- data.frame(height_m = seq(0.1, 20, by=0.1))
#data.plot.sens1 <- data.frame(x=x.pred[,1],y=2^predict(IDA_first_fit_sens1, x.pred))
summary(IDA_quad_fit_sens2) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 5.03 | 4.68 – 5.39 | <0.001 |
| height m [log2] | 1.94 | 1.74 – 2.15 | <0.001 |
| height m | -0.18 | -0.46 – 0.10 | 0.216 |
| Observations | 1009 | ||
| R2 / R2 adjusted | 0.651 / 0.650 | ||
RESULTS:
In the range from 0.4 to 5 m the estimated relationship between height and weight remains similar to the one obtained using all data (in green in the graph). The standard error of the regression coefficients increase substantially.
my.plot.quad <- ggplot(Pokemon_clean, aes(height_m, weight_kg)) + geom_point(pos="jitter", alpha=.7) + scale_x_continuous(trans="log2")+ scale_y_continuous(trans="log2") + labs(x="Height (m), log 2 scale", y="Weight (kg) - log 2 scale") +
geom_smooth(aes(col="GAM smoother")) +
geom_line(data=data.plot, aes(x,y,col="log2(Weight)~log2(Height)+Height"), cex=2) +
scale_colour_manual("Regression model",
#breaks = c("A", "B", "log Weight ~ log Height + Height ","Loess smoother"),
breaks = c("Weight~Height", "log2(Weight)~log2(Height)", "log2(Weight)~log2(Height)+Height","GAM smoother", "Excluding taller than 5 m"),
values=c("Weight~Height"="red", "log2(Weight)~log2(Height)"="yellow", "log2(Weight)~log2(Height)+Height"="purple", "GAM smoother"="blue", "Excluding taller than 5 m"="green")
) + theme_bw() + guides(colour = guide_legend(nrow = 2, byrow = TRUE)) + theme(
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9)
)
data.plot.sens2 <- data.frame(x=IDA_quad_fit_sens2$x[,2], y=IDA_quad_fit_sens2$fitted)
my.plot.quad + geom_label_repel(data=Pokemon_clean %>% filter(BMI<0.1 | BMI>1000), aes(label=name)) + geom_line(data=data.plot.sens2, aes(2^x, 2^y, col="Excluding taller than 5 m"), cex=1) +
geom_label_repel(data=Pokemon_clean[rows_to_label,], aes(label=name))
A2 - Multivariable linear regression models
Log2-transformed weight is the outcome variable. Height is the key explanatory variable, modeled using both log2-transformed height and untransformed height (linear term), as motivated by the univariable analysis above. For simplicity, we refer to these as ‘weight’ and ‘height’ below.
Inclusion of covariates
Based on the initial data analysis, we include the categorical gender variable, base statistics, categorical happiness, special status, and evolutionary phase as covariates in addition to height.
Evolutionary stage is strongly associated with the other covariates. We nonetheless retain it in the model, since it was regarded as a key explanatory variable in the study plan; its association with the other covariates suggests that its independent contribution to explaining the outcome may be limited, a question addressed through variable selection below.
We include all six base statistics rather than their total, since IDA suggested that they carry some information not captured by the total. As the base stats are correlated, some may be redundant in explaining the outcome; this is also addressed through variable selection.
Catch rate is excluded, given its strong correlation with the base stats observed during IDA.
COMMENT: To take the expected correlations between explanatory variables into account, we amended the SAP and decided to use a backward variable selection method based on AIC to eliminate possible explanatory variables that do not contribute to the explanation of the outcome.
Gender is defined as four categories (genderless, equal male and female, mostly male, and mostly female). The reference level is ‘Equal’ to understand the difference in estimated weight compared to categories mostly female, mostly male, or genderless. Happiness is defined as a categorical variable with normal level as reference.
Full model with all variables
Pokemon_clean$gender_cat <- factor(Pokemon_clean$gender_cat )
Pokemon_clean$gender_cat <-relevel(Pokemon_clean$gender_cat, ref = "Equal")
Pokemon_clean$PokemonPhase_5cat <-relevel(Pokemon_clean$PokemonPhase_5cat, ref = "Base")
Pokemon_clean$happiness_cat <-relevel(factor(Pokemon_clean$happiness_cat), ref = "Normal")
mv_fit <- 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) + special_status + happiness_cat + PokemonPhase_5cat , data = Pokemon_clean, x=T, y=T)
summary(mv_fit) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 3.98 | 3.40 – 4.55 | <0.001 |
| height m [log2] | 1.77 | 1.63 – 1.92 | <0.001 |
| height m | -0.32 | -0.42 – -0.21 | <0.001 |
|
gender cat [Mostly female] |
-0.33 | -0.67 – -0.00 | 0.050 |
| gender cat [Mostly male] | 0.24 | 0.03 – 0.46 | 0.027 |
| gender cat [Genderless] | 0.62 | 0.29 – 0.95 | <0.001 |
| hp/10 | 0.13 | 0.09 – 0.17 | <0.001 |
| attack/10 | 0.03 | -0.01 – 0.06 | 0.174 |
| defense/10 | 0.11 | 0.07 – 0.15 | <0.001 |
| sp attack/10 | -0.08 | -0.11 – -0.04 | <0.001 |
| sp defense/10 | 0.01 | -0.03 – 0.05 | 0.571 |
| speed/10 | -0.04 | -0.07 – -0.00 | 0.031 |
|
special status [sub-legendary] |
-0.46 | -0.92 – 0.00 | 0.052 |
|
special status [legendary] |
-0.07 | -0.63 – 0.49 | 0.799 |
| special status [mythical] | -0.72 | -1.36 – -0.08 | 0.027 |
| happiness cat [Happy] | -0.27 | -0.76 – 0.22 | 0.276 |
| happiness cat [Unhappy] | 0.24 | -0.03 – 0.52 | 0.085 |
| PokemonPhase 5cat [Baby] | 0.47 | -0.15 – 1.09 | 0.137 |
|
PokemonPhase 5cat [Phase1] |
-0.11 | -0.38 – 0.16 | 0.427 |
|
PokemonPhase 5cat [Phase2] |
-0.04 | -0.42 – 0.34 | 0.837 |
|
PokemonPhase 5cat [NotEvolving] |
-0.29 | -0.64 – 0.05 | 0.097 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.720 / 0.714 | ||
anova(mv_fit)## Analysis of Variance Table
##
## Response: log2(weight_kg)
## Df Sum Sq Mean Sq F value Pr(>F)
## log2(height_m) 1 3632.4 3632.4 2353.9801 < 2.2e-16 ***
## height_m 1 42.9 42.9 27.7984 1.648e-07 ***
## gender_cat 3 35.2 11.7 7.6010 4.990e-05 ***
## I(hp/10) 1 57.0 57.0 36.9130 1.751e-09 ***
## I(attack/10) 1 7.6 7.6 4.9284 0.026641 *
## I(defense/10) 1 84.6 84.6 54.8390 2.767e-13 ***
## I(sp_attack/10) 1 64.6 64.6 41.8457 1.540e-10 ***
## I(sp_defense/10) 1 0.0 0.0 0.0001 0.994215
## I(speed/10) 1 15.6 15.6 10.0785 0.001546 **
## special_status 3 17.8 5.9 3.8434 0.009438 **
## happiness_cat 2 8.8 4.4 2.8484 0.058405 .
## PokemonPhase_5cat 4 9.2 2.3 1.4932 0.202092
## Residuals 1004 1549.2 1.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#par(mfrow=c(2,2))
#plot(mv_fit)
autoplot(mv_fit, which = 1:6, ncol = 2) #ggfortify package
RESULTS: The use of additional explanatory variables increases the adjusted determination coefficient of the model from 0.66 to 0.71.
Backward selection of variables from the full model
Since we observed large associations between the potential explanatory variables in the IDA report, we decided to use backward variable selection based on AIC to eliminate variables that might not contribute to the explanation of variation in the outcome or inflate the standard errors.
mv_fit_step <- step(mv_fit, trace = FALSE)
summary(mv_fit_step) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.24 | 3.76 – 4.72 | <0.001 |
| height m [log2] | 1.80 | 1.66 – 1.94 | <0.001 |
| height m | -0.33 | -0.43 – -0.22 | <0.001 |
|
gender cat [Mostly female] |
-0.33 | -0.66 – -0.00 | 0.049 |
| gender cat [Mostly male] | 0.29 | 0.08 – 0.50 | 0.008 |
| gender cat [Genderless] | 0.56 | 0.24 – 0.88 | 0.001 |
| hp/10 | 0.13 | 0.09 – 0.17 | <0.001 |
| defense/10 | 0.11 | 0.08 – 0.14 | <0.001 |
| sp attack/10 | -0.08 | -0.11 – -0.05 | <0.001 |
| speed/10 | -0.04 | -0.07 – -0.01 | 0.007 |
|
special status [sub-legendary] |
-0.51 | -0.95 – -0.06 | 0.025 |
|
special status [legendary] |
-0.10 | -0.64 – 0.45 | 0.733 |
| special status [mythical] | -0.72 | -1.35 – -0.10 | 0.024 |
| happiness cat [Happy] | -0.29 | -0.77 – 0.19 | 0.243 |
| happiness cat [Unhappy] | 0.27 | -0.01 – 0.54 | 0.057 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.718 / 0.714 | ||
RESULTS: The backward elimination based on AIC removed evolutionary stage, attack and special defense, selecting a model with the same adjusted determination coefficient as the full model. The estimated association between height and weight is similar in the full and the selected models.
Residuals vs. height (linearity check)
ggplot(data.frame(height_m = mv_fit_step$model$height_m, resid = residuals(mv_fit_step)),
aes(x = height_m, y = resid)) +
geom_point(alpha = 0.5) +
geom_smooth(se = FALSE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(x = "Height (m)", y = "Residuals", title = "Residuals vs. Height (selected multivariable model)") +
theme_bw()
x.pred <- data.frame(height_m = seq(0.1, 20, by=0.1), gender_cat="Equal", hp =100,
attack=100, defense=100, sp_attack=100, sp_defense=100, speed=100, PokemonPhase_5cat="Base", special_status = "normal", happiness_cat="Normal"
)
tmp <- predict(mv_fit, newdata=x.pred)
data.plot2 <- data.frame(x=x.pred[,1], y=2^tmp)
tmp <- predict(mv_fit_step, newdata=x.pred)
data.plot <- data.frame(x=x.pred[,1], y=2^tmp)
ggplot(Pokemon_clean, aes(height_m, weight_kg)) + geom_point(pos="jitter") +
#geom_line(aes(height_m, 2^IDA_quad_fit$fitted.values), col="red", cex=2) +
scale_x_continuous(trans="log2")+ scale_y_continuous(trans="log2") + labs(x="Height (m), log 2 scale", y="Weight (kg) - log 2 scale") +
scale_colour_manual("Regression model",
breaks = c("Selected multivariable model","GAM smoother", "Excluding taller than 5 m","Full model"),
values=c("Selected multivariable model"="purple", "GAM smoother"="blue", "Excluding taller than 5 m"="green", "Full model"="red")
) + theme_bw()+ guides(colour = guide_legend(nrow = 2, byrow = TRUE)) + theme(
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 9),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9)
) + geom_line(data=data.plot, aes(x, y, col="Selected multivariable model"), cex=2) + #mult
geom_line(data=data.plot2, aes(x, y, col="Full model"), cex=2) + #atypical bmi
geom_smooth(aes(col="GAM smoother"))
Sensitivity analyses
Correlation within evolutionary family (SA1)
Observations may not be fully independent, as species within the same evolutionary family are expected to be more similar. This is partially addressed by including evolutionary phase as a covariate; the residual correlation within evolutionary family is assessed directly below (SA1), using cluster-robust standard errors rather than mixed-effects models or GEE, since evolutionary families are typically small (see IDA additional report), with many families including a single species, making reliable estimation of a random-effects variance impractical.
To take into account the correlation within Pokemon evolutionary family, we fit a model that includes the variables selected with the backward procedure and derive the robust standard errors of the estimated regression coefficients.
library(sandwich)
# Cluster-robust standard errors and P-values (evolutionary family),
# using the same format as for the previous models
tab_model(mv_fit_step,
vcov.fun = "vcovCL",
vcov.args = list(cluster = ~PokemonID))| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.24 | 3.48 – 5.01 | <0.001 |
| height m [log2] | 1.80 | 1.56 – 2.04 | <0.001 |
| height m | -0.33 | -0.47 – -0.18 | <0.001 |
|
gender cat [Mostly female] |
-0.33 | -0.76 – 0.10 | 0.134 |
| gender cat [Mostly male] | 0.29 | 0.09 – 0.48 | 0.004 |
| gender cat [Genderless] | 0.56 | 0.12 – 1.01 | 0.014 |
| hp/10 | 0.13 | 0.07 – 0.18 | <0.001 |
| defense/10 | 0.11 | 0.06 – 0.16 | <0.001 |
| sp attack/10 | -0.08 | -0.12 – -0.04 | <0.001 |
| speed/10 | -0.04 | -0.08 – -0.01 | 0.024 |
|
special status [sub-legendary] |
-0.51 | -1.04 – 0.03 | 0.062 |
|
special status [legendary] |
-0.10 | -0.63 – 0.44 | 0.726 |
| special status [mythical] | -0.72 | -1.53 – 0.09 | 0.080 |
| happiness cat [Happy] | -0.29 | -0.86 – 0.28 | 0.324 |
| happiness cat [Unhappy] | 0.27 | -0.13 – 0.66 | 0.182 |
| Observations | 1025 | ||
| R2 / R2 adjusted | 0.718 / 0.714 | ||
# alternative coding
# Cluster-robust standard errors by cluster ID
# cluster_se <- vcovCL(mv_fit_step, cluster = ~PokemonID)
# Robust inference
# ct <- coeftest(mv_fit_step, vcov = cluster_se)
# Display results
#kable(as.table(ct), digits = 3, caption = "Cluster-robust standard errors and P-values") RESULTS: As expected, taking into account the correlation within evolutionary family increased the standard errors of most of the estimated regression coefficients.
In terms of statistical significance, two of the three special status levels (sub-legendary, mythical) and the “Mostly female” level of gender_cat move from p<0.05 to p>0.05 after this correction. Happiness_cat remains non-significant, as it already was in the non-robust model. All other predictors remain significant.
Extreme height (SA2)
We refit the selected model limited to Pokémon shorter than 5 m.
# sensitivity to large heights
mv_fit_step_sens1 <- update(mv_fit_step, data=Pokemon_clean[Pokemon_clean$height_m<5,])
summary(mv_fit_step_sens1) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.53 | 3.96 – 5.09 | <0.001 |
| height m [log2] | 1.96 | 1.75 – 2.17 | <0.001 |
| height m | -0.57 | -0.86 – -0.28 | <0.001 |
|
gender cat [Mostly female] |
-0.33 | -0.66 – -0.00 | 0.047 |
| gender cat [Mostly male] | 0.27 | 0.06 – 0.48 | 0.014 |
| gender cat [Genderless] | 0.60 | 0.28 – 0.92 | <0.001 |
| hp/10 | 0.13 | 0.09 – 0.17 | <0.001 |
| defense/10 | 0.11 | 0.08 – 0.15 | <0.001 |
| sp attack/10 | -0.08 | -0.12 – -0.05 | <0.001 |
| speed/10 | -0.04 | -0.08 – -0.01 | 0.006 |
|
special status [sub-legendary] |
-0.52 | -0.98 – -0.07 | 0.025 |
|
special status [legendary] |
0.10 | -0.50 – 0.70 | 0.737 |
| special status [mythical] | -0.73 | -1.36 – -0.10 | 0.023 |
| happiness cat [Happy] | -0.30 | -0.79 – 0.18 | 0.221 |
| happiness cat [Unhappy] | 0.26 | -0.02 – 0.53 | 0.069 |
| Observations | 1009 | ||
| R2 / R2 adjusted | 0.707 / 0.703 | ||
RESULTS: The estimated association between height and weight is similar as with the full model. See picture below.
Identification of species with large influence (SA3)
Below we display the measures of influence (leverage, Cook’s distance and studentized residuals) obtained based on the selected multivariable model, for the species that have the value of at least one of the measures among the top 40 species; threshold values are also displayed graphically (\(n\) is the sample size, \(p\) is the number of variables in the model; 2p/n for leverage, 4/n for Cook’s distance, 2 for absolute value of studentized residuals)
# ----------------------------
# 1. Compute diagnostic measures and define cut-off values
# ----------------------------
lev <- hatvalues(mv_fit_step) # leverage
cooks <- cooks.distance(mv_fit_step) # Cook's distance
stud <- rstudent(mv_fit_step) # studentized residuals
n <- nrow(mv_fit_step$model)
p <- length(coef(mv_fit_step))
cutoff_lev <- 2 * p / n
mean_lev <- p / n
cutoff_cooks <- 4 / n # common Cook’s distance rule-of-thumb
cutoff_resid <- 2 # studentized residual threshold
# Align rows with model frame
idx <- as.numeric(rownames(mv_fit_step$model))
# ----------------------------
# 2. Combine into a data frame
# ----------------------------
diag_tbl <- Pokemon_clean %>%
mutate(
lev = lev,
cooks = cooks,
stud_res = stud
) %>%
select(name, special_status, type1, type2, lev, cooks, stud_res)
# ----------------------------
# Select top N influential Pokémon
# ----------------------------
top_n <- 40
top_diag <- diag_tbl %>%
arrange(desc(lev)) %>%
slice(1:top_n)
# Pivot into long format for plotting
top_diag_long <- top_diag %>%
pivot_longer(cols = c("lev", "cooks", "stud_res"),
names_to = "measure", values_to = "value")
# Order Pokémon names by leverage (highest first)
top_diag_long$name <- factor(top_diag_long$name, levels = top_diag$name)
# Thresholds as a data frame (so we can facet with them)
thresholds <- data.frame(
measure = c("lev", "lev", "cooks", "stud_res", "stud_res"),
value = c(mean_lev, cutoff_lev, cutoff_cooks, -cutoff_resid, cutoff_resid),
type = c("mean leverage", "cutoff leverage", "Cook's threshold", "±2 resid", "±2 resid")
)
# ----------------------------
# Plot: vertical alignment (Pokémon names on Y-axis, facets for measures)
# ----------------------------
ggplot(top_diag_long, aes(x = value, y = fct_rev(name), color = special_status)) +
geom_point(size = 3) +
geom_vline(data = thresholds, aes(xintercept = value, linetype = type),
color = "black", inherit.aes = FALSE) +
facet_wrap(~measure, scales = "free_x") +
scale_linetype_manual(values = c("dotted", "dashed", "dashed", "dashed", "dashed")) +
theme_minimal() +
labs(
x = "Influence / Residual",
y = NULL,
title = "Top 40 large leverage Pokémon in Weight Regression",
color = "Special Status",
linetype = "Threshold"
)
# ----------------------------
# Select top N Cook distance Pokémon
# ----------------------------
top_n <- 40
top_diag <- diag_tbl %>%
arrange(desc(cooks)) %>%
slice(1:top_n)
# Pivot into long format for plotting
top_diag_long <- top_diag %>%
pivot_longer(cols = c("lev", "cooks", "stud_res"),
names_to = "measure", values_to = "value")
# Order Pokémon names by leverage (highest first)
top_diag_long$name <- factor(top_diag_long$name, levels = top_diag$name)
# ----------------------------
# Define thresholds
# ----------------------------
n <- nrow(mv_fit_step$model)
p <- length(coef(mv_fit_step))
cutoff_lev <- 2 * p / n
mean_lev <- p / n
cutoff_cooks <- 4 / n # common Cook’s distance rule-of-thumb
cutoff_resid <- 2 # studentized residual threshold
# Thresholds as a data frame (so we can facet with them)
thresholds <- data.frame(
measure = c("lev", "lev", "cooks", "stud_res", "stud_res"),
value = c(mean_lev, cutoff_lev, cutoff_cooks, -cutoff_resid, cutoff_resid),
type = c("mean leverage", "cutoff leverage", "Cook's threshold", "±2 resid", "±2 resid")
)
# ----------------------------
# Plot: vertical alignment (Pokémon names on Y-axis, facets for measures)
# ----------------------------
ggplot(top_diag_long, aes(x = value, y = fct_rev(name), color = special_status)) +
geom_point(size = 3) +
geom_vline(data = thresholds, aes(xintercept = value, linetype = type),
color = "black", inherit.aes = FALSE) +
facet_wrap(~measure, scales = "free_x") +
scale_linetype_manual(values = c("dotted", "dashed", "dashed", "dashed", "dashed")) +
theme_minimal() +
labs(
x = "Influence / Residual",
y = NULL,
title = "Top 40 large Cook's distance Pokémon in Weight Regression",
color = "Special Status",
linetype = "Threshold"
)
As an illustrative sensitivity analysis, we refit the selected model after removing all species with Cook’s distance above a commonly used threshold (4/n), rather than examining each flagged species individually; this mechanical removal is a simplification for teaching purposes and is not a general recommendation for handling influential points in practice. Overall, 62 species are removed from the analysis. The excluded species are highlighted in the height-weight scatterplot.
idx_infl <- which(cooks > cutoff_cooks)
Pokemon_clean$excluded <- ifelse(1:nrow(Pokemon_clean) %in% idx_infl,
"Excluded (Cook's D > cutoff)",
"Included")
# ----------------------------
# 1. Scatterplot with jitter, original scales on axes
# ----------------------------
p1 <- ggplot(Pokemon_clean, aes(x = log2(height_m), y = log2(weight_kg), color = excluded)) +
geom_jitter(size = 3, alpha = 0.8, width = 0.05, height = 0.05) +
scale_x_continuous(name = "Height (m)", labels = function(x) round(2^x, 2)) +
scale_y_continuous(name = "Weight (kg)", labels = function(y) round(2^y, 1)) +
theme_minimal(base_size = 12) +
labs(
title = "Height vs. Weight with Jitter",
subtitle = "Points colored by inclusion/exclusion based on Cook's Distance",
color = "Status"
)
p1
We provide a list of the names of the removed species.
# ----------------------------
# 1. Influence measures
# ----------------------------
#cooks <- cooks.distance(mv_fit_step)
#n <- nrow(mv_fit_step$model)
#cutoff_cooks <- 4 / n # standard cutoff
# ----------------------------
# 2. Identify influential Pokémon
# ----------------------------
#idx_infl <- which(cooks > cutoff_cooks)
cat("Influential Pokémon removed (Cook's D > cutoff):\n")## Influential Pokémon removed (Cook's D > cutoff):print(Pokemon_clean$name[idx_infl])## [1] "Ekans" "Gastly" "Haunter" "Chansey" "Dratini"
## [6] "Dragonair" "Misdreavus" "Wobbuffet" "Blissey" "Larvitar"
## [11] "Shedinja" "Aron" "Cacnea" "Beldum" "Jirachi"
## [16] "Deoxys" "Happiny" "Munchlax" "Rotom" "Uxie"
## [21] "Mesprit" "Azelf" "Heatran" "Manaphy" "Cryogonal"
## [26] "Durant" "Tornadus" "Thundurus" "Doublade" "Hoopa"
## [31] "Volcanion" "Pyukumuku" "Minior" "Cosmog" "Cosmoem"
## [36] "Lunala" "Celesteela" "Kartana" "Magearna" "Poipole"
## [41] "Meltan" "Melmetal" "Silicobra" "Polteageist" "Hatterene"
## [46] "Cursola" "Alcremie" "Eternatus" "Regieleki" "Glastrier"
## [51] "Calyrex" "Enamorus" "Tinkatuff" "Tinkaton" "Wiglett"
## [56] "Dondozo" "Scream Tail" "Brute Bonnet" "Flutter Mane" "Iron Treads"
## [61] "Terapagos" "Pecharunt"The figure below shows the impact of the removal of the influential species on the fitted regression coefficients.
# ----------------------------
# 3. Refit model without influential Pokémon
# ----------------------------
Pokemon_no_infl <- Pokemon_clean[-idx_infl, ]
mv_fit_clean <- lm(
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_no_infl
)
# ----------------------------
# 4. Extract coefficients (with 95% CI) from both models
# ----------------------------
coef_orig <- tidy(mv_fit_step, conf.int = TRUE) %>%
mutate(model = "Original")
coef_clean <- tidy(mv_fit_clean, conf.int = TRUE) %>%
mutate(model = "Without influential points")
coef_compare <- bind_rows(coef_orig, coef_clean)
# ----------------------------
# 5. Plot coefficient comparison
# ----------------------------
ggplot(coef_compare, aes(x = estimate, y = term, color = model)) +
geom_point(position = position_dodge(width = 0.6)) +
geom_errorbar(aes(xmin = conf.low, xmax = conf.high),
width = 0.2,
position = position_dodge(width = 0.6)) +
geom_vline(xintercept = 0, linetype = "dotted", color = "grey50") +
theme_minimal(base_size = 12) +
labs(
title = "Comparison of Regression Coefficients",
subtitle = "Original model vs. model without influential Pokémon (Cook's D > cutoff)",
x = "Coefficient estimate (with 95% CI)",
y = "Predictor",
color = "Model"
)
RESULTS: The estimated coefficients with and without influential points are similar. The estimated coefficients for happiness are closer to 0, while the coefficients estimated for special statuses and for genderless species are larger. The adjusted determination coefficient increases substantially from 0.71 to 0.81 and the residual standard error decreases.
Overall, the estimated association between height and weight is similar to the one obtained using all data; the most marked difference is in the large range of heights (see picture below).
# ----------------------------
# 1. Mark excluded vs included Pokémon
# ----------------------------
Pokemon_clean$excluded <- ifelse(1:nrow(Pokemon_clean) %in% idx_infl,
"Excluded",
"Included")
# ----------------------------
# 2. Prepare prediction dataset (typical values for other covariates)
# ----------------------------
num_medians <- sapply(Pokemon_clean %>% select(hp, defense, sp_attack, speed), median)
cat_refs <- sapply(Pokemon_clean %>% select(gender_cat, special_status, happiness_cat), function(x) levels(x)[1])
height_seq <- seq(min(Pokemon_clean$height_m), max(Pokemon_clean$height_m), length.out = 100)
pred_df_clean <- data.frame(
height_m = height_seq,
log2.height_m = log2(height_seq),
hp = num_medians["hp"],
defense = num_medians["defense"],
sp_attack = num_medians["sp_attack"],
speed = num_medians["speed"],
gender_cat = cat_refs["gender_cat"],
special_status = cat_refs["special_status"],
happiness_cat = cat_refs["happiness_cat"]
)
# ----------------------------
# 3. Predictions for cleaned model
# ----------------------------
pred_clean <- predict(mv_fit_clean, newdata = pred_df_clean, se.fit = TRUE)
pred_df_clean$log2_weight_kg <- pred_clean$fit
pred_df_clean$se <- pred_clean$se.fit
pred_df_clean$lower <- pred_df_clean$log2_weight_kg - 1.96 * pred_df_clean$se
pred_df_clean$upper <- pred_df_clean$log2_weight_kg + 1.96 * pred_df_clean$se
pred_df_clean$model <- "Cleaned model"
# ----------------------------
# 4. Predictions for original model
# ----------------------------
pred_orig <- predict(mv_fit_step, newdata = pred_df_clean, se.fit = TRUE)
pred_df_orig <- pred_df_clean
pred_df_orig$log2_weight_kg <- pred_orig$fit
pred_df_orig$se <- pred_orig$se.fit
pred_df_orig$lower <- pred_df_orig$log2_weight_kg - 1.96 * pred_df_orig$se
pred_df_orig$upper <- pred_df_orig$log2_weight_kg + 1.96 * pred_df_orig$se
pred_df_orig$model <- "Original model"
# Combine for plotting
pred_df <- bind_rows(pred_df_clean, pred_df_orig)
# ----------------------------
# 5. Plot on log2 axes, tick labels in original scale
# ----------------------------
ggplot() +
# Model fits
geom_line(data = pred_df, aes(x = log2(height_m), y = log2_weight_kg, color = model), size = 1.2) +
geom_ribbon(data = pred_df, aes(x = log2(height_m), ymin = lower, ymax = upper, fill = model), alpha = 0.2, color = NA) +
# Actual points
geom_point(data = Pokemon_clean, aes(x = log2(height_m), y = log2(weight_kg), color = excluded),
size = 2, alpha = 0.8) +
scale_color_manual(values = c("Included" = "gray50", "Excluded" = "red", "Original model" = "blue", "Cleaned model" = "green")) +
scale_fill_manual(values = c("Original model" = "blue", "Cleaned model" = "green")) +
# Tick labels in original scale
scale_x_continuous(name = "Height (m)", labels = function(x) round(2^x, 2)) +
scale_y_continuous(name = "Weight (kg)", labels = function(y) round(2^y, 1)) +
theme_minimal(base_size = 12) +
labs(
title = "Estimated Weight vs. Height (log2 scale)",
subtitle = "Other covariates at median/reference; comparison of original vs. cleaned model",
color = "Points / Model",
fill = "Model Fit"
)
Comparison of IDA profile imbalance index to MDA influence diagnostics (continuation of SA3)
During IDA, a profile imbalance index was defined to identify Pokémon species with unusual combinations of base-stat values, using only the predictor distributions and without reference to the outcome. Here we compare the ten species with the highest profile imbalance index to the species identified as influential by Cook’s distance in the fitted model, to assess whether these influential points could have been anticipated from the predictor distributions alone, before model fitting
base_stats <- c("hp","attack","defense","sp_attack","sp_defense","speed")
z_stats <- scale(Pokemon_clean[, base_stats])
Pokemon_clean$profile_imbalance <- apply(z_stats, 1, sd)
top10_imbalance <- Pokemon_clean %>%
arrange(desc(profile_imbalance)) %>%
slice(1:10) %>%
pull(name)
influential_cooks <- diag_tbl %>%
filter(cooks > cutoff_cooks) %>%
pull(name)
overlap <- intersect(top10_imbalance, influential_cooks)
data.frame(species = union(top10_imbalance, influential_cooks)) %>%
mutate(
flagged_by_IDA_index = species %in% top10_imbalance,
flagged_by_cooks_distance = species %in% influential_cooks
) %>%
arrange(desc(flagged_by_IDA_index), desc(flagged_by_cooks_distance)) %>%
kbl(booktabs = TRUE,
caption = "Comparison: species flagged by the IDA profile imbalance index (top 10) vs. species exceeding the Cook's distance threshold (4/n) in the selected model") %>%
kable_styling()| species | flagged_by_IDA_index | flagged_by_cooks_distance |
|---|---|---|
| Blissey | TRUE | TRUE |
| Chansey | TRUE | TRUE |
| Wobbuffet | TRUE | TRUE |
| Regieleki | TRUE | TRUE |
| Shuckle | TRUE | FALSE |
| Guzzlord | TRUE | FALSE |
| Stakataka | TRUE | FALSE |
| Regice | TRUE | FALSE |
| Regidrago | TRUE | FALSE |
| Steelix | TRUE | FALSE |
| Ekans | FALSE | TRUE |
| Gastly | FALSE | TRUE |
| Haunter | FALSE | TRUE |
| Dratini | FALSE | TRUE |
| Dragonair | FALSE | TRUE |
| Misdreavus | FALSE | TRUE |
| Larvitar | FALSE | TRUE |
| Shedinja | FALSE | TRUE |
| Aron | FALSE | TRUE |
| Cacnea | FALSE | TRUE |
| Beldum | FALSE | TRUE |
| Jirachi | FALSE | TRUE |
| Deoxys | FALSE | TRUE |
| Happiny | FALSE | TRUE |
| Munchlax | FALSE | TRUE |
| Rotom | FALSE | TRUE |
| Uxie | FALSE | TRUE |
| Mesprit | FALSE | TRUE |
| Azelf | FALSE | TRUE |
| Heatran | FALSE | TRUE |
| Manaphy | FALSE | TRUE |
| Cryogonal | FALSE | TRUE |
| Durant | FALSE | TRUE |
| Tornadus | FALSE | TRUE |
| Thundurus | FALSE | TRUE |
| Doublade | FALSE | TRUE |
| Hoopa | FALSE | TRUE |
| Volcanion | FALSE | TRUE |
| Pyukumuku | FALSE | TRUE |
| Minior | FALSE | TRUE |
| Cosmog | FALSE | TRUE |
| Cosmoem | FALSE | TRUE |
| Lunala | FALSE | TRUE |
| Celesteela | FALSE | TRUE |
| Kartana | FALSE | TRUE |
| Magearna | FALSE | TRUE |
| Poipole | FALSE | TRUE |
| Meltan | FALSE | TRUE |
| Melmetal | FALSE | TRUE |
| Silicobra | FALSE | TRUE |
| Polteageist | FALSE | TRUE |
| Hatterene | FALSE | TRUE |
| Cursola | FALSE | TRUE |
| Alcremie | FALSE | TRUE |
| Eternatus | FALSE | TRUE |
| Glastrier | FALSE | TRUE |
| Calyrex | FALSE | TRUE |
| Enamorus | FALSE | TRUE |
| Tinkatuff | FALSE | TRUE |
| Tinkaton | FALSE | TRUE |
| Wiglett | FALSE | TRUE |
| Dondozo | FALSE | TRUE |
| Scream Tail | FALSE | TRUE |
| Brute Bonnet | FALSE | TRUE |
| Flutter Mane | FALSE | TRUE |
| Iron Treads | FALSE | TRUE |
| Terapagos | FALSE | TRUE |
| Pecharunt | FALSE | TRUE |
cat("Overlap:", length(overlap), "of", length(top10_imbalance), "IDA-flagged species were also flagged by Cook's distance.\n")## Overlap: 4 of 10 IDA-flagged species were also flagged by Cook's distance.RESULTS: TOUPDATE:4 of the 10 species with the highest profile imbalance index (Blissey, Chansey, Wobbuffet, Regieleki) are also among those with Cook’s distance above the 4/n threshold. This overlap is well above what would be expected by chance (about 0.6 species, given that 62 of 1025 species, roughly 6%, exceed the Cook’s distance threshold).
CONSEQUENCES: This indicates that species with unbalanced predictor profiles were, to some extent, anticipated as influential before model fitting, supporting the usefulness of IDA-stage screening. At the same time, most Cook’s-distance-flagged species (58 of 62) were not among the top 10 by profile imbalance, showing that predictor-based screening alone cannot anticipate most outcome-dependent influential points; the two approaches capture different, complementary aspects of the data.
BMI of Pokémon (SA4)
We remove the 5 Pokémon with the most extreme BMI values or identified as influential (same 5 species as in the univariable analysis).
# sensitivity to Pokémon with atypical BMI
mv_fit_step_sens2 <- update(mv_fit_step, data=Pokemon_clean[-c(790,92,93,890,977 ),])
summary(mv_fit_step_sens2) %>% tab_model()| log 2(weight kg) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.61 | 4.14 – 5.07 | <0.001 |
| height m [log2] | 1.97 | 1.83 – 2.11 | <0.001 |
| height m | -0.43 | -0.55 – -0.30 | <0.001 |
|
gender cat [Mostly female] |
-0.34 | -0.64 – -0.04 | 0.028 |
| gender cat [Mostly male] | 0.23 | 0.03 – 0.42 | 0.023 |
| gender cat [Genderless] | 0.55 | 0.25 – 0.85 | <0.001 |
| hp/10 | 0.11 | 0.08 – 0.15 | <0.001 |
| defense/10 | 0.09 | 0.06 – 0.12 | <0.001 |
| sp attack/10 | -0.07 | -0.10 – -0.04 | <0.001 |
| speed/10 | -0.05 | -0.08 – -0.02 | 0.002 |
|
special status [sub-legendary] |
-0.50 | -0.91 – -0.08 | 0.018 |
|
special status [legendary] |
-0.45 | -0.97 – 0.06 | 0.086 |
| special status [mythical] | -0.64 | -1.22 – -0.06 | 0.031 |
| happiness cat [Happy] | -0.29 | -0.73 – 0.16 | 0.206 |
| happiness cat [Unhappy] | 0.20 | -0.05 – 0.46 | 0.119 |
| Observations | 1020 | ||
| R2 / R2 adjusted | 0.751 / 0.748 | ||
RESULTS: As in the univariable case, the removal of the 5 Pokémon increases the adjusted determination coefficient (from 0.71 to 0.75) and decreases the residual standard error. The estimated regression coefficient of the legendary group changes from -0.10 to -0.44. As anticipated in IDA, the effect is due to the large influence of Cosmoem, which is a legendary Pokémon (a rare group, including only 31 Pokémon).
For both sensitivity analyses the estimated association between height and weight is similar as with the full model.
Gender coding (SA5)
We compare the coding of gender used in the main model (gender_cat, 3 categories among Pokémon with a defined gender) to a finer-grained coding using percentage_male as a 7-level factor. This analysis is restricted to the 871 Pokémon with a defined gender; it does not evaluate whether genderless Pokémon should be modeled as a separate category, which is a structural feature of the data rather than an alternative coding choice.
Pokemon_gendered <- droplevels(Pokemon_clean[Pokemon_clean$gender_cat != "Genderless", ])
mv_fit_step_gendercat <- update(mv_fit_step, data = Pokemon_gendered)
mv_fit_step_gendernum <- update(mv_fit_step, . ~ . - gender_cat + factor(percentage_male), data = Pokemon_gendered)
tab_model(mv_fit_step_gendercat, mv_fit_step_gendernum,
dv.labels = c("Gender as 3-category factor", "Percentage male as 7-level factor"))| Gender as 3-category factor | Percentage male as 7-level factor | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 4.10 | 3.60 – 4.60 | <0.001 | 3.66 | 3.01 – 4.31 | <0.001 |
| height m [log2] | 1.78 | 1.63 – 1.94 | <0.001 | 1.77 | 1.62 – 1.93 | <0.001 |
| height m | -0.42 | -0.55 – -0.29 | <0.001 | -0.42 | -0.55 – -0.29 | <0.001 |
|
gender cat [Mostly female] |
-0.45 | -0.75 – -0.15 | 0.004 | |||
| gender cat [Mostly male] | 0.23 | 0.04 – 0.43 | 0.019 | |||
| hp/10 | 0.14 | 0.10 – 0.18 | <0.001 | 0.14 | 0.10 – 0.18 | <0.001 |
| defense/10 | 0.12 | 0.08 – 0.15 | <0.001 | 0.12 | 0.09 – 0.15 | <0.001 |
| sp attack/10 | -0.08 | -0.11 – -0.05 | <0.001 | -0.08 | -0.11 – -0.05 | <0.001 |
| speed/10 | -0.03 | -0.06 – 0.01 | 0.101 | -0.03 | -0.06 – 0.01 | 0.101 |
|
special status [sub-legendary] |
0.51 | -0.31 – 1.33 | 0.222 | 0.59 | -0.25 – 1.43 | 0.167 |
|
special status [legendary] |
-0.29 | -1.32 – 0.74 | 0.583 | -0.20 | -1.26 – 0.86 | 0.710 |
| happiness cat [Happy] | -0.01 | -0.58 – 0.55 | 0.960 | -0.01 | -0.58 – 0.55 | 0.961 |
| happiness cat [Unhappy] | 0.16 | -0.14 – 0.47 | 0.294 | 0.17 | -0.13 – 0.48 | 0.260 |
| percentage male [12.5] | 1.11 | -0.51 – 2.73 | 0.179 | |||
| percentage male [25] | -0.18 | -0.76 – 0.40 | 0.541 | |||
| percentage male [50] | 0.42 | 0.03 – 0.81 | 0.035 | |||
| percentage male [75] | 1.22 | 0.58 – 1.86 | <0.001 | |||
| percentage male [87.5] | 0.60 | 0.18 – 1.03 | 0.006 | |||
| percentage male [100] | 0.48 | -0.11 – 1.07 | 0.108 | |||
| Observations | 871 | 871 | ||||
| R2 / R2 adjusted | 0.711 / 0.707 | 0.714 / 0.709 | ||||
AIC(mv_fit_step_gendercat, mv_fit_step_gendernum)## df AIC
## mv_fit_step_gendercat 14 2709.653
## mv_fit_step_gendernum 18 2709.588RESULTS:
The AIC values for the two models are nearly identical (2709.65 for the 3-category gender_cat model vs. 2709.59 for the 7-level percentage_male model, ΔAIC = 0.07), far below the threshold typically used to prefer one model over another. R² adjusted is also nearly identical (0.707 vs. 0.709). Estimates for all other covariates (height, base statistics, special status, happiness) are essentially unchanged between the two codings.
COMMENT:
The finer 7-level coding does not improve model fit relative to gender_cat, supporting the 3-category coding as an adequate simplification for the main model. The 7-level coding suggests a possible non-monotonic pattern (percentage_male = 75 has the largest, most significant coefficient), but given the wide confidence intervals and the lack of a consistent monotonic trend across levels, this is not sufficient evidence to prefer the more complex coding.
Interpretation of the results
Here we first interpret the results of the selected multivariable model.
Compared to the full model, backward selection based on AIC removed evolutionary phase, attack, and special defense. The selected model has the same adjusted R² as the full model, and the estimated height-weight association is nearly identical between the two, both in point estimates and in confidence intervals. Estimates for the other retained covariates are also similar between the two models, including for the sub-legendary category (-0.46 in the full model vs. -0.51 in the selected model), whose p-value shifts slightly, from 0.052 to 0.025, after removing the collinear variables.
Height remains the strongest predictor of weight. The model retains both the log2-transformed and the linear height terms; together they describe a relationship where weight increases with height, but at a decelerating rate for the tallest Pokémon, consistent with the pattern already observed in the univariable analysis and shown in the graph above.
The coefficient for the genderless category is positive and large, suggesting that Pokémon species who are genderless are heavier than other species, after accounting for height. This follows the expert understanding that many genderless Pokémon species are legendary, and legendary species tend to be larger than species with an equal male/female gender split. The coefficients for mostly male and mostly female are positive and negative, respectively, which follows a real-world interpretation of females of many species being smaller than males.
After accounting for the other variables, sub-legendary and mythical Pokémon are estimated to be lighter than normal; the legendary category itself is not distinguishable from normal in this model.
When interpreting base statistics, we use increases of 10 to interpret coefficients rather than a 1-point change, since domain knowledge and IDA univariable descriptions show that a 1-point change is not meaningful for these variables. HP and defense are positively associated to weight, speed and special attack negatively (after accounting for the other variables).
Unhappy Pokémon are estimated to be heavier than normal, happy lighter (but not statistically significantly different than normal).
Robustness of the results across sensitivity analyses
The estimated association between height and weight is consistent across all sensitivity analyses (SA1-SA5): its direction and approximate magnitude do not change materially when accounting for within-family correlation, restricting to shorter Pokémon, removing influential species, or recoding gender.
Some coefficients are less robust than others, however. Taking into account the correlation within evolutionary family (SA1) increases standard errors throughout, and moves the coefficients for sub-legendary and mythical special status, and for the “mostly female” gender category, from p<0.05 to p>0.05; happiness_cat remains non-significant, as it already was without this correction.
Comparing the species removed based on Cook’s distance (SA3) to those flagged during IDA by the profile imbalance index (a predictor-only measure, computed before model fitting) shows a modest but above-chance overlap: 4 of the 10 species with the highest profile imbalance were also flagged by Cook’s distance, versus about 0.6 expected by chance. Most influential species, however, are not captured by the profile imbalance index, since influence also depends on the relationship with the outcome, which the index does not use; the two diagnostics are complementary rather than interchangeable.
Removing the species with the most extreme height or BMI values (SA2, SA4) increases the adjusted R² substantially and changes the estimated coefficient for the legendary special-status category, driven by the outsized influence of Cosmoem, a legendary species with an extreme weight-for-height ratio; the estimated height-weight association itself is not affected.
The choice of gender coding (SA5) does not materially affect the results: a finer 7-level coding of percentage_male gives essentially the same AIC and R² as the 3-category gender_cat used in the main model (when the analysis is restricted to species with known gender), and the estimates of all other covariates remain stable between the two.