Как выполнить pslda в R, показывая частоту ошибок для нескольких элементов списка?

Question

Я выполняю splsda-модель в R на 10 фреймах данных (данные 10 областей исследования), хранящихся в виде списка (datalist). Все эти фреймы данных похожи, с одинаковыми переменными, но просто разными значениями. Для этого я использую библиотеку micromics.

Это руководитель первого направления обучения. Он сравнивает отсутствие или наличие водно-болотных угодий (факторная переменная - водно-болотные угодья или отсутствие водно-болотных угодий) в зависимости от его значения TPI в различных диапазонах.

> head(datalist[[1]])
  OID POINTID WETLAND TPI200 TPI350 TPI500 TPI700 TPI900 TPI1000 TPI2000 TPI3000 TPI4000 TPI5000 TPI2500
1  -1       1 no wetl     70     67     55     50     48      46      53      47      49      63      48
2  -1       2 no wetl     37     42     35     29     32      16      17      35      49      63      26
3  -1       3 no wetl     45     55     45     39     41      41      53      47      49      63      48
4  -1       4 no wetl     46     58     51     43     46      36      54      47      49      62      49
5  -1       5 no wetl     58     55     53     49     47      46      54      47      49      62      49
6  -1       6 no wetl     56     53     51     49     46      46      54      47      49      61      49

Я выполнил этап перекрестной проверки, используя следующий код:

library(mixOmics)

for (i in 1: length(model_list))
{
  myperf_plsda <- perf(model_list[[i]], validation = "Mfold", folds = 10, 
                            progressBar = FALSE, nrepeat = 10, auc = TRUE)
  save(myperf_plsda, file="performancePLSDA.RData")
}

model_list - это список, полученный из функции spslda.

Но теперь я застрял на следующем шаге, который должен посмотреть на частоту ошибок (в целом и по классам)

Только для одного фрейма данных (область исследования) я могу использовать следующий код :

# cross-validation error in function of nr of PCs
# can see how many PCs is best
plot(myperf_plsda, col = color.mixo(5:7), sd = TRUE, 
     legend.position = "horizontal")

# error rate overall and per class
myperf_plsda$error.rate
myperf_plsda$error.rate.class
myperf_plsda$auc

Итак, первый , я пытаюсь построить график ошибка в функции основных компонентов (= график, первый код здесь выше для одной области исследования). Результат будет примерно таким: I would like to have it in a pdf.

Second, I want to know the overall error rate and error rate per class, from which the code is mentioned above for one study area. The result for one study area is then for example:

• overall error rate:
• error rate per class: частота ошибок на класс

Я пробовал несколько способов для всех этих кодов в for l oop или используя lapply, чтобы получить эти результаты для 10 областей исследования.

, например:

### To see how many PCs is best ###
pdf('overallerrorrate_wetlall_small.pdf')
for (i in 1:length(myperf_plsda))
{
plot(model_list[[i]], col = color.mixo(5:7), sd = TRUE, 
     legend.position = "horizontal")
}
dev.off()

или

for (i in 1:length(myperf_plsda))
  {plot(myperf_plsda, col = color.mixo(5:7), sd = TRUE, 
     legend.position = "horizontal")}

или

for (i in 1:length(myperf_plsda))
{myperf_plsda[[1]]error.rate
  myperf_plsda[[1]]error.rate.class
  myperf_plsda[[i]]auc
}

или

lapply(myperf_plsda, [[, 'error.rate')`

Но все эти коды не работают! Как я могу запустить код для нескольких элементов в списке? Большое спасибо!

Answer 1

На основе ваших результатов вам нужно будет создать новый список и сохранить в нем результаты. Использование только myperf_plsda может привести к перезаписи каждого шага в l oop. Также большинство необходимых мер - это списки, поэтому я добавил некоторые функции обработки для доступа к фреймам данных. Я использовал следующие фиктивные данные:

library(mixOmics)
#Function
custom_splsda <- function(datalist, ncomp, keepX, ..., Xcols, Ycol){
  Y <- datalist[[Ycol]]
  X <- datalist[Xcols]
  res <- splsda(X, Y, ncomp = ncomp, keepX = keepX, ...)
  res
}
#Data
datalist <- list(df1 = structure(list(OID = c(-1, -1, -1, -1, -1, -1), POINTID = c(1, 
2, 3, 4, 5, 6), WETLAND = c("no wetl", "no wetl", "no wetl", 
"wetl", "wetl", "wetl"), TPI200 = c(70, 37, 45, 46, 58, 56), 
    TPI350 = c(67, 42, 55, 58, 55, 53), TPI500 = c(55, 35, 45, 
    51, 53, 51), TPI700 = c(50, 29, 39, 43, 49, 49), TPI900 = c(48, 
    32, 41, 46, 47, 46), TPI1000 = c(46, 16, 41, 36, 46, 46), 
    TPI2000 = c(53, 17, 53, 54, 54, 54), TPI3000 = c(47, 35, 
    47, 47, 47, 47), TPI4000 = c(49, 49, 49, 49, 49, 49), TPI5000 = c(63, 
    63, 63, 62, 62, 61), TPI2500 = c(48, 26, 48, 49, 49, 49)), row.names = c(NA, 
6L), class = "data.frame"), df2 = structure(list(OID = c(-1, 
-1, -1, -1, -1, -1), POINTID = c(1, 2, 3, 4, 5, 6), WETLAND = c("no wetl", 
"no wetl", "no wetl", "wetl", "wetl", "wetl"), TPI200 = c(70, 
37, 45, 46, 58, 56), TPI350 = c(67, 42, 55, 58, 55, 53), TPI500 = c(55, 
35, 45, 51, 53, 51), TPI700 = c(50, 29, 39, 43, 49, 49), TPI900 = c(48, 
32, 41, 46, 47, 46), TPI1000 = c(46, 16, 41, 36, 46, 46), TPI2000 = c(53, 
17, 53, 54, 54, 54), TPI3000 = c(47, 35, 47, 47, 47, 47), TPI4000 = c(49, 
49, 49, 49, 49, 49), TPI5000 = c(63, 63, 63, 62, 62, 61), TPI2500 = c(48, 
26, 48, 49, 49, 49)), row.names = c(NA, 6L), class = "data.frame"))

Теперь код, я создам пустой список myperf_plsda:

#Create model_list, you must have the object created
model_list <- lapply(datalist, custom_splsda,
                     ncomp = 2, keepX = c(5, 5),
                     Xcols = 4:8, Ycol = "WETLAND")
#Create empty list
myperf_plsda <- list()
#Loop for objects and saving
for (i in 1: length(model_list))
{
  myperf_plsda[[i]] <- perf(model_list[[i]], validation = "Mfold", folds = 3, 
                       progressBar = FALSE, nrepeat = 3, auc = TRUE)
  object <- myperf_plsda[[i]]
  save(object,file = paste0("performancePLSDA.",i,".RData"))
}
#Process the object myperf_plsda
#First function to get elements
extract1 <- function(x)
{
  #Object
  error.rate <- x$error.rate
  error.rate <- lapply(error.rate, as.data.frame)
  #Process
  O1 <- do.call(rbind,error.rate)
  #Separate vars
  O1$id <- rownames(O1)
  rownames(O1) <- NULL
  O1$id1 <- gsub("\\..*","", O1$id )
  O1$id2 <- gsub(".*\\.","", O1$id )
  O1$id <- NULL
  return(O1)
}
#Function 2
extract2 <- function(x)
{
  #Object
  error.rate.class <- x$error.rate.class
  names(error.rate.class) <- gsub('.','_',names(error.rate.class),fixed = T)
  error.rate.class <- lapply(error.rate.class, as.data.frame)
  #Process
  O2 <- do.call(rbind,error.rate.class)
  #Separate vars
  O2$id <- rownames(O2)
  rownames(O2) <- NULL
  O2$id1 <- gsub("\\..*","", O2$id )
  O2$id2 <- gsub(".*\\.","", O2$id )
  O2$id <- NULL
  return(O2)
}
#Function 3
extract3 <- function(x)
{
  #Object
  auc <- x$auc
  #Modify for dataframe
  change <- function(x)
  {
    y <- as.data.frame(x)
    y$id1 <- rownames(y)
    rownames(y)<-NULL
    y$id1 <- gsub('.','_',y$id1,fixed = T)
    return(y)
  }
  auc <- lapply(auc, change)
  #Process
  O3 <- do.call(rbind,auc)
  #Separate vars
  O3$id2 <- rownames(O3)
  rownames(O3) <- NULL
  O3$id2 <- gsub("\\..*","", O3$id2 )
  return(O3)
}
#Apply functions and save in lists for late process
L1 <- lapply(myperf_plsda,extract1)
L2 <- lapply(myperf_plsda,extract2)
L3 <- lapply(myperf_plsda,extract3)
#Assign the same names from model_list
names(L1) <- names(model_list)
names(L2) <- names(model_list)
names(L3) <- names(model_list)
#Bind the data
#Error rate
error.rate.df <- do.call(rbind,L1)
error.rate.df$genid <- gsub("\\..*","", rownames(error.rate.df) )
rownames(error.rate.df) <- NULL
#Error rate class
error.rate.class.df <- do.call(rbind,L2)
error.rate.class.df$genid <- gsub("\\..*","", rownames(error.rate.class.df) )
rownames(error.rate.class.df) <- NULL
#Auc
auc.df <- do.call(rbind,L3)
auc.df$genid <- gsub("\\..*","", rownames(auc.df) )
rownames(auc.df) <- NULL

С предыдущим кодом вы получите три фрейма данных, которые содержат значения, которые определены в соответствии с именами model_list, вы можете перемещаться по переменным id1, id2 и genid, чтобы увидеть меры, компоненты и наборы данных:

error.rate.df

   max.dist centroids.dist mahalanobis.dist     id1   id2 genid
1 0.2222222      0.2222222        0.2222222 overall comp1   df1
2 0.2777778      0.3888889        0.2777778 overall comp2   df1
3 0.2222222      0.2222222        0.2222222     BER comp1   df1
4 0.2777778      0.3888889        0.2777778     BER comp2   df1
5 0.2222222      0.2222222        0.2222222 overall comp1   df2
6 0.2777778      0.3333333        0.2777778 overall comp2   df2
7 0.2222222      0.2222222        0.2222222     BER comp1   df2
8 0.2777778      0.3333333        0.2777778     BER comp2   df2

error.rate.class.df

       comp1     comp2              id1     id2 genid
1  0.3333333 0.3333333         max_dist no wetl   df1
2  0.1111111 0.2222222         max_dist    wetl   df1
3  0.3333333 0.6666667   centroids_dist no wetl   df1
4  0.1111111 0.1111111   centroids_dist    wetl   df1
5  0.3333333 0.3333333 mahalanobis_dist no wetl   df1
6  0.1111111 0.2222222 mahalanobis_dist    wetl   df1
7  0.3333333 0.3333333         max_dist no wetl   df2
8  0.1111111 0.2222222         max_dist    wetl   df2
9  0.3333333 0.5555556   centroids_dist no wetl   df2
10 0.1111111 0.1111111   centroids_dist    wetl   df2
11 0.3333333 0.3333333 mahalanobis_dist no wetl   df2
12 0.1111111 0.2222222 mahalanobis_dist    wetl   df2

auc.df

           x      id1   id2 genid
1 0.62966667 AUC_mean comp1   df1
2 0.06414361   AUC_sd comp1   df1
3 0.81483333 AUC_mean comp2   df1
4 0.06414361   AUC_sd comp2   df1
5 0.62966667 AUC_mean comp1   df2
6 0.06414361   AUC_sd comp1   df2
7 0.77780000 AUC_mean comp2   df2
8 0.11110000   AUC_sd comp2   df2

Наконец, для графиков вы можете использовать следующий код (я назначил имя набора данных метке x, чтобы вы могли идентифицировать его на графиках):

#Plot and save
#Assign names
names(myperf_plsda) <- names(model_list)
pdf('example.pdf')
for (i in 1:length(myperf_plsda))
{
  plot(myperf_plsda[[i]], col = color.mixo(5:7), sd = TRUE, 
       legend.position = "horizontal",xlab = paste0(names(myperf_plsda)[i],' (Comp)'))
}
dev.off()

Как примечание, я изменил количество складок, чтобы сделать код работает, но с вашими реальными данными вы можете установить исходные значения, которые у вас есть.