Final Project.Rmd

---
title: "CS555 R Final Project Notebook"
output:
  html_document:
    df_print: paged
  pdf_document:
    latex_engine: xelatex
  html_notebook:
    toc: true 
editor_options: 
  chunk_output_type: inline
---
```{r}
library('tidyverse')
library('reshape2')
library('car')
```
### Research Scenario
We have soccer player dataset of 500 rows. With this research we aim to answer the following questions

#### Research questions
-   Is their a relationship between value and wage of the soccer player ?
-   How strong is the relationship ?
-   Which other factors are associated with player's wage ?
-   How large is the association between each factor and wage ?
-   How accurately can we predict wage of the player ?
-   Is the relationship linear ?

The goal of the research is to predict wage of the soccer player

## 1. Data Source [link](https://raw.githubusercontent.com/kiat/R-Examples/master/projects/datasets/project-10.csv)

This is the FIFA19 dataset consists of 15 variables such as age, nationality, name, club, jersey number, height, weight, value, wage, reaction, balance, strength etc of a list of 500 soccer players exist in the game who are currently playing in different soccer clubs around the world.

In our analysis we plan to use following variables

-   Age (integer): The age of the player.
-   Nationality (character): The nationality of the player.
-   Overall (integer): The overall rating of the player's skills and abilities.
-   Potential (integer): The potential future rating of the player's skills and abilities.
-   Club (character): The football club the player is currently associated with.
-   Value_in_M (double): The market value of the player in millions of dollars.
-   Wage_in_K (integer): The weekly wage of the player in thousands of dollars.
-   Height_in_Inches (integer): The height of the player in inches.
-   Weight_in_lbs (integer): The weight of the player in pounds.
-   Reactions (integer): The player's reaction attribute rating.
-   Balance (integer): The player's balance attribute rating.
-   Strength (integer): The player's strength attribute rating.

```{r}
data <- read.csv('https://raw.githubusercontent.com/kiat/R-Examples/master/projects/datasets/project-10.csv')
data
```
```{r}
glimpse(data)
```

## 2. Data Cleaning

Renaming of the club names with proper names with no special characters is done so that data can be parsed and processed.
Rows with incorrect column values are removed like players with no wage.

```{r}
data$Club[data$Club == "Atl\xe9tico Madrid"] <- "Atletico Madrid"
data$Club[data$Club == "FC Bayern M\xfcnchen"] <- "FC Bayern Munchen"
data$Club[data$Club == "Borussia M\xf6nchengladbach
"] <- "Borussia Monchengladbach"
data$Club[data$Club == "Be?ikta? JK"] <- "Besiktas JK"
data$Club[data$Club == "Fenerbah\xe7e SK"] <- "Fenerbahce SK"
data$Club[data$Club == "Gr\xeamio"] <- "Gremio"
data$Club[data$Club == "AS Saint-\xc9tienne"] <- "AS Saint Etienne"
data$Club[data$Club == "Medipol Ba?ak?ehir FK"] <- "Basaksehir FK"
data$Club[data$Club == "Borussia M\xf6nchengladbach"] <- "Borussia Monchengladbach"
data$Club[data$Club == "1. FC K\xf6ln"] <- "1. FC Koln"
data$Club[data$Club == "Deportivo Alav\xe9s"] <- "Deportivo Alaves"
data$Club[data$Club == "Guangzhou R&F; FC"] <- "Guangzhou City F.C."
data$Club[data$Club == "Atl\xe9tico Mineiro"] <- "Atletico Mineiro"


data <- data[data$Wage_in_K!=0,]
names(data) <- make.names(names(data))

write.csv(data,'clean_data.csv')
```


## 3. Data Analysis

Various statistical methods like **Correlation tests**, **Simple Linear Regression**, **Multiple Linear Regression**, **ANOVA and ANCOVA** are used for analysis

### Correlation test

Check whether value of the player is associated with the his wage ?

```{r}
attach(data)
view(data)
```

```{r}
plot( Value_in_M, Wage_in_K, xlab = 'Value', ylab = 'Wage', main = 'Relation between wage and value of the player')
cor(Value_in_M, Wage_in_K)
```

The correlation coefficient of **0.75** indicates *value* of the player is **strongly** associated with his *wage* and the direction is **positive**.

### 3.1    Simple Least Square Linear Regression


\[ \text{wage} = f(\text{value}) \]
\[Y_{wage} = \beta_0 + \beta_1 \times x_{value}\]


```{r}
# SLR
m <- lm(Wage_in_K ~ Value_in_M)
# plot SLR
plot( Value_in_M, Wage_in_K, xlab = 'Value', ylab = 'Wage', main = 'SLR between wage and value of player')
abline(m, col='red')
# SLR summary
summary(m)
mse <- mean(residuals(m)^2)
print(mse)
# CI of coefficients
upper.c1 <- coef(m)[1]+qt(0.95,500-2)*4.3343
lower.c1 <- coef(m)[1]-qt(0.95,500-2)*4.3343

upper.value <- coef(m)[2]+qt(0.95,500-2)*0.1367
lower.value <- coef(m)[2]-qt(0.95,500-2)*0.1367
```

-   As **p-value** of the model \<\< **0.05** for \\(\\alpha\\)=0.05, model is **significant**. Hence their is good evidence that ***Value*** is significant predictor of ***Wage***.

-   ***Value*** predictor explains **56.8%** of variance in ***Wage***.

-   **95% confidence interval** for \\(\\beta_0\\) is [-8.28, 6] and the 95% confidence interval for \\(\\beta_1\\) is [3.27, 3.72]. Therefore we can conclude that if the player has no value, his wage will fall between 0 to 6K dollars. Futhermore, for each 100 million increase in value of the player, there will be an average increase in wage of between 327K to 372K dollars.

### -   Residual plot analysis

```{r}
plot(Value_in_M, residuals(m), xlab = 'Value of the player (X)', ylab = 'Residual')
abline(h=0)
plot(Wage_in_K, residuals(m), xlab = 'Wage of the player(Y)', ylab = 'Residual')
abline(h=0)
hist(residuals(m))
```

-   Clear **obvious pattern** emerges in residual plot. 
-   As X increase **variability** of outcome increase. 
-   Residuals are normally distributed. 
-   It clearly **violates** the linearity and constant variance assumptions. Hence, non-linear relationship exist between the factors Wage and value.

### -   Transformation

```{r}
new.data <- data[Wage_in_K!=0,]
```

#### -   Take log of response variable

\[log(Y_{wage}) = \beta_0 + \beta_1 x_{value}\]

```{r}
cor(new.data$Value_in_M, log(new.data$Wage_in_K))

m <- lm(log(new.data$Wage_in_K)~new.data$Value_in_M)
summary(m)
```

Strength of linear relationship (correlation) decreases to 0.54. R2 of model 0.30

#### -   Take square root of response variable
\[\sqrt{Y_{wage}} = \beta_0 + \beta_1 x_{value}\]

```{r}
cor(sqrt(data$Wage_in_K), Value_in_M)

m <- lm(sqrt(data$Wage_in_K)~Value_in_M)
summary(m)

```

Correlation is 0.70, R2 is 0.49

#### -   Squaring explanatory variable
\[Y_{wage} = \beta_0 + \beta_1 x_{value}^2\]
```{r}
cor(Wage_in_K, Value_in_M**2)

m <- lm(Wage_in_K~Value_in_M**2)
summary(m)
```

Correlation is 0.72 and R2 of model is 0.57

#### -   Log of explanatory variable

\[Y_{wage} = \beta_0 + log_(\beta_1 x_{value})\]
```{r}
new.data <- data[Value_in_M!=0,]
cor(new.data$Wage_in_K, log(new.data$Value_in_M))

m <- lm(new.data$Wage_in_K~log(new.data$Value_in_M))
summary(m)
```

Correlation is 0.66 and R2 of model is 0.43

#### -   Polynomial regression
\[Y_{wage} = \beta_0 + \beta_1 x_{value} + \beta_2 x_{value}^2\]

```{r}
m <- lm(Wage_in_K ~ Value_in_M + I(Value_in_M**2))
summary(m)
```

-   Model is significant but **R2** is same and the coefficient of newly added quadratic term is insignificant, by t-test.

-   By doing **transformations**, we saw that the **correlation** between Wage and Value of the player **reduces** which in turn **drops** linear model's performance. Even with the transformations we don't have sufficient evidence of linear relationship between wage and value of the player, their may be other factors that cause variation in wage.

### -   Corelation Matrix

```{r}
numeric_cols <- data %>% select_if(is.numeric)

cor_matrix <- cor(numeric_cols)
cor_matrix_melted <- melt(cor_matrix)
ggplot(cor_matrix_melted, aes(Var1, Var2, fill = value, label = round(value, 2))) +
  geom_tile(color = "white") +
  geom_text(size = 3, color = "black") +  # Add text labels for correlation coefficients
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1, 1), 
                       space = "Lab", name = "Correlation") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Correlation Heatmap Matrix", x = "Variables", y = "Variables")
```

-   From the above plot we can see that ***Overall score*** and ***Value*** variables have strong association with the ***wage*** of the player.
-   ***Reactions***,***Potential*** variables are moderately associated.
-   While the ***Balance***, ***Strength*** and ***Age*** are weakly associated.

### 3.2 Multiple Linear Regression


```{r}
mlr.m2 <- lm(Wage_in_K ~ Value_in_M+Overall+Balance+Strength+Potential+Age+Height_in_Inches+Weight_in_lbs+Reactions, data = data)
mean(residuals(mlr.m2)^2)

mlr.m2 <- lm(Wage_in_K ~ Value_in_M+Overall+Balance+Strength+Potential+Age+Height_in_Inches+Reactions, data = data)
mean(residuals(mlr.m2)^2)
mlr.m2 <- lm(Wage_in_K ~ Value_in_M+Overall+Balance+Strength+Potential+Age+Reactions, data = data)
mean(residuals(mlr.m2)^2)
mlr.m2 <- lm(Wage_in_K ~ Value_in_M+Overall+Balance+Strength+Potential+Age, data = data)
mean(residuals(mlr.m2)^2)
mlr.m2 <- lm(Wage_in_K ~ Value_in_M+Overall+Balance+Strength+Age, data = data)
mean(residuals(mlr.m2)^2)
mlr.m2 <- lm(Wage_in_K ~ Value_in_M+Overall+Balance+Strength, data = data)
mean(residuals(mlr.m2)^2)
summary(mlr.m2)

```
#### MLR Model Summary 

|     | MLR Model                                                                           | Global test | Adjusted R2 | Insignificant variable(s)   \\(\\alpha\\)=0.05                    |
|-----|-------------------------------------------------------------------------------------|-------------|-------------|------------------------------------------------|
| 1.   | \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{potential}+\beta_3 x_{value}+\beta_4 x_{reaction}+\beta_5 x_{balance}+\beta_6 x_{strength}+\beta_7 x_{age}+\beta_8 x_{weight}+\beta_9 x_{height}\] | Significant | \\(0.6465\\)      | potential, reactions, strength, weight, height
 |
| 2.   | \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{potential}+\beta_3 x_{value}+\beta_4 x_{reaction}+\beta_5 x_{balance}+\beta_6 x_{strength}+\beta_7 x_{age}+\beta_8 x_{height}\]         | Significant | \\(0.6528\\)      | potential, strength, height, reactions         |
| 3.   | \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{potential}+\beta_3 x_{value}+\beta_4 x_{reaction}+\beta_5 x_{balance}+\beta_6 x_{strength}+\beta_7 x_{age}\]                 | Significant | \\(0.6461\\)      | reactions, potential                           |
| 4.   | \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{potential}+\beta_3 x_{value}+\beta_4 x_{balance}+\beta_5 x_{strength}+\beta_6 x_{age}\]                            | Significant | \\(0.6467\\)      | potential                                      |
| 5.   | \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}+\beta_5 x_{age}\]                                       | Significant | \\(0.647\\)       |    age                                            |
| 6.   | \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}\]                                            | Significant | \\(0.6457\\)      |                                                |


-   All models have almost same **R2** value with **marginal differences**. 
-   All the models explain around **64.5%** variability in player's ***wage***. 
-   The model **\(Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}\)** is the best MLR model as it explains the **same** amount of **variation** with **less variables**. 
-   Hence important variables in determining wage are, player's *value*, *potential*, *balance* and *strength*. 
-   MLR model has better R2 of **64.5%** than that of SLR model of **56%**. Hence it explains **greater** variability in Y.


```{r}
plot(Wage_in_K, residuals(mlr.m2))
```

### -   Adding group variable `Club`

```{r}
length(unique(Club))
length(unique(Nationality))
```

In the dataset their are 92 different Clubs from 59 different Nationality.

### -   Distribution of Wage across Clubs

**Select top 20 Clubs with highest mean wage**

```{r}
# get top 20 clubs by highest mean wage
data %>% group_by(Club) %>% summarise(Mean_Wage = mean(Wage_in_K)) %>% arrange(desc(Mean_Wage)) %>% select(Club) %>% head(20) -> top.clubs

data %>% filter(Club %in% top.clubs$Club) -> top.clubs.data

```

```{r}
options(repr.plot.width = 100, repr.plot.height = 6)
ggplot(data=top.clubs.data, aes(x=Club, y=Wage_in_K)) +
  geom_boxplot()+
  labs(x = "Clubs", y = "Wage (in K)", title = 'Distribution of Wage across Clubs')

```

Club with the mean highest wage is *Real Madrid* followed by FC Barcelona and Manchester City.

### -   Does a player's wage differ across soccer clubs ?

#### ANOVA Test
\[ \text{wage} = f(\text{club}) \]
```{r}
m <- aov(Wage_in_K ~ Club, data=data)

summary(m)
```
As p-value of ***Club*** in anova model \<\< 0.05. we have significant evidence that wage of the player **varies** by Clubs they play for.
Therefore wages change depending on the different Clubs.

```{r}
# pairwise t-test using bonferroni method
pairwise.t.test(Wage_in_K, Club, data=data,p.adj='bonferroni')
```

```{r}
# tukeys pair wise t-test
tukey_result <- TukeyHSD(m)
# convert to dataframe
comparisons <- as.data.frame(tukey_result$Club)
# filter pairs whose p-value<0.05
comparisons %>% filter(comparisons$`p adj`<0.05) %>% select('p adj') -> diff.groups

diff.groups
```

By doing pairwise comparisons using **Tukeys method** the mean wages of the above club combinations are different.\
Out of **4095 pairs**, **280** club pairs are **different**.

### -   Does wage differ by clubs after controlling other variables ?
#### ANCOVA Test
\[ \text{wage} = f(\text{value}+\text{overall}+\text{strength}+\text{club}) \]
```{r}
Anova(lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Club, data=data), type=3)
```
From individual f-test, p-value of ***Club*** << 0.05. ***Club*** is a **significant** predictor of a player's wage after controlling other factors like value, balance, strength and overall.

```{r}
library(emmeans)
m <- lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Club, data=data)
emmeans.result <- emmeans(m, specs = "Club" , contr = "pairwise",  adjust="tukey")$contrast
emmeans.data <- as.data.frame(emmeans.result)
emmeans.data %>% filter(p.value<0.05) %>% view()
```


```{r}
m <- lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Club, data=data) #R2=0.875
mean(residuals(m)^2)
m <- lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Potential+Club, data=data) #R2=0.905
mean(residuals(m)^2)
m <- lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Potential+Age+Club, data=data) #R2=0.909
mean(residuals(m)^2)
summary(m)
```

| | Model | Adjusted R2 | Global Test |
|-|:-------:|:-----------:|:---------|
|1.| \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}+\beta_5 x_{club}\] | \\(0.875\\) | Significant |
|2.| \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}+\beta_5 x_{potential}+\beta_6 x_{club}\] | \\(0.905\\) | Significant |
|3.| \[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}+\beta_5 x_{potential}+ +\beta_6 x_{age} + \beta_7 x_{club}\] | \\(0.909\\) | Significant |


-   From the analysis after including group variable `Club` in the model, it explains much **higher** variability in wage i.e from **65%** to around **87%**.
-   By adding other variables, **R2** of the model increases to **0.90**.
-   **\(Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}+\beta_5 x_{potential}+ +\beta_6 x_{age} + \beta_7 x_{club}\)** is the **best model** as it has highest accuracy.

### -   Adding group variable `'Nationality'`

#### -    Does wage vary by nationalaity ?

#####  ANOVA test
```{r}
m <- aov(Wage_in_K ~ Nationality, data=data)
summary(m)
```
Model is **not** significant as p-value of *nationality* > 0.05. We have significant evidence that wages of player does not vary by their nationality.

### Final model
\[Y_{wage} = \beta_0+\beta_1 x_{overall}+\beta_2 x_{value}+\beta_3 x_{balance}+\beta_4 x_{strength}+\beta_5 x_{potential}+ +\beta_6 x_{age} + \beta_7 x_{club}\]

#### Aceses linear assumption
```{r}
mlr.final = lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Potential+Age+Club, data=data)
plot(mlr.final)
plot(Wage_in_K ,residuals(mlr.final))
```

For extremely high wages ( >400K ), the residual error is very large. Hence could be possible outliers.

### Outlier Removal
```{r}
data %>% filter(Wage_in_K>400) %>% view()
data %>% ggplot(aes(x = Wage_in_K)) +
  geom_histogram(color = 'black') +
  geom_vline(xintercept = median(Wage_in_K), colour = 'red', show.legend = TRUE) +
  geom_vline(xintercept = mean(Wage_in_K), colour = 'blue', show.legend = TRUE) +
  labs(x = 'Wage (in K)', y = 'Frequency', title = 'Histogram of Wage Distribution')
```
Right skewed distribution
Most soccer players have below median wage and very few players have significantly large wage.


```{r}
Q1 <- quantile(Wage_in_K, 0.25)
Q3 <- quantile(Wage_in_K, 0.75)
IQR <- Q3 - Q1

Med <- median(Wage_in_K)

filtered.data<-data[!(Wage_in_K>(Q3+1.5*IQR)|Wage_in_K<(Q1-1.5*IQR)),]

m <- lm(Wage_in_K ~ Overall+Value_in_M+Balance+Strength+Potential+Club, data=filtered.data)
summary(m)
plot(m)
plot(filtered.data$Wage_in_K,residuals(m))
```


## 3.3 Logistic Regression
```{r}
view(data)
data %>% group_by(Club) %>% summarise(Mean_wage = median(Wage_in_K)) %>% view()
data %>% group_by(Club) %>% summarise(Mean_wage = mean(Wage_in_K))  -> tmp

tmp %>% ggplot(aes(x=Mean_wage))+
  geom_histogram(binwidth = 10,fill = "skyblue", color = "black")+
  geom_vline(xintercept = mean(tmp$Mean_wage), color='red')+
  labs(title = 'Distribution of Mean wages of clubs', x='Club Mean wage')
```
Right Skewed distribution
Most of the clubs have very low wages, while very few clubs have significantly high wages 

```{r}
data$High_Wage <- ifelse(Wage_in_K > mean(tmp$Mean_wage), 1, 0)
attach(data)
m <- glm(High_Wage ~ Value_in_M, family = binomial)
summary(m)
```

## Conclusion
**Wage** of soccer player depends on **value, overall rating, balance rating, potential, strength rating, age and club**. Given these factors are known, wage of a player could be calculated with **90%** accuracy. 

Although model explains 90% of wage, its linear assumptions may **not** be met since their is increasing variance of residuals with higher fitted values.