The weight of Pokémon species - Initial data analysis

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library(dplyr)
library(naniar) #missing values
library(Hmisc)  #data summaries
library(kableExtra) #display tables
library(forcats)  #in tidyverse, to work with categorical variables
library(GGally) #includes ggpairs
library(pheatmap) #to plot heatmaps
library(gt) #for tables that use gt

library(dendextend) #coloring dendrograms
library(viridis) #colors
library(RColorBrewer) #colors

library(lme4) #linear mixed models

library(plotly)
library(reshape2)

library(tidyr)
library(ggrepel) # to add labels to busy graphs
library(scales) # to evaluate percentages in graphs

color <- function(x, color="red") {
  if (knitr::is_latex_output()) {
    sprintf("\\textcolor{%s}{%s}", color, x)
  } else if (knitr::is_html_output()) {
    sprintf("<span style='color: %s;'>%s</span>", color, 
      x)
  } else x
}

#can be used inline to change the color of the text, default is red
#`r color("try red color")`, `r color("try blue color", "blue")`

# specification used in Hmisc to print the tables
mu <- markupSpecs$html  

my.urls <- c("https://img.pokemondb.net/artwork/large/mewtwo.jpg", "https://img.pokemondb.net/artwork/large/eevee.jpg", "https://img.pokemondb.net/artwork/large/pikachu.jpg",
             "https://img.pokemondb.net/artwork/large/hypno.jpg",  "https://img.pokemondb.net/artwork/large/gengar.jpg",
               "https://img.pokemondb.net/artwork/large/psyduck.jpg", "https://img.pokemondb.net/artwork/large/mew.jpg",  
             "https://img.pokemondb.net/artwork/large/charizard.jpg",  "https://img.pokemondb.net/artwork/large/palkia.jpg",
             "https://img.pokemondb.net/artwork/large/arceus.jpg"
             #,  "https://img.pokemondb.net/artwork/large/dialga.jpg"
               )

knitr::include_graphics(my.urls)

Disclaimer: Pokémon images are rendered in this document using references to online resourcesPokémon Database, using URLs in the https://img.pokemondb.net/artwork/ directory and are not stored locally.

Comments on the IDA report

COMMENT: In the following of this document we describe the data screening steps of initial data analysis (IDA). The aim is to explore the data characteristics to gain knowledge about data properties, which enables informed statistical analysis and proper interpretation of the results. The IDA data screening plan was developed based on the checklist proposed by Heinze et al. for cross-sectional studies (link to the paper).

An IDA report typically focuses on presenting missing data, univariable, and multivariable descriptions. The research team would then discuss these IDA findings and decide on further actions. In this report a discussion of individual findings and potential consequences are included.

Boxes with a blue left border (DEVIATION FROM SAP v1) flag points where this report goes beyond, or deviates from, what was planned in the initial version of the statistical analysis plan (SAP v1) — for example, categorizations of variables decided during data screening. These deviations have been incorporated into the updated plan (SAP v2).

To reproduce the content of this file on your computer**

Create a project in R (in RStudio, File -> New project), with two subfolders: Data and Report; save the Rmd version of this document in the Report folder and the data sets describing the properties of the Pokémon and the meta data in the Data folder. Set the working directory to project directory (menu: Session), working within a project makes importing data - working with relative paths - easier.

Importing data

Data are saved in a text tab delimited file and are read them using the read.delim function from the Data folder.

# local reading of data

pokemon <- read.delim("../Data/Pokemon_clean.tsv", 
                      strip.white = TRUE, na.string =c(NA, ""), 
                      blank.lines.skip = TRUE, 
                      stringsAsFactors = FALSE)

# Converting some variables into factors (PokemonPhase)
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"))
#arrange by pokedex number
pokemon <- 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"))

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"))

COMMENT:

Information about data sources and data cleaning and curation can be found in the data dictionary and in the statistical analysis plan. Data cleaning and data management was performed on this dataset.

Research aim

The aim of this analysis project is to study the variability in the weight of Pokémon species.

Objective: to describe how differences in weight of the Pokémon species are associated with differences in height and other variables included in the dataset: evolutionary phase, gender proportion within species, attack/defense characteristics (base stats or their total), special status (mythical, legendary, sub-legendary, or normal), using a linear regression model. Base happiness and catch rate can also be considered as exploratory variables.

Hypothesis: There is a positive linear relationship between the height and weight of Pokémon species, even when various other characteristics of the species are taken into account.

Evolution phase and generation are structural variables for IDA purposes: in the multivariable descriptions, the main characteristics of the Pokémon species will be described stratifying the summaries by the levels of these variables. Gender is examined separately, as a dedicated subgroup analysis (see below), motivated by the missingness finding for Pokémon with unknown gender.

Initial data analysis

The dataset includes 1025 different Pokémon, which are part of 541 different evolution families.

Missing values

Prevalence

RESULT: Most of the variables in the dataset are complete.

The variables with missing values are: percentage_male, where missing values indicate that the species is genderless (which represents a category by itself, as some Pokémon do not have a gender), and second type of Pokémon (variable type2, it is possible that a Pokémon has only one type).

tab1 <- pokemon %>%
  select(PokemonPhase_5cat, generation,  height_m, weight_kg, type1, type2, special_status, 
         hp, attack, defense,   sp_attack, sp_defense, speed, total, catch_rate,    percentage_male,  base_happiness_8    ) %>%  miss_var_summary() %>%
  gt::gt() %>%
    gt::cols_label(
    variable = "Variable",
    n_miss = "Missing (count)",
    pct_miss = "Missing (%)"
  ) %>%
  gt::fmt_number(
    columns = c(pct_miss),
    decimals = 2
  )# Item (per variable) missingness

# to get scrollable gt tables 
tab1 %>% gt::tab_options(container.height = px(200))

Variable	Missing (count)	Missing (%)
type2	489	47.7
percentage_male	154	15.0
PokemonPhase_5cat	0	0
generation	0	0
height_m	0	0
weight_kg	0	0
type1	0	0
special_status	0	0
hp	0	0
attack	0	0
defense	0	0
sp_attack	0	0
sp_defense	0	0
speed	0	0
total	0	0
catch_rate	0	0
base_happiness_8	0	0

COMMENT: Genderless Pokémon are a separate category by definition. Missing values do not indicate that the value is unknown to the researcher, but not defined. See also the univariable explorations for further considerations on the use of the variable percentage of male gender in modeling.

CONSEQUENCE: We define an arbitrary code value for missing gender (-1) that will help summarizing the information about genderless Pokémon (with blank value in male percentage). Summaries will be computed as for factors, as the number of different values among the non missing entries is small.

DEVIATION FROM SAP v1: SAP v1 anticipated that missing values in percentage_male would be categorized and their treatment assessed during IDA (Section 4.1), but did not specify the coding scheme. The choice of an arbitrary placeholder (-1) and the subsequent definition of the gender_cat categorical variable (see Univariable descriptions, below) were determined during this IDA and have been incorporated into SAP v2, Section 6.6.

COMMENT: Type 1 and type 2 should not be interpreted hierarchically. Namely, for a dual type Pokémon, excluding type 2 (rather than type 1) would be a completely arbitrary choice, as the order in which the types are listed does not matter. Our current project does not consider type as a potential explanatory variable. For projects that consider type as a potential explanatory variable, given the large proportion of single type Pokémon (=missing value for type 2), we would suggest that the type could be modeled using a binary variable to represent each type, rather than the two separate types. Also, all observed combinations between type 1 and type 2, regardless of order, could be considered as a new variable, but this would require a very large number of degrees of freedom to be used for modelling.

pokemon$percentage_male[is.na(pokemon$percentage_male)] <- -1

CONSEQUENCE: Additional IDA summaries for Pokémon species that are single or dual type; we also provide some summaries by type, using the is_(type), binary variables that indicate if a Pokémon is of that type. Single-type Pokémon are of one (pure) type, dual-types are of two types (the dual_type variable). The summary statistics for type 1 and type 2 are not reported.

DEVIATION FROM SAP v1: SAP v1 discussed 18 binary type indicators only as a possible modeling choice for projects that use type as an explanatory variable (Section 3, Type). Here type is not an explanatory variable; constructing the is_(type) indicators for descriptive IDA purposes (stratified summaries) was decided during IDA and has been incorporated into SAP v2, Table 1 (Type2 row) and Section 4.2.

Complete cases

To evaluate the impact of analyzing only Pokémon with known gender as would be the case in a complete case analysis), here we report tables that compare the main characteristics of Pokémon with known or unknown gender.

Characteristics by known or unknown gender.
	Known gender N=871	Genderless N=154
PokemonPhase_5cat : Baby	0.02 17⁄871	0.00 0⁄154
Base	0.37 326⁄871	0.10 16⁄154
Phase1	0.39 336⁄871	0.12 18⁄154
Phase2	0.12 105⁄871	0.04 6⁄154
NotEvolving	0.10 87⁄871	0.74 114⁄154
weight_kg	8.45 24.50 56.00 48.08 ±  77.33	20.62 93.60 215.00 173.73 ± 225.54
height_m	0.500 0.900 1.400 1.083 ± 0.971	0.600 1.500 2.375 1.949 ± 2.104
special_status : normal	0.98 856⁄871	0.42 64⁄154
sub-legendary	0.01 10⁄871	0.27 41⁄154
legendary	0.01 5⁄871	0.17 26⁄154
mythical	0.00 0⁄871	0.15 23⁄154
hp	50.0 65.0 80.0 67.5 ±  25.1	67.0 88.5 100.0 85.3 ±  29.8
attack	53.5 71.0 95.0 74.5 ±  28.5	73.2 95.0 120.0 94.6 ±  31.2
defense	50.0 65.0 85.0 68.8 ±  27.8	75.0 95.0 109.5 93.2 ±  29.0
sp_attack	45.0 60.0 85.0 65.6 ±  26.8	70.5 95.0 124.8 95.3 ±  32.5
sp_defense	50.0 65.0 80.0 66.5 ±  24.7	75.0 90.0 107.8 91.0 ±  28.0
speed	44.0 61.0 84.0 63.8 ±  26.8	62.0 90.0 105.8 86.3 ±  31.5
total	315 420 490 407 ± 101	506 570 600 546 ± 105
catch_rate	45.0 75.0 180.0 105.1 ±  74.8	3.0 10.0 45.0 37.9 ±  58.1
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.

RESULT: Compared to others, genderless Pokémon are often not evolving, are considerably larger and stronger, are often mythical, sub-legendary or legendary.

CONSEQUENCE: A naive analysis that only uses the Pokémon with complete data would focus on a target population that is considerably different from the complete population. Note that these are not randomly missing data, since the variable is not defined for this category. Thus, imputation methods would not be appropriate.

Univariable descriptions

The univariable descriptive statistics (at Pokémon level) are reported in a table, the variables are also presented graphically.

For easier data management, we define an additional dataset pokemon_by_type that facilitates the summarization of the Pokémon by types (see code below).

pokemon_by_type <- pokemon %>% pivot_longer(cols = c("type1", "type2"), values_drop_na = TRUE, names_to = NULL, values_to = "type")

Summary statistic in table formats

pokemon.to.describe

34 Variables   1025 Observations

---

weight_kg

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	481	1	66.96	89.59	1.2	3.0	8.5	28.0	70.0	162.3	260.0

lowest : 0.1 0.2 0.3 0.4 0.5 , highest: 820 888 920 950 999.9

---

height_m

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	54	0.997	1.213	0.9738	0.3	0.3	0.5	1.0	1.5	2.0	2.8

lowest : 0.1 0.2 0.3 0.4 0.5 , highest: 8.8 9.2 12 14.5 20

---

PokemonPhase_5cat

n	missing	distinct
1025	0	5

 Value             Baby        Base      Phase1      Phase2 NotEvolving
 Frequency           17         342         354         111         201
 Proportion       0.017       0.334       0.345       0.108       0.196 

---

generation

n	missing	distinct	Info	Mean	Gmd
1025	0	9	0.985	4.749	3.023

 Value          1     2     3     4     5     6     7     8     9
 Frequency    151   100   135   107   156    72    80   104   120
 Proportion 0.147 0.098 0.132 0.104 0.152 0.070 0.078 0.101 0.117 

For the frequency table, variable is rounded to the nearest 0

---

percentage_male

n	missing	distinct	Info	Mean	Gmd
1025	0	8	0.761	46.37	27.24

 Value       -1.0   0.0  12.5  25.0  50.0  75.0  87.5 100.0
 Frequency    154    37     2    25   631    19   131    26
 Proportion 0.150 0.036 0.002 0.024 0.616 0.019 0.128 0.025 

For the frequency table, variable is rounded to the nearest 0

---

special_status

n	missing	distinct
1025	0	4

 Value             normal sub-legendary     legendary      mythical
 Frequency            920            51            31            23
 Proportion         0.898         0.050         0.030         0.022 

---

hp

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	107	0.998	70.18	28.12	36.2	40.0	50.0	68.0	85.0	100.0	111.0

lowest : 1 10 20 25 28 , highest: 190 200 223 250 255

---

attack

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	117	0.999	77.52	33.88	30	40	55	75	100	120	130

lowest : 5 10 15 20 22 , highest: 147 150 160 165 181

---

defense

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	108	0.998	72.51	32.02	35	40	50	70	90	110	125

lowest : 5 10 15 20 23 , highest: 180 184 200 211 230

---

sp_attack

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	113	0.999	70.08	33.37	30.0	35.0	47.0	65.0	90.0	111.6	125.0

lowest : 10 15 20 21 23 , highest: 145 150 151 154 173

---

sp_defense

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	103	0.998	70.21	29.49	31.2	40.0	50.0	67.0	86.0	105.0	120.0

lowest : 20 23 25 30 31 , highest: 142 150 154 200 230

---

speed

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	119	0.999	67.19	32.53	25	30	45	65	88	105	115

lowest : 5 10 13 15 20 , highest: 145 150 151 160 200

---

total

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	207	1	427.7	128.6	250	280	325	450	508	570	600

lowest : 175 180 185 190 194 , highest: 660 670 680 690 720

---

catch_rate

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1025	0	37	0.97	95.01	81.87	3	25	45	60	140	225	255

lowest : 3 5 6 10 15 , highest: 205 220 225 235 255

---

base_happiness_8

n	missing	distinct	Info	Mean	Gmd
1025	0	8	0.459	46.89	12.79

 Value          0    20    35    50    70    90   100   140
 Frequency     80     3    73   835     3     5    15    11
 Proportion 0.078 0.003 0.071 0.815 0.003 0.005 0.015 0.011 

For the frequency table, variable is rounded to the nearest 0

---

dual_type

n	missing	distinct
1025	0	2

 Value        Dual Single
 Frequency     536    489
 Proportion  0.523  0.477 

---

is_water

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.383	154	0.1502	0.2556

---

is_normal

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.334	131	0.1278	0.2232

---

is_grass

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.323	126	0.1229	0.2158

---

is_flying

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.285	109	0.1063	0.1903

---

is_psychic

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.271	103	0.1005	0.181

---

is_bug

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.245	92	0.08976	0.1636

---

is_poison

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.226	84	0.08195	0.1506

---

is_fire

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.221	82	0.08	0.1473

---

is_ground

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.203	75	0.07317	0.1358

---

is_rock

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.201	74	0.0722	0.1341

---

is_dark

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.198	73	0.07122	0.1324

---

is_fighting

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.198	73	0.07122	0.1324

---

is_dragon

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.191	70	0.06829	0.1274

---

is_electric

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.188	69	0.06732	0.1257

---

is_steel

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.178	65	0.06341	0.1189

---

is_ghost

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.178	65	0.06341	0.1189

---

is_fairy

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.176	64	0.06244	0.1172

---

is_ice

n	missing	distinct	Info	Sum	Mean	Gmd
1025	0	2	0.144	52	0.05073	0.09641

---

COMMENT: The is_(type)* variables are binary variables that indicate if a Pokémon is of that type. Single-type Pokémon are of one (pure) type, dual-types are of two types (the dual_type variable. The summary statistics for type 1 and type 2 are not reported.

Additionally, we provide a table that summarizes the descriptive statistics, as it would be usually presented in a paper.

Overall characteristics.
	N=1025
weight_kg	8.5 28.0 70.0 67.0 ± 121.2
height_m	0.50 1.00 1.50 1.21 ± 1.25
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
generation : 1	0.15 151⁄1025
2	0.10 100⁄1025
3	0.13 135⁄1025
4	0.10 107⁄1025
5	0.15 156⁄1025
6	0.07 72⁄1025
7	0.08 80⁄1025
8	0.10 104⁄1025
9	0.12 120⁄1025
special_status : normal	0.90 920⁄1025
sub-legendary	0.05 51⁄1025
legendary	0.03 31⁄1025
mythical	0.02 23⁄1025
hp	50.0 68.0 85.0 70.2 ± 26.6
attack	55.0 75.0 100.0 77.5 ±  29.8
defense	50.0 70.0 90.0 72.5 ± 29.3
sp_attack	47.0 65.0 90.0 70.1 ± 29.7
sp_defense	50.0 67.0 86.0 70.2 ± 26.6
speed	45.0 65.0 88.0 67.2 ± 28.7
total	325 450 508 428 ± 113
catch_rate	45.0 60.0 140.0 95.0 ±  76.4
percentage_male : -1	0.15 154⁄1025
0	0.04 37⁄1025
12.5	0.00 2⁄1025
25	0.02 25⁄1025
50	0.62 631⁄1025
75	0.02 19⁄1025
87.5	0.13 131⁄1025
100	0.03 26⁄1025
base_happiness_8 : 0	0.08 80⁄1025
20	0.00 3⁄1025
35	0.07 73⁄1025
50	0.81 835⁄1025
70	0.00 3⁄1025
90	0.00 5⁄1025
100	0.01 15⁄1025
140	0.01 11⁄1025
is_water	0.15 154⁄1025
is_normal	0.13 131⁄1025
is_grass	0.12 126⁄1025
is_flying	0.11 109⁄1025
is_psychic	0.1 103⁄1025
is_bug	0.09 92⁄1025
is_poison	0.08 84⁄1025
is_fire	0.08 82⁄1025
is_ground	0.07 75⁄1025
is_rock	0.07 74⁄1025
is_dark	0.07 73⁄1025
is_fighting	0.07 73⁄1025
is_dragon	0.07 70⁄1025
is_electric	0.07 69⁄1025
is_steel	0.06 65⁄1025
is_ghost	0.06 65⁄1025
is_fairy	0.06 64⁄1025
is_ice	0.05 52⁄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.

RESULTS - NUMERICAL VARIABLES: Most of the numerical variables have positively asymmetric distributions, most notably height and weight, as it can be observed comparing the median and the mean values, or interpreting the interquartile ranges. This aspect is explored further with graphical displays presented below.

Percentage of male gender and base happiness have few distinct values and might be discretized without substantial loss of information. For gender, this is helpful also for taking into account the problem of missing values in percentage male, as the Pokémon with undefined gender can be grouped in their own category.

The range of the base stats is large, possible range from 1 to 255, but the observed ranges are less than 200. Thus, a unit difference is a small change.

CONSEQUENCES: We define a new variable gender_cat that discretizes the percentage of male gender; the categories are: genderless, mostly females (0-49 percentage males), equal, mostly males (51-100 percentage males). We also define a new variable happiness_cat that discretizes the base happiness; the categories are: Unhappy (<50), Normal (50),Happy (>50). (Note that the for convenience the code to generate these variables is provided below, but the dataset used in this project already includes this variables).

In regression modeling it is sensible to present regression coefficients of base stats that refer to a meaningful difference (for example, in units of 10 or 20).

DEVIATION FROM SAP v1: The categorization of percentage_male into gender_cat (4 levels) and of base_happiness_8 into happiness_cat (3 levels) was not specified in SAP v1; both variables and their categorization schemes were defined during IDA, based on the limited number of distinct values observed. This has been incorporated into SAP v2, Section 6.6, and both are used as categorical explanatory variables in the regression model (Section 4.2).

pokemon$gender_cat <- pokemon$gender_unknown
pokemon$gender_cat[pokemon$gender_unknown==1] = "Genderless"
pokemon$gender_cat[pokemon$percentage_male>=0] = "Mostly females"
pokemon$gender_cat[pokemon$percentage_male==50] = "Equal"
pokemon$gender_cat[pokemon$percentage_male>50] = "Mostly males"

pokemon$gender_cat <- relevel(factor(pokemon$gender_cat), ref="Equal")

pokemon$happiness_cat <- pokemon$base_happiness_8
pokemon$happiness_cat[pokemon$base_happiness_8==50] = "Normal"
pokemon$happiness_cat[pokemon$base_happiness_8<50] = "Unhappy"
pokemon$happiness_cat[pokemon$base_happiness_8>50] = "Happy"
pokemon$happiness_cat <- relevel(factor(pokemon$happiness_cat), ref ="Normal")

We summarize the descriptive statistics with the new categorical variables.

Overall characteristics - new categorized variables.
	N=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
happiness_cat : Happy	0.03 34⁄1025
Normal	0.81 835⁄1025
Unhappy	0.15 156⁄1025

RESULTS - CATEGORICAL VARIABLES Generation 1 and 5 included the most Pokémon, 6 the least. Most Pokémon are Base or Phase 1, baby Pokémon is a rare category, most evolutionary chains include 2 species. There are 18 Pokémon types, water is the most common type among all Pokémon, while the most common type for dual-type Pokémon is flying. Approximately 50% of Pokémon are single-type (type 2 missing) and 50% are dual-type. Most Pokémon species have equal gender representations (62%), and about 15% of Pokémon species are genderless.

Only 31 Pokémon species are legendary, 51 Pokémon species are sub-legendary, and only 23 are mythical.

Some categorical variables have categories with few observed Pokémon, like baby phase in the evolution phase, or happy Pokémon for the happiness variable, or mythical Pokémon for the variable indicating the special status of Pokémon.

CONSEQUENCES: In regression modeling we can use reference categories: Base for evolution phase, normal happiness, equal gender, normal for special status variable.

Graphical displays

The numerical variables are displayed using histograms with small classes (to identify specific features of the data, like digit preference or values that are very common), for some variables also larger classes are used (to display more clearly the shape of the distributions). The number (and proportion) of Pokémon classified by any type is presented with a graph (each Pokémon is classified in 1 or 2 types).

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
ggplot(pokemon_by_type, aes(forcats::fct_infreq(type)))+geom_bar(aes(fill = dual_type)) +coord_flip()+scale_fill_manual(values=cbPalette)+labs(x = "Type", fill = "Single/Dual type") + theme_bw()

RESULTS: The graphical displays confirm that most of the numerical variables have positively asymmetric distributions, most notably height and weight. The variables related to attacks are positively asymmetric, values ending with 0 are very common, numbers multiple of 0.1 are common for height; the total number of points has a bimodal/trimodal distribution.

CONSEQUENCES: It is not sensible to assume that the numerical variables follow a Gaussian distribution, log transformation might be advisable, especially for the outcome variable of weight, and height.

We also plot log2-transformed of height and weight with normal qqplots to assess normality after transformation.

RESULTS: We note that both distributions are significantly closer to normality after log transformation.

CONSEQUENCES: The transformed versions of the variables might be used for modeling.

From the qqplot it is clear that relatively few distinct values of height are observed in the sample.

RESULTS: Extreme values (far from the rest of the distribution) are observed for many variables. There are very few data points for extremely large height or weight; extreme values appear also in the lower and upper part of the distribution for HP, in the upper tails of speed, special defense.

Only 18 Pokémon have a weight greater than 500kg and only 5 Pokémon have a weight greater than 900kg. The 10 most extreme values are shown below, a table reports the details of the 5 with weight above 900 kg.

tail(sort(pokemon$weight_kg), 10)

##  [1] 750.0 800.0 800.0 820.0 888.0 920.0 950.0 950.0 999.9 999.9

db  <-  data.frame(name=pokemon$name[which(pokemon$weight_kg>900)],  weight=pokemon$weight_kg[which(pokemon$weight_kg>900)], height = pokemon$height_m[which(pokemon$weight_kg>900)], special = pokemon$special_status[which(pokemon$weight_kg>900)], image = "")

path_images = paste0("https://img.pokemondb.net/artwork/large/", tolower(pokemon$name[which(pokemon$weight_kg>900)]), ".jpg") 



db %>% arrange(desc(weight)) %>% 
  kbl(booktabs = T, align = "c") %>%
  kable_styling() %>%
  kable_paper(full_width = T) %>%
  column_spec(5, image = spec_image(path_images[order(-db$weight)], 280, 200) # dimensions of the images
              ) 

name	weight	height	special	image
Cosmoem	999.9	0.1	legendary
Celesteela	999.9	9.2	sub-legendary
Groudon	950.0	3.5	legendary
Eternatus	950.0	20.0	legendary
Mudsdale	920.0	2.5	normal

The 10 most extreme values of height are shown below, indicating that the highest values are far from the rest of the distribution.

tail(sort(pokemon$height_m), 10)

##  [1]  5.8  6.2  6.5  7.0  8.8  9.2  9.2 12.0 14.5 20.0

Only 3 Pokémon have more than 10 m of height, 15 more than 5 m, 40 more than 3 m, 101 more than 2 m; only 5% of Pokémon are taller than 2.8 m.

A table reports the details of the 6 with height above 7.5 m.

db  <-  data.frame(name=pokemon$name[which(pokemon$height_m>7.5)],  weight=pokemon$weight_kg[which(pokemon$height_m>7.5)], height = pokemon$height_m[which(pokemon$height_m>7.5)], 
special = pokemon$special_status[which(pokemon$height_m>7.5)], 
image = "")

path_images = paste0("https://img.pokemondb.net/artwork/large/", tolower(pokemon$name[which(pokemon$height_m>7.5)]), ".jpg") 



db %>% arrange(desc(height)) %>%
  kbl(booktabs = T, align = "c") %>%
  kable_styling() %>% 
  kable_paper(full_width = T) %>%
  column_spec(5, image = spec_image(path_images[order(-db$height)], 280, 200) # dimensions of the images
              ) 

name	weight	height	special	image
Eternatus	950.0	20.0	legendary
Wailord	398.0	14.5	normal
Dondozo	220.0	12.0	normal
Steelix	400.0	9.2	normal
Celesteela	999.9	9.2	sub-legendary
Onix	210.0	8.8	normal

RESULTS: These data descriptions indicate that the estimation of the relationship between height and weight for large values of height will depend heavily on few species.

RESULTS: The estimates obtained for large values of height might be very imprecise, due to the extreme values. The estimates for height values larger than 10 m will be based only on 3 data points. Some Pokémon with extreme values in some of the explanatory variables might have large influence on the results.

Additional aspects

See the additional IDA report for additional descriptions of aspect related to evolutionary family and type.

Multivariable descriptions

COMMENT: In multivariable descriptions we evaluate the association between exploratory variables. Special attention is devoted to structural variables (in our example: generation and evolution phase), and to other selected explorations (associations with height, associations among base stats, identification of observations that might be influential).

IDA does not directly explore the association between the exploratory variables and the outcome, which will be addressed during modeling. For this reason we do not present results that explore the association between weight and the exploratory variable (such as scatterplots or summary statistics).

Association (between explanatory variables and structural variables)

Structural variables are generation and evolutionary phase. Below we report descriptive statistics for the explanatory variables, stratified by the values of these variables. Gender is examined separately in a dedicated subgroup analysis (see “Descriptive statistics of selected variables, by gender”, below).

Descriptive statistics of selected variables, by generation

Here we describe the characteristic of the Pokémon by generation.

RESULTS:

Several characteristics vary systematically across generations. Baby Pokémon were included only in generations 2, 3 and 4; the proportion of base Pokémon decreased in later generations, while the proportion of non-evolving Pokémon substantially increased. The most recent generations (7 and later) were characterized by a larger proportion of legendary Pokémon, of non-evolving Pokémon, of genderless Pokémon, and of Pokémon with low happiness. Some attack statistics are higher in generation 4 and later.

CONSEQUENCES: These generational differences overlap with variables already included in the model (evolutionary phase, special status, gender category, happiness category). Generation is therefore unlikely to explain additional variation in the outcome once these covariates are included, and is not used as a separate covariate in the model.

Characteristics by generation.
	1 N=151	2 N=100	3 N=135	4 N=107	5 N=156	6 N=72	7 N=80	8 N=104	9 N=120
PokemonPhase_5cat : Baby	0.00 0⁄151	0.08 8⁄100	0.01 2⁄135	0.07 7⁄107	0.00 0⁄156	0.00 0⁄ 72	0.00 0⁄ 80	0.00 0⁄104	0.00 0⁄120
Base	0.45 68⁄151	0.39 39⁄100	0.33 45⁄135	0.24 26⁄107	0.36 56⁄156	0.36 26⁄ 72	0.28 22⁄ 80	0.28 29⁄104	0.26 31⁄120
Phase1	0.37 56⁄151	0.30 30⁄100	0.33 45⁄135	0.40 43⁄107	0.35 55⁄156	0.38 27⁄ 72	0.28 22⁄ 80	0.38 39⁄104	0.31 37⁄120
Phase2	0.11 16⁄151	0.10 10⁄100	0.11 15⁄135	0.11 12⁄107	0.12 19⁄156	0.11 8⁄ 72	0.11 9⁄ 80	0.11 11⁄104	0.09 11⁄120
NotEvolving	0.07 11⁄151	0.13 13⁄100	0.21 28⁄135	0.18 19⁄107	0.17 26⁄156	0.15 11⁄ 72	0.34 27⁄ 80	0.24 25⁄104	0.34 41⁄120
height_m	0.700 1.000 1.500 1.195 ± 0.963	0.600 0.900 1.400 1.163 ± 1.083	0.600 1.000 1.500 1.230 ± 1.478	0.500 1.000 1.400 1.134 ± 0.893	0.575 1.000 1.400 1.032 ± 0.626	0.475 0.800 1.500 1.085 ± 0.975	0.400 0.850 1.800 1.281 ± 1.378	0.475 1.150 1.925 1.490 ± 2.084	0.500 1.200 1.650 1.354 ± 1.356
special_status : normal	0.97 146⁄151	0.94 94⁄100	0.93 125⁄135	0.87 93⁄107	0.92 143⁄156	0.92 66⁄ 72	0.76 61⁄ 80	0.81 84⁄104	0.90 108⁄120
sub-legendary	0.02 3⁄151	0.03 3⁄100	0.04 5⁄135	0.06 6⁄107	0.04 6⁄156	0.00 0⁄ 72	0.16 13⁄ 80	0.11 11⁄104	0.03 4⁄120
legendary	0.01 1⁄151	0.02 2⁄100	0.02 3⁄135	0.03 3⁄107	0.02 3⁄156	0.04 3⁄ 72	0.06 5⁄ 80	0.04 4⁄104	0.06 7⁄120
mythical	0.01 1⁄151	0.01 1⁄100	0.01 2⁄135	0.05 5⁄107	0.03 4⁄156	0.04 3⁄ 72	0.01 1⁄ 80	0.05 5⁄104	0.01 1⁄120
hp	45.0 60.0 80.0 64.2 ± 28.6	50.0 67.5 86.2 71.0 ± 31.2	50.0 61.0 80.0 65.7 ± 25.2	57.5 70.0 85.0 73.1 ± 24.7	55.0 70.0 80.0 70.3 ± 21.6	52.2 65.5 80.5 68.9 ± 21.7	50.8 68.0 80.0 70.5 ± 28.2	56.0 70.0 90.0 73.1 ± 27.4	55.8 77.0 90.0 77.4 ± 28.2
attack	51.0 70.0 92.0 72.9 ±  26.8	49.8 67.5 85.0 68.3 ±  28.4	50.0 70.0 90.0 73.1 ±  30.4	60.0 80.0 100.0 80.2 ±  30.9	55.0 77.5 100.0 81.0 ±  29.4	52.0 69.0 90.5 72.5 ±  25.6	58.8 75.0 107.8 82.6 ±  32.5	60.0 85.0 112.0 84.6 ±  31.3	60.0 79.5 103.5 82.4 ±  29.3
defense	50.0 65.0 84.0 68.2 ± 26.9	45.0 65.0 85.0 69.7 ± 35.2	46.5 63.0 84.0 69.0 ± 31.1	50.0 70.0 97.5 75.1 ± 30.6	50.0 70.0 85.0 71.2 ± 23.0	53.8 67.5 88.0 75.1 ± 31.3	53.8 75.0 95.0 77.0 ± 30.1	55.0 70.0 95.0 75.0 ± 30.0	60.0 75.0 95.5 76.8 ± 27.8
sp_attack	45.0 65.0 87.5 67.1 ± 28.5	40.0 64.0 85.0 64.5 ± 25.6	48.5 65.0 90.0 67.9 ± 28.3	50.0 69.0 93.0 73.3 ± 31.2	45.0 61.5 90.0 69.2 ± 29.8	50.0 66.0 91.8 72.5 ± 28.0	50.0 65.5 95.0 73.3 ± 33.5	50.0 69.0 90.0 73.2 ± 30.1	50.0 65.0 91.2 72.9 ± 31.5
sp_defense	49.0 65.0 80.0 66.1 ± 24.2	50.0 65.0 86.2 72.3 ± 31.5	50.0 60.0 81.0 66.5 ± 28.5	52.0 70.0 95.0 74.4 ± 27.6	50.0 65.0 80.0 67.3 ± 21.9	54.8 69.5 89.2 74.6 ± 30.4	50.0 72.5 95.0 74.7 ± 29.5	50.0 70.0 85.0 70.0 ± 23.6	55.0 70.0 88.2 72.5 ± 24.6
speed	46.5 70.0 90.0 69.1 ± 27.0	40.0 60.0 85.0 61.4 ± 27.2	44.0 60.0 80.0 61.6 ± 26.9	45.5 70.0 90.5 69.5 ± 27.6	45.0 64.5 90.5 66.6 ± 28.2	47.5 60.0 78.5 65.7 ± 25.9	42.0 60.0 83.2 63.8 ± 28.4	44.8 69.5 88.5 69.9 ± 35.2	50.0 75.0 95.8 75.5 ± 29.4
total	320.0 405.0 490.0 407.6 ±  99.9	319.5 415.0 492.5 407.2 ± 112.5	303.5 410.0 475.0 403.7 ± 115.6	337.0 480.0 525.0 445.6 ± 117.2	322.2 445.5 498.0 425.8 ± 102.5	333.8 450.0 501.8 429.3 ± 111.7	333.5 477.0 530.0 441.9 ± 118.1	333.8 485.0 520.0 445.8 ± 117.5	347.5 488.5 555.0 457.4 ± 116.6
catch_rate	45.0 75.0 190.0 106.2 ±  77.1	45.0 60.0 120.0 91.9 ±  71.7	45.0 90.0 190.0 113.4 ±  83.8	45.0 45.0 120.0 78.9 ±  69.5	45.0 75.0 190.0 103.1 ±  76.6	45.0 62.5 165.0 100.4 ±  72.5	45.0 45.0 120.0 84.5 ±  72.2	45.0 45.0 120.0 86.7 ±  80.8	28.8 45.0 120.0 77.8 ±  70.7
gender_cat : Equal	0.60 91⁄151	0.67 67⁄100	0.70 95⁄135	0.55 59⁄107	0.64 100⁄156	0.64 46⁄ 72	0.54 43⁄ 80	0.57 59⁄104	0.59 71⁄120
Mostly female	0.08 12⁄151	0.08 8⁄100	0.04 6⁄135	0.07 7⁄107	0.06 9⁄156	0.07 5⁄ 72	0.07 6⁄ 80	0.07 7⁄104	0.03 4⁄120
Mostly male	0.23 35⁄151	0.17 17⁄100	0.13 18⁄135	0.22 24⁄107	0.20 31⁄156	0.19 14⁄ 72	0.12 10⁄ 80	0.13 14⁄104	0.11 13⁄120
Genderless	0.09 13⁄151	0.08 8⁄100	0.12 16⁄135	0.16 17⁄107	0.10 16⁄156	0.10 7⁄ 72	0.26 21⁄ 80	0.23 24⁄104	0.27 32⁄120
happiness_cat : Happy	0.03 4⁄151	0.03 3⁄100	0.02 3⁄135	0.11 12⁄107	0.04 7⁄156	0.03 2⁄ 72	0.00 0⁄ 80	0.03 3⁄104	0.00 0⁄120
Normal	0.93 140⁄151	0.83 83⁄100	0.76 103⁄135	0.75 80⁄107	0.83 130⁄156	0.89 64⁄ 72	0.81 65⁄ 80	0.82 85⁄104	0.71 85⁄120
Unhappy	0.05 7⁄151	0.14 14⁄100	0.21 29⁄135	0.14 15⁄107	0.12 19⁄156	0.08 6⁄ 72	0.19 15⁄ 80	0.15 16⁄104	0.29 35⁄120
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.

Descriptive statistics of selected variables, by Pokémon evolution phase

RESULTS: The Pokémon classified by the evolution phases differ by most of the variables. As expected, Pokémon in later evolution phases are stronger and heavier. Pokémon that do not evolve are similar to Phase 2 Pokémon (stronger in some aspects, mostly speed and total points), and they show more variability in most of the variables. The Pokémon that have special status generally do not evolve (evolving special Pokémon were first introduced in generation 7).

CONSEQUENCES: The strong association between the evolution phases and the other variables might indicate that the evolution phase might not be needed as an explanatory variable in the regression model.

Characteristics by Pokémon evolution phase.
	Baby N=17	Base N=342	Phase1 N=354	Phase2 N=111	NotEvolving N=201
factor(generation) : 1	0.00 0⁄ 17	0.20 68⁄342	0.16 56⁄354	0.14 16⁄111	0.05 11⁄201
2	0.47 8⁄ 17	0.11 39⁄342	0.08 30⁄354	0.09 10⁄111	0.06 13⁄201
3	0.12 2⁄ 17	0.13 45⁄342	0.13 45⁄354	0.14 15⁄111	0.14 28⁄201
4	0.41 7⁄ 17	0.08 26⁄342	0.12 43⁄354	0.11 12⁄111	0.09 19⁄201
5	0.00 0⁄ 17	0.16 56⁄342	0.16 55⁄354	0.17 19⁄111	0.13 26⁄201
6	0.00 0⁄ 17	0.08 26⁄342	0.08 27⁄354	0.07 8⁄111	0.05 11⁄201
7	0.00 0⁄ 17	0.06 22⁄342	0.06 22⁄354	0.08 9⁄111	0.13 27⁄201
8	0.00 0⁄ 17	0.08 29⁄342	0.11 39⁄354	0.10 11⁄111	0.12 25⁄201
9	0.00 0⁄ 17	0.09 31⁄342	0.10 37⁄354	0.10 11⁄111	0.20 41⁄201
height_m	0.300 0.500 0.600 0.476 ± 0.225	0.400 0.500 0.700 0.633 ± 0.572	0.800 1.100 1.500 1.326 ± 1.094	1.300 1.600 1.900 1.635 ± 0.563	0.700 1.500 2.000 1.829 ± 2.000
special_status : normal	1.00 17⁄ 17	0.99 337⁄342	0.99 349⁄354	0.98 109⁄111	0.54 108⁄201
sub-legendary	0.00 0⁄ 17	0.01 3⁄342	0.01 3⁄354	0.00 0⁄111	0.22 45⁄201
legendary	0.00 0⁄ 17	0.00 1⁄342	0.00 1⁄354	0.02 2⁄111	0.13 27⁄201
mythical	0.00 0⁄ 17	0.00 1⁄342	0.00 1⁄354	0.00 0⁄111	0.10 21⁄201
hp	40.0 45.0 50.0 55.6 ±  30.9	40.0 49.5 60.0 51.5 ±  21.0	60.0 70.0 85.0 75.6 ±  22.8	75.0 80.0 95.0 85.1 ±  18.0	70.0 82.0 100.0 85.4 ±  26.2
attack	23.0 30.0 40.0 37.4 ±  23.5	44.2 55.0 66.8 56.3 ±  19.7	65.0 80.5 100.0 83.6 ±  25.5	80.0 100.0 118.5 98.0 ±  25.9	75.0 96.0 115.0 95.0 ±  28.3
defense	28.0 37.0 48.0 38.5 ±  21.2	40.0 50.0 60.0 53.0 ±  20.5	60.0 70.0 90.0 77.9 ±  26.2	71.5 85.0 96.5 87.3 ±  22.2	71.0 90.0 105.0 90.9 ±  29.9
sp_attack	35.0 40.0 65.0 45.2 ±  21.5	35.0 45.0 60.0 49.4 ±  18.8	55.0 71.0 95.0 75.4 ±  24.6	70.0 90.0 108.5 88.1 ±  25.6	65.0 85.0 115.0 88.0 ±  33.1
sp_defense	45.0 55.0 65.0 58.7 ±  23.8	40.0 48.0 56.0 50.5 ±  18.3	60.0 70.0 82.8 74.2 ±  21.5	75.0 85.0 95.0 85.7 ±  18.8	70.0 89.0 101.0 89.2 ±  28.4
speed	20.0 35.0 60.0 40.4 ±  26.5	35.0 50.0 65.0 51.3 ±  20.6	50.0 67.5 90.0 69.8 ±  27.3	60.0 79.0 96.0 78.1 ±  24.5	70.0 90.0 103.0 86.0 ±  29.5
total	218.0 280.0 310.0 275.8 ±  63.2	280.0 306.5 330.0 311.9 ±  60.5	415.8 474.0 500.0 456.5 ±  62.6	500.0 525.0 534.0 522.3 ±  47.6	475.0 570.0 590.0 534.5 ±  87.3
catch_rate	50.00 130.00 170.00 127.35 ±  72.11	71.25 190.00 232.50 161.60 ±  80.25	45.00 60.00 90.00 70.57 ±  32.30	45.00 45.00 45.00 43.50 ±  8.65	3.00 30.00 50.00 50.46 ±  66.18
gender_cat : Equal	0.41 7⁄ 17	0.72 245⁄342	0.70 249⁄354	0.59 66⁄111	0.32 64⁄201
Mostly female	0.29 5⁄ 17	0.06 20⁄342	0.06 23⁄354	0.05 6⁄111	0.05 10⁄201
Mostly male	0.29 5⁄ 17	0.18 61⁄342	0.18 64⁄354	0.30 33⁄111	0.06 13⁄201
Genderless	0.00 0⁄ 17	0.05 16⁄342	0.05 18⁄354	0.05 6⁄111	0.57 114⁄201
happiness_cat : Happy	0.12 2⁄ 17	0.01 4⁄342	0.02 6⁄354	0.00 0⁄111	0.11 22⁄201
Normal	0.88 15⁄ 17	0.91 312⁄342	0.90 318⁄354	0.88 98⁄111	0.46 92⁄201
Unhappy	0.00 0⁄ 17	0.08 26⁄342	0.08 30⁄354	0.12 13⁄111	0.43 87⁄201
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.

# 1. Evolutionary Phase vs Gender Category
ggplot(pokemon, aes(x = factor(PokemonPhase_5cat), fill = gender_cat)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = percent_format()) +
  labs(
    x = "Evolutionary Phase",
    y = "Proportion",
    fill = "Gender Category",
    title = "Proportion of Gender Category by Evolutionary Phase"
  ) +
  theme_minimal()

# 2. Evolutionary Phase vs Happiness Category
ggplot(pokemon, aes(x = factor(PokemonPhase_5cat), fill = happiness_cat)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = percent_format()) +
  labs(
    x = "Evolutionary Phase",
    y = "Proportion",
    fill = "Happiness Category",
    title = "Proportion of Happiness Category by Evolutionary Phase"
  ) +
  theme_minimal()

# 4. Evolutionary Phase vs Dual Type
ggplot(pokemon, aes(x = factor(PokemonPhase_5cat), fill = as.factor(dual_type))) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = percent_format()) +
  labs(
    x = "Evolutionary Phase",
    y = "Proportion",
    fill = "Dual Type",
    title = "Proportion of Single vs Dual-Type Pokémon by Evolutionary Phase"
  ) +
  theme_minimal()

# 5. Evolutionary Phase vs Special Status
ggplot(pokemon, aes(x = factor(PokemonPhase_5cat), fill = special_status)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = percent_format()) +
  labs(
    x = "Evolutionary Phase",
    y = "Proportion",
    fill = "Special Status",
    title = "Proportion of Pokémon by Special Status and Evolutionary Phase"
  ) +
  theme_minimal()

Descriptive statistics of selected variables, by gender

Here we describe the characteristic of the Pokémon by gender.

COMMENT: A subgroup analysis is planned by gender category.

RESULT: The genderless category differs most from the other categories. The mostly female category are the shortest species.

Characteristics by gender.
	Equal N=631	Mostly female N=64	Mostly male N=176	Genderless N=154
PokemonPhase_5cat : Baby	0.01 7⁄631	0.08 5⁄ 64	0.03 5⁄176	0.00 0⁄154
Base	0.39 245⁄631	0.31 20⁄ 64	0.35 61⁄176	0.10 16⁄154
Phase1	0.39 249⁄631	0.36 23⁄ 64	0.36 64⁄176	0.12 18⁄154
Phase2	0.10 66⁄631	0.09 6⁄ 64	0.19 33⁄176	0.04 6⁄154
NotEvolving	0.10 64⁄631	0.16 10⁄ 64	0.07 13⁄176	0.74 114⁄154
height_m	0.500 0.900 1.400 1.108 ± 1.089	0.500 0.650 1.200 0.833 ± 0.481	0.600 1.000 1.500 1.082 ± 0.553	0.600 1.500 2.375 1.949 ± 2.104
special_status : normal	1.00 629⁄631	0.94 60⁄ 64	0.95 167⁄176	0.42 64⁄154
sub-legendary	0.00 1⁄631	0.05 3⁄ 64	0.03 6⁄176	0.27 41⁄154
legendary	0.00 1⁄631	0.02 1⁄ 64	0.02 3⁄176	0.17 26⁄154
mythical	0.00 0⁄631	0.00 0⁄ 64	0.00 0⁄176	0.15 23⁄154
hp	50.0 65.0 80.0 66.4 ±  24.1	53.5 69.0 80.0 74.7 ±  38.3	54.0 66.5 80.0 68.9 ±  22.0	67.0 88.5 100.0 85.3 ±  29.8
attack	53.0 71.0 95.0 74.2 ±  28.4	40.8 55.0 75.2 58.5 ±  26.6	60.0 80.0 102.5 81.3 ±  27.3	73.2 95.0 120.0 94.6 ±  31.2
defense	50.0 65.0 85.0 69.8 ±  29.0	45.0 62.0 75.0 60.6 ±  25.4	50.0 65.0 80.0 68.5 ±  23.6	75.0 95.0 109.5 93.2 ±  29.0
sp_attack	44.5 60.0 80.0 63.3 ±  26.1	45.0 64.5 82.8 68.3 ±  29.1	53.0 66.5 95.0 73.1 ±  27.0	70.5 95.0 124.8 95.3 ±  32.5
sp_defense	45.5 61.0 80.0 64.8 ±  24.6	55.0 70.0 96.5 75.8 ±  28.6	51.5 65.0 85.0 69.4 ±  22.2	75.0 90.0 107.8 91.0 ±  28.0
speed	40.0 60.0 80.0 62.1 ±  26.9	40.8 63.0 90.0 65.2 ±  28.9	50.0 65.0 85.0 69.5 ±  25.1	62.0 90.0 105.8 86.3 ±  31.5
total	309.5 420.0 485.0 400.6 ± 100.6	300.0 427.0 490.0 403.0 ± 109.1	350.0 420.0 525.0 430.7 ±  94.4	505.8 570.0 600.0 545.6 ± 105.2
catch_rate	45.0 90.0 190.0 117.9 ±  76.5	45.0 75.0 175.0 106.5 ±  74.6	45.0 45.0 45.0 58.6 ±  44.4	3.0 10.0 45.0 37.9 ±  58.1
happiness_cat : Happy	0.01 9⁄631	0.12 8⁄ 64	0.03 5⁄176	0.08 12⁄154
Normal	0.90 570⁄631	0.84 54⁄ 64	0.93 164⁄176	0.31 47⁄154
Unhappy	0.08 52⁄631	0.03 2⁄ 64	0.04 7⁄176	0.62 95⁄154
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.

Association of selected variables with height

RESULTS: Most of the attack/defense variables are strongly and positively associated to height in the middle ranges of the values of height, while the association sometimes is not present for the most extreme values of height (height is log-2 transformed). The estimates are rather imprecise for values of height above 4 m, as few Pokémon are as high.

CONSEQUENCES: The non-linear association between base stats and (log2) height might anticipate non-linear associations also with the outcome variable.

Correlation between explanatory variables

We use the generalized pairs plot to display the associations (summarized also with Pearson’s correlation) among numerical variables.

COMMENT: Total is the sum of hp, attack, defense, sp_defense, sp_attack and speed; therefore, a strong positive correlation of Total with these variables is expected. In regression modeling a choice must be made between the use of total or of the other six base stats.

RESULT: Most of the numerical variables were positively correlated; speed was weakly correlated to defense and height, and in general it was the least correlated variable to the other base stats and to total. Height was very strongly correlated to total, and generally positively associated to the other base stats (but the correlation was weaker than with total).

Spearman’s correlations among variables are graphically displayed using a heatmap. We add some binary variables to the list of the previously explored numerical explanatory variables.

#GGally::ggpairs(pokemon[,c(11:18,26)])

my.data <- cbind.data.frame(log2height= log2(pokemon$height_m),  pokemon %>% select("total", "hp", "attack", "defense", "speed", "sp_attack", "sp_defense"))
# Create plot with larger correlation text in upper panel
pair_plot <- GGally::ggpairs(
  my.data,
  aes(alpha = 0.1),
  lower = list(continuous = wrap("smooth", method = "loess", size = 0.3)),
  upper = list(continuous = wrap("cor", size = 5)),  # <- Increased from 3 to 5
  diag = list(continuous = wrap("barDiag", colour = "black", fill = "grey70"))
) #+my_theme

pair_plot

my.data <- cbind.data.frame(log2height= log2(pokemon$height_m),  pokemon %>% select("total", "hp", "attack", "defense", "speed", "sp_attack", "sp_defense",  "catch_rate" ))



#setting the breaks for the heatmaps, to make colors comparable
Breaks <- seq(-1, 1, by=.01)
# selecting a color palette, setting 1 to gray
my_palette <- c( colorRampPalette( rev(RColorBrewer::brewer.pal(n = 7, name ="RdYlBu")))(length(Breaks)-2) , "grey80", "grey80")


pheatmap(cor(my.data, use="p", method="s"), cluster_cols=TRUE, cluster_rows=TRUE, display_numbers = TRUE, main="Correlation between variables", breaks=Breaks, 
         #labels_row = paste("Wave", 1:7),  labels_col = paste("Wave", 1:7), 
         color= my_palette)

Below we show graphically the association between total and the variables that constitute them, using different colors for the mythical, legendary/sub-legendary and normal groups. We also highlight the names of the Pokémon with unusual combinations of the values variables, which might be multivariable outliers.

RESULTS: The base stats variables are generally strongly associated to total, but not necessarily over the whole range of their values (total does not increase with HP, special defense and defense for their large values, some Pokémon can have large/low values of one variable but not necessarily also in the total). The base stats are less associated to height compared to total. The six base stats seem to carry additional information compared to Total, speed has distinct features from others, for example speed and defense have an extremely low correlation.

Catch rate and happiness (to a lesser extent) were negatively associated to the other variables (as higher values are generally associated to stronger Pokémon).

From the graphs it can also be seen that legendary Pokémon have consistently large values of total (only few legendary Pokémon have low values, Cosmoem, Cosmog and Terapagos) and that their values in the base stats is rather homogeneous, while there is more heterogeneity for mythical and sub-legendary Pokémon.

CONSEQUENCES: The association between the base stats and total is approximated fairly by a linear relationship for attack, speed and special attack, while the association is not visible for large values of HP, defense and special defense; some Pokémon with very large values of these variables do not have also a very large total.

In regression model the use of the six base stats might be more informative than the use of total, as the six base stats provide to some extents some heterogeneous information about the Pokémon strengths. It is however possible that the correlation between the base stats might inflate the standard errors of the estimated coefficient, or cause some instability in the model, or that some of the variables will not contribute to improve prediction (for example special attack and special defense, see also the graphs in the next section).

The very large (negative) correlation between catch rate and the base stats indicates that the use of this variable in modeling might not contribute much to prediction.

pokemon$log2height <- log2(pokemon$height_m)
which.variables <- c("log2height", "total", "hp", "attack", "defense", "speed", "sp_attack", "sp_defense", "catch_rate")

my.mat <- cor(pokemon[, which.variables], use="p", method="s")
diag(my.mat) <- apply(pokemon[, which.variables], 2, sd, na.rm=TRUE)
my.mat[lower.tri(my.mat)] <- cov(pokemon[, which.variables], use="p")[lower.tri(cov(pokemon[, which.variables], use="p"))]

kable(my.mat, digits=2, caption="Correlation (upper triangle, Spearman) / SD (diagonal) / covariance (lower triangle) between variables") %>%
  kable_styling()

Correlation (upper triangle, Spearman) / SD (diagonal) / covariance (lower triangle) between variables

	log2height	total	hp	attack	defense	speed	sp_attack	sp_defense	catch_rate
log2height	1.07	0.72	0.66	0.62	0.51	0.33	0.44	0.48	-0.53
total	83.19	112.77	0.76	0.71	0.69	0.54	0.70	0.73	-0.73
hp	16.70	2003.52	26.63	0.59	0.49	0.27	0.47	0.50	-0.54
attack	18.32	2399.02	377.65	29.77	0.53	0.36	0.30	0.30	-0.52
defense	13.91	2078.94	233.21	405.97	29.29	0.08	0.30	0.58	-0.49
speed	9.54	1793.45	137.56	300.16	6.77	28.72	0.42	0.27	-0.41
sp_attack	13.22	2345.08	284.04	249.11	182.50	360.40	29.66	0.56	-0.53
sp_defense	11.49	2097.24	261.84	179.66	392.76	163.88	389.41	26.64	-0.54
catch_rate	-42.17	-6298.50	-997.54	-1177.22	-1024.29	-865.52	-1182.86	-1051.06	76.36

Identification of species with unusual combinations of base stats

To identify Pokémon species with unusual combinations of base stats in a systematic and reproducible way, we define a profile imbalance index: each of the six base stats is standardized (z-scores, computed across all 1,025 species), and for each Pokémon species we compute the standard deviation of its own six standardized values. A low value indicates a balanced profile (similar standing across all six stats relative to other species); a high value indicates an imbalanced profile, combining an extreme value in one or more stats with comparatively low values in others. We identify the ten species with the highest profile imbalance index.

Pokémon species with the most unbalanced combination of base stats

name	hp	attack	defense	sp_attack	sp_defense	speed	profile_imbalance
Shuckle	20	10	230	10	230	5	4.02
Blissey	255	10	10	75	135	55	3.47
Chansey	250	5	5	35	105	50	3.46
Guzzlord	223	101	53	97	53	43	2.50
Stakataka	61	131	211	53	101	13	2.32
Wobbuffet	190	33	58	33	58	33	2.28
Regice	80	50	100	100	200	50	2.08
Regidrago	200	100	50	100	50	80	2.08
Steelix	75	85	200	55	65	30	1.99
Regieleki	80	100	50	100	50	200	1.99

RESULTS: We defined a profile imbalance index for each species — the standard deviation of its six standardized (z-score) base stats (hp, attack, defense, sp_attack, sp_defense, speed) — to systematically identify species with unusual, imbalanced combinations of base-stat values, rather than relying on ad hoc visual inspection of scatterplots. The 10 species with the highest profile imbalance index are Shuckle, Blissey, Chansey, Guzzlord, Stakataka, Wobbuffet, Regice, Regidrago, Steelix and Regieleki. These species combine extremely high values in some base stats with extremely low values in others (for example, Shuckle has very high defense and special defense but very low attack, special attack and speed; Blissey and Chansey have very high hp but low attack and defense).

CONSEQUENCES: These species are candidates for having large influence in the regression model. Their actual influence will be assessed with model diagnostics (Cook’s distance, leverage) during regression modeling; the species flagged as influential will be compared to those already identified by the profile imbalance index, to assess whether influential points could have been anticipated from the predictor distributions alone, before model fitting.

DEVIATION FROM SAP v1: SAP v1 specified only that outliers and influential points would be examined using residual plots and Cook’s distance during regression modeling (Section 4.2) — an outcome-informed, post-fitting approach. During IDA we additionally defined a profile imbalance index, an outcome-blind, predictor-only screening tool, to flag species with unusual base-stat combinations before any model is fitted. This addition, which replaces an earlier, less systematic description of individual atypical species, has been incorporated into SAP v2, Sections 5.3 (SA3) and 7.2.2.

Visualization of the results

z_df <- as.data.frame(z_stats)
z_df$name <- pokemon$name

z_long <- z_df %>%
  filter(name %in% top_imbalance$name) %>%
  pivot_longer(cols = all_of(base_stats), names_to = "stat", values_to = "z")

ggplot(z_long, aes(x = stat, y = reorder(name, -match(name, top_imbalance$name)), fill = z)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +
  labs(x = "Base stat", y = "Pokémon species", fill = "z-score",
       title = "Standardized base stats for the most unbalanced profiles") +
  theme_bw()

Exploration of the special status of Pokémon

COMMENT: Here we explore the characteristics of the different types of special Pokémon. The aim is to assess if the groups of mythical, legendary and sub-legendary are similar in their characteristics, and if it is sensible to use a single category for legendary and sub-legendary.

RESULTS: The three special types differ among each other in several aspects. In terms of height, mythical are similar to normal and are are considerably shorter than sub-legendary, which are shorter than legendary. Mythical have very homogeneous values of total (almost all have value 600), legendary have larger values, sub-legendary smaller (and rather homogeneous); the difference between the base stats is less pronounced.

CONSEQUENCE: The differences observed between the three groups seem to indicate that it is sensible to maintain the three special groups separate in modeling, if the special status is used as an explanatory variable.

The strong separation of subgroups (especially legendary and sub-legendary) in terms of height might be problematic in modeling, as it is possible that the inclusion of this variable in modeling might inflate the standard errors of the height-related regression coefficient; if the groups are associated to weight beyond the differences due to height, it is also possible that some spurious associations between the groups and the outcomes will be estimated.

Also, the interpretation of the regression coefficients from the multivariable model might be problematic (as they refer to comparisons among species of same height but with different status).

Exploration of the happiness of Pokémon

COMMENT: Here we explore the characteristics of the different happiness categories of Pokémon. The aim is to assess how the groups of unhappy / normal and happy Pokémon differ to gain insight about the potential usefulness of the variable as an explanatory variable.

RESULTS: The three happiness groups differ among each other in several aspects, it might be unexpected that the happy and unhappy group are similar, and differ from the normal group in most characteristics. Normal types are weaker than happy/unhappy in all aspects, however the height of normal and happy is similar, and unhappy are higher.

CONSEQUENCE: Happiness group might contribute some additional information compared to the other explanatory variables. These explorations indicate also that a non-linear relationship might be expected between the outcome and the numerical variable and, due also to the limited number of distinct values, the categorization of the variable might prove beneficial.

Clustering of the Pokémon based on the explanatory variables

The clustering of the Pokémon (based on attack only variables and total points) and a graphical summary of their characteristics is given in the dendrogram below. Shades of red indicate three different groups. For Phases, yellow is for base, orange Phase 1, red Phase 2, Black not evolving, pink for baby Pokémon. For the special categories, blue indicates mythical, red legendary and orange sub-legendary Pokémon. For happiness: white is normal, blue is happy, red is unhappy.

# new version of the graph

  clust.compl <- hclust(dist(pokemon[,12:18]), method="complete")
#  table(cutree(clust.compl, h=150))

  nodePar <- list(lab.cex = 0.6, pch = c(NA, 5), 
                   cex = 0.7, col = "blue")
  
  #plot(dend1, nodePar = nodePar, leaflab = "none")
  
  
  par(mfrow=c(1,1))
  par(mar = c(12,4,1,1))
  
  #save the results in the format dealt with by library library(dendextend)
  dend <- as.dendrogram(hclust(dist(pokemon[,12:18]), method="complete"), nodePar = nodePar, leaflab = "none")
  
  dend1 <- color_branches(dend, h=150)
  
  #library(magrittr)
  #dend1 <- dend1 %>% set("labels_cex", .25) 
  
  
  #divide the values of some of the variables in 3 groups for color coding of the values for the individuals
  #cols3 <- diverging_hcl(3, h = c(240, 0))
  #cols3 <- inferno(3)
  
  cols3 <- brewer.pal(3, name="Reds")
  
  colhp <- (cols3[cut(pokemon$hp, breaks=quantile(pokemon$hp, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colspeed <- (cols3[cut(pokemon$speed, breaks=quantile(pokemon$speed, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colattack <- (cols3[cut(pokemon$attack, breaks=quantile(pokemon$attack, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colspatt <- (cols3[cut(pokemon$sp_attack, breaks=quantile(pokemon$sp_attack, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colspdef <- (cols3[cut(pokemon$sp_defense, breaks=quantile(pokemon$sp_defense, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colstotal <- (cols3[cut(pokemon$total, breaks=quantile(pokemon$total, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colsdef <- (cols3[cut(pokemon$defense, breaks=quantile(pokemon$defense, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colsweight <- (cols3[cut(pokemon$weight_kg, breaks=quantile(pokemon$weight_kg, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colsheight <- (cols3[cut(pokemon$height_m, breaks=quantile(pokemon$height_m, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
  colsev <- c("white", "yellow", "orange", "red", "black")[pokemon$PokemonPhase_5cat] 
  colleg <- c("white", "yellow", "orange", "red")[factor(pokemon$special_status)] 
  colUnknowGender <- c("white", "black")[factor(pokemon$gender_unknown)] 
  colmit <- c("white", "red", "blue", "orange")[factor(pokemon$special_status)] #white: normal, red: legendary, blue" mythical, orange: sub-legendary
  colshappy <-  c("white", "blue", "red")[factor(pokemon$happiness_cat)]
  colcapt <- (cols3[cut(pokemon$catch_rate, breaks=quantile(pokemon$catch_rate, c(0, 0.33, 0.66, 1)), include.lowest=TRUE)])
 # coltype <- factor(pokemon$typeCards1) 
 
  
  plot(dend1, nodePar = nodePar, leaflab = "none")
  
  colored_bars(cbind(colsev, colUnknowGender, colmit, #coltype, 
                     #pokemon$type1, 
                     colsheight, colshappy, colcapt, colhp,colsdef, colattack, colspeed, colspatt,  colspdef, colstotal), dend1, 
               rowLabels=c("Phase", "Genderless", "Mythical/(sub/)Legendary", #"Type",
                           #"Type1", 
                           "Height", "Happiness", "Catch rate", "HP", "Defense", "Attack", "Speed", "Special attack", "Special defense", 
                                                                     "Total"))

RESULTS: This alternative graphical display confirms the findings from previous data explorations. For example, one group includes most of the legendary/mythical Pokémon that tend to cluster together, most of genderless Pokémon, and exhibits large values of the base stats; this group includes a large proportion of happy and unhappy Pokémon. Roughly, the other three groups are mostly expressing the differences in evolution phases. The association of base stats with their total can be also observed. For example, large values of speed can be observed also for Pokémon that have low overall base stats.

IDA Result: Analysis-ready dataset

Save out analysis-ready dataset: here we save a copy of the newly defined data sets, where some additional variables were defined based on the results from IDA.

saveRDS(pokemon,"../Data/Pokemon_afterIDA.Rdata")