Modélisation, prédiction et clustering de données fonctionnelles issues de problématiques sportives

Author

Arthur Leroy

Published

October 15, 2025

Nous utiliserons dans ce document les conventions et les fonctions du tidyverse pour le formattage et la manipulation des données.

library(tidyverse)

I) Traitement basique de données fonctionnelles discrètes : l’exemple de la charge d’entrainement cumulée

Importation des package et données

library(REDI)
data_workload <- read_csv("Data/data_workload_fm.csv")

Mise en pratique

Objectif

Calculer la charge d’entrainement cumulé avec le Robust Exponential Decreasing Index (REDI) et visualiser les résultats.

  1. Convertir le jeu de données dans le format attendu (un tibble/dataframe avec une colonne Input pour la date, et Output pour la charge d’entrainement) par les fonctions du package REDI.

  2. Calculer la valeur de charge d’entrainement cumulée à la dernière date présente dans le jeu de données.

  3. Appliquer la fonction loop_redi() pour calculer le REDI successivement à toutes les dates du jeu de données.

  4. Afficher les résultats dans un graphique sous forme d’une séries temporelle des valeurs succéssives du REDI.

  5. Utiliser la fonction wrapprer redi() pour calculer le REDI pour différentes coefficients de décroissance λ et afficher les résultats.

  6. Analyser les résultats obtenus et discuter de l’impact du choix du coefficient λ, et proposer un critère de choix de ce coefficient dans un contexte sportif.

  7. (Bonus) Introduire des données manquantes dans le jeu de données en retirant aléatoirement 10%, 25% puis 50% des observations. Recalculer les valeurs correspondantes du REDI et afficher les trois courbes sur un même graphique. Discuter de sa robustesse aux données manquantes.

Vous pouvez utiliser les fonctions suivantes :

  1. format_data()
  2. compute_redi()
  3. loop_redi()
  4. plot_redi()
  5. redi()
  6. /
  7. slice_sample()

Assurez-vous de bien comprendre les arguments de chaque fonction en consultant la documentation si nécessaire : https://grenouil.github.io/REDI/

## 1 ##  Convert the dataset to the correct format
db_workload <- format_data(
    data = data_workload,
    input = 'Date',
    output = 'Training load'
    )

## 2 ## Compute REDI for the last date in the dataset
compute_redi(db_workload)

## 3 ## Apply loop_redi() to compute REDI for all dates in the dataset 
db_loop <- loop_redi(data = db_workload)

db_loop

## 4 ## Display results as a time series of REDI values
plot_redi(redi = db_loop,
          x_axis = 'Date',
          y_axis = 'Workload',
          plot_data = TRUE)

## 5 ## Apply redi() on data using a vector of coefficients and display results
db_full_redi <- redi(data = db_workload, coef = c(0.05, 0.1, 0.5), plot = TRUE)

## 6 ## 
# Le choix de lambda peut typiquement venir d'une connaissance expert, en fonction du sport
# étudié, la charge d'activité peut être plus ou moins persistante dans le temps (très
# différent entre un sport d'endurance et un sport de sprint par exemple). Si aucune
# connaissance experte n'est disponible, on peut imaginer choisir lambda en optimisant, à
# partir de données réelles, un critère de performance prédictive (par exemple en
# cross-validation), ou une fonction de coût/perte spécifique au contexte sportif.

## 7 ## Create datasets with 10%, 25%, 50% of missing values 
db_10 <- db_workload %>% slice_sample(prop = 0.9)
db_25 <- db_workload %>% slice_sample(prop = 0.75)
db_50 <- db_workload %>% slice_sample(prop = 0.5)

db_missing <- redi(db_10, coef = 0.1) %>%
  mutate(Missing  = '10') %>%
  bind_rows(
    redi(db_25, coef = 0.1) %>%
      mutate(Missing  = '25')
  ) %>%
  bind_rows(
    redi(db_50, coef = 0.1) %>%
      mutate(Missing  = '50')
  )
 
# Plot results coloured according to the missing data ratio

ggplot(db_missing) +
  geom_line(aes(x = Input, y = REDI, col = Missing)) +
  geom_point(data = db_workload, aes(x = Input, y = Output), 
  col = 'grey', alpha = 0.7) +
  theme_classic()

II) Analyse de données fonctionnelles : B-splines, FPCA et clustering de courbes

Importation des package et données

library(splines)
library(fdapace)
data_swim <- read_csv("Data/db_100m_freestyle.csv")

Mise en pratique

Objectif

  1. Définir une base de splines de degré 3, sur un intervalle couvrant l’ensemble des données entre 10 et 20 ans, avec 4 knots équidistants. Afficher cette base de fonctions sur un graphique.

  2. Extraire les données de performance au cours du temps, spécifiques à un.e athlète, et ajuster un modèle de régression spline sur ses données. Afficher les données et le courbe ajustée.

  3. Répéter l’ajustement par spline pour toutes les athlètes du jeu de données, et afficher les courbes ajustées pour chaque athlète sur un même graphique.

  1. Commencer par extraire un sous jeu de données pour 100 individus pour diminuer le temps de calcul Effectuer une analyse en composantes principales fonctionnelle (FPCA) sur le jeu de données à l’aide du package fdapace. Afficher la proportion de variance expliquée par chaque composante principale.

  2. Choisir le nombre de fonctions principales à retenir et afficher leur courbe correspondante, représentant les modes de variation principaux des données. Analyser ces modes de variation et discuter de leur interprétation dans un contexte sportif.

  3. (Bonus) Effectuer un clustering (via un algorithme des k-means par exemple) des athlètes basé sur les scores des premières composantes principales. Afficher les courbes moyennes ajustées par spline pour chaque cluster et discuter des différences entre les groupes.

  1. La fonction bs() du package splines permet de définir une base de splines. Vous pouvez utiliser la fonction matplot() de base R ou ggplot2 pour afficher les différentes fonctions de la base.

  2. La syntaxe est proche d’une régression linéaire classique, mais en utilisant la fonction bs() pour définir la base de splines dans la formule. Vous pouvez utiliser la fonction predict() pour obtenir les valeurs ajustées sur une grille de points.

  3. La syntaxe n’est pas très intuitive pour faire une boucle sur tous les athlètes, je vous conseille d’aller directement dans l’onglet Solutions pour récupérer et executer le code.

  4. Utiliser la fonction split() pour formater les données par athlète comme attendu par la fonction FPCA.

  5. La part de variance expliquée est contenue dans l’élément $FVE du résultat de la fonction FPCA. Les fonctions principales sont contenues dans l’élément $phi.

  6. Les scores des athlètes sur les fonctions principales sont contenus dans l’élément $xiEst. Vous pouvez utiliser la fonction kmeans() de base R pour effectuer le clustering.

## 1 ## Define a B-spline basis of degree 3 with 4 equidistant knots on the interval [10, 20]

# Define the grid of evaluation and knots
x_grid = seq(10, 20, 0.1)
knots <- seq(12, 18, by = 2)

# Create the B-spline basis
sp_basis = bs(x_grid, degree = 3, knots = c(12, 14, 16, 18), intercept = TRUE)

## Easy (and old school) plot using base R
matplot(x_grid, sp_basis, type = "l", lty = 1, main = "Base de splines cubiques")

## Cool kids version using ggplot2
sp_basis_df = sp_basis %>%
  as_tibble() %>%
  mutate(Input = x_grid) %>%
  pivot_longer(cols = -Input, names_to = 'Basis', values_to = 'Value')

ggplot(sp_basis_df) + geom_line(aes(x = Input, y = Value, col = Basis)) +
  theme_classic() + theme(legend.position = 'none') + ylab('Basis value')

## 2 ## Extract data for a single athlete and fit a spline regression model

# Select a swimmer
swimmer = data_swim %>% filter(ID == 'ID_1')  

# Fit a spline and predict of a grid of values
fit <- lm(Performance ~ bs(Age, degree = 3, knots = knots, Boundary.knots = c(10, 20)), data = swimmer)
pred_spline <- predict(fit, newdata = data.frame(Age = x_grid))

# Display the results
tibble(Age = x_grid, Performance = pred_spline) %>%
  ggplot() + geom_line(aes(x = Age, y = Performance), col = '#DD3497') +
  geom_point(data = swimmer, aes(x = Age, y = Performance)) + theme_classic() +
  xlab('Age') + ylab('Performance')

## 3 ## Fit splines basis to the full dataset

# Fit a spline for each swimmer and predict on a grid of values
pred_all_spline <- data_swim %>% 
  group_by(ID) %>%
  do(data.frame(Age = x_grid,
                y = predict(lm(Performance ~ bs(Age, degree=3, knots=knots, Boundary.knots=c(10,20)),
                data=.),
                newdata = data.frame(Age = x_grid))))

# Display the resulting curves
ggplot(pred_all_spline) + 
  geom_line(aes(Age, y, group = ID)) +
  theme_classic() + xlab('Age') + ylab('Performance')

## 4 ## Perform and FPCA on the full dataset

# Subset the data to 100 swimmers for computational time
sub_db_swim = data_swim %>% filter(ID %in% paste0('ID_', 1:100))

# Format the data and apply FPCA
x_grid <- split(sub_db_swim$Age, sub_db_swim$ID)
y_grid <- split(sub_db_swim$Performance, sub_db_swim$ID)
fpca <- FPCA(y_grid, x_grid)

## 5 ## Display the results of the FPCA

# Display the variance explained by each component
barplot(fpca$FVE, names.arg = 1:length(fpca$FVE), main = "Variance expliquée (%)")

# Display the principal components
matplot(fpca$workGrid, fpca$phi[,1], type="l", lty=1)

## 6 ## Clustering based on the FPCA scores

# Extract the scores on the first component
scores <- fpca$xiEst[,1]

# Perform a k-means clustering on the scores
cl <- kmeans(scores, 3)$cluster

# Add the cluster information to the dataset 
swim_with_clust <- data.frame(ID = names(x_grid), Cluster = cl) %>%
  left_join(sub_db_swim, by = "ID")~

# Display the clusters of curves
ggplot(swim_with_clust) +
  geom_smooth(aes(Age, Performance, col=factor(Cluster))) + 
  theme_classic() + labs(col='Cluster')

III) Formule de Bayes et modélisation probabiliste

Cette partie s’effectue sur feuille, sans code R.

Objectif

Utiliser la formule de Bayes pour calculer des probabilités conditionnelles dans le cas dicret et continu.

Imaginez qu’un.e athlète, à chaque séance d’entrainement, a 1 chance sur 1000 de se rompre le tendon d’Achille. Un nouveau test de prévention, lancé sur le marché, permet de détecter correctement 99% des ruptures du tendon d’Achille à l’avance si le test était positif. Le test en question permet même d’identifier correctement que 99% des séances ne provoqueront pas de blessure au tendon s’il est négatif.

Si cet.te athlète obtient un résultat positif au test avant une séance, pensez vous que c’est complétement inconscient de prendre ce risque ?

Quelle est la probabilité, connaissant ce résultat positif au test, qu’il/elle se rompe le tendon d’Achille lors de son entrainement ?

Rappelez vous de la formule des probabilités totales. La probabilité “totale” d’obtenir un test positif est la somme de la probabilité d’obtenir un test positif puis de se blesser et de la probabilité d’obtenir un test positif sans se blesser ensuite.

Ou, pour les puristes (pour tout évènement T et B):

P(T)=P(T∩B)+P(T∩B¯)

Allez, si vraiment le terme d’évidence ne l’est pas tant que ça, rappelez vous également que l’intersection de deux évènements peut toujours s’écrire comme la probabilité conditionnelle de l’un sachant l’autre fois la probabilité de l’autre. Ou plus formellement : P(T∩X)=P(T|X)×P(X),    ∀X

Allez un dernier indice, pour rappeler que pour tout évènement binaire (par exemple un test positif/négatif), la probabilité de l’évènement est égal à 1 moins la probabilité de son complémentaire. Ou plus formellement :

P(T)=1−P(T¯)

Ce qui reste vrai pour les évènements conditionnels :

P(T|X)=1−P(T¯|X),    ∀X

Commençons par définir mathématiquement les évènements :

  • B = “se rompre le tendon d’Achille”,
  • B¯ = “ne pas se rompre le tendon d’Achille”,
  • T = “test positif”,
  • T¯ = “test négatif”.

A partir de l’énoncé, nous avons les probabilités suivantes :

  • P(B)=0.001, “A chaque séance d’entrainement, 1 chance sur 1000 de se rompre le tendon d’Achille”,
  • P(T|B)=0.99, “permet de détecter correctement 99% des ruptures du tendon d’Achille à l’avance si le test était positif”,
  • P(T¯|B¯)=0.99, “identifier correctement que 99% des séances ne provoqueront pas de blessure au tendon s’il est négatif”.

Nous avons donc à notre disposition les termes de prior P(B) et de vraisemblance P(T|B), il ne reste qu’à calculer le terme d’évidence P(T) pour pouvoir appliquer la formule de Bayes.

Pour cela, nous allons utiliser la formule des probabilités totales : P(T)=P(T∩B)+P(T∩B¯)

Puis en utilisant la formule de l’intersection, nous obtenons : P(T)=P(T|B)×P(B)+P(T|B¯)×P(B¯)=P(T|B)×P(B)+(1−P(T¯|B¯))×(1−P(B))=0.99×0.001+(1−0.99)×(1−0.001)=0.00099+0.00999=0.01098

Nous pouvons maintenant appliquer la formule de Bayes pour obtenir la probabilité conditionnelle recherchée :

P(B|T)=P(T|B)×P(B)P(T)=0.99×0.0010.01098≈0.09016

Donc, même avec un test positif, la probabilité que l’athlète se rompe le tendon d’Achille lors de son entrainement est d’environ 9.016%. Cela signifie qu’il y a environ 90.984% de chances de ne pas se blesser pas malgré le résultat positif du test.

IV) Apprentissage par processus gaussiens :

Importation des package et données

library(MagmaClustR)
data_gps <- read_csv("Data/data_GPS_rugby_with_speed.csv")

Mise en pratique

Objectif

Se familiariser avec l’apprentissage par processus gaussiens pour la modélisation et la prédiction de données fonctionnelles. Faire des prédictions probabilistes pour des courbes avec leur incertitude associée. La documentation du package MagmaClustR, permettant cela, est disponible ici : https://arthurleroy.github.io/MagmaClustR/

  1. Convertir le jeu de données ‘data_swim’ au format attendu par les fonctions du package MagmaClustR. La colonne représentant le temps doit être nommée ‘Input’ et la colonne représentant les performances ‘Output’, tandis que la colonne identifiant les individus doit être nommée ‘ID’.

  2. Filtrer les données d’un.e seul.e individu.e et retirer la colonne ‘Gender’ du jeu de données. Découper aléatoirement le jeu de données en un ensemble d’entrainement (80% des données) et un ensemble de test (20% des données).

  3. Entrainer un modèle de processus gaussien sur ces données en optimisant les hyper-paramètres du noyau de covariance par défaut. Afficher la valeur des hyper-paramètres optimaux obtenus.

  4. Calculer les prédictions sur les données de test en obtenant la loi normale a posteriori, caractérisée par sa moyenne et l’incertitude associée (variance) en tout point de cette grille. Définir ensuite une fine grille de points d’observations (Input) sur un intervalle plus large, typiquement entre 10 et 20 ans (en veillant à ne pas dépasser une vecteur de taille > 2000 maximum environ), puis relancer la prédiction sur cette grille.

  5. Utiliser la fonction plot_gp() les résultats sous forme d’un graphique personnalisé. Commencer par afficher la prédiction minimale avec la moyenne prédite et une bande de confiance à 95%, puis via les arguments supplémentaires, ajouter les données d’entarinement. Ajouter ‘à la main’ au graphique ggplot les points de test, en rouge, pour visualiser et discuter les erreurs de prédiction.

  6. Répéter les étapes 2 à 5 pour un.e nouvel.le d’individu.e du jeu de données. Essayer de définir une moyenne a priori de votre choix (via l’argument prior_mean) et d’utiliser un autre kernel (par exemple linéaire via l’argument kern = ‘LIN’ ou périodique via kern = ‘PER’). Essayer également de sélectionner des points de test de la fin de l’intervalle d’observation pour voir comment le modèle se comporte en extrapolation. Discuter des résultats obtenus, et le rôle de l’incertitude des prédictions.

  7. (Bonus) Utiliser les fonctions pred_gif() et plot_gif() pour créer votre propre animation de la prédiction par processus gaussien et son adaptation aux nouvelles données collectées au fil du temps.

  8. (Bonus) Utiliser le jeu de données data_GPS_rugby_with_speed pour répéter les étapes précédentes sur des données en 2-dimensions. En terme de format, toute colonne supplémentaire (autre que ID, Input, Output) est considéré comme un Input supplémentaire, quelque soit son nom. Vous pouvez ainsi modéliser et prédire la vitesse d’un joueur de rugby en fonction de sa position sur le terrain (coordonnées X et Y).

Vous pouvez utiliser :

  1. La fonction rename() de dplyr pour renommer les colonnes du jeu de données.

  2. Les fonctions filter(), select(), sample() et slice() pour extraire les données d’un individu, supprimer la colonne Gender et découper le jeu de données en entrainement/test.

  3. La fonction train_gp() directement sur le dataset pour entrainer le modèle de processus gaussien.

  4. La fonction pred_gp() pour calculer les prédictions sur les données de test et sur une grille d’inputs. Les arguments importants sont data (les données d’entrainement), hp (les hyper-paramètres optimaux obtenus à l’étape précédente), grid_inputs (une grille d’inputs où évaluer la prédiction).

  5. La fonction plot_gp() pour afficher les résultats. Les arguments importants sont pred (les résultats de la prédiction), data (les données d’entrainement, optionnel).

  6. Les fonction head() et tail() pour choisir des points de train/test en début/fin d’intervalle, et les arguments pour la moyenne a priori (prior_mean), le kernel (kern, options possibles 'PERIO' ou 'LIN').

  7. La fonction pred_gif() pour créer une liste de prédictions successives, et la fonction plot_gif() pour afficher l’animation. L’argument important est export_gif = TRUE pour sauvegarder l’animation au format GIF.

# Define a seed for reproducibility
set.seed(42)

## 1 ## Format the dataset 
db_swim = data_swim %>% 
  rename(Input = Age, Output = Performance)

## 2 ## Extract data of one swimmer and split data into train and test sets
swimmer = db_swim %>% 
  filter(ID == "ID_42") %>%
  select(- Gender) 

train_indices = sample(1:nrow(swimmer), size = 0.8*nrow(swimmer))

swimmer_train = swimmer %>%
  slice(train_indices)

swimmer_test = swimmer %>%
  slice(-train_indices)

## 3 ## Train a default GP model on the training set and display hyper-parameters values

mod_gp = train_gp(data = swimmer_train)

print(mod_gp)

## 4 ## Predict on the test set and display results

pred_gp = pred_gp(data = swimmer_train, hp = mod_gp)

grid_inputs = seq(10, 20, 0.1)

pred_gp_grid = pred_gp(data = swimmer_train, hp = mod_gp, grid_inputs = grid_inputs)

## 5 ## Plot results with the dedicated function

# Minimal visualisation
plot_gp(pred_gp_grid)

# Add training points
plot_gp(pred_gp_grid, data = swimmer_train)

# Add test points
plot_gp(pred_gp_grid, data = swimmer_train) + 
  geom_point(data = swimmer_test, aes(x = Input, y = Output), col = 'red')

## 6 ## Try with another swimmer in a forecating context

# Format the train/test sets
swimmer2 = db_swim %>% 
  filter(ID == "ID_17") %>%
  select(- Gender) 

swimmer_train2 = swimmer2 %>%
  head(round(0.8*nrow(swimmer2)))
swimmer_test2 = swimmer2 %>%
  tail(-round(0.8*nrow(swimmer2)))

# Train a periodic kernel and forecast
mod_gp_perio = train_gp(data = swimmer_train2, prior_mean = 50,  kern = 'PERIO')
pred_gp_perio = pred_gp(data = swimmer_train2, hp = mod_gp_perio, kern = 'PERIO', grid_inputs = grid_inputs, plot = FALSE)
plot_gp(pred_gp_perio, data = swimmer_train2) + 
  geom_point(data = swimmer_test2, aes(x = Input, y = Output), col = 'red')

# Train a linear kernel and forecast
mod_gp_lin = train_gp(data = swimmer_train2, prior_mean = 50, kern = 'LIN')
pred_gp_lin = pred_gp(data = swimmer_train2, hp = mod_gp_lin, kern = 'LIN', grid_inputs = grid_inputs, plot = FALSE)
plot_gp(pred_gp_lin, data = swimmer_train2) + 
  geom_point(data = swimmer_test2, aes(x = Input, y = Output), col = 'red')

## 7 ## Create a sequence of predictions to display an animation
 
loop_pred = pred_gif(data = swimmer, hp = mod_gp, grid_inputs = grid_inputs)

gif = plot_gif(loop_pred, data = swimmer, export_gif = TRUE)

## 8 ## Try with a 2D example

# Import the dataset and format the train/test sets
db_gps = data_gps %>%
  rename(ID = player, Input = x, Input2 = y, Time = frame_id, Output = speed) %>%
  select(ID, Input, Input2, Output)

player = db_gps %>%
  filter(ID == 9) 

player_train = player %>%
  head(round(0.8*nrow(player)))
player_test = player %>%
  tail(-round(0.8*nrow(player)))

# Train a 2D-Input GP model for this player
mod_gp_2D = train_gp(data = player_train)

# Create a 2D grid of inputs for predictions with a dedicated function
grid_inputs_2D = expand_grid_inputs(
  Input = seq(60, 95, 1),
  Input2 = seq(0, 50, 1)
)

# Make predictions on the grid
pred_gp_2D = pred_gp(data = player_train, hp = mod_gp_2D, grid_inputs = grid_inputs_2D, plot = FALSE)

# Plot results and add 2D testing points with labels of true Output values
plot_gp(pred_gp_2D, data = player_train) + 
  geom_text(data = player_test,
     aes(x = Input, y = Input2, label = round(Output, 1)), col = "#16e30fff", check_overlap = TRUE)

V) Apprentissage par processus gaussiens multi-tâches :

Objectif

Découvrir le paradime multi-tâches pour l’apprentissage par processus gaussiens, pour la modéliser simultanément de multiples données fonctionnelles issues d’un même phénomène. L’implémentation repose toujours sur le package MagmaClustR dont la documentation est disponible ici : https://arthurleroy.github.io/MagmaClustR/

  1. Reprendre le format attendu par les fonctions du package MagmaClustR (ID, Input, Ouput), en conservant cette fois plusieurs individus pour l’entrainement. Pour conserver des temps de calculs raisonnables, conserver 30 individus. Choisir et extraire un individu de l’ensemble d’entrainement, et séparer ses données en 50% pour la prédiction et 50% pour le testing.

  2. Entrainer un modèle de processus gaussien multi-tâche via l’algorithme Magma sur ces données. Afficher le tableau des hyper-paramètres optimaux obtenus pour les individus, puis la distribution a posteriori du processus moyen (stockée dans l’élément $hyperpost$pred du modèle entrainé) a posteriori via la fonction plot_gp().

  3. Entrainer un nouveau modèle en optimisant des hyper-paramètres individu-spécifiques grâce à l’argument common_hp_i = FALSE. Afficher le nouveau tableau des hyper-paramètres optimaux des individus, ainsi que le graphique du processus moyen a posteriori.

  4. Utiliser à nouveau la grille d’inputs pour calculer la loi a posteriori de l’individu de test en fournissant ses données partiellement observées à la fonction de prédiction. Afficher le tableau des prédictions (moyenne et variance) pour chaque point de la grille, puis relancer la fonction de prédiction en spécifiant l’argument get_hyperpost = TRUE pour obtenir également la loi hyper-posterior du processus moyen.

  5. Utiliser la fonction plot_magma() pour personaliser la visualisation des résultats. Commencer par afficher la prédiction minimale, puis via les arguments supplémentaires, ajouter au graphique les données d’entrainement, les données de prédiction. Constater que le graphique inclus en ligne pointillée le processus moyen issu de l’entrainement multi-tâches (conséquence de l’argument get_hyperpost = TRUE pendant la prédiction). Ajouter au graphique obtenu les points de test, en rouge, pour visualiser et discuter les erreurs de prédiction.

  6. Calculer une nouvelle prédiction en utilisant le modèle entrainé avec l’argument common_hp_i = FALSE et comparer les résultats obtenus dans le cas d’hyper-paramètres partagés et individu-spécifiques. Que constatez-vous ?

  1. (Bonus) Utiliser les fonctions pred_gif() et plot_gif() pour créer votre propre animation de la prédiction par processus gaussien multi-tâches au fil du temps. Comparer cette animation à celle précédemment obtenue pour un processus gaussien simple.

  2. (Bonus) Utiliser le jeu de données data_GPS_rugby_with_speed pour répéter les étapes précédentes sur des données en 2-dimensions. Les données observées en 2D sur une grille très fine peuvent rapidement devenir coûteuses à entrainer dans un cadre multi-tâches. Il peut alors être utile de sous-échantillonner en utilisant la fonction dédiée du package regularise_data(). Vous pouvez ainsi modéliser et prédire la vitesse d’un joueur de rugby en fonction de sa position sur le terrain (coordonnées X et Y). Comparer à nouveau les résultats obtenus avec un processus gaussien simple et multi-tâches.

Vous pouvez utiliser :

  1. La fonction train_magma() pour entraîner un GP multi-tâches sur un sous-ensemble d’individus (data = db).

    • Comparer hyper-paramètres partagés vs spécifiques avec l’argument common_hp = TRUE/FALSE.
    • Inspecter MODEL$hp_i (tibble des HP par individu) et MODEL$hyperpost$pred (processus moyen a posteriori).

  2. La création d’une grille d’inputs avec seq() (1D) ou expand_grid_inputs() (2D) pour les prédictions (grid_inputs = ...).

    • Conserver une taille raisonnable (< ~2000 points) pour limiter le temps de calcul.

  3. La fonction pred_magma() pour la prédiction d’un individu partiellement observé (data = swimmer_pred, trained_model = mod_magma, grid_inputs = ...).

    • Ajouter get_hyperpost = TRUE pour récupérer la loi hyper-postérieure du processus moyen ; utiliser plot = FALSE pour obtenir le tibble (moyenne/variance).

  4. La fonction plot_magma() pour visualiser la distribution prédictive.

    • Paramètres utiles : pred (objet de pred_magma()), data (données de l’individu), data_train (pour afficher les autres individus).
    • Ajouter au besoin les points de test : + geom_point(..., col = 'red').

  5. La comparaison des réglages via un second entraînement avec train_magma(..., common_hp = FALSE) puis pred_magma(...) sur le même individu/grille pour comparer les prédictions.

  6. L’animation avec pred_gif() puis plot_gif(..., export_gif = TRUE, path = "magma.gif") pour sauvegarder un GIF (fournir un chemin .gif valide afin d’éviter la création de PNG temporaires).

# Set a seed for reproducibility
set.seed(123)

## 1 ## Subset 100 swimmers, format, extract a testing individual and split into prediction and testing sets
db = db_swim %>%
  filter(ID %in% paste0('ID_', 1:30)) %>%
  select(- Gender)

swimmer = db_swim %>%
  filter(ID == 'ID_101') %>%
  select(- Gender)

swimmer_pred = swimmer %>%
  head(nrow(swimmer)/2)

swimmer_test = swimmer %>%
  tail(nrow(swimmer)/2)

## 2 ## Train a Multi-Task GP model on the dataset with the algorithm Magma
mod_magma = train_magma(data = db)

# Display the tibble of all individual hyperparameters of the kernels
mod_magma$hp_i

# Display the hyperposterior distribution of the mean process
plot_gp(mod_magma$hyperpost$pred)

## 3 ## Retrain a Multi-Task GP model with individual-specific hyperparameters
mod_magma_dif_hp = train_magma(data = db, common_hp = FALSE)

mod_magma_dif_hp$hp_i
plot_gp(mod_magma_dif_hp$hyperpost$pred)

## 4 ## Make predictions on a grid of inputs for the testing individual
pred = pred_magma(data = swimmer_pred, trained_model = mod_magma, grid_inputs = grid_inputs)

plot_magma(pred, data = swimmer_pred, data_train = db) +
  geom_point(data = swimmer_test, aes(x = Input, y = Output), col = 'red') 

# Compute the predictive distribution and display it as a tibble
pred_magma(
  data = swimmer_pred,
  trained_model = mod_magma,
  grid_inputs = grid_inputs, 
  plot = FALSE)

# Compute the same distribution but keep the hyperposterior distribution of the mean process
pred = pred_magma(
  data = swimmer_pred,
  trained_model = mod_magma,
  grid_inputs = grid_inputs,
  plot = FALSE,
  get_hyperpost = TRUE)

## 5 ## Plot the predictive distribution along with the training and testing data

# Minimal plot with hyperposterior of the mean process (dashed line)
plot_magma(pred)

# Add data from training individuals
plot_magma(pred, data_train = db)

# Add prediction data from the testing individual
plot_magma(pred, data = swimmer_pred, data_train = db)

# Add testing data from the testing individual
plot_magma(pred, data = swimmer_pred, data_train = db) +
  geom_point(data = swimmer_test, aes(x = Input, y = Output), col = 'red') 

## 6 ## Compute predictions with individual-specific hyperparameters

# This will optimise new hyperparameters for the testing individual
pred_dif_hp = pred_magma(
  data = swimmer_pred,
  trained_model = mod_magma_dif_hp,
  grid_inputs = grid_inputs,
  plot = FALSE,
  get_hyperpost = TRUE)

# Display the prediction with individual-specific hyperparameters and compare with the previous one
plot_magma(pred_dif_hp, data = swimmer_pred, data_train = db) +
  geom_point(data = swimmer_test, aes(x = Input, y = Output), col = 'red') 

## 8 ## Create and plot an animation of the prediction over time

loop_pred_magma = pred_gif(
  data = swimmer,
  trained_model = mod_magma,
  grid_inputs = grid_inputs,
  )

gif_magma = plot_gif(loop_pred_magma, data = swimmer, export_gif = TRUE, path = "gif_magma.gif")

## 9 ## Repeat the previous steps with the 2D rugby data 

# Sub-sample the data set by regularising the input space on a grid of size 50x50
sub_db_gps = db_gps %>%
  regularise_data(size_grid = 50)

# Visualise the difference between original and sub-sampled data
ggplot(db_gps) + geom_point(aes(x = Input, y = Input2, col = Output, group = ID)) + theme_classic()
ggplot(sub_db_gps) + geom_point(aes(x = Input, y = Input2, col = Output, group = ID)) + theme_classic()

# Select a player (sub-sampling is not required because prediction is cheap), split into train/test sets
player = db_gps %>%
  filter(ID == 9)

player_train = player %>%
  head(round(0.8*nrow(player)))
player_test = player %>%
  tail(-round(0.8*nrow(player)))

# Train a 2D-Input GP model on the sub-sampled dataset (without the testing player)
mod_magma_2D = train_magma(data = sub_db_gps %>% filter(ID != 9))

# Make predictions on the previous 2D grid of inputs
pred_magma_2D = pred_magma(
  data = player_train,
  trained_model = mod_magma_2D,
  grid_inputs = grid_inputs_2D,
  plot = FALSE)

# Plot results and add 2D testing points with labels of true Output values
plot_gp(pred_magma_2D, data = player_train) + 
  geom_text(data = player_test,
     aes(x = Input, y = Input2, label = round(Output, 1)), col = "#16e30fff", check_overlap = TRUE)

VI) Modélisation, clustering et prédiction simultanés :

Objectif

Etendre le cadre de l’apprentissage GPs multi-tâches pour introduire la prise en compte d’un clustering (via un modèle de mélange de GPs), en complément de la modélisation et de la prédiction de données fonctionnelles. Rappel de la documentation : https://arthurleroy.github.io/MagmaClustR/

  1. Réutiliser la base de données d’entrainement au format attendu, ainsi que l’individu de test, partiellement observé (50% des données pour la prédiction, 50% pour le test). Entrainer un modèle de mélange de GPs multi-tâches via l’algorithme MagmaClust sur ces données, en spécifiant un nombre de clusters égal à 2 et une moyenne a priori égale à la moyenne empirique des Ouputs dans le dataset (via l’argument prior_mean_k).

  2. Afficher le tableau des probabilités d’appartenance à chaque cluster pour chaque individu (via l’élément $hyperpost$mixture), ainsi que la distribution a posteriori chaque processus moyen cluster-spécifique (stockée dans un élément du type $hyperpost$pred$K1 du modèle entrainé) a posteriori via la fonction plot_gp().

  3. Calculer une prédiction ‘générique’ pour un nouvel individu jamais observé en laissant l’argument data=NULL dans la fonction de prédiction. A quoi correspond cette prédiction selon vous, et en quoi est-elle différente des processus moyens cluster-spécifiques ?

  4. Calculer maintenant la prédiction pour l’individu de test partiellement observées sur une grille d’inputs, en spécifiant les arguments get_hyperpost = TRUE pour obtenir la loi hyper-posteriori du processus moyen, ainsi que get_full_cov = TRUE pour la matrice de covariance prédictive complète (et pas seulement la variance ponctuelle). Afficher la probabilité d’appartenance à chaque cluster pour l’individu de test.

  5. Calculer le cluster d’appartenance le plus probable pour tous les individus d’entrainement à l’aide de la fonction data_allocate_cluster(), et afficher le nouveau tibble de données obtenu.

  6. Utiliser la fonction plot_magmaclust() pour personnaliser la visualisation des résultats. Ajouter les données d’entrainement coloriées par cluster d’appartenance (avec l’argument col_clust = TRUE), les données de prédiction, les processus moyens cluster-spécifiques (via l’argument prior_mean = MODEL$hyperpost$mean), puis représenter l’incertitude avec des échantillons plutôt qu’une bande de crédibilité (grâce à l’argument samples = TRUE). Ajouter au graphique les points de test, en vert, pour visualiser et discuter les erreurs de prédiction.

  1. (Bonus) Chercher sur internet la valeur actuelle du record de France du 100m nage libre masculin et féminin. Utiliser la fonction sample_magmaclust() pour générer 1000 échantillons à partir de la prédiction d’un.e nageur.euse utilisant ses données de performance jusqu’à 16 ans. En calculant la proportion d’échantillons dont la trajectoire descend sous le temps du record de France, et déterminer la probabilité qu’il/elle batte ce record entre 16 et 20 ans ?

Vous pouvez utiliser :

  1. La fonction train_magmaclust() pour entraîner un mélange de GPs multi-tâches (data = db, nb_cluster = K). Explorer MODEL$hyperpost$mixture (probas d’appartenance) et MODEL$hyperpost$pred$K1/K2/... (processus moyens cluster-spécifiques).

  2. La fonction pred_magmaclust() pour la prédiction :

    • Générique (nouvel individu sans données) avec data = NULL.
    • Individu partiel avec data = swimmer_pred, grid_inputs = ..., et options get_hyperpost = TRUE (processus moyen) et get_full_cov = TRUE (matrice de covariance complète).

  3. La grille d’inputs comme en section V (seq() / expand_grid_inputs()), en veillant à des tailles raisonnables pour la stabilité et la vitesse.

  4. L’attribution de cluster la plus probable pour les individus d’entraînement avec data_allocate_cluster(MODEL) afin d’obtenir un tibble enrichi (ID + cluster).

  5. La fonction plot_magmaclust() pour visualiser les prédictions du mélange :

    • Paramètres utiles : pred, data (individu), data_train (avec col_clust = TRUE), prior_mean = MODEL$hyperpost$mean, samples = TRUE (incertitude par échantillons).
    • Ajouter les points de test en vert : + geom_point(..., col = 'green').

  6. La comparaison entre prédiction du mélange et prédictions cluster-spécifiques (en changeant cluster ou en traçant MODEL$hyperpost$pred$K1/K2/...) pour discuter des effets d’appartenance.

  7. La simulation avec sample_magmaclust(pred_obj, nb_samples = 1000) pour estimer des probabilités d’événements (ex. franchir un seuil) ; visualiser le seuil avec geom_hline() et fixer une seed pour la reproductibilité.

# Set a seed for reproducibility
set.seed(42)

## 1 ## Train a Magmaclust model on a subset of swimmers

mod_clust = train_magmaclust(data = db, nb_cluster = 2)

## 2 ## Display the probability of cluster memberships for all swimmers

mod_clust$hyperpost$mixture %>% print(n=30)

# Plot the cluster-specific mean process for each cluster
plot_gp(mod_clust$hyperpost$pred$K1)
plot_gp(mod_clust$hyperpost$pred$K2)

## 3 ## Predict on a grid of inputs for a new swimmer with no data point

pred_swim = pred_magmaclust(data = NULL, trained_model = mod_clust)
## This prediction is a mixture of cluster-specific predictions weighted by the membership proportions
## It does not use any individual-specific information, and stands as a 'generic' prediction

## 4 ## Compute predictions for the testing swimmer

pred_clust = pred_magmaclust(
  data = swimmer_pred, 
  trained_model = mod_clust,
  grid_inputs = grid_inputs,
  get_hyperpost = TRUE, 
  get_full_cov = TRUE)

# Display the cluster membership probabilities for the testing swimmer
pred_clust$mixture

## 5 ## Compute the most probable cluster for each individual for the next plot

db_with_clust = data_allocate_cluster(mod_clust)

db_with_clust %>% print(n=30)

## 6 ## Display the mixture of GPs prediction for the testing swimmer

plot_magmaclust(
  pred = pred_clust,
  data = swimmer_pred, 
  data_train = db_with_clust, 
  prior_mean = mod_clust$hyperpost$mean,
  col_clust = TRUE, 
  samples = TRUE) + 
  geom_point(data = swimmer_test, aes(x = Input, y = Output), col = 'green', size = 3)

## 7 ## Compute the probability to break the French record of 100m freestyle before age 20

# French record for 100m freestyle : 46s94 (Male) - 52s74 (Female)
# Search for the best swimmer in the dataset and compte prediction from data up to 16 years
best_swimmer = db %>%
  filter(Output == min(Output)) %>%
  pull(ID)

best_swimmer_pred = db %>%
  filter(ID == best_swimmer, Input <= 16)

best_swimmer_test = db %>%
  filter(ID == best_swimmer, Input > 16)

# Compute the MagmaClust prediction for the best swimmer
pred_best = pred_magmaclust(
  data = best_swimmer_pred,
  trained_model = mod_clust,
  grid_inputs = grid_inputs,
  get_full_cov = TRUE) 

# Visual illustration of the prediction with the record line
plot_magmaclust(
  pred = pred_best,
  data = best_swimmer_pred, 
  data_train = db_with_clust, 
  prior_mean = mod_clust$hyperpost$mean,
  col_clust = TRUE, 
  samples = TRUE) + 
  geom_point(data = best_swimmer_test, aes(x = Input, y = Output), col = 'green') + 
  geom_hline(yintercept = 46.94, lty = 2) +
  scale_y_continuous(breaks = sort(c(seq(40, 80, 10), 46.94)))

# Generate 10 000 sample curves from the posterior predictive distribution
samples_pred = sample_magmaclust(pred_best, nb_samples = 1000)

# Compute the proportion of samples breaking the record after age 16
prop_record = samples_pred %>%
  filter(Output < 46.94, Input >= 16) %>% 
  pull(Sample) %>%
  n_distinct() / 1000 * 100

cat(paste0('Probability to break the French record before age 20 : ', round(prop_record,2), '%\n'))