Problème pour le calcul du Q2 et Q2cum d'une PLS avec plsr()

Bonjour à tous,

Je suis nouveau sur R et cherche actuellement à mettre en place une régression PLS sous R. Pour le moment j'utilise plsr() de la librairie(pls).

Pour valider mon approche, j'ai choisi de travailler sur un jeu de données externe à R (exemple XLSTAT téléchargé sur le lien : https://help.xlstat.com/customer/fr/...stat?b_id=9283).

Mon objectif est de faire la prédiction de la variable expliquée Y "J1" en fonction de toutes les variables explicatives XX présentent dans ce jeu de données. J'ai centré puis réduit les variables XX, j'arrive à retrouver les valeurs de R2X et R2Y mais je n'arrive pas a reproduire le Q2cum (Table 1). J'ai essayé de recalculer manuellement le Q2 et je retrouve exactement les mêmes résultats que R, mais ceci ne concorde pas avec les résultats proposés par XLSTAT (Q2 R toujours > Q2cum proposé par XLSTAT). Je sais également que la formule du Q2cum diffère et prend en compte les produit des rapport (PRESS/SS) pour chacune des variables, mais je ne tombe pas sur les mêmes résultats. J'ai également testé la fonction plsreg1() de la librarie(plsdepot) qui me donne a peu de chose prêt la même chose que plsr().

Je bloque sur le calcul du Q2 et Q2 cum, pouvez-vous me donner un petit coup de pouce ?

Merci par avance,

Baptiste

Le code que j'utilise pour faire ma PLS et déterminer les critères de qualité de ma PLS est le suivant :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

# input m <- read.delim(File, sep=",",header=TRUE) rownames(m) <- m$Jus_orange m$Jus_orange <- NULL   # PLS train ncomp = 4 train_pls <- plsr(m_pls$J1~.,data=m_pls,ncomp=ncomp,scale=TRUE) R2(train_pls) # Valeur de R2Y summary(train_pls) # Valeur de R2X   # PLS CV CV_pls <- plsr(m_pls$J1~.,data=m_pls,ncomp=4,scale=TRUE,validation="LOO") # Cross-Validation Leave-One-Out MSEP(CV_pls) RMSEP(CV_pls) R2(CV_pls,estimate="all") # R2 cross-validated (Q2)[/I][/I]       Table 1 : résultats proposés sur XLSTAT   Comp1 Comp2 Comp3 Comp4 Q² cum 0,334 0,765 0,936 0,994 R²Y cum 0,615 0,949 0,997 1,000 R²X cum 0,574 0,763 0,851 0,939

---

Fichier csv en entrée (dataset) :

Citation:

---

Jus_orange,Glucose,Fructose,Saccharose,Pouvoir_sucrant,pH_brut,pH_apres_centrifuga,Titre,Acide_citrique,Vitamine_C,intensite_odeur,typicite_odeur,caractere_pulpeux,intensite_gout,caractere_acide,caractere_amer,caractere_sucre,J1
pampryl_amb,25.32,27.36,36.45,89.95,3.59,3.55,13.98,0.84,43.44,2.82,2.53,1.66,3.46,3.15,2.97,2.6,2
tropicana_amb,17.33,20,44.15,82.55,3.89,3.84,11.14,0.67,32.7,2.76,2.82,1.91,3.23,2.55,2.08,3.32,2
fruvita_fr,23.65,25.65,52.12,102.22,3.85,3.81,11.51,0.69,37,2.83,2.88,4,3.45,2.42,1.76,3.38,3
joker_amb,32.42,34.54,22.92,90.71,3.6,3.58,15.75,0.95,36.6,2.76,2.59,1.66,3.37,3.05,2.56,2.8,2
tropicana_fr,22.7,25.32,45.8,94.87,3.82,3.78,11.8,0.71,39.5,3.2,3.02,3.69,3.12,2.33,1.97,3.34,4
pampryl_fr,27.16,29.48,38.94,96.51,3.68,3.66,12.21,0.74,27,3.07,2.73,3.34,3.54,3.31,2.63,2.9,3

---

Une proposition de solution malgré le retard.
D'abord, la formule est erronée. Avec

Code:

---

plsr(m_pls$J1~., data = m_pls, ncomp = ncomp, scale = TRUE)

---

ce qui sera évalué est ceci : la variable J1 du jeu de données m_pls va être expliquée en fonction de toutes les autres variables restantes dans m_pls (111 variables). Or, on voulait plutôt l'expliquer en fonctions de seulement 16 variables. Voir help("formula") pour plus de détails sur l'écriture d'une formule.
Ensuite, une seule formule suffit à chaque fois qu'on veut expliquer une seule variable. Mais dès qu'on passe à l'explication de plusieurs variables à la fois, utiliser plutôt une matrice contenant des modèles. Voir help("model.matrix") pour les explications détaillées.
Enfin, la page sur XLSTAT n'indique pas l'algorithme utilisé ; peut-être que c'est celui qui est dans le papier d'où est tiré le jeu de données. La bibliothèque pls a des paramètres par défaut qui indiquent quels algorithmes seront appliqués suivant les contextes ; actuellement, il y a des choses du genre kernelpls, cppls et svdpc comme valeurs par défaut.
Tout cela réuni fait qu'il n'y aura pas à la virgule près les mêmes résultats entre XLSTAT et pls (standard). Voici un exemple de comment aller au bout sans chercher à reproduire nécessairement l'exemple de XLSTAT.

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

library("pls") pls.options() # Consulter les paramètres par défaut; ils sont modifiés avec la même méthode.   orange_dataset <- read.csv("orange_dataset.csv", quote = "", row.names = 1L, stringsAsFactors = FALSE) characteristics_ind <- seq_len(16L)   our_formula_responses <- paste(colnames(orange_dataset[, -characteristics_ind]), collapse = " + ")   our_formula_predictors <- paste(colnames(orange_dataset[, characteristics_ind]), collapse = " + ")   orange_dataset$all_responses <- model.matrix( as.formula(paste0(" ~ ", our_formula_responses)), data = orange_dataset)   one_response_formula <- as.formula(paste("J1 ~", our_formula_predictors))   our_formula <- as.formula(paste0("all_responses ~ ", our_formula_predictors))   xxx <- plsr(formula = one_response_formula, data = orange_dataset, method = "simpls")   yyy <- plsr(formula = our_formula, data = orange_dataset, method = "kernelpls")   zzz <- pcr(formula = our_formula, data = orange_dataset)

---

De multiples autres options sont disponibles et documentées dans la bibliothèque. Pour traiter spécifiquement du jeu de données de Michel TENENHAUS et alii, la deuxième bibliothèque que vous aviez mentionnée clame être la mieux indiquée. Voir ce qu'en dit son auteur. Malheureusement, comme il n'y a que 6 observations dans le jeu de données, il n'y aura pas de validation croisée des modèles, donc pas de calcul des quantités que vous cherchiez. Voir help("plsreg1") et help("plsreg2). On peut quand même faire quelque chose du genre :

Code:

---

1 2 3

library("plsdepot") vvv <- plsreg2(orange_dataset[, characteristics_ind], orange_dataset[, -characteristics_ind], comps = 4L, crosval = FALSE)

---

Une autre possibilité : les auteurs de la bibliothèque plsRglm disent qu'ils améliorent pls. Voici ce que donne un essai assez rapide:

Code:

---

1 2 3 4 5 6 7

library("plsRglm") library("plsdof") aaa <- plsR(formula = one_response_formula, data = orange_dataset, nt = 4, typeVC = "standard") print(aaa$Q2) print(aaa$R2) print(aaa$Q2cum)

---

Avant de réutiliser our_formula de la même façons que one_response_formula, quelques retouches sont nécessaires afin d'éviter les valeurs manquantes.
Pour explorer d'autres possibilités, voir le point « In views: » sur la page CRAN de pls.

Bonjour juliatheric,

Tout d'abord merci d'avoir pris du temps pour m'aider et pour votre réponse même si elle est tardive. J'ai fait les constatations que les votres, et en effet, je me suis rabatu sur plsRglm qui utilise là même méthodologie que SIMCA-Plus et par conséquent XLSTAT. Les résultats en utilisant la fonction plsR me donne donc les mêmes résultats à présents.

J'ai voulu aller encore plus loin pour cette pls avec plsRglm en essayant de l'utiliser dans les fonctions disponibles sous caret, mais pour le moment mes résultats ne sont pas concluants, car la répétition du job me donne des résultats non répétables même en fixant le set.seed() dans R.

Par conséquent, j'utilise plsR() uniquement en ajoutant les fonctionalités non disponibles de base par moi-même comme par exemple le calcul des VIPs. Par ailleur, ceci peut-faire l'objet d'une autre question, car je retrouve à peu de chose prêt les mêmes valeurs de VIPs en les calculant par moi-même à partir des résultats de plsRglm. Mais, la différence des VIPs observés entre mon modèle et les résultats XLSTAT augmente lorsque que le nombre de composantes étudiées augmente. Ceci peut provenir de l'algorithme plsRglm qui ne donne pas exactement les mêmes valeurs que XLSTAT pour les scores de X-Space et Y-Space, mais je sais aussi qu'un facteur prenant en compte le risque alpha de 5% est injecté dans le caclul des VIPs ... Chose que je n'ai pas pris en compte pour le moment, car je ne sais pas comment ce facteur est intégré par XLSTAT (non renseigné dans la documentation technique).

Bien cordialement,

Baptiste

Bonjour Baptiste,
Sympa de savoir que vous aviez résolu le problème plus tôt.

Citation:

---

Envoyé par baptiste45000

[...] J'ai voulu aller encore plus loin pour cette pls avec plsRglm en essayant de l'utiliser dans les fonctions disponibles sous caret, mais pour le moment mes résultats ne sont pas concluants, car la répétition du job me donne des résultats non répétables même en fixant le set.seed() dans R.

Par conséquent, [...]

---

Je suis curieux de voir comment vous avez géré l'entropie avec les nombres pseudo(quasi)-aléatoires (PQA), i.e. set.seed(). J'ai l'impression qu'il y avait des itérations à faire et que des nombres PQA étaient générés à chaque nouvelle itération, ce qui faussera bien entendu les résultats si les états précédents ne sont pas sauvegardés puis restaurés si de besoin. Ce papier^(*) (et la bibliothèque qui va avec) constituent un bon point d'entrée en la matière. Alors, pourriez-vous poster le code non-reproductible ?

___
(*) Certaines choses dites ont évolué depuis la date de sa rédaction, par exemple les paramètres par défaut de RNGkind() ne sont plus les mêmes qu'à l'époque. Cela dit, la partie scientifique reste valide.

Bonjour juliatheric,

Je suis dans l'incapacité de vous transmettre le script puisque les résultats n'étant pas assez performants pour mon étude, il a été écarté et donc supprimé de mon système.

Par ailleur caret propose une fonction permettant de calculer les "Variable Importance in Projection" globales du modèle et non composantes par composantes (grosse limitation dans le développement d'un modèle QSAR/QSPR). Ajouté au problème de répétabilité observé, j'ai donc développé ma propre méthodologie sur la base de plsRglm. Cette fois encore le script que j'utilise actuellement ne peut pas être transmi, mais il me permet à présent de connaitre toutes les VIPs par composantes de mes variables explicatives et cela sur les résultats en sortie de plsRglm.

Bien cordialement,

Baptiste

Bonjour Baptiste,
Parfait, l'essentiel était que vous dépassiez la difficulté initiale.
Le code fini pourrait intéresser des gens passant sur le forum, l'historique du forum montre qu'à des intervalles un peu longs, des questions sur PLS surgissent.
Bonne continuation.

En tout cas je tiens a vous remercier pour le temps que vous avez pris pour me répondre. Cela m'a énormément aidé pour mieux comprendre les causes de mon problème.

Je serai présent à ce sujet si mon expertise permet de répondre aux futures problématiques posées ;-).

Bien cordialement,

Baptiste