Tune sPLSDA: Error in solve.default(t(Pmat) %*% Wmat)

[Kim-Anh Le Cao]

See email from Carine

> rm(list=ls())
> workDir <- "C:/Users/cagenet/Documents/A/AMHAROC/10_Human_data/Analyse_R"
> setwd(workDir)
> library(RColorBrewer)
> library(mixOmics)
> q<-read.csv("qpcr_med_ano.csv", sep=";", header=TRUE)
> rownames(q)<-q[,1]
> E<-read.csv("E_med_ano.csv", sep=";", header = TRUE, stringsAsFactors = FALSE)
> rownames(E)<-E$sample
> E$treatment<-(c(rep("CTRL",14),rep("PCOS",16)))
> X<-q[,-c(1,56)]
> Y<-E[,-c(1,10)]
> cbind(rownames(X),rownames(Y))
      [,1]   [,2] 
 [1,] "C105" "C105"
[2,] "C110" "C110"
[3,] "C111" "C111"
[4,] "C112" "C112"
[5,] "C115" "C115"
[6,] "C121" "C121"
[7,] "C27"  "C27"
 [8,] "C28"  "C28"
 [9,] "C52"  "C52"
[10,] "C53"  "C53"
[11,] "C63"  "C63"
[12,] "C84"  "C84"
[13,] "C85"  "C85"
[14,] "C89"  "C89"
[15,] "P107" "P107"
[16,] "P11"  "P11"
[17,] "P114" "P114"
[18,] "P116" "P116"
[19,] "P118" "P118"
[20,] "P13"  "P13"
[21,] "P22"  "P22"
[22,] "P30"  "P30"
[23,] "P33"  "P33"
[24,] "P34"  "P34"
[25,] "P70"  "P70"
[26,] "P74"  "P74"
[27,] "P83"  "P83"
[28,] "P96"  "P96"
[29,] "P98"  "P98"
[30,] "P99"  "P99"
> result<-rcc(X,Y, ncomp = 3, method = 'shrinkage' )
> plot(result, scree.type = "barplot")
> res.plsda<-plsda(X, as.factor(q$Treatment), ncomp=7)
> perf.plsda <- perf(res.plsda, validation = "Mfold", folds = 4,
+                           progressBar = TRUE, auc = TRUE, nrepeat = 100)

comp 1
  |===================================================================================| 100%
comp 2
  |===================================================================================| 100%
comp 3
  |===================================================================================| 100%
comp 4
  |===================================================================================| 100%
comp 5
  |===================================================================================| 100%
comp 6
  |===================================================================================| 100%
comp 7
  |===================================================================================| 100%
> perf.plsda$error.rate
$overall
        max.dist centroids.dist mahalanobis.dist
comp 1 0.3110000      0.3080000        0.3080000
comp 2 0.3243333      0.3193333        0.3236667
comp 3 0.3160000      0.3070000        0.3176667
comp 4 0.3406667      0.3176667        0.3400000
comp 5 0.3260000      0.3093333        0.3263333
comp 6 0.3233333      0.2933333        0.3226667
comp 7 0.3460000      0.3113333        0.3460000

$BER
        max.dist centroids.dist mahalanobis.dist
comp 1 0.3156250      0.3115179        0.3115179
comp 2 0.3260714      0.3203125        0.3247321
comp 3 0.3196875      0.3098214        0.3211161
comp 4 0.3426339      0.3196429        0.3418750
comp 5 0.3269196      0.3113393        0.3272321
comp 6 0.3246875      0.2955804        0.3239732
comp 7 0.3470982      0.3139286        0.3470982

> plot(perf.plsda, overlay='measure', sd=TRUE)
> plot(perf.plsda, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal")
> list.keepX <- c(1:10,  seq(10, 40, 10))
> res.plsda<-plsda(X, as.factor(q$Treatment), ncomp=2)
> plotIndiv(res.plsda , comp = 1:2,
+            group = X4$treatment, ind.names = FALSE,
+            ellipse = TRUE, legend = TRUE, title = 'PLSDA on QPCR HUMAN')
Error in is.factor(group) : object 'X4' not found
> plotIndiv(res.plsda , comp = 1:2,
+            group = X$treatment, ind.names = FALSE,
+            ellipse = TRUE, legend = TRUE, title = 'PLSDA on QPCR HUMAN')
Error in `[<-`(`*tmp*`, classification == groups[j], j, value = 1) :
  subscript out of bounds
> plotIndiv(res.plsda , comp = 1:2,
+            group = q$Treatment, ind.names = FALSE,
+            ellipse = TRUE, legend = TRUE, title = 'PLSDA on QPCR HUMAN')
> tune.splsda <- tune.splsda(X, as.factor(q$Treatment), ncomp = 1, validation = 'Mfold', folds = 10, dist = 'centroids.dist', test.keepX = list.keepX, nrepeat=50)

comp 1
  |                                                                                   |   0%Error in solve.default(t(Pmat) %*% Wmat) :
  Lapack routine dgesv: system is exactly singular: U[2,2] = 0
> tune.splsda <- tune.splsda(X, as.factor(q$Treatment), ncomp = 1, validation = 'loo', folds = 10, dist = 'centroids.dist', test.keepX = list.keepX)

comp 1
  |                                                                                   |   0%Error in solve.default(t(Pmat) %*% Wmat) :
  system is computationally singular: reciprocal condition number = 3.18528e-17
> tune.splsda <- tune.splsda(X, as.factor(q$Treatment), ncomp = 1, validation = 'loo', folds = 5, dist = 'centroids.dist', test.keepX = list.keepX)

comp 1
  |                                                                                   |   0%Error in solve.default(t(Pmat) %*% Wmat) :
  system is computationally singular: reciprocal condition number = 3.18528e-17
> tune.splsda <- tune.splsda(X, as.factor(q$Treatment), ncomp = 1, validation = 'Mfold', folds = 5, dist = 'centroids.dist', test.keepX = list.keepX, nrepeat=50)

comp 1
  |                                                                                   |   0%Error in solve.default(t(Pmat) %*% Wmat) :
  Lapack routine dgesv: system is exactly singular: U[2,2] = 0
> tune.splsda <- tune.splsda(X, as.factor(q$Treatment), ncomp = 1, validation = 'loo', folds = 3, dist = 'centroids.dist', test.keepX = list.keepX)

comp 1
  |                                                                                   |   0%Error in solve.default(t(Pmat) %*% Wmat) :
  system is computationally singular: reciprocal condition number = 3.18528e-17
> tune.splsda <- tune.splsda(X, as.factor(q$Treatment), ncomp = 2, validation = 'loo', folds = 3, dist = 'centroids.dist', test.keepX = list.keepX)

comp 1
  |                                                                                   |   0%Error in solve.default(t(Pmat) %*% Wmat) :
  system is computationally singular: reciprocal condition number = 3.18528e-17

[GENET Carine]

Hi all, I tried to use splsda but with another dataset than above. and I also obtain an error.. So I certainly missed something important about MixoMics, right ? Or is it a bug in the latest version ? How should I treated my data before splsda ? For my X dataset here, I used RNAseq norm counts. I filtered my data by discarding genes with variance <27. Should I have to proceed differently ? Thanks in advance... below is my code.

dim(X) = 8 12343 dim(Y)=8 148 dim(dead)=8 19 result.rcc<-rcc(X, dead, ncomp = 7, method = 'shrinkage' )

plot(result.rcc, scree.type = "barplot") par(mfrow=c(1,1)) plotIndiv(result.rcc, comp = 1:2, ind.names = phys2$sample, + group = phys2$treatment, rep.space = "XY-variate", + legend = TRUE, title = 'Mouton, rCCA XY-space') plotIndiv(result.rcc, comp = 1:2, ind.names = phys2$sample, + group = phys2$treatment, legend = TRUE, title = 'Mouton rCCA, Each subspace') par(mfrow=c(1,1)) col.nutri <- color.mixo(as.factor(phys2$treatment)) plotArrow(result.rcc, col = col.nutri, title = 'Mouton, arrow plot') plotVar(result.rcc, comp=1:2, var.names = c(TRUE, TRUE),cex = c(4,4), + title='AMHAROC sheep, rCCA comp 1-2', cutoff=0.5) res.plsda<-plsda(X, as.factor(phys2$treatment), ncomp=7) perf.plsda <- perf(res.plsda, validation = "loo", + progressBar = TRUE, auc = TRUE)

comp 1 |==================================================================================| 100% comp 2 |==================================================================================| 100% comp 3 |==================================================================================| 100% comp 4 |==================================================================================| 100% comp 5 |==================================================================================| 100% comp 6 | | 0%Error in solve.default(Sr) : system is computationally singular: reciprocal condition number = 2.02407e-17

res.plsda<-plsda(X, as.factor(phys2$treatment), ncomp=5) perf.plsda <- perf(res.plsda, validation = "loo", + progressBar = TRUE, auc = TRUE)

comp 1 |==================================================================================| 100% comp 2 |==================================================================================| 100% comp 3 |==================================================================================| 100% comp 4 |==================================================================================| 100% comp 5 |==================================================================================| 100%

perf.plsda$error.rate $overall max.dist centroids.dist mahalanobis.dist comp 1 0.125 0.125 0.125 comp 2 0.000 0.000 0.000 comp 3 0.000 0.000 0.000 comp 4 0.000 0.000 0.000 comp 5 0.000 0.000 0.000

$BER max.dist centroids.dist mahalanobis.dist comp 1 0.125 0.125 0.125 comp 2 0.000 0.000 0.000 comp 3 0.000 0.000 0.000 comp 4 0.000 0.000 0.000 comp 5 0.000 0.000 0.000

plot(perf.plsda, overlay='measure', sd=TRUE) plot(perf.plsda, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal") list.keepX <- c(1:20, seq(20, 50, 10)) #je choisis ncomp=2 res.plsda<-plsda(X, as.factor(phys2$treatment), ncomp=2) plotIndiv(res.plsda , comp = 1:2, + group = phys2$treatment, ind.names = FALSE, + ellipse = TRUE, legend = TRUE, title = 'PLSDA on RNAseq sheep') tune.splsda <- tune.splsda(X, as.factor(phys2$treatment), ncomp = 2, validation = 'loo', dist = 'centroids.dist', test.keepX = list.keepX, nrepeat=100)

comp 1 | | 0%Error in solve.default(t(Pmat) %*% Wmat) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0 In addition: Warning message: In tune.splsda(X, as.factor(phys2$treatment), ncomp = 2, validation = "loo", : Leave-One-Out validation does not need to be repeated: 'nrepeat' is set to '1'.

tune.splsda <- tune.splsda(X, as.factor(phys2$treatment), ncomp = 2, validation = 'loo', dist = 'centroids.dist', test.keepX = list.keepX)

comp 1 | | 0%Error in solve.default(t(Pmat) %*% Wmat) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0

[Florian Rohart]

solved in commit 6e0cbef will be available with the next release.

Reason: two keepX are the same in the list, easy fix from the user side: use only unique keepX