---
title: "The weight of Pokémon species - Initial data analysis"
author: "Lara Lusa, Sam Manski, Marianne Huebner"
date: "July 2026"
output:
  html_document:
    code_folding: "hide"
    toc: true
    toc_float: TRUE
---

```{r setup}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```


```{css, echo=FALSE}
.infobox {
  padding: 1em 1em 1em 4em;
  margin-bottom: 10px;
  border: 2px solid orange;
  border-radius: 10px;
  background: #d3d3d3 5px center/3em no-repeat;
  background-image: url("https://img.pokemondb.net/artwork/large/psyduck.jpg");
  }

.caution {
  background-image: url("https://img.pokemondb.net/artwork/large/psyduck.jpg");
}

.resultbox {
  border-color: #2f855a !important;
  background-image: url("https://img.pokemondb.net/artwork/large/alakazam.jpg") !important;
}

.consequencebox {
  border-color: #6b46c1 !important;
  background-image: url("https://img.pokemondb.net/artwork/large/hypno.jpg") !important;
}

.sapdev {
  padding: 1em 1em 1em 4em;
  margin-bottom: 10px;
  border: 2px solid #2b6cb0;
  border-radius: 10px;
  background: #eaf2fb 5px center/3em no-repeat;
  background-image: url("https://img.pokemondb.net/artwork/large/eevee.jpg");
}
```


<script>
document.addEventListener("DOMContentLoaded", function() {
  document.querySelectorAll('.infobox').forEach(function(box) {
    var txt = box.textContent.trim();
    if (/^RESULT/i.test(txt)) { box.classList.add('resultbox'); }
    else if (/^CONSEQUENCE/i.test(txt)) { box.classList.add('consequencebox'); }
  });
});
</script>

```{r packages, message=FALSE, warning=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

```




```{r extrafunctions0}
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  

```


```{r, out.width="150px"}
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 resources[Pokémon Database](https://pokemondb.net), using URLs in the `https://img.pokemondb.net/artwork/` directory and are not stored locally.


## Comments on the IDA report

:::: {.infobox data-latex=""}

**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](https://www.researchsquare.com/article/rs-3580334/v1)). 

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. 


```{r, eval=TRUE}
# local reading of data

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


```{r}

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

```{r, eval = FALSE}
# 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"))


```



:::: {.infobox data-latex=""}

**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



:::: {.infobox data-latex=""}

The dataset includes `r length(unique(pokemon$pokedex_number))` different Pokémon, which are part of `r length(unique(pokemon$PokemonID))` different evolution families. 

::::


### Missing values

#### Prevalence

:::: {.infobox data-latex=""}

**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).


```{r ItemMissingness}

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

```



:::: {.infobox data-latex=""}

**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.
::::


:::: {.infobox data-latex=""}

**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. 
::::


:::: {.sapdev data-latex=""}

**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.
::::


:::: {.infobox data-latex=""}

**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. 
::::


```{r chImputeMissing}
pokemon$percentage_male[is.na(pokemon$percentage_male)] <- -1
```



:::: {.infobox data-latex=""}

**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.
::::

:::: {.sapdev data-latex=""}

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

```{r chCompleteCases1, results="asis", message=FALSE, warning =FALSE ,eval=FALSE, echo=FALSE}


s <-
     Hmisc::summaryM(
         PokemonPhase_5cat+  weight_kg + height_m+  special_status + 
         hp + attack+defense+  sp_attack+ sp_defense+ speed+ total+ gender_unknown ~ factor(dual_type),
         overall = FALSE,
         test = FALSE, 
         data = pokemon
     )


Hmisc::html(
    s,
    caption = 'Characteristics by single or dual type',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)


#:::: {.infobox data-latex=""}

#**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. 
#::::

```





```{r chCompleteCases2, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


s <-
     Hmisc::summaryM(
         PokemonPhase_5cat+  weight_kg + height_m+  special_status + 
         hp+ attack+defense+  sp_attack+ sp_defense+ speed+ total+ catch_rate~ factor(gender_unknown, labels=c("Known gender", "Genderless")),
         data = pokemon,
         overall = FALSE,
         test = FALSE 
         
     )


Hmisc::html(
    s,
    caption = 'Characteristics by known or unknown gender',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```



:::: {.infobox data-latex=""}
**RESULT:**
Compared to others, genderless Pokémon are often not evolving, are considerably larger and stronger, are often mythical, sub-legendary or legendary.
::::



:::: {.infobox data-latex=""}
**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).


```{r, eval = TRUE}
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

```{r UnivariableTable2, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


#pokemon.to.describe <- pokemon[,-c(1:3, 6:7, 8:9)]
pokemon.to.describe <- pokemon %>% select( "weight_kg", "height_m",   "PokemonPhase_5cat", "generation",  "percentage_male",   "special_status" ,      "hp"  ,  "attack", "defense", "sp_attack",         "sp_defense",    "speed" ,            "total" ,            "catch_rate",        "base_happiness_8" , "dual_type", (starts_with("is_") & !starts_with("is_spec"))  )       


# re-ordering the is_type columns by frequency

#retrieve the names
names.is_types <- pokemon %>% select(starts_with("is_") & !starts_with("is_spec")  ) %>% names()

#retrieve the absolute frequencies
freq.is_types <- pokemon %>% select(starts_with("is_") & !starts_with("is_spec")  ) %>% colSums()

#re-order the dataset, works as the is_type variables are at the end of the selected variables
pokemon.to.describe <- pokemon.to.describe[,c(
  (1:(ncol(pokemon.to.describe)-length(freq.is_types))), ncol(pokemon.to.describe)-length(freq.is_types)+order(-freq.is_types)) ]


s <-      Hmisc::describe(pokemon.to.describe)
#       height_m + PokemonPhase_5cat +  generation + weight_kg + type1 + type2 + is_mythical + is_legendary + 
        # hp+ attack+defense+  sp_attack+ sp_defense+ speed+ total+ catch_rate+    percentage_male+ base_happiness_8+  IQ_group        )


Hmisc::html(
    s,
    caption = 'Overall characteristics',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```

:::: {.infobox data-latex=""}
**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. 


```{r UnivariableTable, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


s <-      Hmisc::summaryM(
       weight_kg  + height_m + PokemonPhase_5cat +  generation + special_status + 
         hp+ attack+defense+  sp_attack+ sp_defense+ speed+ total+ catch_rate+    percentage_male+ base_happiness_8 + 
is_water+is_normal+is_grass+is_flying+is_psychic+is_bug+is_poison+is_fire+is_ground+is_rock+is_dark+is_fighting+is_dragon+is_electric+is_steel+ is_ghost + is_fairy + +is_ice

~ 1,
         data =  pokemon %>% mutate(percentage_male=factor(percentage_male)),
         overall = FALSE,
         test = FALSE
     )


Hmisc::html(
    s,
    caption = 'Overall characteristics',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```

:::: {.infobox data-latex=""}
**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.
::::



:::: {.infobox data-latex=""}
**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).
::::


:::: {.sapdev data-latex=""}

**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).
::::




```{r, eval = FALSE}
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")

```

```{r, eval = FALSE}
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.


```{r UnivariableTableNew, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


s <-      Hmisc::summaryM(
     gender_cat + happiness_cat   ~ 1,
         data =  pokemon,
         overall = FALSE,
         test = FALSE
     )


Hmisc::html(
    s,
    caption = 'Overall characteristics -  new categorized variables',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```


:::: {.infobox data-latex=""}
**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 `r sum(pokemon$special_status=="legendary")` Pokémon species are legendary, 
`r sum(pokemon$special_status=="sub-legendary")` Pokémon species are sub-legendary, 
and only `r sum(pokemon$special_status=="mythical")` 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.  

::::


:::: {.infobox data-latex=""}
**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).


```{r UnivariableGraphs, out.width="45%", echo=FALSE, warning=FALSE, message = FALSE}


ggplot(pokemon, aes(weight_kg)) + geom_histogram(binwidth = 5, color="black", fill="white") + labs(x="Weight (kg)") + theme_bw()
ggplot(pokemon, aes(weight_kg)) + geom_histogram(binwidth = 20, color="black", fill="white") + labs(x="Weight (kg)") + theme_bw()


ggplot(pokemon, aes(height_m)) + geom_histogram(binwidth = .1, color="black", fill="white") + labs(x="Height (m)") + theme_bw()
ggplot(pokemon, aes(height_m)) + geom_histogram(binwidth = .5, color="black", fill="white") + labs(x="Height (m)") + theme_bw()


ggplot(pokemon, aes(hp)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Hit points (HP)") + theme_bw()+lims(x=c(0, 255))
ggplot(pokemon, aes(hp)) + geom_histogram(binwidth = 10, color="black", fill="white") + labs(x="Hit points (HP)") + theme_bw()+lims(x=c(0, 255))


ggplot(pokemon, aes(attack)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Attack") + theme_bw()+lims(x=c(0, 255))
ggplot(pokemon, aes(attack)) + geom_histogram(binwidth = 10, color="black", fill="white") + labs(x="Attack") + theme_bw()+lims(x=c(0, 255))


ggplot(pokemon, aes(defense)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Defense") + theme_bw()+lims(x=c(0, 255))
ggplot(pokemon, aes(defense)) + geom_histogram(binwidth = 10, color="black", fill="white") + labs(x="Defense") + theme_bw()+lims(x=c(0, 255))


ggplot(pokemon, aes(sp_attack)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Special attack") + theme_bw()+lims(x=c(0, 255))
ggplot(pokemon, aes(sp_attack)) + geom_histogram(binwidth = 10, color="black", fill="white") + labs(x="Special attack") + theme_bw()+lims(x=c(0, 255))


ggplot(pokemon, aes(sp_defense)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Special defense") + theme_bw()+lims(x=c(0, 255))
ggplot(pokemon, aes(sp_defense)) + geom_histogram(binwidth = 10, fill="white", color="black") + labs(x="Special defense") + theme_bw()+lims(x=c(0, 255))


ggplot(pokemon, aes(speed)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Speed") + theme_bw()+lims(x=c(0, 255))
ggplot(pokemon, aes(speed)) + geom_histogram(binwidth = 10, fill="white", color="black") + labs(x="Speed") + theme_bw()+lims(x=c(0, 255))


ggplot(pokemon, aes(total)) + geom_histogram(binwidth = 2, color="black", fill="white") + labs(x="Total") + theme_bw()
ggplot(pokemon, aes(total)) + geom_histogram(binwidth = 20, fill="white", color="black") + labs(x="Total") + theme_bw()


ggplot(pokemon, aes(catch_rate)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Catch rate") + theme_bw()
ggplot(pokemon, aes(percentage_male)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Percentage male") + theme_bw()

ggplot(pokemon, aes(base_happiness_8)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Base happiness") + theme_bw()



```



```{r, out.width="50%"}
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()


```






:::: {.infobox data-latex=""}
**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. 
::::


:::: {.infobox data-latex=""}
**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. 


```{r UnivariableGraphsLog, out.width="45%", echo=FALSE, warning=FALSE, message = FALSE}

ggplot(pokemon, aes(weight_kg)) + geom_histogram(binwidth = 1, color="black", fill="white") + labs(x="Weight (kg) - log2") + theme_bw() + scale_x_continuous(trans = 'log2') 


ggplot(pokemon, aes(height_m)) + geom_histogram(binwidth = 2, color="black", fill="white") + labs(x="Height (m) - log2") + 
theme_bw() + scale_x_continuous(trans = 'log2') 


 qqnorm(log2(pokemon$weight_kg), main="log2(Weight)")
 qqline(log2(pokemon$weight_kg))

 
 qqnorm(log2(pokemon$height_m), main="log2(Height)")
 qqline(log2(pokemon$height_m))




```

:::: {.infobox data-latex=""}
**RESULTS: **
We note that both distributions are significantly closer to normality after log transformation.
:::: 

:::: {.infobox data-latex=""}
**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. 




:::: {.infobox data-latex=""}
**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 `r sum(pokemon$weight_kg>500)` Pokémon have a weight greater than 500kg and only `r  sum(pokemon$weight_kg>900)` Pokémon have a weight greater than 900kg. The 10 most extreme values are shown below, a table reports the details of the `r  sum(pokemon$weight_kg>900)` with weight above 900 kg. 

```{r extremesWeight}
tail(sort(pokemon$weight_kg), 10)

```


```{r tabwithFigs1}
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
              ) 


```

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


```{r extremesHeight}
tail(sort(pokemon$height_m), 10)

```

 Only `r sum(pokemon$height_m>10)` Pokémon have more than 10 m of height, `r sum(pokemon$height_m>5)` more than 5 m, `r sum(pokemon$height_m>3)` more than 3 m,  `r sum(pokemon$height_m>2)` more than 2 m; only 5% of Pokémon are taller than `r as.numeric(quantile(pokemon$height_m, 0.95))` m. 
 
A table reports the details of the `r  sum(pokemon$height_m>7.5)` with height above 7.5 m. 



```{r tabwithFigs}
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
              ) 

```



:::: {.infobox data-latex=""}
**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.
:::: 

:::: {.infobox data-latex=""}
**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


:::: {.infobox data-latex=""}
**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.


:::: {.infobox data-latex=""}
**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.
:::: 

:::: {.infobox data-latex=""}
**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. 
:::: 


```{r ch153, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


s <-
     Hmisc::summaryM(
         PokemonPhase_5cat+  height_m+  special_status + 
         hp+ attack+defense+  sp_attack+ sp_defense+ speed+ total+ catch_rate+    gender_cat + happiness_cat  ~ factor(generation),
         data = pokemon %>% mutate(percentage_male=factor(percentage_male)),
         overall = FALSE,
         test = FALSE
     )


Hmisc::html(
    s,
    caption = 'Characteristics by generation',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```

```{r chGraphsByGeneration, out.width="49%", echo= FALSE}

ggplot(pokemon, aes(x = factor(generation), fill = gender_cat)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Generation",
    y = "Proportion",
    fill = "Gender Category",
    title = "Proportion of Gender Category by Generation"
  ) +
  theme_minimal()

ggplot(pokemon, aes(x = factor(generation), fill = happiness_cat)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Generation",
    y = "Proportion",
    fill = "Happiness Category",
    title = "Proportion of Happiness Category by Generation"
  ) +
  theme_minimal()


ggplot(pokemon, aes(x = factor(generation), fill = PokemonPhase_5cat)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Generation",
    y = "Proportion",
    fill = "Evolutionary Phase",
    title = "Proportion of Evolutionary Phases by Generation"
  ) +
  theme_minimal()




ggplot(pokemon, aes(x = factor(generation), fill = as.factor(dual_type))) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Generation",
    y = "Proportion",
    fill = "Dual Type",
    title = "Proportion of Single vs Dual-Type Pokémon by Generation"
  ) +
  theme_minimal()

ggplot(pokemon, aes(x = factor(generation), fill = special_status)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Generation",
    y = "Proportion",
    fill = "Special Status",
    title = "Proportion of Pokémon by Special Status and Generation"
  ) +
  theme_minimal()
#mosaicplot(table(pokemon$generation, pokemon$is_mythical), xlab = "Generation", ylab = "Mythical Status", las = 1, main = "Generation vs Mythical Status Mosaic")


```


```{r, echo=FALSE, out.width="49%", eval =FALSE}

mosaicplot(table(pokemon$generation, pokemon$gender_cat), xlab = "Generation", ylab = "Gender", las = 1, main = "Generation vs Categorical Gender Mosaic")

mosaicplot(table(pokemon$generation, pokemon$happiness_cat), xlab = "Generation", ylab = "Happiness", las = 1, main = "Generation vs Categorical Happiness Mosaic")

mosaicplot(table(pokemon$generation, pokemon$PokemonPhase_5cat), xlab = "Generation", ylab = "Evolutionary Phase", las = 1, main = "Generation vs Evolutionary Phase Mosaic")


mosaicplot(table(pokemon$generation, pokemon$dual_type), xlab = "Generation", ylab = "Dual Type", las = 1, main = "Generation vs Dual Type Mosaic")

mosaicplot(table(pokemon$generation, pokemon$special_status), xlab = "Generation", ylab = "Legendary Status", las = 1, main = "Generation vs Legendary Status Mosaic")

#mosaicplot(table(pokemon$generation, pokemon$is_mythical), xlab = "Generation", ylab = "Mythical Status", las = 1, main = "Generation vs Mythical Status Mosaic")


```



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



:::: {.infobox data-latex=""}
**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). 
:::: 

:::: {.infobox data-latex=""}
**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.
:::: 



```{r ch259, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


s <-
     Hmisc::summaryM(
         # Wave+ 
          #country+ 
          factor(generation) + height_m+  special_status + 
         hp+ attack+defense+  sp_attack+ sp_defense+ speed+ total+ catch_rate+    gender_cat + happiness_cat ~ factor( PokemonPhase_5cat),
         data =  pokemon %>% mutate(percentage_male=factor(percentage_male)),
         overall = FALSE,
         test = FALSE
     )


Hmisc::html(
    s,
    caption = 'Characteristics by Pokémon evolution phase',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```



```{r chGraphsByPhase, out.width="49%", echo=FALSE}


 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),height_m)) + geom_boxplot() + labs(y="Height (m)", x="Evolution Phase") + theme_bw() + scale_y_continuous(trans = "log2")
 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),hp)) + geom_boxplot() + labs(y="Hit points (HP)", x="Evolution Phase") + theme_bw() 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),attack)) + geom_boxplot() + labs(y="Attack", x="Evolution Phase") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),defense)) + geom_boxplot() + labs(y="Defense", x="Evolution Phase") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),sp_attack)) + geom_boxplot() + labs(y="Special attack", x="Evolution Phase") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),sp_defense)) + geom_boxplot() + labs(y="Special defense", x="Evolution Phase") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),speed)) + geom_boxplot() + labs(y="Speed", x="Evolution Phase") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(PokemonPhase_5cat),total)) + geom_boxplot() + labs(y="Total", x="Evolution Phase") + theme_bw() 
 


```

```{r graphsByPhase, out.width="45%"}

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


```


```{r, echo=FALSE, out.width="49%", eval= FALSE}
mosaicplot(table(pokemon$PokemonPhase_5cat, pokemon$gender_cat), xlab = "Evolutionary Phase", ylab = "Gender", las = 1, main = "Evolutionary Phase vs Categorical Gender Mosaic")

mosaicplot(table(pokemon$PokemonPhase_5cat, pokemon$happiness_cat), xlab = "Evolutionary Phase", ylab = "Happiness", las = 1, main = "Evolutionary Phase vs Categorical Happiness Mosaic")

mosaicplot(table(pokemon$PokemonPhase_5cat, pokemon$generation), xlab = "Evolutionary Phase", ylab = "Generation", las = 1, main = "Evolutionary Phase vs Generation Mosaic")

mosaicplot(table(pokemon$PokemonPhase_5cat, pokemon$chain_length), xlab = "Evolutionary Phase", ylab = "Chain Length", las = 1, main = "Evolutionary Phase vs Chain Length Mosaic")

mosaicplot(table(pokemon$PokemonPhase_5cat, pokemon$dual_type), xlab = "Evolutionary Phase", ylab = "Dual Type", las = 1, main = "Evolutionary Phase vs Dual Type Mosaic")

mosaicplot(table(pokemon$PokemonPhase_5cat, pokemon$special_status), xlab = "Evolutionary Phase", ylab = "Legendary Status", las = 1, main = "Evolutionary Phase vs Legendary Status Mosaic")





```

##### Descriptive statistics of selected variables, by gender

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


:::: {.infobox data-latex=""}
**COMMENT: **
A subgroup analysis is planned by gender category.
:::: 


:::: {.infobox data-latex=""}
**RESULT: **
The genderless category differs most from the other categories. The mostly female category are the shortest species.
:::: 


```{r multi_gender, results="asis", message=FALSE, warning =FALSE , echo=FALSE}


s <-
     Hmisc::summaryM(
         PokemonPhase_5cat+  height_m+  special_status + 
         hp+ attack+defense+  sp_attack+ sp_defense+ speed+ total+ catch_rate+     happiness_cat  ~ factor(gender_cat),
         data = pokemon,
         overall = FALSE,
         test = FALSE
     )


Hmisc::html(
    s,
    caption = 'Characteristics by gender',
    exclude1 = TRUE,
    npct = 'both',
    digits = 3,
    prmsd = TRUE,
    brmsd = TRUE,
    msdsize = mu$smaller2
)
```


```{r chGraphsByGender, out.width="90%", echo= FALSE}

ggplot(pokemon, aes(x = log2(height_m), group = gender_cat)) +
  geom_histogram() +
  labs(
    x = "log2(height)",
    y = "Frequency",
    title = "Distributions of Pokémon height by gender"
  ) + facet_grid(.~gender_cat)+
  theme_minimal()


```


##### Association of selected variables with height

:::: {.infobox data-latex=""}
**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. 

:::: 

:::: {.infobox data-latex=""}
**CONSEQUENCES: **
The non-linear association between base stats and (log2) height might anticipate non-linear associations also with the outcome variable.
:::: 





```{r chGraphsWithWeight, out.width="49%", echo=FALSE}



 ggplot(pokemon, aes(height_m, hp )) + geom_point(pos="jitter") + geom_smooth(pos="jitter") + labs(x="Height (m) - log2", y="HP") + theme_bw()+scale_x_continuous(trans="log2")
 
 ggplot(pokemon, aes(height_m, attack )) +geom_point(pos="jitter") + geom_smooth(pos="jitter") + labs(y="Attack", x="Height (m) - log2") + theme_bw() +scale_x_continuous(trans="log2")
 
 ggplot(pokemon, aes(height_m, defense)) + geom_point(pos="jitter") + geom_smooth() + labs(y="Defense", x="Height (m) - log2") + theme_bw() +scale_x_continuous(trans="log2")
 
 ggplot(pokemon, aes(height_m, sp_attack)) + geom_point(pos="jitter") +  geom_smooth() + labs(y="Special attack", x="Height (m) - log2") + theme_bw() +scale_x_continuous(trans="log2")
 
 ggplot(pokemon, aes(height_m, sp_defense)) + geom_point(pos="jitter") +  geom_smooth() + labs(y="Special defense", x="Height (m) - log2") + theme_bw() +scale_x_continuous(trans="log2")
 
 ggplot(pokemon, aes(height_m, speed)) + geom_point(pos="jitter") +  geom_smooth() + labs(y="Speed", x="Height (m) - log2") + theme_bw() +scale_x_continuous(trans="log2")
 
 ggplot(pokemon, aes(height_m, total)) + geom_point(pos="jitter") +  geom_smooth() + labs(y="Total", x="Height (m) - log2") + theme_bw() +scale_x_continuous(trans="log2")
 


```




#### Correlation between explanatory variables 

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


:::: {.infobox data-latex=""}
**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.
:::: 




:::: {.infobox data-latex=""}
**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. 



```{r MultivariableCorrelation2}
#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
```





```{r MultivariableHeatmap}

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.





:::: {.infobox data-latex=""}
**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.
::::



:::: {.infobox data-latex=""}
**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. 
::::

```{r}
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()
```





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




```{r, out.width="45%", echo=FALSE}
base_stats <- c("hp","attack","defense","sp_attack","sp_defense","speed")

# standardize the 6 base stats (z-scores)
z_stats <- scale(pokemon[, base_stats])

# within-profile imbalance: SD of each Pokemon's own 6 standardized values
pokemon$profile_imbalance <- apply(z_stats, 1, sd)

top_imbalance <- pokemon %>%
  arrange(desc(profile_imbalance)) %>%
  select(name, all_of(base_stats), profile_imbalance) %>%
  slice(1:10)

top_imbalance %>%
  kbl(booktabs = TRUE, digits = 2, align = "c",
      caption = "Pokémon species with the most unbalanced combination of base stats") %>%
  kable_styling() %>%
  kable_paper(full_width = TRUE)
```

:::: {.infobox data-latex=""}
**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).
::::


:::: {.infobox data-latex=""}
**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.
::::


:::: {.sapdev data-latex=""}

**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**

```{r, out.width="45%", echo=FALSE}


ggplot(pokemon, aes(defense, hp)) + geom_point(pos="jitter", alpha=.5, aes(color=special_status)) +
  geom_label_repel(data=pokemon[pokemon$name %in% top_imbalance$name,], aes(label=name, color=special_status)) +
  labs(x="Defense", y="HP") + theme_bw() + geom_smooth()

ggplot(pokemon, aes(attack, hp)) + geom_point(pos="jitter", alpha=.5, aes(color=special_status)) +
  geom_label_repel(data=pokemon[pokemon$name %in% top_imbalance$name,], aes(label=name, color=special_status)) +
  labs(x="Attack", y="HP") + theme_bw() + geom_smooth()

ggplot(pokemon, aes(defense, speed)) + geom_point(pos="jitter", alpha=.5, aes(color=special_status)) +
  geom_label_repel(data=pokemon[pokemon$name %in% top_imbalance$name,], aes(label=name, color=special_status)) +
  labs(x="Defense", y="Speed") + theme_bw() + geom_smooth()

ggplot(pokemon, aes(speed, sp_defense)) + geom_point(pos="jitter", alpha=.5, aes(color=special_status)) +
  geom_label_repel(data=pokemon[pokemon$name %in% top_imbalance$name,], aes(label=name, color=special_status)) +
  labs(x="Speed", y="Special defense") + theme_bw() + geom_smooth()

ggplot(pokemon, aes(sp_attack, sp_defense)) + geom_point(pos="jitter", alpha=.5, aes(color=special_status)) +
  geom_label_repel(data=pokemon[pokemon$name %in% top_imbalance$name,], aes(label=name, color=special_status)) +
  labs(x="Special attack", y="Special defense") + theme_bw() + geom_smooth()
```



```{r}
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


:::: {.infobox data-latex=""}
**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.

:::: 
```{r chGraphsBySpecial, out.width="49%", echo=FALSE}


 
 ggplot(pokemon, aes(as.factor(special_status),height_m)) + geom_boxplot() + labs(y="Height (m)", x="Special Pokémon") + theme_bw() + scale_y_continuous(trans = "log2")
 

 ggplot(pokemon, aes(as.factor(special_status),hp)) + geom_boxplot() + labs(y="Hit points (HP)", x="Special Pokémon") + theme_bw() 
 ggplot(pokemon, aes(as.factor(special_status),attack)) + geom_boxplot() + labs(y="Attack", x="Special Pokémon") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(special_status),defense)) + geom_boxplot() + labs(y="Defense", x="Special Pokémon") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(special_status),sp_attack)) + geom_boxplot() + labs(y="Special attack", x="Special Pokémon") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(special_status),sp_defense)) + geom_boxplot() + labs(y="Special defense", x="Special Pokémon") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(special_status),speed)) + geom_boxplot() + labs(y="Speed", x="Special Pokémon") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(special_status),total)) + geom_boxplot() + labs(y="Total", x="Special Pokémon") + theme_bw() 
 


```

:::: {.infobox data-latex=""}
**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. 


::::

:::: {.infobox data-latex=""}
**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

:::: {.infobox data-latex=""}
**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.

:::: 


```{r chGraphsByHappiness, out.width="49%", echo=FALSE}


 ggplot(pokemon, aes(as.factor(happiness_cat),height_m)) + geom_boxplot() + labs(y="Height (m)", x="Happiness group") + theme_bw() + scale_y_continuous(trans = "log2")
 

 ggplot(pokemon, aes(as.factor(happiness_cat),hp)) + geom_boxplot() + labs(y="Hit points (HP)", x="Special Pokémon") + theme_bw() 
 ggplot(pokemon, aes(as.factor(happiness_cat),attack)) + geom_boxplot() + labs(y="Attack", x="Happiness group") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(happiness_cat),defense)) + geom_boxplot() + labs(y="Defense", x="Happiness group") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(happiness_cat),sp_attack)) + geom_boxplot() + labs(y="Special attack", x="Happiness group") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(happiness_cat),sp_defense)) + geom_boxplot() + labs(y="Special defense", x="Happiness group") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(happiness_cat),speed)) + geom_boxplot() + labs(y="Speed", x="Special Pokémon") + theme_bw() 
 
 ggplot(pokemon, aes(as.factor(happiness_cat),total)) + geom_boxplot() + labs(y="Total", x="Happiness group") + theme_bw() 
 


```



:::: {.infobox data-latex=""}
**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.  
::::

:::: {.infobox data-latex=""}
**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.

```{r dendrogram1, fig.height=4}

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

```

:::: {.infobox data-latex=""}
**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. 

```{r saveFiles, eval=TRUE}

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


```