Additional aspects related to initial data analysis of Pokémon data

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 percetages 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.

Content of the report

In this report we present two aspects that are not addressed in depth in the main IDA report

• summary statistics by type of Pokémon species

• summary statistics relative to the evolution families

These are aspects that can be useful to address in projects that aim to use type as explanatory variable in a regression model (as described in the additional regression report), or that aim at taking into account the possible correlation between species that belong to the same evolutionary family with modelling (for example, by using a mixed effect model or robust standard errors, as described in the sensitivity analysis in the main regression model).

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.

Import data

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

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

Data can be read also online - use the code below.

# read data online 
pokemon <- read.delim("https://www.dropbox.com/scl/fi/ccdidhioa2qqajrg8g133/Pokemon_clean.tsv?rlkey=yiwoou90ict4stsk0vyy6xte0&dl=1", 
                      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 <- arrange(pokemon, pokedex_number)

NOTE:

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.

Initial data analysis

NOTE: In the following of this document we describe data screening steps of initial data analysis (IDA) that focus on type and evolutionary families. The other IDA explorations are reported in the IDA main report.

RESULT: 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 gender is unknown (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: Expert domain knowledge suggests that 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 (with 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.

NOTE: Considerations about missing values in gender are included in the main IDA report, here we focus on type only.

CONSEQUENCE: in our further IDA summaries we summarize if Pokémon species 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.

Complete cases

To evaluate the impact of analyzing only Pokémon with single type (as in a complete case analysis), here we report tables that compare the main characteristics of Pokémon with single and dual type.

Characteristics by single or dual type.
	Dual N=536	Single N=489
PokemonPhase_5cat : Baby	0.01 6⁄536	0.02 11⁄489
Base	0.26 139⁄536	0.42 203⁄489
Phase1	0.35 188⁄536	0.34 166⁄489
Phase2	0.15 78⁄536	0.07 33⁄489
NotEvolving	0.23 125⁄536	0.16 76⁄489
weight_kg	11.0 35.0 82.6 76.7 ± 124.3	7.6 21.4 55.0 56.2 ± 117.0
height_m	0.60 1.20 1.60 1.35 ± 1.36	0.50 0.80 1.30 1.07 ± 1.09
special_status : normal	0.88 472⁄536	0.92 448⁄489
sub-legendary	0.06 30⁄536	0.04 21⁄489
legendary	0.04 21⁄536	0.02 10⁄489
mythical	0.02 13⁄536	0.02 10⁄489
hp	55.0 70.0 88.2 72.4 ± 24.7	50.0 65.0 80.0 67.8 ± 28.5
attack	59.0 80.0 101.0 81.0 ±  29.8	52.0 70.0 95.0 73.7 ±  29.3
defense	55.0 75.0 95.0 77.3 ± 30.4	48.0 63.0 82.0 67.3 ± 27.1
sp_attack	50.0 70.0 95.0 74.6 ± 30.2	43.0 60.0 83.0 65.1 ± 28.3
sp_defense	55.0 71.0 90.0 74.0 ± 26.6	48.0 60.0 81.0 66.0 ± 26.0
speed	50.0 70.0 90.2 70.6 ± 28.1	42.0 60.0 85.0 63.4 ± 28.9
total	355 478 525 450 ± 110	309 409 490 403 ± 111
gender_unknown	0.19 100⁄536	0.11 54⁄489
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: Single type Pokémon are smaller, weaker, more often of base type and with known gender and not with a special status, compared to dual type Pokémon.

CONSEQUENCE: A naive analysis that included only Pokémon with a defined second type would focus on a target population considerably different from the complete population. Note that this is not properly missing data, as it is not due to unavailable information but represents a case where the value is not defined (a single-type Pokémon simply does not have a second type); imputation methods would therefore not be appropriate.

Univariate descriptions

The univariate 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

19 Variables   1025 Observations

---

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

---

NOTE: 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
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

RESULTS - TYPE 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.

Graphical displays

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

Evolutionary families

Pokémon can be classified in evolution families (they evolve in different Pokémon); the information is summarized in the variable PokemonID, which indicates the Pokémon that belong to the same evolution family.

For example, the evolution family of Charmander includes three Pokémon and it is displayed below; the evolution of Pikachu also includes 3 phases, but starting from Baby phase, which is uncommon.

NOTE: Below we show some examples of (famous) evolution families. It is expected that within evolution families the outcome (weight) is associated, therefore in regression modeling it might be sensible to consider this association, taking the correlation within the chain into account, for example using mixed models.

which.pokemon <- which(pokemon$name=="Pikachu")
which.use <- which(pokemon$PokemonID==pokemon$PokemonID[which.pokemon])
which.use <- which.use[c(3,1,2)] # change the order to have the baby as the first

my.urls <- paste0("https://img.pokemondb.net/artwork/large/", tolower(pokemon[which.use,]$name), ".jpg")


#my.tab <- rbind(Image = sprintf("![](%s)", paste0(pokemon$name[which.use], ".jpg")), Name=pokemon$name[which.use], Phase=pokemon$PokemonPhase[which.use])

#kable(my.tab) %>% kable_styling()

my.tab <- rbind(Image = sprintf("![](%s)", my.urls), Name=pokemon$name[which.use], Phase=as.character(pokemon$PokemonPhase[which.use]), 
                Total = pokemon$total[which.use],
                Height =pokemon$height_m[which.use]
                
                )

kable(my.tab) %>% kable_styling()

Image
Name	Pichu	Pikachu	Raichu
Phase	Baby	Base	Phase1
Total	205	320	485
Height	0.3	0.4	0.8

The Pokémon that do not evolve are usually stronger than the Base Pokémon.

which.pokemon <- which(is.element(pokemon$name, c("Arceus","Mew","Mewtwo" )))
which.use <- which.pokemon

my.urls <- paste0("https://img.pokemondb.net/artwork/large/", tolower(pokemon[which.use,]$name), ".jpg")


#my.tab <- rbind(Image = sprintf("![](%s)", paste0(pokemon$name[which.use], ".jpg")), Name=pokemon$name[which.use], Phase=pokemon$PokemonPhase[which.use])

#kable(my.tab) %>% kable_styling()

my.tab <- rbind(Image = sprintf("![](%s)", my.urls), Name=pokemon$name[which.use], Phase=as.character(pokemon$PokemonPhase[which.use]), 
                Total = pokemon$total[which.use],
                Height =pokemon$height_m[which.use]
                
                )

kable(my.tab) %>% kable_styling()

Image
Name	Mewtwo	Mew	Arceus
Phase	NotEvolving	NotEvolving	NotEvolving
Total	680	600	720
Height	2	0.4	3.2

The table shows the number of Pokémon within each evolution family; most evolution families include 2 Pokémon (n=211), followed by 1 (n=201) (Pokémon that do not evolve) and 3 (n=201).

kable(table(table(pokemon$PokemonID)), caption ="Number of Pokémon in the evolution family", col.names=c("Number in the evolution family", "Frequency" )) %>% kable_styling()

Number of Pokémon in the evolution family

Number in the evolution family	Frequency
1	201
2	211
3	121
4	5
5	2
9	1

Few evolution families (8) include more than 3 Pokémon. We display the evolution family of the Pokémon Eevee (Base) that includes 8 Pokémon of Phase 1, besides Eevee.

RESULTS: Evolutionary families are typically small: across the 541 families, the median size is 2 species (mean 1.9), and 37% of families (201 of 541) include only a single species, with no other member to compare it to; about 76% of families include at most 2 species. Only a handful of families are larger (up to 9 species, the Eevee family shown above).

CONSEQUENCE: With most families this small, a random-effects variance for evolutionary family could not be estimated reliably in a mixed-effects model, since the majority of clusters contain too few (or only one) observations to inform the within-family variance component. This supports the use of cluster-robust standard errors, rather than a mixed-effects model or GEE, in the main regression analysis (see SA1 in the main regression report).

db <- cbind.data.frame(pokemon[,c("PokemonPhase_5cat", "name", "generation", "height_m", "weight_kg", "type1", "PokemonID")] ) %>% filter(PokemonID==133) %>% cbind.data.frame(path_images=my.urls) %>% cbind.data.frame(image="")


db %>%
  kbl(booktabs = T, align = "c", digits=2, caption="Evolution family of Eevee") %>%
  kable_styling() %>% 
  kable_paper(full_width = T) %>%
  column_spec(9, image = spec_image(db$path_images, 280, 200) # dimensions of the images
              )  %>% remove_column(c(7,8)) %>% collapse_rows(columns = 1, valign = "top")

Evolution family of Eevee

PokemonPhase_5cat	name	generation	height_m	weight_kg	type1	image
Base	Eevee	1	0.3	6.5	normal
Phase1	Vaporeon	1	1.0	29.0	water
	Jolteon	1	0.8	24.5	electric
	Flareon	1	0.9	25.0	fire
	Espeon	2	0.9	26.5	psychic
	Umbreon	2	1.0	27.0	dark
	Leafeon	4	1.0	25.5	grass
	Glaceon	4	0.8	25.9	ice
	Sylveon	6	1.0	23.5	fairy

Multivariate descriptions

NOTE: In multivariate descriptions we evaluate the association between explanatory variables. Special attention is devoted to the association between type and other variables. Dual vs single type Pokemon were described in the missing value domain.

Association (between explanatory variables and structural variables)

Structural variable is type. Below we report descriptive statistics for each type. Note that the species that belong to two types are included in two groups in the description.

Descriptive statistics of selected variables, by type

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

RESULTS: The Pokémon classified by the type differ by many of the variables. Dragon, steel, rock, ground and fighting types tend to be strong and large; bug, fairy, grass, poison and normal tend to be small and weak. Electric, dragon and fighting are the fastest. Dragon, fighting and psychic have most species with special status, but the least. Dark, dragon and steel have most unhappy species, psychic, electric and normal most happy species. Fairy is the type with most females, while fighting, fire, rock and water are the types with most males; dragon, electric, psychic and steel have most species that are genderless.

NOTE: We examine here the association between type and the other explanatory variables considered in the main regression model (gender, happiness, evolutionary phase, dual type, special status), since a strong association between type and these variables could explain why their estimated coefficients change when type is added to the model (see additional regression report).

Type and height: joint distribution relevant to the planned interaction

NOTE: Type is considered as a moderator of the height-weight association in the additional regression report (type×height interaction). Interactions cannot be evaluated using the outcome during IDA; instead, we examine here the joint distribution of type and height only, following the recommendation that “if interactions of predictors have been pre-specified, IDA may evaluate the joint distribution of these predictors” (Heinze et al., 2024, “Regression without regrets”).

RESULTS: The height boxplot by type shown above indicates that height varies systematically by type (e.g., dragon and steel tend to be tall, bug and fairy tend to be short), confirming that type and height are associated predictors. The number of species per type ranges from 52 (ice) to 154 (water), counting dual-type species under both types; the interquartile range of height per type varies from 0.8 m (bug, dark) to 1.9 m (dragon). No type is so sparse or so narrowly distributed in height as to make its interaction term inestimable, but types with fewer species or a narrower height range (e.g., ice, fairy) will have less precise type-specific height slopes than more numerous or more height-variable types (e.g., water, dragon).

CONSEQUENCE: The association between type and height means that the type×height interaction terms in the additional regression report come with some estimation uncertainty by design, more so for less numerous or less height-variable types. This should be kept in mind when interpreting the interaction coefficients, particularly for types with small sample size.