Ensemble filters do not select the top ranked features

[annette987]

The code below shows two methods of getting the top-ranked features from ensemble filter E-mean. Both should return the same result, if the base filters are the same, but they do not.

library(mlr)
#> Loading required package: ParamHelpers
library(survival)
library(testthat)
library(checkmate)

data(veteran)
vet.task <- makeSurvTask(id = "VET", data = veteran, target = c("time", "status"))
vet.task <- createDummyFeatures(vet.task)

# Below are two methods for getting the top-ranked features from E-mean
# Both should return the same result, if the base filters are the same

# Get the 5 top-ranked features generated by ensemble filter E-mean, 
# by calling the function for that filter directly
set.seed(24601)
filt = mlr:::.FilterEnsembleRegister[["E-mean"]]
cox.lrn <- makeLearner(cl="surv.coxph", id = "coxph", predict.type="response")
dat = filt$fun(vet.task, 
               base.methods = c("univariate.model.score", "randomForestSRC_importance"),
               nselect = 5,
               more.args = list("univariate.model.score"=list(perf.learner=cox.lrn))
)
ordered_dat = dat[with(dat, order(method, -value)), ]
top_features1 = ordered_dat[ordered_dat$method == "E-mean", ]$name[1:5]

# Get the 5 top-ranked features generated by ensemble filter E-mean, 
# by caling filterFeatures with filter E-mean
set.seed(24601)
new.task = filterFeatures(vet.task,
                          method="E-mean",
                          abs = 5,
                          base.methods = c("univariate.model.score", "randomForestSRC_importance"))
top_features2 = getTaskFeatureNames(new.task)
expect_equal(sort(top_features1), sort(top_features2))
#> Error: sort(top_features1) not equal to sort(top_features2).
#> 3/5 mismatches
#> x[1]: "age"
#> y[1]: "celltype.adeno"
#> 
#> x[2]: "celltype.adeno"
#> y[2]: "celltype.large"
#> 
#> x[5]: "prior"
#> y[5]: "karno"

^{Created on 2019-11-25 by the reprex package (v0.3.0)}

There are several issues that are preventing ensemble filters from selecting the top-ranked features.

1. Each ensemble filter calls calcBaseFilters which gets the values produced by the base filters and ranks the features in order of the value they produce. However, the method used to rank the features in this function does not work. Here is evidence:

library(testthat)

toy.data= data.frame(
  name = rep(c("a", "b", "c", "d", "e"), 2),
  method = c(rep("m1", 5), rep("m2", 5)),
  value = rep(c(0.3, 0.5, 0.1, 0.2, 0.4), 2)
)
toy.data
#>    name method value
#> 1     a     m1   0.3
#> 2     b     m1   0.5
#> 3     c     m1   0.1
#> 4     d     m1   0.2
#> 5     e     m1   0.4
#> 6     a     m2   0.3
#> 7     b     m2   0.5
#> 8     c     m2   0.1
#> 9     d     m2   0.2
#> 10    e     m2   0.4

toy.data.ranked = transform(toy.data,
                            rank = ave(1:nrow(toy.data), method,
                                       FUN = function(x) order(toy.data$value[x])))
toy.data.ranked
#>    name method value rank
#> 1     a     m1   0.3    3
#> 2     b     m1   0.5    4
#> 3     c     m1   0.1    1
#> 4     d     m1   0.2    5
#> 5     e     m1   0.4    2
#> 6     a     m2   0.3    3
#> 7     b     m2   0.5    4
#> 8     c     m2   0.1    1
#> 9     d     m2   0.2    5
#> 10    e     m2   0.4    2
expect_equal(toy.data.ranked$rank, rep(c(3, 1, 5, 4, 2), 2))
#> Error: toy.data.ranked$rank not equal to rep(c(3, 1, 5, 4, 2), 2).
#> 6/10 mismatches (average diff: 2.67)
#> [2] 4 - 1 ==  3
#> [3] 1 - 5 == -4
#> [4] 5 - 4 ==  1
#> [7] 4 - 1 ==  3
#> [8] 1 - 5 == -4
#> [9] 5 - 4 ==  1

^{Created on 2019-11-25 by the reprex package (v0.3.0)}

2. In filterFeatures.R, the method to use for the final ranking must be selected. This is done by the line below, but this has no effect:
   fval = fval[fval$method == fval$method, ]
   As a result fval still contains the data for all methods, not just the ensemble method.

3.When the top features are selected by setting a threshold, this must be done after the subsetting referred to in point 2. If not nselect can be greater than the total number of features.

  if (select == "threshold") {
    nselect = sum(fval[["value"]] >= threshold, na.rm = TRUE)
  }

If fval contains data for all methods, there will be multiple rows for each feature - one per method. So the calculation above can include the same feature multiple times and thus nselect could be greater than the number of features.

[annette987]

Please note that the issue of nselect being able to be greater than the number of features when a threshold is used does not show up in practice because getTaskData, called from subsetTask, at the end of filterFeatures, corrects this problem when it uses set operations like setdiff and union. But if you examine the value of nselect in filterFeatures you will see it.

[pat-s]

@annette987 Thanks for the effort in this issue. Even though I do not fully understand what you want to do in your "mlr:::.FilterEnsembleRegister" example, I think there is an issue in the aggregation of ensemble filters right now.
I've created a small reprex of the part which is in my view causing the issue right now - the creation of the "rank" column in

mlr/R/FilterEnsemble.R

Lines 275 to 277 in 27f0ba4

	fval.all.calced.simple = transform(fval.calc$data,
	rank = ave(1:nrow(fval.calc$data), method,
	FUN = function(x) order(fval.calc$data$value[x])))

does not work as intended.

There might be a regression in some commits in the past that caused this.
I have an idea but need to confirm.
Also we need some tests to make this more safe.

Please note that the issue of nselect being able to be greater than the number of features when a threshold is used does not show up in practice because getTaskData, called from subsetTask, at the end of filterFeatures, corrects this problem when it uses set operations like setdiff and union. But if you examine the value of nselect in filterFeatures you will see it.

Is this a separate issue? I've read it twice now and do not yet think that I understand the problem.
For now, let's focus on getting the right ordering of ranks again in this issue.

library(mlr)
#> Loading required package: ParamHelpers

data = structure(list(task.desc = structure(list(id = "VET", type = "surv",
  target = c("time", "status"), size = 137L, n.feat = c(numerics = 9L,
    factors = 0L, ordered = 0L, functionals = 0L), has.missings = FALSE,
  has.weights = FALSE, has.blocking = FALSE, has.coordinates = FALSE), class = c("SurvTaskDesc",
  "SupervisedTaskDesc", "TaskDesc")), data = structure(list(name = c("trt",
  "karno", "diagtime", "age", "prior", "celltype.squamous", "celltype.smallcell",
  "celltype.adeno", "celltype.large", "trt", "karno", "diagtime",
  "age", "prior", "celltype.squamous", "celltype.smallcell", "celltype.adeno",
  "celltype.large"), type = c("numeric", "numeric", "numeric",
  "numeric", "numeric", "numeric", "numeric", "numeric", "numeric",



  "numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
  "numeric", "numeric", "numeric"), method = structure(c(1L, 1L,
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
), .Label = c("univariate.model.score", "variance"), class = "factor"),
value = c(0.397014925373134, 0.620398009950249, 0.551741293532338,
  0.46865671641791, 0.5, 0.530845771144279, 0.580597014925373,
  0.561194029850746, 0.5, 0.251824817518248, 401.585229712323,
  112.617539716617, 111.125912408759, 20.8243881494204, 0.191605839416058,
  0.229282954057535, 0.159403177329326, 0.159403177329326)), row.names = c(NA,
  -18L), class = c("data.table", "data.frame"))))

# rank base filters by method
transform(data$data,
  rank = ave(1:nrow(data$data), method,
    FUN = function(x) order(data$data$value[x])))
#>                   name    type                 method       value rank
#>  1:                trt numeric univariate.model.score   0.3970149    1
#>  2:              karno numeric univariate.model.score   0.6203980    4
#>  3:           diagtime numeric univariate.model.score   0.5517413    5
#>  4:                age numeric univariate.model.score   0.4686567    9
#>  5:              prior numeric univariate.model.score   0.5000000    6
#>  6:  celltype.squamous numeric univariate.model.score   0.5308458    3
#>  7: celltype.smallcell numeric univariate.model.score   0.5805970    8
#>  8:     celltype.adeno numeric univariate.model.score   0.5611940    7
#>  9:     celltype.large numeric univariate.model.score   0.5000000    2
#> 10:                trt numeric               variance   0.2518248    8
#> 11:              karno numeric               variance 401.5852297    9
#> 12:           diagtime numeric               variance 112.6175397    6
#> 13:                age numeric               variance 111.1259124    7
#> 14:              prior numeric               variance  20.8243881    1
#> 15:  celltype.squamous numeric               variance   0.1916058    5
#> 16: celltype.smallcell numeric               variance   0.2292830    4
#> 17:     celltype.adeno numeric               variance   0.1594032    3
#> 18:     celltype.large numeric               variance   0.1594032    2

^{Created on 2019-12-03 by the reprex package (v0.3.0)}

[pat-s]

The part shown above should behave as follows: It should rank the values by method and then order by rank (actually in reverse because high values of values should get a higher rank value).

I have not much time tomorrow but will focus on getting this fixed in the next days including a fast release to CRAN afterwards.

library(mlr)
#> Loading required package: ParamHelpers
library(data.table)

data = structure(list(task.desc = structure(list(id = "VET", type = "surv",
  target = c("time", "status"), size = 137L, n.feat = c(numerics = 9L,
    factors = 0L, ordered = 0L, functionals = 0L), has.missings = FALSE,
  has.weights = FALSE, has.blocking = FALSE, has.coordinates = FALSE), class = c("SurvTaskDesc",
  "SupervisedTaskDesc", "TaskDesc")), data = structure(list(name = c("trt",
  "karno", "diagtime", "age", "prior", "celltype.squamous", "celltype.smallcell",
  "celltype.adeno", "celltype.large", "trt", "karno", "diagtime",
  "age", "prior", "celltype.squamous", "celltype.smallcell", "celltype.adeno",
  "celltype.large"), type = c("numeric", "numeric", "numeric",
  "numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
  "numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
  "numeric", "numeric", "numeric"), method = structure(c(1L, 1L,
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
), .Label = c("univariate.model.score", "variance"), class = "factor"),
value = c(0.397014925373134, 0.620398009950249, 0.551741293532338,
  0.46865671641791, 0.5, 0.530845771144279, 0.580597014925373,
  0.561194029850746, 0.5, 0.251824817518248, 401.585229712323,
  112.617539716617, 111.125912408759, 20.8243881494204, 0.191605839416058,
  0.229282954057535, 0.159403177329326, 0.159403177329326)), row.names = c(NA,
  -18L))))

data$data = as.data.table(data$data)


setkey(data$data, method)

data$data[, rank := frank(-value, ties.method = "first"), by = method]
setorderv(data$data, c("method", "rank"))

data$data
#>                   name    type                 method       value rank
#>  1:              karno numeric univariate.model.score   0.6203980    1
#>  2: celltype.smallcell numeric univariate.model.score   0.5805970    2
#>  3:     celltype.adeno numeric univariate.model.score   0.5611940    3
#>  4:           diagtime numeric univariate.model.score   0.5517413    4
#>  5:  celltype.squamous numeric univariate.model.score   0.5308458    5
#>  6:              prior numeric univariate.model.score   0.5000000    6
#>  7:     celltype.large numeric univariate.model.score   0.5000000    7
#>  8:                age numeric univariate.model.score   0.4686567    8
#>  9:                trt numeric univariate.model.score   0.3970149    9
#> 10:              karno numeric               variance 401.5852297    1
#> 11:           diagtime numeric               variance 112.6175397    2
#> 12:                age numeric               variance 111.1259124    3
#> 13:              prior numeric               variance  20.8243881    4
#> 14:                trt numeric               variance   0.2518248    5
#> 15: celltype.smallcell numeric               variance   0.2292830    6
#> 16:  celltype.squamous numeric               variance   0.1916058    7
#> 17:     celltype.adeno numeric               variance   0.1594032    8
#> 18:     celltype.large numeric               variance   0.1594032    9

^{Created on 2019-12-03 by the reprex package (v0.3.0)}

[annette987]

Be careful about ranking in descending order of value in calcBaseFilters. filterFeatures also ranks in descending order of value:

  if (nselect > 0L) {

    # order by method and (desc(value))
    features = fval[with(fval, order(method, -value)), ]
    # select names of top n
    features = features[1:nselect, ]$name

  } else {
    features = NULL
  }

So if the code that does the ranking (in calcBaseFilters/FilterEnsemble.R) is called via filterFeatures then it will be ranked in descending order twice - which will mean it will be ranked in ascending order.
Apart from that your code will solve the ranking problem. Just curious as to why you prefer this to my one line solution. Easier to understand perhaps?

[annette987]

Regarding the issue with the timing of the calculation of nselect:

if (select == "threshold") {
    nselect = sum(fval[["value"]] >= threshold, na.rm = TRUE)
  }

In the situation where the user has passed a threshold, nselect is currently calculated on the full data frame returned from generateFilterValuesData. Hopefully this reprex will explain what is happening:

library(survival)
library(mlr)
#> Loading required package: ParamHelpers

data(veteran)
set.seed(24601)
vet.task <- makeSurvTask(id = "VET", data = veteran, target = c("time", "status"))
vet.task <- createDummyFeatures(vet.task)

cox.lrn <- makeLearner(cl="surv.coxph", id = "coxph", predict.type="response")
fval =generateFilterValuesData(vet.task, 
                             method = list("E-mean", c("univariate.model.score", "randomForestSRC_importance")),
                             more.args=list("univariate.model.score"=list(perf.learner=cox.lrn), "randomForestSRC_importance"=list(ntree=100))
)
fval = fval$data
fval
#>                   name    type                     method         value
#>  1:              prior numeric randomForestSRC_importance -0.0034653967
#>  2:                trt numeric randomForestSRC_importance -0.0002310563
#>  3:           diagtime numeric randomForestSRC_importance  0.0007997677
#>  4:                age numeric randomForestSRC_importance  0.0020504755
#>  5:     celltype.large numeric randomForestSRC_importance  0.0084291392
#>  6:     celltype.adeno numeric randomForestSRC_importance  0.0088886429
#>  7:  celltype.squamous numeric randomForestSRC_importance  0.0111077710
#>  8: celltype.smallcell numeric randomForestSRC_importance  0.0137876876
#>  9:              karno numeric randomForestSRC_importance  0.1285040870
#> 10:              prior numeric     univariate.model.score  0.4085623679
#> 11:                trt numeric     univariate.model.score  0.4474747475
#> 12:                age numeric     univariate.model.score  0.4882100750
#> 13:     celltype.large numeric     univariate.model.score  0.5371655104
#> 14:           diagtime numeric     univariate.model.score  0.5669050051
#> 15:  celltype.squamous numeric     univariate.model.score  0.5669882101
#> 16:     celltype.adeno numeric     univariate.model.score  0.5731948566
#> 17: celltype.smallcell numeric     univariate.model.score  0.5972083749
#> 18:              karno numeric     univariate.model.score  0.6612903226
#> 19:                age numeric                     E-mean  6.5000000000
#> 20:     celltype.adeno numeric                     E-mean  7.0000000000
#> 21:     celltype.large numeric                     E-mean  2.0000000000
#> 22: celltype.smallcell numeric                     E-mean  7.0000000000
#> 23:  celltype.squamous numeric                     E-mean  7.0000000000
#> 24:           diagtime numeric                     E-mean  3.5000000000
#> 25:              karno numeric                     E-mean  1.0000000000
#> 26:              prior numeric                     E-mean  6.0000000000
#> 27:                trt numeric                     E-mean  5.0000000000
#>                   name    type                     method         value

threshold = 0.5
nselect = sum(fval[["value"]] >= threshold, na.rm = TRUE)
nselect
#> [1] 15

^{Created on 2019-12-04 by the reprex package (v0.3.0)}

In this case nselect is set to 15, greater than the number of features, because the threshold is evaluated across the values returned by all 3 filters.

Instead nselect needs to be calculated after the filter to be used for the final ranking is determined i.e. after this block of code:

  # in case multiple filters have been calculated, choose which ranking is used
  # for the final subsetting
  if (length(levels(as.factor(fval$method))) >= 2) {
    # if 'method' is an ensemble method, we always choose the ensemble method
    # unless select.method is specified specifically. Method[[1]] should usually
    # be the ensemble method
    if (is.null(select.method) && !method[[1]] %in% ls(.FilterEnsembleRegister)) {
      stopf("You supplied multiple filters. Please choose which should be used for the final subsetting of the features.")
    }
    if (is.null(select.method)) {
      fval = fval[fval$method == fval$method, ]
    } else {
      assertSubset(select.method, choices = unique(fval$method))
      fval = fval[fval$method == select.method, ]
    }
  }

I hope that makes more sense.

[pat-s]

Be careful about ranking in descending order of value in calcBaseFilters. filterFeatures also ranks in descending order of value:

Thanks. I've added two tests to ensure that both the ranking and subsetting is being done correctly. Tests check whether the variables "value" and "rank" are sorted and that the corresponding maxima of both have the same row.

Just curious as to why you prefer this to my one line solution. Easier to understand perhaps?

We are enforcing more and more data.table syntax throughout the whole mlr-org projects and since we're more often operating internally on data.table's than on data.frame's, I'd like to go with this solution.

Would you mind having a close look at #2698 and report shortly if you agree that now everything is being returned as expected?

---

I'll move the thresholding part to a new issue: #2699

[annette987]

Yes I will have a look and get back to you.

[annette987]

Sorry I couldn't get back to this sooner. I have had a good look at #2698 and it is passing all my tests. Thanks.