3 Example
In this section, we show the R script necessary to conduct the significant feature selection with “sigFeature” package on the example dataset. First, we load the “sigFeature” package and the example dataset (as mentioned in the data section) into the memory of the machine. After that, step by step execution takes place. Backward feature elimination procedure is a time-consuming procedure. If we are dealing with the full dataset, we may have to wait for a long time to process the dataset. Alternatively, we can take advantage of parallel processing using the R package “parallel” to reduce the processing time. To reduce the build-up time of this package sometime, we use to process data in this vignette.
## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: GenomicRanges## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, basename, cbind, colnames, dirname, do.call,
## duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
## lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:base':
##
## I, expand.grid, unname## Loading required package: IRanges## Loading required package: GenomeInfoDb## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians## class: SummarizedExperiment
## dim: 2204 27
## metadata(0):
## assays(1): counts
## rownames(2204): 1494_f_at 179_at ... AFFX-M27830_M_at
## AFFX-r2-Hs28SrRNA-5_at
## rowData names(0):
## colnames(27): GSM42246 GSM42248 ... GSM42270 GSM42271
## colData names(1): sampleLabelsIn this example, the variable “x” is containing the microarray expression data, and the “y” contains the sample labels. The variable “x” is a two-dimensional matrix in which the row contains the sample names, and the column contains the gene/feature names. The variable “y” contains the sample labels as -1 and 1.
Before going into details of the vignette, let us plot the unadjusted p-value for each features which are estimated by using the expression values of the two classes. The function named “sigFeaturePvalue()” is used here to calculate the p-value using t-statistic. In this case, we assume unequal variances and using the Welch approximation to estimate the appropriate degrees of freedom.
pvals <- sigFeaturePvalue(x,y)
hist(unlist(pvals),breaks=seq(0,0.08,0.0015),col="skyblue",
xlab="p value",ylab="Frequency",main="")
The idea of “sigFeature()” (Significant Feature Selection): The recursive feature elimination that uses support vector machine can rank better classifiers, but those classifiers may or may not be differently significant between the classes. Features with notable classification accuracy and also differentially significant between the classes have a predominant role in a biological aspect. The algorithm, introduced here is very efficient in ranking classifiers, as well as in determining differently significant classes among them.
In data mining and optimization, feature selection is a very active field of research where the SVM-RFE is distinguished as one of the most effective methods. SVM-RFE is a greedy method that only hopes to find the best possible combination for classification. In this proposed algorithm, we have included a t-statistics along with SVM-RFE. The primary objective of this algorithm is to enumerate the ranking weights for all features and sort the features according to weight vectors as the classification basis. The coefficient and the expression mean difference between two compared groups of each individual feature is used to calculate the weight value of that particular gene or the feature. The iteration process is followed by backward removal of the feature. The iteration process is continued until there is only one feature remaining in the dataset. The smallest ranking weight will be removed by the algorithm and stored in a stack, while remaining the feature variables have a significant impact for producing weight. Finally, the feature variables which are stored in the stack will be listed in the descending order of descriptive different degree.
The variable “sigfeatureRankedList” will return the address of the genes/features after execution. To visualize the output, print command can be used that will print the address of top 10 ranked genes/features.
## [1] 2064 370 2032 2035 1519 1573 1446 2105 997 611For producing the model vector for training the data set by using the top 1000 genes/features, the R script is given as below. In this case the package “e1071” is used. For producing the model vector, linear kernel is chosen to get better performance.
library(e1071)
sigFeature.model=svm(x[ ,sigfeatureRankedList[1:1000]], y,
type="C-classification", kernel="linear")
summary(sigFeature.model)##
## Call:
## svm.default(x = x[, sigfeatureRankedList[1:1000]], y = y, type = "C-classification",
## kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1
##
## Number of Support Vectors: 16
##
## ( 10 6 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1The R script below shows the prediction result of the same dataset that was trained before. The user can also separate the dataset into two parts e.g. training dataset and testing dataset and can implement the algorithm.
## y
## pred -1 1
## -1 8 0
## 1 0 19To solve the classification problem with the help of ranking the features, another algorithm was proposed by Guyon et. al 2002 named as SVM-RFE. In this algorithm, the dataset has been trained with SVM linear kernel model to produce a weight vector for each feature. After that, the feature with smallest rank is recursively removed and stored in a stack data structure. The iteration process continues till the last feature remains. This criterion is the w value of the decision hyperplane given by the SVM. The weight vector is calculated by squiring the variable w. Finally the features are selected in a list in the descending order. The function “svmrfeFeatureRanking()” is included in this package to compare the performance with “sigFeature()”.
#system.time(featureRankedList <- svmrfeFeatureRanking(x, y))
data(featureRankedList)
print("Top 10 features are printed below:")## [1] "Top 10 features are printed below:"## [1] 1073 1404 1152 5 1253 1557 105 1207 792 57For producing the model vector of the training dataset by using the top 1000 genes/features, the R script is given below. The linear kernel is chosen for producing the model vector.
RFE.model=svm(x[ ,featureRankedList[1:1000]], y,
type="C-classification", kernel="linear")
summary(RFE.model)##
## Call:
## svm.default(x = x[, featureRankedList[1:1000]], y = y, type = "C-classification",
## kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1
##
## Number of Support Vectors: 17
##
## ( 10 7 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1The R script below will show the prediction result of the same dataset which is trained before. The user can also separate the dataset into two parts e.g. training dataset and testing dataset and can implement the algorithm.
## y
## pred -1 1
## -1 8 0
## 1 0 19In this example section we will observe the p value of the selected top 100 features. The function named “sigFeaturePvalue()” will calculate the p value using t-statistic. The expression value will be taken from the dataset of top 100 features.
pvalsigFe <- sigFeaturePvalue(x, y, 100, sigfeatureRankedList)
pvalRFE <- sigFeaturePvalue(x, y, 100, featureRankedList)
par(mfrow=c(1,2))
hist(unlist(pvalsigFe),breaks=50, col="skyblue", main=paste("sigFeature"),
xlab="p value")
hist(unlist(pvalRFE),breaks=50, col="skyblue",
main=paste("SVM-RFE"), xlab="p value")
In this section a comparison of p values produced by using top 100 features (using “sigFeature” and “SVM-RFE” algorithm) is shown by using scatter box plot.
mytitle<-'Box Plot'
boxplot(unlist(pvalsigFe), unlist(pvalRFE), main=mytitle,
names=c("sigFeature", "SVM-RFE"),
ylab="p value", ylim=c(min(unlist(pvalsigFe)), max(unlist(pvalRFE))))
stripchart(unlist(pvalsigFe), vertical=TRUE, method="jitter", add=TRUE, pch=16,
col=c('green'))
stripchart(unlist(pvalRFE), vertical=TRUE, at=2, method="jitter", add=TRUE,
pch=16, col=c('blue'))
grid(nx=NULL, ny=NULL, col="black", lty="dotted")
In this section a heatmap is shown below with the expression value of top 20 features generated by using “sigFeature()” function.
library("pheatmap")
library("RColorBrewer")
pheatmap(x[ ,sigfeatureRankedList[1:20]], scale="row",
clustering_distance_rows="correlation")
In this section a heatmap is shown below with the expression value of top 20 features generated by using “svmrfeFeatureRanking()” function.
For estimation of feature selection algorithms, Ambroise and McLachlan (2002) suggested external cross-validation, in which the feature selection and validation are performed on different parts of the sample set, to obtain an unbiased performance estimate. The stratified cross-validation (Breiman et al.,1984) is slightly different from regular cross-validation. In k-fold stratified cross-validation, a data set is randomly partitioned into k equally sized folds such that the class distribution in each fold is approximately the same as that in the entire data set. In contrast, regular cross-validation randomly partitions the data set into k-folds without considering class distributions. Kohavi (1995) reported that stratified cross-validation has smaller bias and variance than regular cross-validation. The “sigFeature()” function thus farther enhances by incorporating external cross-validation method. In this external k-fold cross-validation procedure, k-1 folds are used for selecting the feature, and one fold remains untouched which will later use as a test sample set.
Detail description of the input parameters is given in the help file. The set.seed function is used here only for the purpose of getting the exact results when regenerating this vignette from its source files. The function is given below.
## List of 1
## $ :List of 5
## ..$ feature.ids : int [1:2204] 2064 2031 2035 1573 370 2023 1605 2032 997 1069 ...
## ..$ train.data.ids : chr [1:24] "GSM42246" "GSM42248" "GSM42250" "GSM42252" ...
## ..$ test.data.ids : chr [1:3] "GSM42263" "GSM42251" "GSM42260"
## ..$ train.data.level: Named num [1:24] 1 1 1 1 1 -1 -1 -1 -1 -1 ...
## .. ..- attr(*, "names")= chr [1:24] "GSM42246" "GSM42248" "GSM42250" "GSM42252" ...
## ..$ test.data.level : Named num [1:3] -1 1 1
## .. ..- attr(*, "names")= chr [1:3] "GSM42263" "GSM42251" "GSM42260"In this example section a new function is introduced named as “sigFeatureFrequency()”. The main purpose of this function is to arrange the feature on the basis of its frequency. In the previous section the dataset is divided into k-folds. Out of which k-1 folds are randomly taken and used for selecting the feature. This procedure is repeated k times and the frequency values are calculated for top 400 genes from the k number of lists. The “sigFeatureFrequency()” function is used here to rank further all the k’th feature list according to its frequency. Details description of this function is given in the help file.
FeatureBasedonFrequency <- sigFeatureFrequency(x, results, 400, 400, pf=FALSE)
str(FeatureBasedonFrequency[1])## List of 1
## $ :List of 5
## ..$ feature.ids : num [1:400] 187 225 246 303 313 370 394 469 479 492 ...
## ..$ train.data.ids : chr [1:24] "GSM42246" "GSM42248" "GSM42250" "GSM42252" ...
## ..$ test.data.ids : chr [1:3] "GSM42263" "GSM42251" "GSM42260"
## ..$ train.data.level: Named num [1:24] 1 1 1 1 1 -1 -1 -1 -1 -1 ...
## .. ..- attr(*, "names")= chr [1:24] "GSM42246" "GSM42248" "GSM42250" "GSM42252" ...
## ..$ test.data.level : Named num [1:3] -1 1 1
## .. ..- attr(*, "names")= chr [1:3] "GSM42263" "GSM42251" "GSM42260"In this section external cross-validation is computed and finally mean cross validation errors and the standard deviation of the errors is returnd. The features, which are produced on the basis of frequency is used for producing the mean cross validation error. Training and test samples are same which were initially split from the main sample set in to k-fold. In each iteration k-1 fold are considered as training sample and remaining one fold is considered as testing sample. In this external cross validation procedure feature number is increased one by one and using that feature training sample is trained. The number of miss classification is averaged out for each fold change which finally represented as error rate.
#inputdata <- data.frame(y=as.factor(y) ,x=x)
#To run the code given bellow will take huge time. Thus the process
#data is given below.
#featsweepSigFe = lapply(1:400, sigCVError, FeatureBasedonFrequency, inputdata)
data("featsweepSigFe")
str(featsweepSigFe[1])## List of 1
## $ :List of 3
## ..$ svm.list:List of 10
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.333
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.333
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.667
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.333
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.333
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 1
## .. .. ..$ cost : num 1.78
## .. .. ..$ error : num 0.667
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.5
## .. .. ..$ dispersion: num NA
## .. ..$ :'data.frame': 1 obs. of 4 variables:
## .. .. ..$ gamma : num 0.000244
## .. .. ..$ cost : num 0.001
## .. .. ..$ error : num 0.5
## .. .. ..$ dispersion: num NA
## ..$ error : num 0.367
## ..$ errorSD : num 0.233In this section the “PlotErrors()” function draw two plots. The first plot is representing the average external cross validation error and the second plot is representing the standard deviation of the of the miss classifications.
Writing the results can be achieved with the “WritesigFeature()” function. Please refer to the help file of “WritesigFeature()” for more details on how to save to results to a specific file.