PGA Tour Data Analysis 2010-2020

Analyse statistique de données

Quentin Deschamps

Avril-Mai 2021

1 Introduction

Durant cette analyse de données, nous allons étudier les statistiques des golfeurs du PGA Tour entre 2010 et 2020.

Cette partie introductive a pour but de présenter le PGA Tour et de se familiariser avec le vocabulaire du golf. Les objectifs et les librairies utilisées sont ensuite exposés.

1.1 PGA Tour

Le PGA Tour (Wikipédia, s. d.b) (Professional Golfers’ Association) est un circuit de golf professionnel qui se déroule aux États-Unis. Il a été créé en 1968 et est considéré comme le meilleur circuit au monde. Il compte aujourd’hui plus de 50 tournois officiels par saison.

La FedExCup (Wikipédia, s. d.a) est un championnat du PGA Tour qui dure toute la saison. Durant chaque tournoi, les golfeurs remportent des points qui comptent pour ce championnat en fonction de leurs résultats à l’évènement. À la fin de la saison, les 30 meilleurs golfeurs du classement s’affrontent pour remporter le trophée.

Un tournoi du PGA Tour se déroule sur 4 jours en 4 rounds. Les joueurs jouent le parcours une fois par jour. Après les 2 premiers jours, seuls les meilleurs se qualifient pour les 2 derniers rounds. On dit qu’ils ont passé le cut. Ceux qui n’arrivent pas à se qualifier dans le cut ne touchent pas d’argent pour leur participation au tournoi.

1.2 Initiation au golf

Voici la définition du golf d’après le Larousse :

Sport consistant à envoyer une balle à l’aide de clubs dans les 18 trous successifs d’un vaste terrain avec le minimum de coups.

Pour les non-initiés, commençons par donner quelques définitions nécessaires à la compréhension de l’analyse (Wikipédia, s. d.c; D. Chessel 2003) :

  • Par : nombre de coups théorique fixé pour un trou (ou un parcours). Le par d’un trou (3, 4 ou 5) dépend de sa longueur. Le par d’un parcours est égal à la somme des pars des trous du parcours (généralement entre 70 et 73).
  • Birdie : trou joué 1 coup sous le par (par 4 en 3 coups par exemple).
  • Drive : coup joué avec le driver. Le driver est le club qui permet de faire les coups les plus longs en distance. Les joueurs utilisent le driver pour jouer le premier coup des pars 5 et souvent celui des pars 4. Il n’est pas utilisé sur les pars 3 car la balle dépasserait le trou.
  • Fairway : zone herbeuse et bien tondue séparant le départ du green, où il est préférable d’envoyer la balle.
  • Rough : partie du parcours de golf longeant les trous et placée sur les côtés du fairway. L’herbe y est plus haute. Elle sanctionne un mauvais drive car il est plus difficile de jouer dans le rough que sur le fairway.
  • Green : le green est la zone du terrain qui entoure le drapeau. Elle est parfaitement tondue pour permettre au joueur de putter.
  • Putt : coup exécuté sur ou très proche du green en faisant rouler la balle. C’est le coup de golf qui nécessite le plus d’application et de précision car il faut étudier la pente du green et contrôler précisément la force du coup.
  • Bunker : un bunker est une zone de sable censée gêner l’approche du green. Il existe aussi des bunkers de parcours.

La figure 1.1 montre les différentes aires de jeu décrites.

Schéma des aires de jeu au golf

Figure 1.1: Schéma des aires de jeu au golf

1.3 Objectifs

L’analyse de données va répondre aux questions suivantes :

  1. Comment les différentes statistiques sont-elles distribuées ?
  2. Y a-t-il une évolution des performances des golfeurs sur ces 10 dernières années ?
  3. Quels sont les pays les plus représentés au PGA Tour ?
  4. Quels sont les liens entre les variables ?
  5. Quelles sont les statistiques qui permettent d’expliquer le score moyen d’un golfeur ?
  6. Y a-t-il différents profils de golfeurs ?
  7. Peut-on déterminer si un golfeur a gagné un tournoi ou non en fonction des autres statistiques ?

Pour répondre à ces questions, les analyses et méthodes suivantes sont mises en oeuvre :

  • Analyse descriptive des données
  • Modèles de régression linéaire
  • Analyse en composantes principales
  • Classification non supervisée
  • Classification supervisée

1.4 Librairies

L’analyse utilise les librairies suivantes :

  • ggplot2 (STHDA, s. d.g, s. d.d, s. d.f; Chang, s. d.) : génération d’histogrammes, de boxplots et diagrammes en bâtons
  • gridExtra (STHDA, s. d.e) : mise en grille des graphes générés avec ggplot2
  • ggcorrplot (STHDA, s. d.b) : matrice de corrélation et heatmap
  • ggfortify (STHDA, s. d.c) : extension de ggplot2 pour l’analyse des résidus en regréssion linéaire
  • car : calcul du VIF pour l’analyse des résidus en régression linéaire
  • mclust : comparaison des clustering
  • cluster : clustering
  • FactoMineR (Husson 2008) : ACP et clustering
  • factoextra (STHDA, s. d.a, s. d.h) : extension de FactoMineR pour la visualisation des analyses
  • MASS : méthodes LDA et QDA
  • pROC : courbes ROC
  • klaR : sélection de modèle pour les méthodes LDA et QDA
  • rpart : arbres de classification
  • randomForest : algorithme de Random forest

2 Jeu de données

Afin de travailler avec des données récentes et de taille raisonnable, j’ai créé mon propre jeu de données grâce à un script en Python effectuant du web scraping sur le site officiel du PGA Tour (« PGA Tour - Official Home of Golf », s. d.). Ce site renseigne les statistiques des golfeurs de façon très détaillée.

Le jeu de données est constitué de 1100 lignes et 15 colonnes. Chaque ligne correspond aux statistiques d’un golfeur du PGA Tour sur une année. Dans ce jeu de données, on retrouve les 100 meilleurs golfeurs au classement de la FedExCup pour toutes les saisons entre 2010 et 2020. On peut alors trouver le même golfeur sur plusieurs lignes différentes avec ses statistiques correspondantes à ses différentes saisons.

On importe les données contenues dans le fichier pgaTour.csv.

pga <- read.csv("../data/pgaTour.csv")
head(pga)
##      PLAYER.NAME       COUNTRY POINTS WINS ROUNDS SCORING    DD    DA   GIR
## 1      Ernie Els  South Africa   1846    2     72  69.843 288.4 60.16 67.86
## 2 Steve Stricker United States   1697    2     73  69.660 282.9 68.50 68.29
## 3      Jim Furyk United States   1691    2     76  69.828 276.0 71.01 67.12
## 4 Phil Mickelson United States   1629    1     76  69.966 299.1 52.66 65.13
## 5    Justin Rose       England   1593    2     78  69.885 287.8 65.17 66.31
## 6   Jeff Overton United States   1536    0     89  70.497 297.3 55.44 66.86
##      SS SCRAMBLING    PA   BA   MONEY YEAR
## 1 48.44      61.23 1.768 3.79 4558861 2010
## 2 54.84      63.75 1.746 3.89 4190235 2010
## 3 48.04      63.47 1.767 3.64 4809622 2010
## 4 53.62      61.84 1.762 3.89 3821733 2010
## 5 58.78      63.85 1.760 3.87 3603331 2010
## 6 51.09      56.80 1.753 3.76 3456356 2010

Voici la liste des variables étudiées :

  • PLAYER NAME : nom du joueur
  • COUNTRY : nationalité du joueur
  • POINTS : les points cumulés par le joueur pour la FedExCup pour l’année indiquée
  • WINS : nombre de victoires
  • ROUNDS : nombre de tours joués
  • SCORING : score moyen par tour joué
  • DD (Driving Distance) : distance moyenne des drives (en yards)
  • DA (Driving Accuracy) : pourcentage de coups joués au départ qui ont atteint le fairway
  • GIR (Greens In Regulation) : pourcentage de greens touchés en régulation, c’est-à-dire atteints au bout de 1 coup pour les pars 3, 2 pour les pars 4, 3 pour les pars 5
  • SS (Sand Saves) : pourcentage de pars effectués alors que le joueur était dans le bunker
  • SCRAMBLING : pourcentage de fois où un joueur manque le GIR mais effectue quand même le par ou mieux
  • PA (Putting Average) : nombre moyen de putts joués par trou
  • BA (Birdie Average) : nombre de birdies divisé par le nombre de rounds joués
  • MONEY : argent officiel remporté durant les tournois (en $)
  • YEAR : année

On peut classer ces variables en trois groupes :

  • Les variables qualitatives : PLAYER NAME, COUNTRY, YEAR
  • Les variables quantitatives mesurant les performances : ROUNDS, SCORING, DD, DA, GIR, SS, SCRAMBLING, PA, BA
  • Les variables quantitatives indiquant les récompenses : POINTS, WINS, MONEY

On transforme la variable quantitative YEAR en variable qualitative. Cela permet d’étudier les statistiques des golfeurs en fonction de l’année.

pga$PLAYER.NAME <- as.character(pga$PLAYER.NAME)
pga$YEAR <- as.factor(pga$YEAR)
str(pga)
## 'data.frame':    1100 obs. of  15 variables:
##  $ PLAYER.NAME: chr  "Ernie Els" "Steve Stricker" "Jim Furyk" "Phil Mickelson" ...
##  $ COUNTRY    : Factor w/ 29 levels "Argentina","Australia",..: 23 27 27 27 9 27 27 27 27 23 ...
##  $ POINTS     : num  1846 1697 1691 1629 1593 ...
##  $ WINS       : int  2 2 2 1 2 0 2 1 0 1 ...
##  $ ROUNDS     : int  72 73 76 76 78 89 85 77 97 88 ...
##  $ SCORING    : num  69.8 69.7 69.8 70 69.9 ...
##  $ DD         : num  288 283 276 299 288 ...
##  $ DA         : num  60.2 68.5 71 52.7 65.2 ...
##  $ GIR        : num  67.9 68.3 67.1 65.1 66.3 ...
##  $ SS         : num  48.4 54.8 48 53.6 58.8 ...
##  $ SCRAMBLING : num  61.2 63.8 63.5 61.8 63.9 ...
##  $ PA         : num  1.77 1.75 1.77 1.76 1.76 ...
##  $ BA         : num  3.79 3.89 3.64 3.89 3.87 3.76 3.54 3.94 3.97 3.52 ...
##  $ MONEY      : num  4558861 4190235 4809622 3821733 3603331 ...
##  $ YEAR       : Factor w/ 11 levels "2010","2011",..: 1 1 1 1 1 1 1 1 1 1 ...

Voici un résumé des variables :

summary(pga)
##  PLAYER.NAME                      COUNTRY        POINTS            WINS       
##  Length:1100        United States     :733   Min.   : 352.0   Min.   :0.0000  
##  Class :character   Australia         : 70   1st Qu.: 584.8   1st Qu.:0.0000  
##  Mode  :character   England           : 48   Median : 745.5   Median :0.0000  
##                     South Africa      : 33   Mean   : 889.2   Mean   :0.3764  
##                     Korea, Republic of: 32   3rd Qu.:1068.8   3rd Qu.:1.0000  
##                     Sweden            : 28   Max.   :4169.0   Max.   :5.0000  
##                     (Other)           :156                                    
##      ROUNDS         SCORING            DD              DA       
##  Min.   : 38.0   Min.   :68.70   Min.   :266.4   Min.   :46.99  
##  1st Qu.: 73.0   1st Qu.:70.21   1st Qu.:287.4   1st Qu.:58.49  
##  Median : 83.0   Median :70.58   Median :293.1   Median :61.93  
##  Mean   : 81.8   Mean   :70.52   Mean   :293.4   Mean   :61.92  
##  3rd Qu.: 91.0   3rd Qu.:70.87   3rd Qu.:299.2   3rd Qu.:65.20  
##  Max.   :122.0   Max.   :72.42   Max.   :322.1   Max.   :76.88  
##                                                                 
##       GIR              SS          SCRAMBLING          PA              BA      
##  Min.   :56.68   Min.   :31.48   Min.   :48.30   Min.   :1.691   Min.   :2.89  
##  1st Qu.:64.89   1st Qu.:46.58   1st Qu.:56.98   1st Qu.:1.749   1st Qu.:3.50  
##  Median :66.52   Median :50.55   Median :59.10   Median :1.766   Median :3.68  
##  Mean   :66.43   Mean   :50.69   Mean   :59.15   Mean   :1.765   Mean   :3.69  
##  3rd Qu.:68.16   3rd Qu.:54.61   3rd Qu.:61.36   3rd Qu.:1.782   3rd Qu.:3.85  
##  Max.   :73.52   Max.   :68.55   Max.   :69.33   Max.   :1.841   Max.   :4.71  
##                                                                                
##      MONEY               YEAR    
##  Min.   :  423907   2010   :100  
##  1st Qu.: 1277580   2011   :100  
##  Median : 1799712   2012   :100  
##  Mean   : 2269694   2013   :100  
##  3rd Qu.: 2841505   2014   :100  
##  Max.   :12030465   2015   :100  
##                     (Other):500

Voici le nombre de joueurs pour chaque année :

table(pga$YEAR)
## 
## 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 
##  100  100  100  100  100  100  100  100  100  100  100

3 Analyse descriptive

Pour commencer, on effectue une analyse descriptive afin de visualiser les données par des diagrammes et des histogrammes. Dans cette partie, les trois premières questions sont abordées :

  1. Comment les différentes statistiques sont-elles distribuées ?
  2. Y a-t-il une évolution des performances des golfeurs sur ces 10 dernières années ?
  3. Quels sont les pays les plus représentés au PGA Tour ?

3.1 Distribution des variables

On étudie d’abord les variables quantitatives du jeu de données.

3.1.1 Histogrammes

On réalise un histogramme (STHDA, s. d.g) pour chacune de ces variables, en superposant l’estimation de la densité afin de voir comment elles sont distribuées.

bins <- 30

h1 <- ggplot(pga, aes(x = POINTS)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of POINTS") +
        theme_minimal()
h1

h2 <- ggplot(pga, aes(x = ROUNDS)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of ROUNDS") +
        theme_minimal()
h2

h3 <- ggplot(pga, aes(x = SCORING)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of SCORING") +
        theme_minimal()
h3

h4 <- ggplot(pga, aes(x = DD)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of DD") +
        theme_minimal()
h4

h5 <- ggplot(pga, aes(x = DA)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of DA") +
        theme_minimal()
h5

h6 <- ggplot(pga, aes(x = GIR)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of GIR") +
        theme_minimal()
h6

h7 <- ggplot(pga, aes(x = SS)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of SS") +
        theme_minimal()
h7

h8 <- ggplot(pga, aes(x = SCRAMBLING)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of SCRAMBLING") +
        theme_minimal()
h8

h9 <- ggplot(pga, aes(x = PA)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of PA") +
        theme_minimal()
h9

h10 <- ggplot(pga, aes(x = BA)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of BA") +
        theme_minimal()
h10

h11 <- ggplot(pga, aes(x = MONEY)) +
        geom_histogram(aes(y = ..density..), bins = bins, color = "darkblue",
                       fill = "lightblue") +
        geom_density(alpha = .2, fill = "#FF6666") +
        ggtitle("Histogram of MONEY") +
        theme_minimal()
h11
HistogrammesHistogrammesHistogrammesHistogrammesHistogrammesHistogrammesHistogrammesHistogrammesHistogrammesHistogrammesHistogrammes

Figure 3.1: Histogrammes

D’après les histogrammes de la figure 3.1, il est clair que les variables POINTS et MONEY ne sont pas distribuées normalement.

On trace le Q-Q plot des variables et on effectue des tests de normalité de Shapiro-Wilk.

qq1 <- ggplot(pga, aes(sample = POINTS)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of POINTS")
qq1

qq2 <- ggplot(pga, aes(sample = ROUNDS)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of ROUNDS")
qq2

qq3 <- ggplot(pga, aes(sample = SCORING)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of SCORING")
qq3

qq4 <- ggplot(pga, aes(sample = DD)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of DD")
qq4

qq5 <- ggplot(pga, aes(sample = DA)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of DA")
qq5

qq6 <- ggplot(pga, aes(sample = GIR)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of GIR")
qq6

qq7 <- ggplot(pga, aes(sample = SS)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of SS")
qq7

qq8 <- ggplot(pga, aes(sample = SCRAMBLING)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of SCRAMBLING")
qq8

qq9 <- ggplot(pga, aes(sample = PA)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of PA")
qq9

qq10 <- ggplot(pga, aes(sample = BA)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of BA")
qq10

qq11 <- ggplot(pga, aes(sample = MONEY)) +
        stat_qq() + stat_qq_line(color = "red") +
        ggtitle("Q-Q Plot of MONEY")
qq11
Q-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plotsQ-Q plots

Figure 3.2: Q-Q plots

shapiro.test(pga$POINTS)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$POINTS
## W = 0.83667, p-value < 2.2e-16
shapiro.test(pga$ROUNDS)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$ROUNDS
## W = 0.99156, p-value = 6.123e-06
shapiro.test(pga$SCORING)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$SCORING
## W = 0.98503, p-value = 3.359e-09
shapiro.test(pga$DD)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$DD
## W = 0.99654, p-value = 0.01564
shapiro.test(pga$DA)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$DA
## W = 0.9974, p-value = 0.07436
shapiro.test(pga$GIR)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$GIR
## W = 0.99638, p-value = 0.01164
shapiro.test(pga$SS)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$SS
## W = 0.9987, p-value = 0.6125
shapiro.test(pga$SCRAMBLING)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$SCRAMBLING
## W = 0.99853, p-value = 0.4888
shapiro.test(pga$PA)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$PA
## W = 0.99767, p-value = 0.1194
shapiro.test(pga$BA)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$BA
## W = 0.98937, p-value = 3.808e-07
shapiro.test(pga$MONEY)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga$MONEY
## W = 0.80708, p-value < 2.2e-16

On conclut que les variables DA, SS, SCRAMBLING et PA semblent être distribuées normalement. En revanche, les variables POINTS, ROUNDS, SCORING, DD, GIR, BA et MONEY ne sont pas distribuées normalement (au risque de se tromper de 5%).

Interprétations :

  • On observe que les histogrammes des variables POINTS et MONEY penchent fortement vers la gauche. Cela signifie que les meilleurs golfeurs gagnent beaucoup plus de points et de gains que les autres. En effet, les récompenses des tournois sont beaucoup plus généreuses pour les golfeurs du top 10 que pour les autres, ce qui explique ce phénomène.
  • Les histogrammes des autres variables témoignent d’un niveau assez similaire entre les golfeurs avec quelques disparités.

3.1.2 Boxplots

On regarde maintenant les boxplots de ces variables par année, en superposant en rouge la moyenne. Cela permet d’étudier l’évolution des performances des golfeurs sur ces 10 dernières années.

b1 <- ggplot(pga, aes(x = YEAR, y = POINTS)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of POINTS") +
        theme_minimal() +
        theme(legend.position = "none")
b1

b2 <- ggplot(pga, aes(x = YEAR, y = ROUNDS)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of ROUNDS") +
        theme_minimal() +
        theme(legend.position = "none")
b2

b3 <- ggplot(pga, aes(x = YEAR, y = SCORING)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of SCORING") +
        theme_minimal() +
        theme(legend.position = "none")
b3

b4 <- ggplot(pga, aes(x = YEAR, y = DD)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of DD") +
        theme_minimal() +
        theme(legend.position = "none")
b4

b5 <- ggplot(pga, aes(x = YEAR, y = DA)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of DA") +
        theme_minimal() +
        theme(legend.position = "none")
b5

b6 <- ggplot(pga, aes(x = YEAR, y = GIR)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of GIR") +
        theme_minimal() +
        theme(legend.position = "none")
b6

b7 <- ggplot(pga, aes(x = YEAR, y = SS)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of SS") +
        theme_minimal() +
        theme(legend.position = "none")
b7

b8 <- ggplot(pga, aes(x = YEAR, y = SCRAMBLING)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of SCRAMBLING") +
        theme_minimal() +
        theme(legend.position = "none")
b8

b9 <- ggplot(pga, aes(x = YEAR, y = PA)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of PA") +
        theme_minimal() +
        theme(legend.position = "none")
b9

b10 <- ggplot(pga, aes(x = YEAR, y = BA)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of BA") +
        theme_minimal() +
        theme(legend.position = "none")
b10

b11 <- ggplot(pga, aes(x = YEAR, y = MONEY)) +
        geom_boxplot(aes(fill = YEAR)) +
        stat_summary(fun = mean, geom = "point", size = 2,
                     color = "red", fill = "red") +
        ggtitle("Boxplot of MONEY") +
        theme_minimal() +
        theme(legend.position = "none")
b11
BoxplotsBoxplotsBoxplotsBoxplotsBoxplotsBoxplotsBoxplotsBoxplotsBoxplotsBoxplotsBoxplots

Figure 3.3: Boxplots

Étudions les boxplots où on observe des phénomènes intéressants.

3.1.2.1 Points à la FedExCup

On regarde les boxplots de la variable POINTS et on affiche les 10 meilleurs golfeurs en fonction de leur nombre de points au classement de la FedExCup pour les saisons 2010, 2015, 2018 et 2020.

b1

# Top 10 de la FedExCup pour chaque année
for (year in c("2010", "2015", "2018", "2020")) {
  df.points <- data.frame(NAMES = pga$PLAYER.NAME[pga$YEAR == year],
                          POINTS = pga$POINTS[pga$YEAR == year])
  df.points$NAMES <- df.points$NAMES[order(df.points$POINTS,
                                           decreasing = TRUE)]
  df.points$POINTS <- df.points$POINTS[order(df.points$POINTS,
                                             decreasing = TRUE)]
  print(head(df.points, 10))
}
##             NAMES POINTS
## 1       Ernie Els   1846
## 2  Steve Stricker   1697
## 3       Jim Furyk   1691
## 4  Phil Mickelson   1629
## 5     Justin Rose   1593
## 6    Jeff Overton   1536
## 7    Hunter Mahan   1528
## 8    Bubba Watson   1498
## 9     Matt Kuchar   1437
## 10      Tim Clark   1409
##             NAMES POINTS
## 1   Jordan Spieth   4169
## 2       Jason Day   2459
## 3    Bubba Watson   2407
## 4    Jimmy Walker   2014
## 5     Justin Rose   1742
## 6    Robert Streb   1720
## 7  Dustin Johnson   1718
## 8    Patrick Reed   1593
## 9       Danny Lee   1561
## 10   Zach Johnson   1559
##                 NAMES POINTS
## 1      Dustin Johnson   2717
## 2       Justin Thomas   2634
## 3       Brooks Koepka   2012
## 4         Justin Rose   1991
## 5        Bubba Watson   1879
## 6           Jason Day   1771
## 7        Webb Simpson   1710
## 8  Francesco Molinari   1682
## 9   Bryson DeChambeau   1617
## 10       Patrick Reed   1555
##                NAMES POINTS
## 1      Justin Thomas   2458
## 2    Collin Morikawa   1902
## 3       Webb Simpson   1878
## 4  Bryson DeChambeau   1657
## 5         Sungjae Im   1633
## 6       Patrick Reed   1426
## 7      Daniel Berger   1347
## 8       Rory McIlroy   1327
## 9       Brendon Todd   1316
## 10          Jon Rahm   1295

On remarque que la moyenne (en rouge) est toujours au dessus de la médiane. En effet, comme le montre les top 10, on voit que les meilleurs golfeurs ont une marge de points assez importante sur les autres. Par exemple, en 2015, Jordan Spieth a terminé en tête de ce classement avec 4169 points, soit 1710 points de plus que le deuxième, Jason Day, et 2155 points de plus que le quatrième.

On peut de plus remarquer que le nombre de points est moins élevé en 2020. Ceci est expliqué par la crise du COVID-19 qui a fait annuler des tournois.

3.1.2.2 Rounds

On affiche le nombre moyen de rounds joués par les golfeurs pour chaque année étudiée.

b2

tapply(pga$ROUNDS, pga$YEAR, mean)
##  2010  2011  2012  2013  2014  2015  2016  2017  2018  2019  2020 
## 84.59 82.44 81.74 77.44 83.94 87.46 82.84 84.80 86.04 80.53 68.00

On voit que le nombre de rounds disputés par les golfeurs est beaucoup moins élevé en 2020 que les autres années. En effet, les joueurs ont disputé en moyenne 68 rounds en 2020 contre plus de 80 sur les autres années. Ceci est encore expliqué par les annulations de tournois dues au COVID-19.

3.1.2.3 Driving distance

On affiche la distance moyenne des drives pour chaque année étudiée.

b4

tapply(pga$DD, pga$YEAR, mean)
##    2010    2011    2012    2013    2014    2015    2016    2017    2018    2019 
## 287.880 293.002 291.611 288.409 291.476 292.614 292.527 296.308 299.018 296.547 
##    2020 
## 298.265

Il est intéressant de voir que la distance moyenne des drives a nettement augmenté au fil des années. En effet, les drives ont gagné en moyenne 10 yards sur ces 10 dernières années.

3.1.2.3.1 Putting average

On affiche le nombre moyen de putts joués par trou pour chaque année étudiée.

b9

tapply(pga$PA, pga$YEAR, mean)
##    2010    2011    2012    2013    2014    2015    2016    2017    2018    2019 
## 1.77278 1.77005 1.77093 1.77343 1.76471 1.76116 1.76741 1.76808 1.76546 1.75537 
##    2020 
## 1.75071

On voit que le nombre moyen de putts joués par trou a diminué. Il est passé de 1.77 en 2010 à 1.75 en 2020.

3.1.2.4 Birdie average

On affiche le nombre moyen de birdies divisé par le nombre de rounds joués pour chaque année étudiée.

b10

tapply(pga$BA, pga$YEAR, mean)
##   2010   2011   2012   2013   2014   2015   2016   2017   2018   2019   2020 
## 3.6083 3.6629 3.5771 3.5619 3.6072 3.7450 3.6485 3.6907 3.7612 3.8639 3.8669

On remarque que le nombre moyen de birdies divisé par le nombre de rounds joués a bien augmenté. Il est passé de 3.6 en 2010 à 3.9 en 2020.

3.1.2.5 Money

On considère les années 2010, 2015, 2018 et 2020. On affiche :

  • Les boxplots des gains remportés durant les tournois par année
  • Les écarts-type et les moyennes de gains par année
  • Les classements des salaires par année
b11

sqrt(tapply(pga$MONEY, pga$YEAR, var))
##    2010    2011    2012    2013    2014    2015    2016    2017    2018    2019 
## 1062532 1240897 1300093 1363190 1406883 1737810 1427342 1693366 1693924 1560821 
##    2020 
## 1295483
tapply(pga$MONEY, pga$YEAR, mean)
##    2010    2011    2012    2013    2014    2015    2016    2017    2018    2019 
## 1893241 2031651 2117782 2044258 2309034 2369100 2304901 2499168 2808228 2502853 
##    2020 
## 2086422
# Classement des salaires par année
for (year in c("2010", "2015", "2018", "2020")) {
  df.money <- data.frame(NAMES = pga$PLAYER.NAME[pga$YEAR == year],
                         MONEY = pga$MONEY[pga$YEAR == year])
  df.money$NAMES <- df.money$NAMES[order(df.money$MONEY,
                                         decreasing = TRUE)]
  df.money$MONEY <- df.money$MONEY[order(df.money$MONEY,
                                         decreasing = TRUE)]
  print(head(df.money))
}
##            NAMES   MONEY
## 1    Matt Kuchar 4910477
## 2      Jim Furyk 4809622
## 3      Ernie Els 4558861
## 4 Dustin Johnson 4473122
## 5 Steve Stricker 4190235
## 6 Phil Mickelson 3821733
##            NAMES    MONEY
## 1  Jordan Spieth 12030465
## 2      Jason Day  9403330
## 3   Bubba Watson  6876797
## 4  Rickie Fowler  5773430
## 5 Dustin Johnson  5509467
## 6    Justin Rose  5462677
##               NAMES   MONEY
## 1     Justin Thomas 8694821
## 2    Dustin Johnson 8457352
## 3       Justin Rose 8130678
## 4 Bryson DeChambeau 8094489
## 5     Brooks Koepka 7094047
## 6      Bubba Watson 5793748
##               NAMES   MONEY
## 1     Justin Thomas 7344040
## 2          Jon Rahm 5959819
## 3    Dustin Johnson 5837267
## 4   Collin Morikawa 5250868
## 5      Webb Simpson 5097742
## 6 Bryson DeChambeau 4998495

On observe que les gains moyens remportés par un golfeur sur une saison ont augmenté jusqu’en 2018 et ont ensuite diminué, ce qui est encore une des conséquences de la crise du COVID-19. En effet, moins de tournois entraîne moins de gains.

On note aussi de grandes différences de gains entre les meilleurs golfeurs et les autres. Par exemple, en 2015, les gains de Jordan Spieth s’élevaient à 12030465 dollars, ce qui représente plus du double des gains du quatrième, Rickie Fowler, avec 5773430 dollars.

3.1.2.6 Conclusions

Pour conclure sur l’interprétation des boxplots, on peut dire que :

  • Le niveau des golfeurs a augmenté durant ces dernières années.
  • La crise du COVID-19 a eu des conséquences sur le nombre de tournois disputés et les gains remportés par les golfeurs en 2020.
  • Il y a de grands écarts au niveau des points au classement de la FedExCup et au niveau des gains remportés entre les meilleurs golfeurs (top 5 par exemple) et les autres.

3.2 Nationalités des golfeurs

Pour terminer l’analyse descriptive, on étudie la variable qualitative COUNTRY afin de déterminer les pays les plus représentés au PGA Tour. Pour les années 2010, 2015, 2018 et 2020, on trie les pays du plus au moins représenté et on montre les 5 premiers pays (STHDA, s. d.d).

for (year in c("2010", "2015", "2018", "2020")) {
  # On trie les pays du plus représenté au moins représenté
  countries <- sort(table(pga$COUNTRY[pga$YEAR == year]), decreasing = TRUE)
  # On crée un dataframe en gardant les 5 pays les plus représentés
  df.countries <- data.frame(countries[1:5])
  # On renomme la colonne 1
  names(df.countries)[1] <- "Country"

  b <- ggplot(data = df.countries, aes(x = Country, y = Freq, fill = Country)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = Freq), hjust = 1.6, color = "white", size = 3.5) +
    ggtitle(paste0("Countries represented in top 100 players in PGA Tour ", year)) +
    coord_flip()
  print(b)
}
Pays représentés dans le top 100 du PGA TourPays représentés dans le top 100 du PGA TourPays représentés dans le top 100 du PGA TourPays représentés dans le top 100 du PGA Tour

Figure 3.4: Pays représentés dans le top 100 du PGA Tour

Les diagrammes en bâtons de la figure 3.4 montrent clairement que les États-Unis sont bien plus représentés que les autres pays. Après les États-Unis, l’Australie et l’Angleterre sont les deux pays qui comptent le plus de golfeurs dans le top 100 du PGA Tour pour chaque saison.

4 Relations entre les variables

Dans cette partie, nous allons étudier les relations entre les variables quantitatives du jeu de données. Nous allons donner de premiers éléments de réponse aux questions suivantes :

  1. Quels sont les liens entre les variables ?
  2. Quelles sont les statistiques qui permettent d’expliquer le score moyen d’un golfeur ?

On définit un dataframe avec seulement les variables quantitatives.

# Variables quantitatives
pga.num <- pga[3:14][-2]
str(pga.num)
## 'data.frame':    1100 obs. of  11 variables:
##  $ POINTS    : num  1846 1697 1691 1629 1593 ...
##  $ ROUNDS    : int  72 73 76 76 78 89 85 77 97 88 ...
##  $ SCORING   : num  69.8 69.7 69.8 70 69.9 ...
##  $ DD        : num  288 283 276 299 288 ...
##  $ DA        : num  60.2 68.5 71 52.7 65.2 ...
##  $ GIR       : num  67.9 68.3 67.1 65.1 66.3 ...
##  $ SS        : num  48.4 54.8 48 53.6 58.8 ...
##  $ SCRAMBLING: num  61.2 63.8 63.5 61.8 63.9 ...
##  $ PA        : num  1.77 1.75 1.77 1.76 1.76 ...
##  $ BA        : num  3.79 3.89 3.64 3.89 3.87 3.76 3.54 3.94 3.97 3.52 ...
##  $ MONEY     : num  4558861 4190235 4809622 3821733 3603331 ...

4.1 Matrice de corrélation

Pour identifier les corrélations entre les variables quantitatives, nous étudions la matrice de corrélation entre ces variables. On construit deux matrices (STHDA, s. d.b) :

  • La matrice de corrélation avec les coefficients de corrélation
  • La matrice de corrélation avec les p-values
# Compute a correlation matrix
pga.corr <- round(cor(pga.num), 2)
pga.corr
##            POINTS ROUNDS SCORING    DD    DA   GIR    SS SCRAMBLING    PA    BA
## POINTS       1.00   0.09   -0.66  0.22 -0.05  0.22  0.20       0.27 -0.37  0.49
## ROUNDS       0.09   1.00    0.11 -0.12  0.08  0.09 -0.04       0.07  0.01 -0.03
## SCORING     -0.66   0.11    1.00 -0.24 -0.05 -0.41 -0.31      -0.46  0.43 -0.58
## DD           0.22  -0.12   -0.24  1.00 -0.61  0.30 -0.17      -0.23 -0.03  0.48
## DA          -0.05   0.08   -0.05 -0.61  1.00  0.26  0.04       0.20  0.09 -0.16
## GIR          0.22   0.09   -0.41  0.30  0.26  1.00 -0.14      -0.02  0.15  0.38
## SS           0.20  -0.04   -0.31 -0.17  0.04 -0.14  1.00       0.58 -0.32  0.14
## SCRAMBLING   0.27   0.07   -0.46 -0.23  0.20 -0.02  0.58       1.00 -0.35  0.14
## PA          -0.37   0.01    0.43 -0.03  0.09  0.15 -0.32      -0.35  1.00 -0.68
## BA           0.49  -0.03   -0.58  0.48 -0.16  0.38  0.14       0.14 -0.68  1.00
## MONEY        0.93  -0.01   -0.71  0.31 -0.06  0.27  0.22       0.28 -0.40  0.57
##            MONEY
## POINTS      0.93
## ROUNDS     -0.01
## SCORING    -0.71
## DD          0.31
## DA         -0.06
## GIR         0.27
## SS          0.22
## SCRAMBLING  0.28
## PA         -0.40
## BA          0.57
## MONEY       1.00
# Compute a matrix of correlation p-values
pga.corr.pmat <- cor_pmat(pga.num)
pga.corr.pmat
##                   POINTS       ROUNDS       SCORING            DD            DA
## POINTS      0.000000e+00 2.421364e-03 4.486332e-138  2.514529e-13  8.348631e-02
## ROUNDS      2.421364e-03 0.000000e+00  1.929372e-04  8.231272e-05  7.265872e-03
## SCORING    4.486332e-138 1.929372e-04  0.000000e+00  5.930619e-16  7.938390e-02
## DD          2.514529e-13 8.231272e-05  5.930619e-16  0.000000e+00 1.963491e-111
## DA          8.348631e-02 7.265872e-03  7.938390e-02 1.963491e-111  0.000000e+00
## GIR         4.472659e-13 3.610021e-03  2.776172e-46  5.257056e-24  7.054084e-19
## SS          9.692371e-12 2.364224e-01  1.710130e-25  2.517885e-08  2.308852e-01
## SCRAMBLING  2.713168e-19 1.744482e-02  7.809477e-59  4.593441e-15  3.420000e-11
## PA          2.331310e-36 7.333580e-01  2.140805e-50  2.724850e-01  2.695786e-03
## BA          2.012598e-68 2.528169e-01 5.172684e-100  7.024429e-66  5.663322e-08
## MONEY       0.000000e+00 6.629974e-01 1.257115e-172  5.149349e-26  5.994701e-02
##                     GIR           SS   SCRAMBLING            PA            BA
## POINTS     4.472659e-13 9.692371e-12 2.713168e-19  2.331310e-36  2.012598e-68
## ROUNDS     3.610021e-03 2.364224e-01 1.744482e-02  7.333580e-01  2.528169e-01
## SCORING    2.776172e-46 1.710130e-25 7.809477e-59  2.140805e-50 5.172684e-100
## DD         5.257056e-24 2.517885e-08 4.593441e-15  2.724850e-01  7.024429e-66
## DA         7.054084e-19 2.308852e-01 3.420000e-11  2.695786e-03  5.663322e-08
## GIR        0.000000e+00 6.066047e-06 4.404226e-01  4.384295e-07  1.392748e-38
## SS         6.066047e-06 0.000000e+00 5.932967e-98  2.191031e-28  2.457553e-06
## SCRAMBLING 4.404226e-01 5.932967e-98 0.000000e+00  7.199359e-33  1.877694e-06
## PA         4.384295e-07 2.191031e-28 7.199359e-33  0.000000e+00 2.577555e-147
## BA         1.392748e-38 2.457553e-06 1.877694e-06 2.577555e-147  0.000000e+00
## MONEY      4.249159e-20 5.461786e-14 2.049932e-21  2.676516e-43  1.862827e-97
##                    MONEY
## POINTS      0.000000e+00
## ROUNDS      6.629974e-01
## SCORING    1.257115e-172
## DD          5.149349e-26
## DA          5.994701e-02
## GIR         4.249159e-20
## SS          5.461786e-14
## SCRAMBLING  2.049932e-21
## PA          2.676516e-43
## BA          1.862827e-97
## MONEY       0.000000e+00

Afin de visualiser facilement les liens entre les variables, on réalise un heatmap (STHDA, s. d.b) de la matrice de corrélation montré par la figure 4.1.

ggcorrplot(pga.corr, hc.order = TRUE, type = "lower", lab = TRUE,
           lab_size = 2.5, tl.cex = 8,
           title = "PGA Tour Correlation matrix heatmap")
Heatmap de la matrice de corrélation

Figure 4.1: Heatmap de la matrice de corrélation

Le heatmap montre les corrélations entre les variables numériques du jeu de données. Les cases en rouge indiquent une corrélation positive, tandis que les cases en bleu indiquent une corrélation négative.

Voici les liens importants détectés :

  • POINTS et MONEY (0.93) : plus un golfeur gagne de points à la FedExCup, plus il gagne d’argent.
  • SCRAMBING et SS (0.58) : les Sand Saves et les Scrambling sont deux situations difficiles pour un golfeur. On en déduit que si un golfeur arrive à gérer l’une de ces deux situations, alors il gère bien l’autre aussi.
  • BA et MONEY (0.57) : les birdies font rapporter plus de gains aux golfeurs, ce qui est logique car cela améliore leurs scores.
  • SCORING et MONEY (-0.71) : plus le score moyen est faible, plus le golfeur gagne d’argent. En effet, quand un golfeur a un score plus faible que les autres, il se situe dans les premières places du classement d’un tournoi, ce qui lui rapporte plus de gains.
  • PA et BA (-0.68) : un golfeur qui réalise moins de putts par trou réalise plus de birdies. En effet, en cas d’échec à un putt sur un trou, il est plus difficile de réaliser un birdie. Dans le cas où le golfeur est en GIR, il est obligé de réussir le putt du premier coup s’il veut faire un birdie.
  • SCORING et POINTS (-0.66) : plus le score moyen d’un golfeur est faible, plus il remporte de points pour la FedExCup. Cela est logique car un golfeur est mieux classé dans les tournois s’il réalise des scores plus faibles que ses adversaires.
  • DA et DD (-0.61) : la balle arrive moins souvent sur le fairway lorsque les drives sont plus longs en moyenne.
  • SCORING et BA (-0.58) : plus le score est faible et plus le nombre moyen de birdies est élevé. En effet, un birdie diminue le score d’un golfeur, ce qui explique cette corrélation. On rappelle qu’au golf, l’objectif est d’obtenir le score le plus faible.

4.2 Régression linéaire

À partir des liaisons entre variables identifiées, nous allons maintenant mettre en oeuvre quelques modèles de régressions.

4.2.1 Régression linéaire simple

Nous allons d’abord montrer des modèles de régressions linéaires simples (STHDA, s. d.f) adaptés au jeu de données.

4.2.1.1 Explication des gains remportés

On peut expliquer les gains remportés durant les tournois en fonction du nombre de points au classement de la FedExCup.

# Linear regression
pga.lm.res <- lm(pga$MONEY ~ pga$POINTS)
summary(pga.lm.res)
## 
## Call:
## lm(formula = pga$MONEY ~ pga$POINTS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1448429  -321854   -73607   217055  3571338 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -483448.00   36689.92  -13.18   <2e-16 ***
## pga$POINTS     3096.23      36.99   83.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 539500 on 1098 degrees of freedom
## Multiple R-squared:  0.8645, Adjusted R-squared:  0.8644 
## F-statistic:  7008 on 1 and 1098 DF,  p-value: < 2.2e-16
# Scatter plot with regression line
ggplot(pga, aes(x = POINTS, y = MONEY)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = lm) +
  ggtitle("Gains remportés en fonction du nombre de points à la FedExCup")
Gains remportés en fonction du nombre de points à la FedExCup

Figure 4.2: Gains remportés en fonction du nombre de points à la FedExCup

Comme le montre la figure 4.2, les meilleurs golfeurs, c’est-à-dire ceux qui ont un nombre de points élevé, sont ceux qui gagnent le plus d’argent, ce qui paraît juste.

4.2.1.2 Explication du pourcentage de birdies

On peut expliquer le nombre moyen de birdies (BA) par le nombre moyen de putts (PA).

# Linear regression
pga.lm2.res <- lm(pga$BA ~ pga$PA)
summary(pga.lm2.res)
## 
## Call:
## lm(formula = pga$BA ~ pga$PA)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69272 -0.13995 -0.00095  0.13482  0.63900 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  17.5337     0.4564   38.42   <2e-16 ***
## pga$PA       -7.8412     0.2585  -30.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2041 on 1098 degrees of freedom
## Multiple R-squared:  0.456,  Adjusted R-squared:  0.4555 
## F-statistic: 920.2 on 1 and 1098 DF,  p-value: < 2.2e-16
# Scatter plot with regression line
ggplot(pga, aes(x = PA, y = BA)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = lm) +
  ggtitle("Nombre moyen de birdies en fonction du nombre moyen de putts")
Nombre moyen de birdies en fonction du nombre moyen de putts

Figure 4.3: Nombre moyen de birdies en fonction du nombre moyen de putts

En effet, d’après la figure 4.3, on en déduit que moins un golfeur joue de putts par trou, plus il a de chance de réaliser un birdie.

4.2.1.3 Explication de la précision des drives

On peut expliquer la précision des drives (DA) par la longueur des drives (DD).

# Linear regression
pga.lm3.res <- lm(pga$DA ~ pga$DD)
summary(pga.lm3.res)
## 
## Call:
## lm(formula = pga$DA ~ pga$DD)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6859  -2.4589   0.1084   2.6847  12.6497 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 155.86676    3.71925   41.91   <2e-16 ***
## pga$DD       -0.32016    0.01267  -25.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.841 on 1098 degrees of freedom
## Multiple R-squared:  0.3677, Adjusted R-squared:  0.3672 
## F-statistic: 638.6 on 1 and 1098 DF,  p-value: < 2.2e-16
# Scatter plot with regression line
ggplot(pga, aes(x = DD, y = DA)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = lm) +
  ggtitle("Précision moyenne des drives en fonction de leur longueur moyenne")
Précision moyenne des drives en fonction de leur longueur moyenne

Figure 4.4: Précision moyenne des drives en fonction de leur longueur moyenne

En effet, d’après la figure 4.4, on en déduit que les golfeurs qui réalisent de longs drives ont plus de chance de rater le fairway que ceux qui préfèrent jouer plus court pour assurer leur coup.

4.2.2 Régression linéaire multiple

Dans cette partie, on veut comprendre quelles variables influencent le score moyen d’un golfeur. On cherche à expliquer la variable SCORING en fonction d’autres variables quantitatives grâce à une régression linéaire multiple.

On considère uniquement les variables en rapport avec les performances d’un golfeur sur le terrain : ROUNDS, DD, DA, GIR, SS, SCRAMBLING, PA et BA.

pga.lm.scoring <- lm(SCORING ~ ROUNDS + DD + DA + GIR + SS + SCRAMBLING + PA + BA, pga)

On vérifie la colinéarité entre les variables.

vif(pga.lm.scoring)
##     ROUNDS         DD         DA        GIR         SS SCRAMBLING         PA 
##   1.050025   3.274633   2.452101   2.473770   1.594171   1.704819   3.898243 
##         BA 
##   5.064998

Les VIF sont inférieurs à 10 donc il n’y a pas de problème de colinéarité entre les variables. On effectue une analyse des résidus (STHDA, s. d.c).

# Diagnostic plots for Linear Model
autoplot(pga.lm.scoring)
## Warning: `arrange_()` was deprecated in dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
Analyse des résidus

Figure 4.5: Analyse des résidus

Aucune observation ne semble perturber le modèle. On vérifie la normalité des résidus avec un test de Shapiro-Wilk.

shapiro.test(pga.lm.scoring$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  pga.lm.scoring$residuals
## W = 0.99888, p-value = 0.7385

La p-value est supérieure à 5% donc il semble que les résidus sont gaussiens. Voici le résumé de la régression linéaire.

summary(pga.lm.scoring)
## 
## Call:
## lm(formula = SCORING ~ ROUNDS + DD + DA + GIR + SS + SCRAMBLING + 
##     PA + BA, data = pga)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06782 -0.24015 -0.00005  0.23005  1.08998 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 74.4985266  1.6070841  46.356  < 2e-16 ***
## ROUNDS       0.0059704  0.0008223   7.260 7.32e-13 ***
## DD          -0.0084875  0.0021390  -3.968 7.72e-05 ***
## DA          -0.0027733  0.0035058  -0.791 0.429087    
## GIR         -0.0838955  0.0070439 -11.910  < 2e-16 ***
## SS          -0.0063163  0.0022826  -2.767 0.005752 ** 
## SCRAMBLING  -0.0637172  0.0044302 -14.383  < 2e-16 ***
## PA           5.1214353  0.8958820   5.717 1.40e-08 ***
## BA          -0.3207765  0.0879403  -3.648 0.000277 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3583 on 1091 degrees of freedom
## Multiple R-squared:  0.5812, Adjusted R-squared:  0.5782 
## F-statistic: 189.3 on 8 and 1091 DF,  p-value: < 2.2e-16

On utilise une procédure de sélection de modèle pour obtenir un modèle de dimension réduite.

pga.lm.scoring0 <- lm(SCORING ~ 1, pga)
pga.lm.scoring.step <- step(pga.lm.scoring0, scope = formula(pga.lm.scoring))
## Start:  AIC=-1307.41
## SCORING ~ 1
## 
##              Df Sum of Sq    RSS     AIC
## + BA          1   112.653 221.86 -1757.1
## + SCRAMBLING  1    70.937 263.58 -1567.6
## + PA          1    61.477 273.04 -1528.8
## + GIR         1    56.746 277.77 -1509.9
## + SS          1    31.598 302.92 -1414.6
## + DD          1    19.369 315.15 -1371.0
## + ROUNDS      1     4.210 330.31 -1319.3
## + DA          1     0.937 333.58 -1308.5
## <none>                    334.52 -1307.4
## 
## Step:  AIC=-1757.1
## SCORING ~ BA
## 
##              Df Sum of Sq    RSS     AIC
## + SCRAMBLING  1    48.652 173.21 -2027.4
## + SS          1    17.315 204.55 -1844.5
## + GIR         1    14.500 207.36 -1829.5
## + DA          1     7.464 214.40 -1792.7
## + ROUNDS      1     2.844 219.02 -1769.3
## + PA          1     0.834 221.03 -1759.2
## + DD          1     0.722 221.14 -1758.7
## <none>                    221.86 -1757.1
## - BA          1   112.653 334.52 -1307.4
## 
## Step:  AIC=-2027.4
## SCORING ~ BA + SCRAMBLING
## 
##              Df Sum of Sq    RSS     AIC
## + GIR         1    19.467 153.75 -2156.5
## + ROUNDS      1     4.988 168.22 -2057.5
## + DD          1     2.876 170.34 -2043.8
## + PA          1     2.550 170.66 -2041.7
## + DA          1     1.396 171.82 -2034.3
## <none>                    173.21 -2027.4
## + SS          1     0.065 173.15 -2025.8
## - SCRAMBLING  1    48.652 221.87 -1757.1
## - BA          1    90.368 263.58 -1567.6
## 
## Step:  AIC=-2156.54
## SCORING ~ BA + SCRAMBLING + GIR
## 
##              Df Sum of Sq    RSS     AIC
## + ROUNDS      1     7.642 146.10 -2210.6
## + PA          1     1.821 151.93 -2167.7
## + DD          1     1.374 152.37 -2164.4
## + SS          1     1.273 152.47 -2163.7
## + DA          1     0.318 153.43 -2156.8
## <none>                    153.75 -2156.5
## - GIR         1    19.467 173.21 -2027.4
## - BA          1    50.063 203.81 -1848.5
## - SCRAMBLING  1    53.618 207.36 -1829.5
## 
## Step:  AIC=-2210.62
## SCORING ~ BA + SCRAMBLING + GIR + ROUNDS
## 
##              Df Sum of Sq    RSS     AIC
## + PA          1     2.452 143.65 -2227.2
## + SS          1     0.875 145.23 -2215.2
## + DD          1     0.758 145.35 -2214.3
## + DA          1     0.267 145.84 -2210.6
## <none>                    146.10 -2210.6
## - ROUNDS      1     7.642 153.75 -2156.5
## - GIR         1    22.121 168.22 -2057.5
## - BA          1    46.358 192.46 -1909.5
## - SCRAMBLING  1    56.793 202.90 -1851.4
## 
## Step:  AIC=-2227.23
## SCORING ~ BA + SCRAMBLING + GIR + ROUNDS + PA
## 
##              Df Sum of Sq    RSS     AIC
## + DD          1     2.548 141.10 -2244.9
## + DA          1     0.847 142.81 -2231.7
## + SS          1     0.735 142.92 -2230.9
## <none>                    143.65 -2227.2
## - PA          1     2.452 146.10 -2210.6
## - BA          1     7.723 151.38 -2171.6
## - ROUNDS      1     8.272 151.93 -2167.7
## - GIR         1    21.952 165.60 -2072.8
## - SCRAMBLING  1    41.626 185.28 -1949.3
## 
## Step:  AIC=-2244.92
## SCORING ~ BA + SCRAMBLING + GIR + ROUNDS + PA + DD
## 
##              Df Sum of Sq    RSS     AIC
## + SS          1     0.940 140.16 -2250.3
## <none>                    141.10 -2244.9
## + DA          1     0.037 141.07 -2243.2
## - BA          1     1.949 143.05 -2231.8
## - DD          1     2.548 143.65 -2227.2
## - PA          1     4.242 145.35 -2214.3
## - ROUNDS      1     7.418 148.52 -2190.6
## - GIR         1    23.617 164.72 -2076.7
## - SCRAMBLING  1    44.110 185.21 -1947.7
## 
## Step:  AIC=-2250.28
## SCORING ~ BA + SCRAMBLING + GIR + ROUNDS + PA + DD + SS
## 
##              Df Sum of Sq    RSS     AIC
## <none>                    140.16 -2250.3
## + DA          1    0.0803 140.08 -2248.9
## - SS          1    0.9399 141.10 -2244.9
## - BA          1    1.7466 141.91 -2238.7
## - DD          1    2.7537 142.92 -2230.9
## - PA          1    4.1601 144.32 -2220.1
## - ROUNDS      1    6.9529 147.12 -2199.0
## - GIR         1   24.4357 164.60 -2075.5
## - SCRAMBLING  1   27.1027 167.27 -2057.8
summary(pga.lm.scoring.step)
## 
## Call:
## lm(formula = SCORING ~ BA + SCRAMBLING + GIR + ROUNDS + PA + 
##     DD + SS, data = pga)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07189 -0.24036 -0.00072  0.23178  1.08779 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 74.226703   1.569649  47.289  < 2e-16 ***
## BA          -0.323994   0.087831  -3.689 0.000236 ***
## SCRAMBLING  -0.064059   0.004408 -14.531  < 2e-16 ***
## GIR         -0.086441   0.006265 -13.798  < 2e-16 ***
## ROUNDS       0.006028   0.000819   7.360 3.61e-13 ***
## PA           5.096219   0.895161   5.693 1.60e-08 ***
## DD          -0.007354   0.001588  -4.632 4.06e-06 ***
## SS          -0.006150   0.002273  -2.706 0.006914 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3583 on 1092 degrees of freedom
## Multiple R-squared:  0.581,  Adjusted R-squared:  0.5783 
## F-statistic: 216.3 on 7 and 1092 DF,  p-value: < 2.2e-16

On en déduit que les variables ROUNDS, DD, GIR, SS, SCRAMBLING, PA et BA sont significatives pour expliquer la variable SCORING. En revanche, Driving Accuracy ne permet pas d’expliquer le score moyen obtenu.

Donc le score ne se déduit pas de la précision des drives. On voit notamment que les Sand Saves et les Scrambling, qui sont des situations difficiles pour les golfeurs, permettent d’expliquer le score moyen. C’est pour cela que les golfeurs qui obtiennent les meilleurs scores moyens sont ceux qui sont complets et qui savent gérer les situations délicates.

Notons tout de même que la variable SS a moins d’influence que les autres variables pour expliquer le score. En effet, les Sand Saves sont des situations assez rares en comparaison avec les Green in Regulation par exemple.

5 Analyse en composantes principales

Dans cette partie, on réalise une analyse en composantes principales (ACP) sur le jeu de données avec le package FactoMineR (Husson 2008; STHDA, s. d.a). Cela permet d’extraire l’essentiel de l’information contenue dans le jeu de données afin de pouvoir l’interpréter plus aisément. Nous allons donc approfondir les réponses aux questions 4 et 5 précédentes.

On considère :

  • WINS comme une variable quantitative supplémentaire
  • YEAR comme une variable qualitative supplémentaire

On affiche le résumé de l’ACP.

pga.pca <- PCA(pga[3:15], quanti.sup = 2, quali.sup = 13, graph = FALSE)
summary(pga.pca)
## 
## Call:
## PCA(X = pga[3:15], quanti.sup = 2, quali.sup = 13, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
## Variance               3.866   2.042   1.487   1.013   0.837   0.772   0.418
## % of var.             35.144  18.567  13.519   9.205   7.609   7.016   3.801
## Cumulative % of var.  35.144  53.711  67.230  76.435  84.043  91.059  94.860
##                        Dim.8   Dim.9  Dim.10  Dim.11
## Variance               0.232   0.179   0.099   0.055
## % of var.              2.107   1.631   0.902   0.500
## Cumulative % of var.  96.967  98.599  99.500 100.000
## 
## Individuals (the 10 first)
##                Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1          |  3.274 |  2.373  0.132  0.526 |  0.022  0.000  0.000 |  0.662
## 2          |  3.937 |  2.936  0.203  0.556 |  1.973  0.173  0.251 |  1.092
## 3          |  4.176 |  1.837  0.079  0.194 |  2.225  0.220  0.284 |  1.874
## 4          |  3.318 |  2.373  0.132  0.512 | -0.602  0.016  0.033 | -1.441
## 5          |  3.182 |  2.379  0.133  0.559 |  1.794  0.143  0.318 |  0.133
## 6          |  2.451 |  1.283  0.039  0.274 | -1.126  0.056  0.211 | -0.470
## 7          |  2.856 |  0.270  0.002  0.009 | -0.094  0.000  0.001 |  2.291
## 8          |  3.301 |  1.672  0.066  0.257 | -2.726  0.331  0.682 |  0.006
## 9          |  4.299 |  3.421  0.275  0.633 |  1.974  0.173  0.211 |  1.468
## 10         |  4.367 |  1.056  0.026  0.058 |  3.783  0.637  0.750 |  1.336
##               ctr   cos2  
## 1           0.027  0.041 |
## 2           0.073  0.077 |
## 3           0.215  0.201 |
## 4           0.127  0.189 |
## 5           0.001  0.002 |
## 6           0.014  0.037 |
## 7           0.321  0.643 |
## 8           0.000  0.000 |
## 9           0.132  0.117 |
## 10          0.109  0.094 |
## 
## Variables (the 10 first)
##               Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr   cos2
## POINTS     |  0.828 17.715  0.685 | -0.015  0.011  0.000 |  0.121  0.992  0.015
## ROUNDS     | -0.024  0.015  0.001 |  0.139  0.947  0.019 |  0.288  5.567  0.083
## SCORING    | -0.861 19.196  0.742 | -0.078  0.300  0.006 | -0.165  1.838  0.027
## DD         |  0.371  3.554  0.137 | -0.809 32.045  0.654 | -0.141  1.334  0.020
## DA         | -0.093  0.222  0.009 |  0.594 17.250  0.352 |  0.662 29.460  0.438
## GIR        |  0.349  3.149  0.122 | -0.277  3.760  0.077 |  0.766 39.453  0.587
## SS         |  0.384  3.808  0.147 |  0.592 17.182  0.351 | -0.341  7.833  0.116
## SCRAMBLING |  0.456  5.373  0.208 |  0.682 22.795  0.466 | -0.094  0.600  0.009
## PA         | -0.633 10.379  0.401 | -0.205  2.060  0.042 |  0.422 11.961  0.178
## BA         |  0.799 16.494  0.638 | -0.266  3.452  0.071 | -0.054  0.194  0.003
##             
## POINTS     |
## ROUNDS     |
## SCORING    |
## DD         |
## DA         |
## GIR        |
## SS         |
## SCRAMBLING |
## PA         |
## BA         |
## 
## Supplementary continuous variable
##               Dim.1   cos2    Dim.2   cos2    Dim.3   cos2  
## WINS       |  0.526  0.276 | -0.096  0.009 |  0.009  0.000 |
## 
## Supplementary categories (the 10 first)
##                Dist    Dim.1   cos2 v.test    Dim.2   cos2 v.test    Dim.3
## 2010       |  0.978 | -0.532  0.296 -2.839 |  0.584  0.356  4.282 |  0.505
## 2011       |  0.553 | -0.321  0.336 -1.709 | -0.264  0.229 -1.940 |  0.052
## 2012       |  0.761 | -0.472  0.385 -2.519 | -0.014  0.000 -0.102 | -0.074
## 2013       |  0.976 | -0.618  0.401 -3.294 |  0.276  0.080  2.024 | -0.065
## 2014       |  0.615 | -0.102  0.027 -0.542 |  0.248  0.163  1.820 | -0.181
## 2015       |  0.636 |  0.247  0.151  1.316 |  0.143  0.051  1.050 |  0.193
## 2016       |  0.370 | -0.086  0.054 -0.459 | -0.038  0.011 -0.281 | -0.002
## 2017       |  0.528 |  0.084  0.025  0.446 | -0.292  0.306 -2.140 | -0.141
## 2018       |  0.952 |  0.591  0.385  3.152 | -0.408  0.184 -2.996 |  0.380
## 2019       |  1.076 |  0.651  0.366  3.471 |  0.105  0.009  0.767 | -0.094
##              cos2 v.test  
## 2010        0.266  4.338 |
## 2011        0.009  0.449 |
## 2012        0.010 -0.638 |
## 2013        0.004 -0.560 |
## 2014        0.087 -1.560 |
## 2015        0.093  1.663 |
## 2016        0.000 -0.019 |
## 2017        0.071 -1.212 |
## 2018        0.159  3.267 |
## 2019        0.008 -0.810 |

On regarde les valeurs propres et parts d’inertie expliquées par chaque axe.

# Extract eigenvalues/variances
get_eig(pga.pca)
##        eigenvalue variance.percent cumulative.variance.percent
## Dim.1  3.86581862       35.1438056                    35.14381
## Dim.2  2.04239475       18.5672250                    53.71103
## Dim.3  1.48707155       13.5188323                    67.22986
## Dim.4  1.01251114        9.2046467                    76.43451
## Dim.5  0.83693588        7.6085080                    84.04302
## Dim.6  0.77181114        7.0164649                    91.05948
## Dim.7  0.41810942        3.8009947                    94.86048
## Dim.8  0.23174137        2.1067397                    96.96722
## Dim.9  0.17944710        1.6313373                    98.59855
## Dim.10 0.09917233        0.9015666                    99.50012
## Dim.11 0.05498672        0.4998792                   100.00000
# Visualize eigenvalues/variances
fviz_screeplot(pga.pca, addlabels = TRUE, ylim = c(0, 50))
Part d'inertie expliquée pour chaque axe

Figure 5.1: Part d’inertie expliquée pour chaque axe

D’après la figure 5.1, on peut garder 4 composantes, ce qui permet d’expliquer environ 76% de l’inertie totale.

La figure 5.2 montre le graphe des variables pour chaque axe, en indiquant la contribution de chaque variable.

for (i in 2:4) {
    # Graph of variables: control variable colors using their contributions
    print(fviz_pca_var(pga.pca, axes = c(1, i), col.var = "contrib",
                       gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                       repel = TRUE))
}
Graphe des variables pour chaque axeGraphe des variables pour chaque axeGraphe des variables pour chaque axe

Figure 5.2: Graphe des variables pour chaque axe

La figure 5.3 montre les contributions des variables pour chaque axe.

for (i in 1:4) {
  # Contributions of variables
  print(fviz_contrib(pga.pca, choice = "var", axes = i, top = 10))
}
Graphe des contributions pour chaque axeGraphe des contributions pour chaque axeGraphe des contributions pour chaque axeGraphe des contributions pour chaque axe

Figure 5.3: Graphe des contributions pour chaque axe

La figure 5.4 montre le graphe des individus pour chaque axe, en habillant avec la variable WINS.

for (i in 2:4) {
    # Graph of individuals
    print(fviz_pca_ind(pga.pca, axes = c(1, i), col.ind = pga$WINS,
                       gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                       label = "none"))
}
Graphe des individus pour chaque axeGraphe des individus pour chaque axeGraphe des individus pour chaque axe

Figure 5.4: Graphe des individus pour chaque axe

Interprétations :

  • L’axe 1 oppose les très bons golfeurs aux golfeurs moyens. En effet, les variables en rapport avec les scores (SCORING, POINTS, BA, PA) contribuent beaucoup à cet axe. On peut voir aussi qu’il oppose les golfeurs gagnant beaucoup de gains (MONEY) et de tournois (WINS) à ceux en gagnant moins.
  • L’axe 2 oppose les joueurs puissants (DD) aux joueurs précis (DA) et qui savent gérer les situations difficiles (SCRAMBLING, SS).
  • L’axe 3 oppose les joueurs bons en GIR et précis aux autres joueurs.
  • L’axe 4 oppose les joueurs disputant beaucoup de rounds (ROUNDS) aux joueurs en disputant moins.

6 Classification non supervisée

Dans cette partie, le but est d’identifier des groupes de golfeurs homogènes, partageant des caractéristiques similaires. Pour cela, on applique différentes méthodes de classification non supervisée sur le jeu de données.

On se restreint d’abord aux statistiques des golfeurs de l’année 2020. On étudie alors 100 golfeurs.

On crée le jeu de données avec uniquement les variables quantitatives.

pga.2020 <- pga[which(pga$YEAR == 2020),]

# Set row names
row.names(pga.2020) <- pga.2020$PLAYER.NAME

# Remove columns PLAYER.NAME, COUNTRY, WINS and YEAR
pga.2020 <- pga.2020[-c(1, 2, 4, 15)]

head(pga.2020)
##                   POINTS ROUNDS SCORING    DD    DA   GIR    SS SCRAMBLING
## Justin Thomas       2458     66  69.128 304.2 57.25 69.61 62.50      63.16
## Collin Morikawa     1902     80  69.856 297.3 65.54 69.03 53.97      59.42
## Webb Simpson        1878     52  68.978 296.2 67.31 70.83 56.34      64.47
## Bryson DeChambeau   1657     62  69.241 322.1 57.37 68.46 62.89      63.92
## Sungjae Im          1633     94  69.987 300.5 64.48 67.08 58.82      62.12
## Patrick Reed        1426     74  69.586 296.0 56.08 65.32 60.61      63.64
##                      PA   BA   MONEY
## Justin Thomas     1.722 4.55 7344040
## Collin Morikawa   1.749 4.05 5250868
## Webb Simpson      1.691 4.67 5097742
## Bryson DeChambeau 1.726 4.42 4998495
## Sungjae Im        1.734 4.15 4337811
## Patrick Reed      1.700 4.23 4250060
str(pga.2020)
## 'data.frame':    100 obs. of  11 variables:
##  $ POINTS    : num  2458 1902 1878 1657 1633 ...
##  $ ROUNDS    : int  66 80 52 62 94 74 64 60 86 58 ...
##  $ SCORING   : num  69.1 69.9 69 69.2 70 ...
##  $ DD        : num  304 297 296 322 300 ...
##  $ DA        : num  57.2 65.5 67.3 57.4 64.5 ...
##  $ GIR       : num  69.6 69 70.8 68.5 67.1 ...
##  $ SS        : num  62.5 54 56.3 62.9 58.8 ...
##  $ SCRAMBLING: num  63.2 59.4 64.5 63.9 62.1 ...
##  $ PA        : num  1.72 1.75 1.69 1.73 1.73 ...
##  $ BA        : num  4.55 4.05 4.67 4.42 4.15 4.23 4.38 4.2 3.91 4.22 ...
##  $ MONEY     : num  7344040 5250868 5097742 4998495 4337811 ...

Voici un résumé du jeu de données étudié.

summary(pga.2020)
##      POINTS           ROUNDS      SCORING            DD              DA       
##  Min.   : 352.0   Min.   :38   Min.   :68.98   Min.   :281.5   Min.   :51.79  
##  1st Qu.: 436.2   1st Qu.:60   1st Qu.:70.15   1st Qu.:292.1   1st Qu.:57.57  
##  Median : 599.5   Median :70   Median :70.50   Median :298.4   Median :60.98  
##  Mean   : 717.5   Mean   :68   Mean   :70.44   Mean   :298.3   Mean   :61.14  
##  3rd Qu.: 881.8   3rd Qu.:76   3rd Qu.:70.82   3rd Qu.:302.8   3rd Qu.:64.49  
##  Max.   :2458.0   Max.   :96   Max.   :72.02   Max.   :322.1   Max.   :73.86  
##       GIR              SS          SCRAMBLING          PA       
##  Min.   :61.48   Min.   :39.02   Min.   :51.46   Min.   :1.691  
##  1st Qu.:65.52   1st Qu.:47.24   1st Qu.:58.00   1st Qu.:1.736  
##  Median :66.92   Median :50.24   Median :60.00   Median :1.750  
##  Mean   :67.02   Mean   :51.57   Mean   :59.99   Mean   :1.751  
##  3rd Qu.:68.55   3rd Qu.:56.37   3rd Qu.:62.24   3rd Qu.:1.765  
##  Max.   :71.69   Max.   :66.99   Max.   :67.45   Max.   :1.824  
##        BA            MONEY        
##  Min.   :3.200   Min.   : 783232  
##  1st Qu.:3.708   1st Qu.:1146694  
##  Median :3.835   Median :1604653  
##  Mean   :3.867   Mean   :2086422  
##  3rd Qu.:4.015   3rd Qu.:2418142  
##  Max.   :4.670   Max.   :7344040

6.1 Classification ascendante hiérarchique

On réalise une classification ascendante hiérarchique (CAH) des golfeurs de l’année 2020.

Les variables étant dans des unités différentes, on travaille sur les données centrées réduites pour accorder la même importance à chaque variable et éviter que les variables à forte variance pèsent indûment sur les résultats.

pga.scale <- scale(pga.2020)
summary(pga.scale)
##      POINTS            ROUNDS           SCORING              DD          
##  Min.   :-0.9352   Min.   :-2.4029   Min.   :-2.5023   Min.   :-2.01552  
##  1st Qu.:-0.7196   1st Qu.:-0.6408   1st Qu.:-0.4959   1st Qu.:-0.74417  
##  Median :-0.3018   Median : 0.1602   Median : 0.1004   Median : 0.02224  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.4205   3rd Qu.: 0.6408   3rd Qu.: 0.6550   3rd Qu.: 0.54521  
##  Max.   : 4.4543   Max.   : 2.2427   Max.   : 2.7167   Max.   : 2.86548  
##        DA                GIR                 SS            SCRAMBLING       
##  Min.   :-1.95202   Min.   :-2.45139   Min.   :-2.0914   Min.   :-2.663182  
##  1st Qu.:-0.74429   1st Qu.:-0.66488   1st Qu.:-0.7216   1st Qu.:-0.622468  
##  Median :-0.03156   Median :-0.04503   Median :-0.2211   Median : 0.002622  
##  Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.000000  
##  3rd Qu.: 0.70049   3rd Qu.: 0.67666   3rd Qu.: 0.8000   3rd Qu.: 0.702630  
##  Max.   : 2.65750   Max.   : 2.06912   Max.   : 2.5711   Max.   : 2.328177  
##        PA                 BA              MONEY        
##  Min.   :-2.52739   Min.   :-2.5424   Min.   :-1.0059  
##  1st Qu.:-0.62264   1st Qu.:-0.6077   1st Qu.:-0.7254  
##  Median :-0.03005   Median :-0.1216   Median :-0.3719  
##  Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.61545   3rd Qu.: 0.5646   3rd Qu.: 0.2561  
##  Max.   : 3.10220   Max.   : 3.0616   Max.   : 4.0584

On calcule la matrice des distances entre individus avec la distance euclidienne, puis on réalise la CAH avec la mesure de Ward. La figure 6.1 montre le dendrogramme obtenu.

# Compute distances and hierarchical clustering
pga.dist <- dist(pga.scale, method = "euclidean")
pga.hc <- hclust(pga.dist, method = "ward.D2")

# Plot dendrogram
fviz_dend(pga.hc, cex = 0.3, labels_track_height = 5,
          main = "Cluster Dendrogram PGA Tour 2020")
Dendrogramme de la CAH

Figure 6.1: Dendrogramme de la CAH

La figure 6.2 montre les hauteurs des branches pour déterminer le nombre de classes à garder.

# Dataframe of heights
df.height <- data.frame(Class = 1:length(pga.hc$height),
                        Height = sort(pga.hc$height, decreasing = TRUE))

# Barplot
ggplot(df.height, aes(x = Class, y = Height, fill = Height)) +
  geom_bar(stat = "identity", width = 0.5) +
  ggtitle("Dendrogram Height") +
  theme_minimal()
Hauteurs des branches du dendrogram de la CAH

Figure 6.2: Hauteurs des branches du dendrogram de la CAH

La perte d’inertie devient plus importante quand on passe de 4 à 3 classes, donc on peut garder 4 classes.

On affiche alors les classes sur le dendrogramme comme le montre la figure 6.3.

k <- 4
# Compute hierarchical clustering and cut into 4 clusters
pga.hcut <- hcut(pga.2020, k = k, stand = TRUE)

# Visualize
fviz_dend(pga.hcut, rect = TRUE, cex = 0.3, labels_track_height = 5,
          main = "Cluster Dendrogram PGA Tour 2020")
Dendrogramme de la CAH avec clusters

Figure 6.3: Dendrogramme de la CAH avec clusters

On forme ensuite les groupes : on liste les golfeurs de chaque groupe.

pga.groups <- cutree(pga.hc, k)

for (i in 1:k)
{
  cat("Group", i,"\n")
  print(rownames(pga.2020)[which(pga.groups == i)])
}
## Group 1 
##  [1] "Justin Thomas"     "Webb Simpson"      "Bryson DeChambeau"
##  [4] "Patrick Reed"      "Daniel Berger"     "Rory McIlroy"     
##  [7] "Jon Rahm"          "Xander Schauffele" "Dustin Johnson"   
## [10] "Hideki Matsuyama"  "Tyrrell Hatton"    "Tony Finau"       
## [13] "Scottie Scheffler" "Harris English"    "Patrick Cantlay"  
## [16] "Adam Scott"       
## Group 2 
##  [1] "Collin Morikawa"  "Sungjae Im"       "Brendon Todd"     "Lanto Griffin"   
##  [5] "Abraham Ancer"    "Sebastián Muñoz"  "Kevin Na"         "Adam Long"       
##  [9] "Joel Dahmen"      "Tom Hoge"         "Kevin Kisner"     "Mark Hubbard"    
## [13] "Jason Day"        "Dylan Frittelli"  "Maverick McNealy" "J.T. Poston"     
## [17] "Matt Kuchar"      "Brian Harman"     "Bud Cauley"       "Alex Noren"      
## [21] "Ian Poulter"      "Brandt Snedeker" 
## Group 3 
##  [1] "Marc Leishman"      "Cameron Champ"      "Kevin Streelman"   
##  [4] "Billy Horschel"     "Joaquin Niemann"    "Viktor Hovland"    
##  [7] "Ryan Palmer"        "Cameron Smith"      "Byeong Hun An"     
## [10] "Gary Woodland"      "Matthew Wolff"      "Richy Werenski"    
## [13] "Brendan Steele"     "Carlos Ortiz"       "Mackenzie Hughes"  
## [16] "Danny Lee"          "Paul Casey"         "Max Homa"          
## [19] "Sung Kang"          "Talor Gooch"        "Charles Howell III"
## [22] "Bubba Watson"       "Denny McCarthy"     "Phil Mickelson"    
## [25] "Xinjun Zhang"       "Sepp Straka"        "Harry Higgs"       
## [28] "Patrick Rodgers"    "Robby Shelton"      "Si Woo Kim"        
## [31] "Matt Jones"         "Cameron Tringale"   "Tommy Fleetwood"   
## [34] "Jason Kokrak"       "Cameron Davis"      "Scott Harrington"  
## [37] "Brooks Koepka"      "Jordan Spieth"     
## Group 4 
##  [1] "Tyler Duncan"      "Nick Taylor"       "Adam Hadwin"      
##  [4] "Michael Thompson"  "Andrew Landry"     "Matt Fitzpatrick" 
##  [7] "Jim Herman"        "Corey Conners"     "Doc Redman"       
## [10] "Henrik Norlander"  "Harold Varner III" "Vaughn Taylor"    
## [13] "Brian Stuard"      "Pat Perez"         "Troy Merritt"     
## [16] "Chez Reavie"       "Nate Lashley"      "Rickie Fowler"    
## [19] "Emiliano Grillo"   "Matthew NeSmith"   "Ryan Armour"      
## [22] "Ryan Moore"        "Louis Oosthuizen"  "Russell Henley"

On affiche le nombre de golfeurs par groupe.

table(pga.groups)
## pga.groups
##  1  2  3  4 
## 16 22 38 24

On regarde maintenant les caractéristiques moyennes de chaque groupe.

# Create matrix
pga.groups.means <- matrix(NA, nrow = k, ncol = ncol(pga.2020))

# Set rownames and colnames
rownames(pga.groups.means) <- 1:k
colnames(pga.groups.means) <- colnames(pga.2020)

# Compute means
for (i in 1:k)
{
  pga.groups.means[i,] <- colMeans(pga.2020[pga.groups == i,])
}

round(pga.groups.means, 2)
##    POINTS ROUNDS SCORING     DD    DA   GIR    SS SCRAMBLING   PA   BA   MONEY
## 1 1244.81  61.56   69.45 305.84 59.31 68.91 54.59      62.73 1.73 4.27 4169090
## 2  811.27  74.45   70.32 292.83 62.19 65.62 55.02      62.61 1.74 3.83 2194328
## 3  588.00  68.82   70.69 302.60 57.99 66.59 50.05      58.23 1.76 3.83 1672304
## 4  484.83  65.08   70.80 291.32 66.37 67.72 48.78      58.56 1.76 3.69 1254750

Interprétations : en 2020, le top 100 des golfeurs au classement de la FedExCup peut être divisé en 4 groupes. On peut même identifier 2 paires de 2 groupes. La première paire est constituée des golfeurs complets, prétendant aux victoires :

  • Le groupe 1 est composé de l’élite des golfeurs. C’est le plus petit groupe. On retrouve notamment Dustin Johnson, le vainqueur de la FedExCup 2020. Ce groupe a les meilleures caractéristiques moyennes pour presque toutes les statistiques. On voit notamment un grand écart au niveau du nombre de points et des gains remportés par rapport aux autres groupes. En revanche, les golfeurs de ce groupe participent à moins de rounds que les autres. On peut alors nommer ce groupe celui des favoris.
  • Le groupe 2 constitue le groupe des outsiders. Les caractéristiques moyennes sont un peu moins bonnes que celles du groupe 1, mais restent néanmoins assez proches pour espérer des podiums et des victoires lors des tournois.

La seconde paire, constituée des groupes 3 et 4, est plus dense et est composée des autres golfeurs du top 100. Les scores moyens sont très proches (70.69 et 70.80) et les précisions des putts sont même égales (1.76). Ces golfeurs ont plus de mal à gérer les situations délicates que les deux premiers groupes (SS et SCRAMBLING). On peut identifier deux styles de jeu différents :

  • Le groupe 3 est composé des golfeurs effectuant des drives plus longs, mais moins précis. En effet :
    • La distance moyenne des drives de ce groupe est de 10 yards supérieure à celle du groupe 4.
    • C’est le groupe ayant la moins bonne précision moyenne aux drives (57.99%).
  • Le groupe 4 est composé des golfeurs effectuant des drives plus courts, mais plus précis. En effet :
    • La distance moyenne des drives est plus de 10 yards inférieure à celle du groupe 3.
    • C’est le groupe ayant nettement la meilleure précision de drives, avec 7% de plus que le groupe 1 constitué des meilleurs golfeurs.

Pour conclure :

  • Les groupes 1 et 2 sont constitués des golfeurs complets qui peuvent prétendre aux victoires.
  • Les groupes 3 et 4 sont les autres golfeurs qui ont des scores moyens similaires mais qui s’appuient sur leurs points forts. Ceux du groupe 3 préfèrent allonger leurs drives pour progresser plus vite au risque d’être moins précis, tandis que ceux du groupe 4 les raccourcissent pour être plus précis.

6.2 Méthode des k-means

On met ensuite en oeuvre la méthode des k-means avec les mêmes données. On utilise comme initialisation les centres de gravité des classes obtenues avec la CAH.

# Create init matrix
init <- matrix(NA, nrow = k, ncol = ncol(pga.2020))
colnames(init) <- colnames(pga.2020)
for (i in 1:k)
{
  init[i,] <- colMeans(pga.scale[pga.groups == i,])
}

# Compute k-means
pga.kmeans <- kmeans(pga.scale, centers = init)
pga.kmeans
## K-means clustering with 4 clusters of sizes 17, 21, 37, 25
## 
## Cluster means:
##        POINTS      ROUNDS    SCORING         DD          DA         GIR
## 1  1.58504832 -0.16490188 -1.6432936  0.7969638 -0.19354702  0.81248118
## 2 -0.03941647  0.27842617 -0.0878207 -0.6238075  0.06058796 -0.78917003
## 3 -0.37312248 -0.04978891  0.3480919  0.5984928 -0.66612251 -0.08864522
## 4 -0.49250175 -0.04805712  0.6760331 -0.9037065  1.06657941  0.24161054
##           SS SCRAMBLING         PA         BA      MONEY
## 1  0.5351646  0.8298330 -1.0409399  1.5770629  1.7505405
## 2  0.8185057  0.7880637 -0.4513150 -0.2586669 -0.1106063
## 3 -0.3379735 -0.4985977  0.3394574 -0.1334580 -0.3590634
## 4 -0.5512559 -0.4883354  0.5845467 -0.6576048 -0.5660444
## 
## Clustering vector:
##      Justin Thomas    Collin Morikawa       Webb Simpson  Bryson DeChambeau 
##                  1                  1                  1                  1 
##         Sungjae Im       Patrick Reed      Daniel Berger       Rory McIlroy 
##                  1                  1                  1                  1 
##       Brendon Todd           Jon Rahm  Xander Schauffele      Lanto Griffin 
##                  2                  1                  1                  2 
##      Abraham Ancer      Marc Leishman     Dustin Johnson    Sebastián Muñoz 
##                  2                  3                  1                  2 
##           Kevin Na   Hideki Matsuyama     Tyrrell Hatton      Cameron Champ 
##                  2                  1                  1                  3 
##          Adam Long    Kevin Streelman         Tony Finau  Scottie Scheffler 
##                  4                  4                  1                  1 
##     Billy Horschel    Joaquin Niemann     Harris English     Viktor Hovland 
##                  3                  3                  1                  3 
##        Ryan Palmer      Cameron Smith      Byeong Hun An    Patrick Cantlay 
##                  3                  2                  3                  1 
##      Gary Woodland      Matthew Wolff       Tyler Duncan         Adam Scott 
##                  3                  3                  4                  3 
##        Nick Taylor        Joel Dahmen           Tom Hoge       Kevin Kisner 
##                  4                  2                  4                  2 
##     Richy Werenski       Mark Hubbard     Brendan Steele        Adam Hadwin 
##                  3                  2                  3                  2 
##          Jason Day   Michael Thompson       Carlos Ortiz      Andrew Landry 
##                  2                  4                  3                  4 
##    Dylan Frittelli   Matt Fitzpatrick   Mackenzie Hughes          Danny Lee 
##                  2                  4                  2                  3 
##         Jim Herman         Paul Casey      Corey Conners           Max Homa 
##                  4                  3                  4                  3 
##   Maverick McNealy        J.T. Poston         Doc Redman          Sung Kang 
##                  2                  2                  4                  3 
##        Talor Gooch        Matt Kuchar Charles Howell III       Bubba Watson 
##                  3                  2                  3                  3 
##     Denny McCarthy     Phil Mickelson   Henrik Norlander       Brian Harman 
##                  3                  3                  4                  2 
##       Xinjun Zhang        Sepp Straka        Harry Higgs  Harold Varner III 
##                  3                  3                  3                  3 
##         Bud Cauley      Vaughn Taylor       Brian Stuard    Patrick Rodgers 
##                  2                  4                  4                  3 
##         Alex Noren          Pat Perez       Troy Merritt      Robby Shelton 
##                  2                  4                  4                  4 
##         Si Woo Kim        Chez Reavie       Nate Lashley        Ian Poulter 
##                  3                  4                  4                  2 
##         Matt Jones   Cameron Tringale      Rickie Fowler    Tommy Fleetwood 
##                  3                  3                  4                  3 
##       Jason Kokrak      Cameron Davis    Emiliano Grillo    Matthew NeSmith 
##                  3                  3                  4                  4 
##   Scott Harrington        Ryan Armour         Ryan Moore      Brooks Koepka 
##                  3                  4                  4                  3 
##    Brandt Snedeker      Jordan Spieth   Louis Oosthuizen     Russell Henley 
##                  2                  3                  3                  4 
## 
## Within cluster sum of squares by cluster:
## [1] 135.3548 114.4067 208.8775 161.4616
##  (between_SS / total_SS =  43.1 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

On compare ce clustering avec celui obtenu avec la CAH.

table(pga.groups, pga.kmeans$cluster)
##           
## pga.groups  1  2  3  4
##          1 15  0  1  0
##          2  2 18  0  2
##          3  0  2 34  2
##          4  0  1  2 21
adjustedRandIndex(pga.groups, pga.kmeans$cluster)
## [1] 0.7036602

Seulement 10 individus sur les 100 sont classés différemment. Donc, les partitionnements sont très proches. De plus, l’index d’ajustement (0.70) est proche de 1 ce qui confirme ce constat.

On lance maintenant une ACP pour interpréter les groupes. On crée un jeu de données des golfeurs de 2020, en ajoutant le numéro de groupe des golfeurs.

pga.2020.pca <- PCA(cbind.data.frame(
  pga.2020, Group = factor(pga.kmeans$cluster)),
  scale.unit = TRUE, quali.sup = 12, graph = FALSE)
summary(pga.2020.pca)
## 
## Call:
## PCA(X = cbind.data.frame(pga.2020, Group = factor(pga.kmeans$cluster)),  
##      scale.unit = TRUE, quali.sup = 12, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
## Variance               4.354   1.761   1.391   1.067   0.891   0.646   0.418
## % of var.             39.583  16.009  12.649   9.702   8.097   5.870   3.801
## Cumulative % of var.  39.583  55.592  68.241  77.943  86.039  91.910  95.711
##                        Dim.8   Dim.9  Dim.10  Dim.11
## Variance               0.192   0.164   0.074   0.042
## % of var.              1.743   1.490   0.677   0.380
## Cumulative % of var.  97.454  98.943  99.620 100.000
## 
## Individuals (the 10 first)
##                       Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## Justin Thomas     |  7.549 |  7.026 11.338  0.866 | -0.267  0.040  0.001 |
## Collin Morikawa   |  4.428 |  3.070  2.164  0.481 |  0.454  0.117  0.011 |
## Webb Simpson      |  6.746 |  6.057  8.425  0.806 |  1.066  0.645  0.025 |
## Bryson DeChambeau |  5.956 |  5.516  6.987  0.858 | -1.422  1.149  0.057 |
## Sungjae Im        |  4.213 |  3.134  2.256  0.553 |  1.188  0.802  0.080 |
## Patrick Reed      |  4.527 |  3.926  3.539  0.752 |  1.150  0.750  0.064 |
## Daniel Berger     |  4.785 |  4.417  4.481  0.852 |  1.306  0.969  0.075 |
## Rory McIlroy      |  4.288 |  3.323  2.536  0.601 | -2.173  2.682  0.257 |
## Brendon Todd      |  4.610 |  1.816  0.757  0.155 |  3.958  8.898  0.737 |
## Jon Rahm          |  5.287 |  4.964  5.660  0.881 | -0.034  0.001  0.000 |
##                    Dim.3    ctr   cos2  
## Justin Thomas      0.433  0.135  0.003 |
## Collin Morikawa    1.895  2.580  0.183 |
## Webb Simpson       1.634  1.919  0.059 |
## Bryson DeChambeau -0.660  0.313  0.012 |
## Sungjae Im         0.764  0.420  0.033 |
## Patrick Reed      -1.474  1.562  0.106 |
## Daniel Berger     -0.089  0.006  0.000 |
## Rory McIlroy       0.158  0.018  0.001 |
## Brendon Todd       0.646  0.300  0.020 |
## Jon Rahm           0.671  0.323  0.016 |
## 
## Variables (the 10 first)
##                      Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## POINTS            |  0.834 15.992  0.696 |  0.042  0.102  0.002 |  0.138  1.369
## ROUNDS            | -0.016  0.006  0.000 |  0.245  3.400  0.060 |  0.245  4.298
## SCORING           | -0.919 19.384  0.844 |  0.010  0.006  0.000 | -0.096  0.661
## DD                |  0.400  3.681  0.160 | -0.834 39.521  0.696 | -0.130  1.219
## DA                | -0.166  0.631  0.027 |  0.611 21.233  0.374 |  0.683 33.483
## GIR               |  0.321  2.368  0.103 | -0.320  5.831  0.103 |  0.728 38.142
## SS                |  0.459  4.843  0.211 |  0.403  9.238  0.163 | -0.405 11.779
## SCRAMBLING        |  0.612  8.598  0.374 |  0.493 13.816  0.243 | -0.183  2.414
## PA                | -0.682 10.678  0.465 | -0.302  5.187  0.091 |  0.252  4.552
## BA                |  0.836 16.070  0.700 | -0.166  1.556  0.027 |  0.141  1.435
##                     cos2  
## POINTS             0.019 |
## ROUNDS             0.060 |
## SCORING            0.009 |
## DD                 0.017 |
## DA                 0.466 |
## GIR                0.531 |
## SS                 0.164 |
## SCRAMBLING         0.034 |
## PA                 0.063 |
## BA                 0.020 |
## 
## Supplementary categories
##                       Dist    Dim.1   cos2 v.test    Dim.2   cos2 v.test  
## Group_1           |  3.785 |  3.742  0.977  8.074 | -0.327  0.007 -1.109 |
## Group_2           |  1.644 |  0.184  0.013  0.453 |  1.347  0.671  5.208 |
## Group_3           |  1.309 | -0.688  0.276 -2.513 | -1.022  0.609 -5.873 |
## Group_4           |  2.097 | -1.681  0.643 -4.628 |  0.603  0.083  2.612 |
##                    Dim.3   cos2 v.test  
## Group_1            0.384  0.010  1.464 |
## Group_2           -0.867  0.278 -3.771 |
## Group_3           -0.369  0.080 -2.388 |
## Group_4            1.014  0.234  4.939 |

On regarde les valeurs propres et parts d’inertie expliquées par chaque axe.

# Extract eigenvalues/variances
get_eig(pga.2020.pca)
##        eigenvalue variance.percent cumulative.variance.percent
## Dim.1  4.35408684       39.5826076                    39.58261
## Dim.2  1.76103252       16.0093865                    55.59199
## Dim.3  1.39137836       12.6488942                    68.24089
## Dim.4  1.06718410        9.7016737                    77.94256
## Dim.5  0.89066213        8.0969285                    86.03949
## Dim.6  0.64574977        5.8704524                    91.90994
## Dim.7  0.41808727        3.8007933                    95.71074
## Dim.8  0.19171554        1.7428685                    97.45360
## Dim.9  0.16386890        1.4897173                    98.94332
## Dim.10 0.07445594        0.6768722                    99.62019
## Dim.11 0.04177863        0.3798058                   100.00000
# Visualize eigenvalues/variances
fviz_screeplot(pga.2020.pca, addlabels = TRUE, ylim = c(0, 50))
Part d'inertie expliquée pour chaque axe

Figure 6.4: Part d’inertie expliquée pour chaque axe

On garde les deux premières composantes. On regarde le cercle de corrélation entre elles.

# Graph of variables: control variable colors using their contributions
fviz_pca_var(pga.2020.pca, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping
             )
Graphe des variables

Figure 6.5: Graphe des variables

On étudie les contributions des deux axes.

# Contributions of variables to PC1
fviz_contrib(pga.2020.pca, choice = "var", axes = 1, top = 10)
# Contributions of variables to PC2
fviz_contrib(pga.2020.pca, choice = "var", axes = 2, top = 10)
Graphe des contributionsGraphe des contributions

Figure 6.6: Graphe des contributions

Interprétation des axes :

  • L’axe 1 oppose les meilleurs golfeurs au moins bons. Importance de SCORING, MONEY, BA, POINTS et PA. On rappelle qu’il est mieux d’avoir PA et SCORING faibles.
  • L’axe 2 oppose les golfeurs faisant de longs drives (en bas) et ceux préférant la précision (en haut).

On regarde maintenant le graphe des individus en montrant les groupes identifiés par la méthode des k-means sur les deux premières composantes (figure 6.7).

# Graph of individuals: use habillage to specify groups for coloring
fviz_pca_ind(pga.2020.pca,
             habillage = factor(pga.kmeans$cluster), # color by groups
             addEllipses = TRUE, # Concentration ellipses
             ellipse.type = "convex",
             repel = TRUE,
             labelsize = 2,
             title = "Cluster plot with k-means, PGA Tour 2020"
             )
Clustering avec la méthode des k-means

Figure 6.7: Clustering avec la méthode des k-means

Interprétations :

Le graphe des individus confirme l’interprétation donnée à la suite de la CAH. En effet :

  • On retrouve le groupe des favoris (le 1) le plus à droite, et celui des outsiders (le 2) juste à sa gauche.
  • Les golfeurs du groupe 3 se trouvent en dessous de l’axe des abscisses ce qui montre leur préférence pour les drives longs mais moins précis.
  • Les golfeurs du groupe 4 se trouvent au dessus de l’axe des abscisses ce qui montre leur préférence pour les drives courts mais plus précis.

On représente ensuite la silhouette des individus (figure 6.8).

pga.sil <- silhouette(pga.kmeans$cluster, pga.dist)

fviz_cluster(pga.kmeans, pga.2020, ellipse.type = "convex", repel = TRUE,
             labelsize = 7, main = "Cluster plot with k-means, PGA Tour 2020")
fviz_silhouette(pga.sil)
##   cluster size ave.sil.width
## 1       1   17          0.21
## 2       2   21          0.16
## 3       3   37          0.15
## 4       4   25          0.10
Clusters et silhouettes des individus avec k-meansClusters et silhouettes des individus avec k-means

Figure 6.8: Clusters et silhouettes des individus avec k-means

Interprétation :

D’après la silhouette, on voit qu’il est assez difficile de classer les golfeurs. En effet, les silhouettes sont globalement faibles. Cela se voit aussi sur le graphe des clusters où on voit les ellipses se chevaucher. On en déduit que les golfeurs ont finalement des niveaux de jeu assez proches.

6.3 Algorithme des k-médoïdes

Enfin, on met en oeuvre l’algorithme des k-médoïdes sur le même jeu de données avec le même nombre de classes. La figure 6.10 montre le clustering obtenu.

# k-medoids
pga.pam <- pam(pga.scale, k)
plot(pga.pam)

# Comparison with k-means
table(pga.pam$clustering, pga.kmeans$cluster)
##    
##      1  2  3  4
##   1 12  0  0  0
##   2  0 13  2  0
##   3  5  4 29  6
##   4  0  4  6 19
Graphes des clusters et silhouettes avec PAMGraphes des clusters et silhouettes avec PAM

Figure 6.9: Graphes des clusters et silhouettes avec PAM

fviz_cluster(pga.pam, pga.2020, ellipse.type = "convex", repel = TRUE,
             labelsize = 7, main = "Cluster plot with PAM, PGA Tour 2020")
Clustering avec l'algorithme des k-médoïdes

Figure 6.10: Clustering avec l’algorithme des k-médoïdes

Sur la figure 6.10, l’axe des abscisses est inversé : les meilleurs golfeurs sont à gauche. En comparant avec la méthode des k-means, on voit que les résultats sont similaires.

6.4 HCPC

On considère désormais le jeu de données complet. On réalise un clustering des golfeurs à partir de l’ACP réalisée dans la partie 5.

On lance donc la commande HCPC (STHDA, s. d.h) de FactoMineR sur l’ACP en laissant choisir le nombre de classes.

# Compute hierarchical clustering on principal components
pga.hcpc <- HCPC(pga.pca, nb.clust = -1, graph = FALSE)

# Dendrogram
# fviz_dend(pga.hcpc, rect = TRUE, rect_fill = TRUE, show_labels = FALSE,
#           main = "Cluster Dendrogram with HCPC, PGA 2010-2020")

# Individuals on the principal component with clustering
fviz_cluster(pga.hcpc, show.clust.cent = TRUE, geom = "point",
             main = "Cluster plot with HCPC, PGA Tour 2010-2020")

# Principal components + tree
plot(pga.hcpc, choice = "3D.map")
Clustering avec HCPCClustering avec HCPC

Figure 6.11: Clustering avec HCPC

Interprétation : avec cette méthode, on ne trouve que 3 classes. Ici, le groupe des favoris et celui des outsiders ne forment qu’un seul groupe (en bleu). Sinon, on retrouve les deux autres groupes : celui des golfeurs puissants (en rouge) et celui des golfeurs précis (en vert). Les résultats trouvés sur l’année 2020 sont donc cohérents avec ceux trouvés sur la période 2010-2020.

7 Classification supervisée

Dans cette partie, le but est de prédire si un golfeur a gagné un tournoi ou non en fonction des autres statistiques. Pour cela, on transforme la variable quantitative WINS, indiquant le nombre de victoires, en une variable qualitative à 2 modalités :

  • 0 indique que le golfeur n’a pas gagné de tournoi cette saison
  • 1 indique que le golfeur a gagné au moins un tournoi cette saison

On garde seulement les variables quantitatives pour mettre en oeuvre le plus de méthodes de classification possible.

# Remove PLAYER.NAME, COUNTRY, YEAR
pga.classif <- pga[3:14]
# WINS variable qualitative à 2 modalités
pga.classif$WINS <- as.factor(as.integer(pga.classif$WINS >= 1))
str(pga.classif)
## 'data.frame':    1100 obs. of  12 variables:
##  $ POINTS    : num  1846 1697 1691 1629 1593 ...
##  $ WINS      : Factor w/ 2 levels "0","1": 2 2 2 2 2 1 2 2 1 2 ...
##  $ ROUNDS    : int  72 73 76 76 78 89 85 77 97 88 ...
##  $ SCORING   : num  69.8 69.7 69.8 70 69.9 ...
##  $ DD        : num  288 283 276 299 288 ...
##  $ DA        : num  60.2 68.5 71 52.7 65.2 ...
##  $ GIR       : num  67.9 68.3 67.1 65.1 66.3 ...
##  $ SS        : num  48.4 54.8 48 53.6 58.8 ...
##  $ SCRAMBLING: num  61.2 63.8 63.5 61.8 63.9 ...
##  $ PA        : num  1.77 1.75 1.77 1.76 1.76 ...
##  $ BA        : num  3.79 3.89 3.64 3.89 3.87 3.76 3.54 3.94 3.97 3.52 ...
##  $ MONEY     : num  4558861 4190235 4809622 3821733 3603331 ...
# Distribution de la variable WINS
table(pga.classif$WINS)
## 
##   0   1 
## 759 341

On crée un échantillon d’apprentissage et un échantillon test à partir des données. L’échantillon test contient 20% de données, le reste constituant l’échantillon d’apprentissage.

# Création d'un échantillon d'apprentissage et d'un échantillon test
data <- pga.classif
set.seed(1)
n <- nrow(data)
p <- ncol(data) - 1
test.ratio <- .2  # ratio of test/train samples
n.test <- round(n * test.ratio)
tr <- sample(1:n, n.test)
data.test <- data[tr,]
data.train <- data[-tr,]
# Distribution de la variable WINS dans l'échantillon d'apprentissage
table(data.train$WINS)
## 
##   0   1 
## 616 264

On compare plusieurs méthodes de classification supervisée en donnant pour chacune :

  • La table de confusion
  • Le taux de bons classements (accuracy)
  • L’aire sous la courbe ROC

7.1 LDA et QDA

On met en oeuvre les méthodes LDA et QDA sur l’échantillon d’apprentissage.

# Modèle
lda.res <- lda(WINS ~ ., data.train)
qda.res <- qda(WINS ~ ., data.train)

# Prédiction
lda.class <- predict(lda.res, newdata = data.test)$class
qda.class <- predict(qda.res, newdata = data.test)$class

# Table de confusion
table(lda.class, data.test$WINS)
##          
## lda.class   0   1
##         0 139  33
##         1   4  44
table(qda.class, data.test$WINS)
##          
## qda.class   0   1
##         0 132  27
##         1  11  50
# Accuracy
lda.accuracy <- mean(lda.class == data.test$WINS)
lda.accuracy
## [1] 0.8318182
qda.accuracy <- mean(qda.class == data.test$WINS)
qda.accuracy
## [1] 0.8272727
# Aire sous la courbe ROC
lda.pred <- predict(lda.res, newdata = data.test)$posterior[, 2]
lda.roc <- roc(data.test$WINS, lda.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
lda.roc$auc
## Area under the curve: 0.9417
qda.pred <- predict(qda.res, newdata = data.test)$posterior[, 2]
qda.roc <- roc(data.test$WINS, qda.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
qda.roc$auc
## Area under the curve: 0.9107

On utilise maintenant une procédure de sélection de modèle pour la LDA et la QDA.

# Sélection de modèle
lda.stepwise <- stepclass(WINS ~ ., data = data.train, method = "lda",
                          direction = "backward")
##  `stepwise classification', using 10-fold cross-validated correctness rate of method lda'.
## 880 observations of 11 variables in 2 classes; direction: backward
## stop criterion: improvement less than 5%.
## correctness rate: 0.86591;  starting variables (11): POINTS, ROUNDS, SCORING, DD, DA, GIR, SS, SCRAMBLING, PA, BA, MONEY 
## correctness rate: 0.87386;  out: "PA";  variables (10): POINTS, ROUNDS, SCORING, DD, DA, GIR, SS, SCRAMBLING, BA, MONEY 
## correctness rate: 0.87614;  out: "SCRAMBLING";  variables (9): POINTS, ROUNDS, SCORING, DD, DA, GIR, SS, BA, MONEY 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##       0.000       0.000       1.107
lda.stepwise
## method      : lda 
## final model : WINS ~ POINTS + ROUNDS + SCORING + DD + DA + GIR + SS + BA + 
##     MONEY
## <environment: 0x56242d6e5878>
## 
## correctness rate = 0.8761
qda.stepwise <- stepclass(WINS ~ ., data = data.train, method = "qda",
                          direction = "backward")
##  `stepwise classification', using 10-fold cross-validated correctness rate of method qda'.
## 880 observations of 11 variables in 2 classes; direction: backward
## stop criterion: improvement less than 5%.
## correctness rate: 0.85455;  starting variables (11): POINTS, ROUNDS, SCORING, DD, DA, GIR, SS, SCRAMBLING, PA, BA, MONEY 
## correctness rate: 0.86364;  out: "MONEY";  variables (10): POINTS, ROUNDS, SCORING, DD, DA, GIR, SS, SCRAMBLING, PA, BA 
## correctness rate: 0.86818;  out: "SS";  variables (9): POINTS, ROUNDS, SCORING, DD, DA, GIR, SCRAMBLING, PA, BA 
## correctness rate: 0.87273;  out: "GIR";  variables (8): POINTS, ROUNDS, SCORING, DD, DA, SCRAMBLING, PA, BA 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##       0.000       0.000       1.267
qda.stepwise
## method      : qda 
## final model : WINS ~ POINTS + ROUNDS + SCORING + DD + DA + SCRAMBLING + PA + 
##     BA
## <environment: 0x56243a840c50>
## 
## correctness rate = 0.8727
# Modèle
lda.stepwise.res <- lda(lda.stepwise$formula, data = data.train)
qda.stepwise.res <- lda(qda.stepwise$formula, data = data.train)

# Prédiction
lda.stepwise.class <- predict(lda.stepwise.res, newdata = data.test)$class
qda.stepwise.class <- predict(qda.stepwise.res, newdata = data.test)$class

# Accuracy
lda.stepwise.accuracy <- mean(lda.stepwise.class == data.test$WINS)
lda.stepwise.accuracy
## [1] 0.8227273
qda.stepwise.accuracy <- mean(qda.stepwise.class == data.test$WINS)
qda.stepwise.accuracy
## [1] 0.8181818
# Aire sous la courbe ROC
lda.stepwise.pred <- predict(lda.stepwise.res,
                             newdata = data.test)$posterior[, 2]
lda.stepwise.roc <- roc(data.test$WINS, lda.stepwise.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
lda.stepwise.roc$auc
## Area under the curve: 0.9392
qda.stepwise.pred <- predict(qda.stepwise.res,
                             newdata = data.test)$posterior[, 2]
qda.stepwise.roc <- roc(data.test$WINS, qda.stepwise.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
qda.stepwise.roc$auc
## Area under the curve: 0.9431

7.2 CART

On lance l’algorithme CART sur les données d’apprentissage pour trouver un arbre optimal.

set.seed(1)

# Tree
tree <- rpart(WINS ~ ., data.train,
              control = rpart.control(minsplit = 5, cp = 0))
printcp(tree)
## 
## Classification tree:
## rpart(formula = WINS ~ ., data = data.train, control = rpart.control(minsplit = 5, 
##     cp = 0))
## 
## Variables actually used in tree construction:
##  [1] BA         DA         DD         GIR        MONEY      PA        
##  [7] POINTS     ROUNDS     SCORING    SCRAMBLING SS        
## 
## Root node error: 264/880 = 0.3
## 
## n= 880 
## 
##           CP nsplit rel error  xerror     xstd
## 1  0.4053030      0  1.000000 1.00000 0.051493
## 2  0.0757576      1  0.594697 0.61364 0.043549
## 3  0.0265152      3  0.443182 0.53788 0.041336
## 4  0.0151515      4  0.416667 0.54545 0.041570
## 5  0.0094697      7  0.371212 0.54167 0.041453
## 6  0.0075758     11  0.333333 0.54167 0.041453
## 7  0.0063131     20  0.265152 0.56439 0.042141
## 8  0.0056818     24  0.238636 0.54167 0.041453
## 9  0.0050505     33  0.185606 0.56061 0.042028
## 10 0.0037879     41  0.143939 0.56439 0.042141
## 11 0.0018939     53  0.098485 0.64773 0.044461
## 12 0.0012626     58  0.087121 0.65909 0.044754
## 13 0.0000000     61  0.083333 0.65909 0.044754
plotcp(tree)

# Optimal cp
cp.opt <- tree$cptable[which.min(tree$cptable[, "xerror"]), "CP"]

# Prune tree
tree.opt <- prune(tree, cp.opt)
rpart.plot(tree.opt, type = 4)
Arbre CART optimalArbre CART optimal

Figure 7.1: Arbre CART optimal

On prédit les classes des individus de l’échantillon test, et on évalue la méthode.

# Prédiction
cart.class <- predict(tree.opt, newdata = data.test, type = "class")

# Table de confusion
table(cart.class, data.test$WINS)
##           
## cart.class   0   1
##          0 130  33
##          1  13  44
# Accuracy
cart.accuracy <- mean(cart.class == data.test$WINS)
cart.accuracy
## [1] 0.7909091
# Aire sous la courbe ROC
cart.pred <- predict(tree.opt, newdata = data.test, type = "prob")[, 2]
cart.roc <- roc(data.test$WINS, cart.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
cart.roc$auc
## Area under the curve: 0.7768

7.3 Random forest

On lance l’algorithme de Random forest sur l’échantillon d’apprentissage.

# Modèle
rf <- randomForest(WINS ~ ., data.train)

# Prédiction
rf.class <- predict(rf, data.test, type = "class")

# Table de confusion
table(rf.class, data.test$WINS)
##         
## rf.class   0   1
##        0 132  23
##        1  11  54
# Accuracy
rf.accuracy <- mean(rf.class == data.test$WINS)
rf.accuracy
## [1] 0.8454545
# Aire sous la courbe ROC
rf.pred <- predict(rf, newdata = data.test, type = "prob")[, 2]
rf.roc <- roc(data.test$WINS, rf.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
rf.roc$auc
## Area under the curve: 0.901

7.4 Régression logistique

On lance une régression logistique sur l’échantillon d’apprentissage.

# Modèle
GLM <- glm(WINS ~ ., family = binomial, data.train)

# Prédiction
glm.pred <- predict(GLM, data.test, type = "response")
# Si pred > 0.5 class 1; Sinon class 0
glm.class <- ifelse(glm.pred > 0.5, 1, 0)

# Table de confusion
table(glm.class, data.test$WINS)
##          
## glm.class   0   1
##         0 132  28
##         1  11  49
# Accuracy
glm.accuracy <- mean(glm.class == data.test$WINS)
glm.accuracy
## [1] 0.8227273
# Aire sous la courbe ROC
glm.roc <- roc(data.test$WINS, glm.pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
glm.roc$auc
## Area under the curve: 0.9431

7.5 Comparaison

Pour conclure, on compare les méthodes précédentes avec les taux de bons classements (Accuracy) et l’aire sous la courbe ROC (AUC).

results <- matrix(NA, ncol = 7, nrow = 2)
rownames(results) <- c('Accuracy', 'AUC')
colnames(results) <- c('LDA', 'LDA.STEP', 'QDA', 'QDA.STEP',
                       'CART', 'RF', 'GLM')

# 1ère ligne : Accuracy
results[1,] <- c(lda.accuracy, lda.stepwise.accuracy, qda.accuracy,
                 qda.stepwise.accuracy, cart.accuracy, rf.accuracy,
                 glm.accuracy)

# 2nde ligne : AUC
results[2,] <- c(lda.roc$auc, lda.stepwise.roc$auc, qda.roc$auc,
                 qda.stepwise.roc$auc, cart.roc$auc, rf.roc$auc, glm.roc$auc)
results
##                LDA  LDA.STEP       QDA  QDA.STEP      CART        RF       GLM
## Accuracy 0.8318182 0.8227273 0.8272727 0.8181818 0.7909091 0.8454545 0.8227273
## AUC      0.9416947 0.9392426 0.9107256 0.9431478 0.7767687 0.9009627 0.9431478
# Meilleure méthode
apply(results, 1, which.max)
## Accuracy      AUC 
##        6        7
# ROC plot
ggroc(list(LDA = lda.roc, LDA.STEP = lda.stepwise.roc, QDA = qda.roc,
           QDA.STEP = qda.stepwise.roc, CART = cart.roc, RF = rf.roc,
           GLM = glm.roc)) +
  ggtitle('ROC Curves') +
  theme_minimal()
Courbes ROC

Figure 7.2: Courbes ROC

D’après le tableau et les courbes ROC de la figure 7.2, les méthodes sont plutôt efficaces pour déterminer si un golfeur a gagné un tournoi ou non pendant une saison : les taux de bons classements et les aires sous la courbe ROC sont proches de 1. Les meilleures méthodes semblent être la régression logistique, la LDA, la QDA et Random forest. La méthode CART donne de moins bons résultats.

7.6 Importance des variables

Pour finir, on regarde les scores d’importance de chaque variable pour l’algorithme de Random forest et pour la méthode CART.

# On crée un dataframe en triant par importance
rf.ord <- order(rf$importance, decreasing = TRUE)
rf.importance.df <- data.frame(Variable = rownames(rf$importance)[rf.ord],
                               Importance = rf$importance[rf.ord])

cart.importance.df <- data.frame(tree.opt$variable.importance)
names(cart.importance.df)[1] <- "Importance"
cart.importance.df$Variable <- row.names(cart.importance.df)

# Rend la colonne Variable ordonnée
rf.importance.df$Variable <- factor(rf.importance.df$Variable,
                                    levels = rf.importance.df$Variable)
cart.importance.df$Variable <- factor(cart.importance.df$Variable,
                                      levels = cart.importance.df$Variable)

# Barplot
ggplot(data = rf.importance.df, aes(x = Variable,
                                    y = Importance, fill = Variable)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = round(Importance)), hjust = 1.6, color = "white", size = 3.5) +
    ggtitle("Importance of variables for Random forest") +
    coord_flip()

ggplot(data = cart.importance.df, aes(x = Variable,
                                      y = Importance, fill = Variable)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = round(Importance)), hjust = 1.6, color = "white", size = 3.5) +
    ggtitle("Importance of variables for CART") +
    coord_flip()
Importances des variablesImportances des variables

Figure 7.3: Importances des variables

D’après la figure 7.3, les variables les plus importantes pour expliquer WINS sont POINTS, MONEY et SCORING. En effet, un golfeur ayant gagné beaucoup de points pour la FedExCup et ayant remporté beaucoup de gains a plus de chance d’avoir gagné un tournoi au cours de la saison.

8 Conclusion

Au cours de cette analyse, nous avons apporté des réponses aux questions posées. Voici un récapitulatif des principaux résultats :

  1. Les statistiques étudiées ne sont pas toutes distribuées normalement. Cela témoigne des disparités entre l’élite des golfeurs et les autres, en particulier pour le nombre de points et les gains remportés.
  2. Le niveau global des golfeurs a progressé sur ces 10 dernières années.
  3. Les États-Unis sont toujours largement plus représentés que les autres pays.
  4. Plusieurs liens existent entre les variables du jeu de données.
  5. Les variables ROUNDS, DD, GIR, SS, SCRAMBLING, PA et BA permettent d’expliquer le score moyen d’un golfeur.
  6. On peut classer les golfeurs en 3 ou 4 groupes : les golfeurs complets (favoris et outsiders), puissants et précis.
  7. Il est possible de déterminer si un golfeur a gagné un tournoi ou non en fonction des autres statistiques. Les méthodes essayées donnent des résultats assez proches. Les variables les plus importantes pour expliquer WINS sont POINTS, MONEY et SCORING.

Références

Chang, Winston. s. d. « ggplot2 colors ». http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/.

D. Chessel, E. Madigou. 2003. « Typologie de golfeurs ». https://pbil.univ-lyon1.fr/R/pdf/pps007.pdf.

Husson, François. 2008. « FactoMineR : Classification Hiérarchique sur Composantes Principales ». http://factominer.free.fr/factomethods/classification-hierarchique-sur-composantes-principales.html.

« PGA Tour - Official Home of Golf ». s. d. https://www.pgatour.com/.

———. s. d.b. « ggcorrplot: Visualization of a correlation matrix using ggplot2 ». http://www.sthda.com/english/wiki/ggcorrplot-visualization-of-a-correlation-matrix-using-ggplot2.

———. s. d.c. « ggfortify : Extension to ggplot2 to handle some popular packages ». http://www.sthda.com/english/wiki/ggfortify-extension-to-ggplot2-to-handle-some-popular-packages-r-software-and-data-visualization.

———. s. d.e. « ggplot2 - Easy way to mix multiple graphs on the same page ». http://www.sthda.com/english/wiki/wiki.php?id_contents=7930.

———. s. d.f. « ggplot2 graphique linéaire : Guide de démarrage rapide ». http://www.sthda.com/french/wiki/ggplot2-graphique-lineaire-guide-de-demarrage-rapide-logiciel-r-et-visualisation-de-donnees.

Wikipédia. s. d.a. « FedEx Cup ». https://en.wikipedia.org/wiki/FedEx_Cup.

———. s. d.b. « PGA Tour ». https://en.wikipedia.org/wiki/PGA_Tour.

———. s. d.c. « Vocabulaire du golf ». https://fr.wikipedia.org/wiki/Vocabulaire_du_golf.

Yihui Xie, Garrett Grolemund, J.J. Allaire. 2021. R Markdown: The Definitive Guide. https://bookdown.org/yihui/rmarkdown/.