P5 : Classification sur données métabolomiques
Antoine Pissoort
March 28, 2017,22:44
data<-read.csv("CarreLatinMBX.csv", sep=';',header=T)1 Préparation et présentation des données
# Définition des matrices de données et tailles associées
n=dim(data)[1]
m=dim(data)[2]-1
X=apply(data[,1+(1:m)],2,as.numeric)
TypeY=rep(c("A","B","C"),c(27,27,27))
classY=as.numeric(as.factor(TypeY)) # classe 1, 2 et 3 utile pour les graphiques
# On extrait les noms des variables, on enlève le "X" et on les remets comme noms de variables
ppm=as.numeric(substr(dimnames(X)[[2]],2,6))
dimnames(X)[[2]]=ppm
# On mets les noms des spectres comme non des lignes de X
dimnames(X)[[1]]=data[,1]
# Centrage des spectres
XC=scale(X,center=T,scale=F)
# Définition de variables pour les graphiques
Xt <- t(X)
library(reshape)
library(reshape2)
library(tidyr)
Xt <- cbind.data.frame(Xt, spectres = as.factor(1:nrow(Xt)))
Xt <- melt(Xt, id.vars = "spectres" )
Xt <- cbind.data.frame(Xt, ppm = ppm)1.1 Dessins des spectres
Il est possible de sélectionner des spectres (droite) si l’ont veut se concentrer sur certains spectres seulement.
library(plotly)
gg_labs <- labs(title = "Représentation des spectres", x = "Intensité", y = "Patient" )
dimnames(Xt)[[2]][2] <- "Spectre"
gg_spectre <- ggplot(Xt) + geom_line(aes(x = ppm, y = value, col = Spectre), size = 0.25) + gg_labs + scale_color_discrete(guide = F) + scale_color_hue(l = 40, c = 200) + theme_piss()
ggplotly(gg_spectre)
De cette manière “dynamique”, il est plus aisément possible de se faire rapidement une idée des données ainsi que de permettre à vérifier visuellement si Pour ce genre de graphes où il y a be beaucoup de lignes représentées avec des différences minimes, la possibilité de zoomer (e.g. en sélectionnant avec la souris une aire sur le graphe) est également intéressante (…)
2 Analyse en composantes principales des données spectrales
2.1 PCA des spectres centrées avec la librairie MBX
ncomp = 8
PCA.res = SVDforPCA(XC,ncomp = ncomp)2.1.1 Valeurs propres
% de la variance expliquée par chaque variable
eig.res=rbind(PCA.res$var,PCA.res$var*100/sum(PCA.res$var),PCA.res$cumvar)
rownames(eig.res)=c("Variances","Prop Var","Cum Eigen Values")
pandoc.table(eig.res,width = 1000)| PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 | |
|---|---|---|---|---|---|---|---|---|
| Variances | 67.2 | 13.4 | 7 | 4.8 | 2.6 | 2.1 | 0.8 | 0.7 |
| Prop Var | 68.15 | 13.59 | 7.099 | 4.868 | 2.637 | 2.13 | 0.8114 | 0.7099 |
| Cum Eigen Values | 67.2 | 80.6 | 87.6 | 92.4 | 95 | 97.1 | 97.9 | 98.6 |
#kable(eig.res)2.1.2 Scores
Représentation graphique des scores pour les 4 premières composantes
gg_pca1 <- DrawSL(PCA.res, type.obj = "PCA", drawNames=TRUE,createWindow=FALSE, main = paste("PCA score plots for PC1 vs PC2, and for PC3 vs PC4"), axes =c(1,2), type.graph="scores",class=classY)
gg_pca2 <- DrawSL(PCA.res, type.obj = "PCA", drawNames=TRUE,createWindow=FALSE, main = paste("PCA score plots for PC1 vs PC2, and for PC3 vs PC4"), axes =c(3,4), type.graph="scores",class=classY)
subplot(hide_legend(ggplotly(gg_pca1)), ggplotly(gg_pca2), shareX = T, shareY= T, which_layout = 2)On remarque que les patients ont tendance à se regrouper entre eux, ce qui est intuitif. En plus de ça, on remarque aussi une tendance intra-pantients à grouper par 3 en fonction de leur moment de prélevement de l’échantillon, ceci étant plus visible sur les graphes des 2 premières composantes (qui expliquent bien mieux les données !)
2.1.3 Graphe des loadings
gg_load <- DrawSL(PCA.res, type.obj = "PCA", drawNames=F,createWindow=FALSE,
axes = c(1:4), type.graph ="loadings", loadingstype="l",xlab = "ppm",ang=90,xaxis="character",nxaxis=10)
ggplotly(gg_load[[1]])2.1.4 Choix des 2 patients
Sur base de cette ACP, nous décidons de choisir les patients A et B afin de les comparer. Nous justifions ce choix sur base de cet ACP seulement, donc assez arbitrairement, car ils semblent être les plus écartés sur les graphes des 2 premières composantes et donc leur “différences” devraient être plus soulignées et donc l’analyse devrait être plus intéressante..
Xdata <- XC[1:(2*27),] # Subset only our 2 types of patients
Xdata <- cbind.data.frame( Xdata, class = as.numeric( c(rep(0, nrow(Xdata)/2), rep(1, nrow(Xdata)/2) ) ) )3 Ajustement de modèles de classification
Notre but ici est de trouver les marqueurs qui différencient le plus possible les spectres.
3.1 Preliminaries and first features selection
df_mean <- t(data.frame(sort(apply(Xdata, 2, mean))))
df_sd <- t(data.frame(sort(apply(Xdata, 2, sd))))
# names(df_mean) <-
# df <- rbind.data.frame(df_mean, df_sd)
rownames(df_mean) <- "mean"
rownames(df_sd) <- "variance"
pander(df_mean[,1:6]) ; pander(df_sd[,1:6])| 3.259 | 0.859 | 3.739 | 1.259 | 3.499 | 3.859 |
|---|---|---|---|---|---|
| -0.001821 | -0.001567 | -0.0008222 | -0.0007686 | -0.0007265 | -0.0003518 |
| 4.539 | 4.579 | 4.619 | 4.659 | 4.699 | 4.739 |
|---|---|---|---|---|---|
| 0 | 0 | 0 | 0 | 0 | 0 |
Après avoir retiré ces variables inutiles, i.e. qui ne varient pas, on va alors regarder d’abord les corrélations linéaires des variables entre elles, qui s’avèrent d’être toutes assez fortes. De plus, on remarque grâce à ça, que l’odre (i.e. nom) des variables n’est pas anodin, et est lié en fait à la corrélation,
## Remove these features
Xdata <- Xdata[ ,colnames(Xdata) %in% names(sort(apply(Xdata, 2, sd))[1:14]) == F]
corrX <- c()
for (i in 0:(ncol(Xdata)-2)){
corrX[i] <- cor(Xdata[,i+1], Xdata[,i+2])
names(corrX)[i] <- paste0(names(Xdata)[i+1],"&", names(Xdata)[i+2])
}
#sort(corrX)
dfcor1 <- t(as.matrix(sort(abs(corrX))[(length(corrX)-4):length(corrX)]))
rownames(dfcor1) <- "Largest corr"
# library(xtable)
# print(xtable(dfcor1), type = "html")
pander(dfcor1)| 0.739&0.779 | 4.459&4.499 | 1.019&1.059 | 4.379&4.419 | |
|---|---|---|---|---|
| Largest corr | 0.9757 | 0.9804 | 0.9877 | 0.9888 |
| 4.419&4.459 | |
|---|---|
| Largest corr | 0.9899 |
dfcor2 <- t(as.data.frame(sort(abs(corrX))[1:4]))
rownames(dfcor2) <- "Lowest corr"
pander(dfcor2)| 4.499&5.099 | 0.939&0.979 | 1.499&1.539 | 0.899&0.939 | |
|---|---|---|---|---|
| Lowest corr | 0.001844 | 0.001911 | 0.01135 | 0.02956 |
library(mRMRe)library(corrplot)
M <- cor(Xdata)
# #corrplot(unname(M), method="color")
#
library(d3heatmap)
# d3heatmap(M, scale = "col", dendrogram = "none",
# #color = scales::col_quantile("Blues", NULL, 5),
# labRow = substr(rownames(M),1,2), labCol = substr(colnames(M),1,2))Il est intéressant aussi surtout de regarder les corrélations linéaires avec la variable d’intérêt
corr_target <- c()
for (i in 1:(ncol(Xdata)-1)){
corr_target[i] <- cor(Xdata[,i], Xdata[,"class"])
names(corr_target)[i] <- names(Xdata)[i]
}
# sort(corr_target)Ou alors de regarder des indicateurs d’importance des variables par rapport à la variable d’intérêt. Ces indicateurs se basent sur l’entropie et sont non-linéaires.
#library(FSelector)
#
# linear.correlation(as.simple.formula(colnames(Xdata),"Heating.Load"), data_X1)
# rank.correlation(as.simple.formula(colnames(Xdata[,]),"Heating.Load"), data_X1)
library(FSelectorRcpp)
info_gain <- information_gain(x = Xdata[,-ncol(Xdata)], y = as.factor(Xdata$class) , type = c("infogain", "gainratio", "symuncert"), threads = 1)
# info_gain[order(info_gain$importance),]## remove features based on info gain (entropy) --> nonlinear AND on linear correlations
# Here we decide to keep features witn linear corr larger than 5% and that have amount of importance
kept_features <- names(Xdata[, -ncol(Xdata)][(info_gain$importance > 4e-1) & (abs(corr_target) > 0.15)])# Enable parallelization in examples
# library(doSNOW) # doSNOW has an option for progress bar
# cl <- makeCluster(2)
# registerDoSNOW(cl)
#
# evaluator_BIC_glm <- function(attributes, data, dependent = "class") {
# extractAIC(
# fit = glm(to_formula(attributes, dependent), family = binomial,
# data = data),
# k = log(nrow(data))
# )[2]
# }
# feature_search(attributes = names(Xdata[, -ncol(Xdata)]),
# fun = evaluator_BIC_glm,
# data = Xdata, parallel = T,
# mode = "greedy")
#
# stopCluster(cl)
library(mlbench)
library(caret)
# prepare training scheme
# control <- trainControl(method="repeatedcv", number=5, repeats=3)
# # train the model
# model <- train(class~., data=Xdata, method="logreg", preProcess="scale", trControl=control)
# # estimate variable importance
# importance <- varImp(model, scale=FALSE)
# # summarize importance
# print(importance)
# # plot importance
# plot(importance)
# define the control using a random forest selection function
control <- rfeControl(functions=rfFuncs, method="cv", number=5)
# run the RFE algorithm
results <- rfe(Xdata[,-(ncol(Xdata))], Xdata$class, sizes= c(1:5, seq(10, 35, by = 5)), rfeControl=control, metric = "accuracy")
# summarize the results
print(results)
Recursive feature selection
Outer resampling method: Cross-Validated (5 fold)
Resampling performance over subset size:
Variables RMSE Rsquared RMSESD RsquaredSD Selected
1 0.0002971 1.0000 0.0004141 8.508e-07 *
2 0.0498470 0.9772 0.0693183 3.329e-02
3 0.0605808 0.9789 0.0602193 2.421e-02
4 0.0640011 0.9660 0.0892768 6.629e-02
5 0.0572107 0.9767 0.0755158 4.568e-02
10 0.0695924 0.9744 0.0711915 4.502e-02
15 0.0740498 0.9738 0.0705607 4.490e-02
20 0.0753214 0.9764 0.0625696 3.839e-02
25 0.0807647 0.9742 0.0644562 4.194e-02
30 0.0797879 0.9756 0.0647704 3.879e-02
35 0.0794946 0.9750 0.0637038 4.085e-02
224 0.0819109 0.9771 0.0574587 3.212e-02
The top 1 variables (out of 1):
1.059# list the chosen features
print(kept_features2 <- predictors(results))[1] "1.059"# plot the results
plot(results, type=c("g", "o"))
#devtools::install_github('mlampros/FeatureSelection')
library(FeatureSelection)
#source("Wrapper_featureselection_package_fun.R")
params_glmnet = list(alpha = 1, family = 'gaussian', nfolds = 5, parallel = TRUE)
params_xgboost = list( params = list("objective" = "reg:linear", "bst:eta" = 0.001, "subsample" = 0.75, "max_depth" = 5, "colsample_bytree" = 0.75, "nthread" = 6),
nrounds = 1000, print.every.n = 1000, maximize = FALSE)
# params_ranger1 = list(probability = FALSE, num.trees = 1000, verbose = TRUE, mtry = 5, min.node.size = 10, num.threads = 6, classifiction = FALSE, importance = 'permutation')
params_features = list(keep_number_feat = 15, union = TRUE) # say we want 15 features max
#Xdata$class <- as.factor(Xdata$class)
feat = wrapper_feat_select(X = Xdata[,-(ncol(Xdata))], y = Xdata$class, params_glmnet = params_glmnet, params_xgboost = params_xgboost, xgb_sort = 'Gain', CV_folds = 5, stratified_regr = FALSE, scale_coefs_glmnet = FALSE, cores_glmnet = 8, params_features = params_features, verbose = F)[1] train-rmse:0.358855
[1000] train-rmse:0.000280
[1] train-rmse:0.357917
[1000] train-rmse:0.000121
[1] train-rmse:0.357739
[1000] train-rmse:0.000266
[1] train-rmse:0.358949
[1000] train-rmse:0.000166
[1] train-rmse:0.358359
[1000] train-rmse:0.000225 #feat$all_feat
feat$union_feat feature importance Frequency
12 3.379 1.000000000 2
21 7.219 0.940823561 2
2 0.899 0.916825743 2
1 0.859 0.914098718 2
5 2.259 0.752113444 2
7 2.739 0.661576220 2
16 4.139 0.541314426 2
33 0.979 0.501772566 1
35 1.059 0.495500409 1
17 4.259 0.490591764 2
61 3.699 0.482956095 1
27 0.579 0.476683938 1
31 0.779 0.470411781 1
3 0.939 0.469593673 1
37 1.139 0.464139624 1
10 3.259 0.458140169 2
38 1.179 0.457867467 1
30 0.699 0.451595310 1
32 0.819 0.445323152 1
60 3.419 0.432778838 1
28 0.619 0.420234524 1
25 0.499 0.413962367 1
29 0.659 0.407690210 1
59 3.299 0.401418053 1
26 0.539 0.395145896 1
74 6.339 0.388873739 1
43 1.419 0.382601582 1
6 2.419 0.380147259 1
41 1.339 0.376329425 1
34 1.019 0.370057268 1
39 1.259 0.363785110 1
53 2.339 0.357512953 1
73 6.259 0.351240796 1
48 1.899 0.344968639 1
71 6.019 0.338696482 1
51 2.139 0.332424325 1
83 7.459 0.319880011 1
58 3.219 0.313607854 1
13 3.459 0.309790019 2
45 1.499 0.307335697 1
62 3.819 0.301063540 1
50 2.019 0.294791383 1
15 4.059 0.290700845 1
78 7.099 0.288519226 1
44 1.459 0.282247068 1
52 2.299 0.275974911 1
18 5.779 0.270520862 2
36 1.099 0.269702754 1
19 5.819 0.268339242 1
94 9.979 0.263430597 1
68 5.459 0.257158440 1
40 1.299 0.250886283 1
49 1.979 0.238341969 1
63 3.979 0.232069812 1
14 3.499 0.223616035 1
90 9.219 0.219525498 1
11 3.339 0.217616580 2
42 1.379 0.213253341 1
85 7.899 0.206981184 1
8 2.939 0.201254431 1
66 4.459 0.200709026 1
69 5.699 0.194436869 1
87 8.539 0.188164712 1
55 2.539 0.181892555 1
9 2.979 0.178892828 2
88 8.659 0.169348241 1
79 7.139 0.163076084 1
67 5.219 0.156803927 1
4 2.059 0.156531224 1
56 2.579 0.144259613 1
47 1.659 0.137987456 1
91 9.299 0.131715299 1
46 1.579 0.125443142 1
77 7.019 0.119170984 1
70 5.979 0.112898827 1
20 5.859 0.111808017 1
80 7.179 0.106626670 1
22 7.499 0.089446414 1
92 9.739 0.087810199 1
82 7.419 0.081538042 1
75 6.619 0.075265885 1
89 8.979 0.068993728 1
93 9.819 0.062721571 1
57 2.779 0.056449414 1
81 7.299 0.050177257 1
54 2.379 0.043905100 1
84 7.659 0.037632942 1
72 6.139 0.031360785 1
76 6.939 0.025088628 1
24 9.939 0.022361603 1
86 8.019 0.018816471 1
64 4.219 0.012544314 1
65 4.419 0.006272157 1
23 9.579 0.000000000 1Xdata[, -ncol(Xdata)] <- as.data.frame(scale(Xdata[, -ncol(Xdata)])) # Scale ? scale les variables (qui sont sur une très petite échelle à la base) ?
Dénition des critères de classification
modnames <- c('Logistique-SW','PLS-DA','OPLS-DA','Logistique-Pen')
# Preparation des matrices de résultats
RMSE <- matrix(nrow = 4,ncol=3)
colnames(RMSE) = c('RMSE','RMSE-CV','Taille du modèle') ; rownames(RMSE) = modnames
coefficient = matrix(0, ncol = ncol(Xdata)-1, nrow=4)
colnames(coefficient) = colnames(XC)[1:(ncol(Xdata)-1)] ; rownames(coefficient) =modnames
PredMat=matrix(nrow = nrow(Xdata),ncol= 4) ; colnames(PredMat)=modnames
listr=list(coefficient=coefficient,RMSE=RMSE,PredMat=PredMat)
# fonction pour stocker et imprimer les résultats de chaque modèle
'printresmod' <- function(i, YC, yfit, yfitCV, RMSECV, coef, parmod, listr){
g1 <- ggplot(data.frame(YC = YC, yfit = yfit)) + geom_point(aes(x = YC, y = yfit), size = .5) + theme_piss()
g2 <- ggplot(data.frame(YC = YC, yfitCV = yfitCV)) + geom_point(aes(x = YC, y = yfitCV), size = .5) + theme_piss()
subplot(ggplotly(g1), ggplotly(g2), shareX = T, shareY= F, titleY=T) %>% layout(showlegend = F) %>% print()
# calcul des RMSE
n=length(YC)
listr$RMSE[i,1] <- sqrt(sum((YC-yfit)^2)/n) ; listr$RMSE[i,2] = RMSECV ; listr$RMSE[i,3] = parmod
# Récupération des coefficients et des réponses
listr$coefficient[i,]=coef ; listr$PredMat[,i]=yfit
# affiichage de résultats
pander(listr$RMSE[i,])
par(mfrow=c(1,1))
dotchart(listr$coefficient[i,],labels=colnames(listr$coefficient),xlab="Coefficients",main=paste("Coefficients ",dimnames(listr$coefficient)[[1]][i]))
return(listr)
}
# Definition d'une fonction qui fait de la k-fold CV sur pour une fonction de régression prédéfinie fctreg
'kfoldCV' <- function(X ,Y, k, argrm = 0, seed = 100){
n=dim(X)[1] ; set.seed=seed ; nbseg = k ; sn=ceiling(n/k)
seg_i <- sample(rep(1:k, length.out = n)) ; yfit.cv = numeric()
for (i in 1:nbseg){
# browser()
val_i <- which(seg_i == i)
Y.train = Y[-val_i] ; X.train =X[-val_i,] ; X.valid=X[val_i,]
#browser()
yfitcv <- fctreg(X.train, Y.train, X.valid, argrm = argrm) ; yfit.cv[val_i] = yfitcv
}
RMSE.cv<- sqrt( (1/n)*sum((Y-yfit.cv)^2) )
return(list(yfit.cv=yfit.cv,RMSE.cv=RMSE.cv))
}3.2 Modèle de régression Logistique Stepwise
Comme premier modèle, nous considérons un modèle logistique de type stepwise forward. Le problème avec ce genre de modèle de régression est lorsque le nombre de de variables est grand par rapport au nombre de points (typiquement \(p>n\)) est qu’on ne saura pas estimer de paramètres, i.e. nous ne bénéficions pas d’assez de points par rapport au nombre de degrés de libertés.
On va donc ici appliquer notre feature selection préliminaire.
Xdata_reg <- Xdata[, kept_features]
Xdata_reg <- cbind.data.frame(Xdata_reg, class = Xdata$class)On remarque ici directement, que ce soit selon le critère AIC ou BIC, c’est le modèle avec comme seule variable 0.899 qui est choisi. Celui-ci permet de prédire parfaitement (sur le training set !) à quel patient appartient l’échantillon de sang.
On va alors plutôt On va plutôt faire du 5-fold car la taille de l’échantillon est déjà très petite (5). L’erreur de validation serait alors trop peu précise car basée sur trop peu de données.
yfit = fitted.values(fit.sw)
coef=coefficients(fit.sw)
# cross validation k-fold pour vérifier le pouvoir prédictif du modèle
'fctreg' <- function(X, Y, XV, argrm){
fit.sw1 <- step(null, scope = list(lower=null, upper=full), direction="both", trace = 0)
coefglm <- coef(fit.sw1)
# browser()
logistyfit <- cbind(1, XV[, colnames(XV) %in% names(fit.sw1$model)[2] ]) %*% as.matrix(coefglm)
exp(logistyfit)/(1+exp(logistyfit))
}
res5cv <- kfoldCV(Xdata_reg[,-ncol(Xdata_reg)], Xdata$class, k = 5)
# Sauvegardes et impression de résultats
listr <- printresmod(1, Xdata$class, yfit, res5cv$yfit.cv, res5cv$RMSE.cv, coef = coef, length(names(coef)), listr)
tc <- trainControl("cv", 5, savePredictions = T) #"cv" = cross-validation, 10-fold
fit <- train(class ~. ,
data = Xdata ,
method = "glm" ,
family = binomial ,
trControl = tc)
# fit$results
# fit$modelType
# fit$modelInfo
# fit$finalModel
#plot(fit)On doit choisir un seuil. On pourrait prendre intuitivement 0.5, mais ce n’est probablement pas optimale.
3.3 PLS-DA
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
'fctreg' <- function(X, Y, XV, argrm){
coef <- coef(pls::plsr(Y ~ ., data = X, ncomp = ncomp))
# browser()
as.matrix(XV) %*% as.matrix(as.numeric(coef[,,]))
}
## Boucle sur le nombre de composantes possibles
RMSECV <- numeric()
yfitCV <- list() ; NCmax <- 20 # tune it
for(i in 1:NCmax){
rescv <- kfoldCV(Xdata[,-ncol(Xdata)], Xdata$class, k = 5, argrm = i)
RMSECV[i] <- rescv$RMSE.cv
yfitCV[[i]] <- rescv$yfit.cv
}
# Recherche du nombre de composantes à garder
df_plsda <- data.frame(N_comp = 1:NCmax, RMSEcv = RMSECV)
ncopt <- which.min(RMSECV)
gg <- ggplot(df_plsda, aes(x = N_comp, y = RMSEcv)) + geom_point() + geom_line() + geom_vline(xintercept = ncopt, col = 2, linetype = 2) +theme_piss(size_c = 12) + labs(title = "RMSE en cross-validation vs complexité", x = "Nombre de composantes", y = "RMSE en cv")
ggplotly(gg)# Ajustement du modèle avec le nombre optimal de composantes
fit.plsr <- plsr(class ~., data = Xdata, ncomp=max(ncopt,2), validation = "none")
print(summary(fit.plsr))Data: X dimension: 54 224 Y dimension: 54 1 Fit method: kernelpls Number of components considered: 7 TRAINING: % variance explained 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps X 28.66 50.51 67.41 75.63 83.99 85.86 88.58 class 73.16 89.09 94.28 98.58 99.36 99.78 99.84 NULL
yfit <- fit.plsr$fitted.values[,,ncopt]
# Sauvegardes et impression de résultats
MY <- mean(Xdata$class)
YC <- scale(Xdata$class, center=T, scale=F)
listr <- printresmod(2, YC+MY, yfit+MY, yfitCV[[ncopt]]+MY, RMSECV[ncopt], coef(fit.plsr), ncopt, listr)
Représentation graphique des scores pour les 4 premières composantes
On refait tout d’abord la PCA avec seulement les individus des 2 groupes choisis.
PCA.res = SVDforPCA(Xdata,ncomp = ncomp)
### Scores
par(mfrow=c(1,1))
PCA.res$scores[,1:2] <- fit.plsr$scores[,1:2]
PCA.res$eigval=rep(0,2)
DrawSL(PCA.res, type.obj = "PCA", drawNames=TRUE,createWindow=FALSE, main = paste("PLS score plot for PC1 and PC2"), axes =c(1,2), type.graph = "scores", class = Xdata$class)
### Graphe des loadings
PCA.res$loadings=fit.plsr$loadings
DrawSL(PCA.res, type.obj = "PCA", drawNames=F,createWindow=FALSE,
axes = c(1:2), type.graph ="loadings", loadingstype="p",ang=90,xaxis="character",nxaxis=10)[[1]]
3.4 OPLS-DA
# Recherche du nombre de composantes optimales par CV
## définition de la fonction d'estimation du modèle
fctreg=function(X,Y,XV,argrm) { as.matrix(XV) %*% OPLSDA(X, Y, no = argrm-1)$b }
## Boucle sur le nombre de composantes possibles
RMSECV=numeric()
yfitCV=list()
indi=2:20
for(i in 1:length(indi)){
rescv <- kfoldCV(as.matrix(Xdata[,-ncol(Xdata)]), Xdata$class, k = 5, argrm = indi[i])
RMSECV[i] <- rescv$RMSE.cv
yfitCV[[i]] <- rescv$yfit.cv
}
# Recherche du nombre de composantes à garder
inc=which.min(RMSECV) ; ncopt=indi[inc]
df_oplsda <- data.frame(indi = indi, RMSEcv = RMSECV)
ncopt <- which.min(RMSECV)
gg <- ggplot(df_oplsda, aes(x = indi, y = RMSEcv)) + geom_point() + geom_line() + geom_vline(xintercept = ncopt, col = 2, linetype = 2) + theme_piss(size_c = 12) + labs(title = "RMSE en cross-validation vs complexité", x = "Nombre de composantes", y = "RMSE en cv")
ggplotly(gg)# Ajustement du modèle avec le nombre optimal de composantes
fit.opls <- OPLSDA(as.matrix(Xdata[,-ncol(Xdata)]), Xdata$class, no = ncopt-1)
coefopls=fit.opls$b
yfit <- as.matrix(Xdata[,-ncol(Xdata)]) %*% coefopls
# Sauvegardes et impression de résultats
MY <- mean(Xdata$class)
listr=printresmod(4,Xdata$class+MY,yfit+MY,yfitCV[[ncopt]]+MY,RMSECV[ncopt],coefopls,ncopt,listr)
### Scores
par(mfrow=c(1,1))
DrawSL(fit.opls, type.obj = "OPLSDA", drawNames=TRUE,createWindow=FALSE, main = paste("OPLS score plot for PC and PC1-Ortho"), axes =c(1,2), type.graph="scores",class=as.factor(Xdata$class))
### Graphe des loadings
#PCA.res$loadings=fit.plsr$loadings
DrawSL(fit.opls, type.obj = "OPLSDA", drawNames=F,createWindow=FALSE,
axes = c(1:2), type.graph ="loadings", loadingstype="l",xlab= "ppm??",ang=90,xaxis="character",nxaxis=10)[[1]]
3.5 Régression Logistique pénalisée
Comme annoncé lors du cours, on laisse tomber le Ridge; le LASSO fera fera, basé sur le même “principe”, fera l’affaire.
Néanmoins, on regardera aussi à
3.5.1 LASSO
library(glmnet)
fit <- glmnet(as.matrix(Xdata[,-ncol(Xdata)]), as.factor(Xdata$class), nlambda = 50, alpha = 1, family = "binomial")
# alpha = 1 for lasso penalty !
plot(fit ,xvar="lambda", main="") # Plot the paths for the fit
fit # look at the fit for each coefficient
Call: glmnet(x = as.matrix(Xdata[, -ncol(Xdata)]), y = as.factor(Xdata$class), family = "binomial", alpha = 1, nlambda = 50)
Df %Dev Lambda
[1,] 0 9.610e-16 0.456900
[2,] 1 1.036e-01 0.415900
[3,] 1 1.915e-01 0.378600
[4,] 1 2.674e-01 0.344600
[5,] 2 3.345e-01 0.313700
[6,] 3 3.971e-01 0.285600
[7,] 4 4.593e-01 0.260000
[8,] 4 5.154e-01 0.236600
[9,] 4 5.647e-01 0.215400
[10,] 5 6.083e-01 0.196100
[11,] 5 6.472e-01 0.178500
[12,] 5 6.816e-01 0.162500
[13,] 5 7.124e-01 0.147900
[14,] 5 7.399e-01 0.134600
[15,] 5 7.645e-01 0.122600
[16,] 5 7.866e-01 0.111600
[17,] 5 8.065e-01 0.101600
[18,] 5 8.244e-01 0.092460
[19,] 5 8.406e-01 0.084160
[20,] 5 8.552e-01 0.076610
[21,] 5 8.684e-01 0.069740
[22,] 5 8.804e-01 0.063490
[23,] 5 8.912e-01 0.057790
[24,] 5 9.010e-01 0.052610
[25,] 5 9.099e-01 0.047890
[26,] 6 9.180e-01 0.043590
[27,] 6 9.254e-01 0.039680
[28,] 6 9.321e-01 0.036120
[29,] 6 9.382e-01 0.032880
[30,] 6 9.437e-01 0.029930
[31,] 7 9.488e-01 0.027250
[32,] 7 9.534e-01 0.024800
[33,] 7 9.576e-01 0.022580
[34,] 7 9.614e-01 0.020550
[35,] 7 9.649e-01 0.018710
[36,] 7 9.681e-01 0.017030
[37,] 7 9.709e-01 0.015500
[38,] 7 9.735e-01 0.014110
[39,] 7 9.759e-01 0.012850
[40,] 8 9.781e-01 0.011690
[41,] 8 9.801e-01 0.010650
[42,] 9 9.819e-01 0.009691
[43,] 9 9.835e-01 0.008821
[44,] 9 9.850e-01 0.008030
[45,] 9 9.863e-01 0.007310
[46,] 9 9.876e-01 0.006654
[47,] 9 9.887e-01 0.006057
[48,] 9 9.897e-01 0.005514
[49,] 9 9.906e-01 0.005019
[50,] 9 9.915e-01 0.004569require(doMC)
registerDoMC(cores = 8) # Go parallel with the method : gain computation time
cv.fit <- cv.glmnet(as.matrix(Xdata[,-ncol(Xdata)]), parallel = T, Xdata$class, nfolds = 5) # Perform cross validation on the fited model
plot(cv.fit) # Plot the mean sq error for the cross validated fit as a function
tpred <- predict(fit, as.matrix(Xdata[,-ncol(Xdata)])) # Predictions on the test data
mse <- apply((tpred - Xdata$class)^2, 2, mean) # Compute mse for the predictions
plot(log(fit$lambda), mse, col="blue", pch="*") # overlay the mse predictions on the plot
legend("topleft",legend=c("10 fold CV","Test"),pch="*",col=c("red","blue"))
fit$df [1] 0 1 1 1 2 3 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 7 7 7
[39] 7 8 8 9 9 9 9 9 9 9 9 9cv.fit$lambda.min[1] 0.005014443.5.2 Elastic Nets
L’avantage de cette méthode comparée au LASSO (voir (Zou and Hastie 2005) ) est :
- Lorsque \(p>n\), comme c’est bien ici le cas (239>54) le LASSO ne sélectionera qu’au plus n variables, ce qui n’est pas forcément optimal.
- Les variables qui apparaissent par “groupes” sont sélectionnées par groupes, et n’e pas ignorer n’ignorent pas les autres variables d’un même groupe.
\[\text{argmin}_{\beta} \Big\{ \sum_{i=1}^n(y_i-x_i'\boldsymbol{\beta})^2 +\lambda_1\sum_{i=1}^n\beta_i^2 + \lambda_2\sum_{j=1}^n|\beta_i|\Big\} \] On voit directement que les 2 termes de régularisation apparaissent ici dans la fonction objective à minimiser.
fit <- glmnet(as.matrix(Xdata[,-ncol(Xdata)]), as.factor(Xdata$class), nlambda = 50, alpha = .5, family = "binomial") # Recherche le paramètre de lissage qui min RMSE
## définition de la fonction d'estimation du modèle
#lambda=0.0001*10^(seq(0.5,2.5,0.5))
lambda <- 1e-4 * 10^(seq(1,2,0.5))
'fctreg' <- function(X, Y, XV, argrm){
fit.lasso <- penalized::penalized(response = Y, penalized = ~ -1 + X, lambda1 = argrm, lambda2 = 0 , trace=F)
attr(fit.lasso,'penalized')
}
fit.lasso <- penalized::penalized(response = Xdata$class, penalized = ~ as.matrix(-1 + Xdata[,-ncol(Xdata)]), lambda1 = lambda[1], lambda2 = 0 , trace=F)
## Boucle sur le nombre de composantes possibles
RMSECV=numeric()
yfitCV=list()
for(i in 1:length(lambda)){
cat("\n ", i)
rescv <- kfoldCV(Xdata[,-ncol(Xdata)], Xdata$class, k = 5, argrm = lambda[1])
RMSECV[i] <- rescv$RMSE.cv
yfitCV[[i]] <- rescv$yfit.cv
}
# Recherche du nombre de composantes à garder
par(mfrow=c(1,1))
plot(lambda,RMSECV,type="o",main="RMSE cross-validation",log="x")
il=which.min(RMSECV)
llassoopt=lambda[il]
abline(v=llassoopt,col=2,lty=2)
# Ajustement du modèle avec le nombre optimal de composantes
fit.lasso <- penalized(response=YC,penalized=~-1+XC,lambda1=llassoopt, lambda2=0 ,trace=F )
coeflasso=attr(fit.lasso,'penalized')
yfit=XC%*%coeflasso
# Sauvegardes et impression de résultats
listr=printresmod(7,YC,yfit,yfitCV[[il]],RMSECV[il],coeflasso,lambda[il],listr)4 Comparaison des performances des modèles
5 Identification des variables ou pics de spectre discriminants
6 Remarques
On voit clairement que les
References
Zou, H., and T. Hastie. 2005. “Regularization and Variable Selection via the Elastic Net.” Journal of the Royal Statistical Society Series B 67. http://dx.doi.org/10.1111/j.1467-9868.2005.00503.x.