特征选择

1
2
3
4
5
6

# Load Packages and prepare dataset
library(TH.data)
library(caret)
data("GlaucomaM", package = "TH.data")
trainData <- GlaucomaM
head(trainData)

1
2

# install.packages('Boruta')
library(Boruta)

1
2

# Perform Boruta search
boruta_output <- Boruta(Class ~ ., data=na.omit(trainData), doTrace=0)  

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

names(boruta_output)
1. ‘finalDecision’
2. ‘ImpHistory’
3. ‘pValue’
4. ‘maxRuns’
5. ‘light’
6. ‘mcAdj’
7. ‘timeTaken’
8. ‘roughfixed’
9. ‘call’
10. ‘impSource’

1
2
3
4
5
6
7

# 得到显著的或者潜在的变量
boruta_signif <- getSelectedAttributes(boruta_output, withTentative = TRUE)
print(boruta_signif)  
 [1] "as"   "ean"  "abrg" "abrs" "abrn" "abri" "hic"  "mhcg" "mhcn" "mhci"
[11] "phcg" "phcn" "phci" "hvc"  "vbss" "vbsn" "vbsi" "vasg" "vass" "vasi"
[21] "vbrg" "vbrs" "vbrn" "vbri" "varg" "vart" "vars" "varn" "vari" "mdn" 
[31] "tmg"  "tmt"  "tms"  "tmn"  "tmi"  "rnf"  "mdic" "emd" 

1
2
3
4
5
6
7

# Do a tentative rough fix
roughFixMod <- TentativeRoughFix(boruta_output)
boruta_signif <- getSelectedAttributes(roughFixMod)
print(boruta_signif)
 [1] "abrg" "abrs" "abrn" "abri" "hic"  "mhcg" "mhcn" "mhci" "phcg" "phcn"
[11] "phci" "hvc"  "vbsn" "vbsi" "vasg" "vbrg" "vbrs" "vbrn" "vbri" "varg"
[21] "vart" "vars" "varn" "vari" "tmg"  "tms"  "tmi"  "rnf"  "mdic" "emd" 

1
2
3
4

# Variable Importance Scores
imps <- attStats(roughFixMod)
imps2 = imps[imps$decision != 'Rejected', c('meanImp', 'decision')]
head(imps2[order(-imps2$meanImp), ])  # descending sort

1
2

# Plot variable importance
plot(boruta_output, cex.axis=.7, las=2, xlab="", main="Variable Importance")  

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# 训练rpart模型，计算变量重要性
library(caret)
set.seed(100)
rPartMod <- train(Class ~ ., data=trainData, method="rpart")
rpartImp <- varImp(rPartMod)
print(rpartImp)
rpart variable importance

  only 20 most important variables shown (out of 62)

     Overall
varg  100.00
vari   93.19
vars   85.20
varn   76.86
tmi    72.31
vbss    0.00
eai     0.00
tmg     0.00
tmt     0.00
vbst    0.00
vasg    0.00
at      0.00
abrg    0.00
vbsg    0.00
eag     0.00
phcs    0.00
abrs    0.00
mdic    0.00
abrt    0.00
ean     0.00

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# Train an RRF model and compute variable importance.
set.seed(100)
rrfMod <- train(Class ~ ., data=trainData, method="RRF")
rrfImp <- varImp(rrfMod, scale=F)
rrfImp
RRF variable importance

  only 20 most important variables shown (out of 62)

     Overall
varg 24.0013
vari 18.5349
vars  6.0483
tmi   3.8699
hic   3.3926
mhci  3.1856
mhcg  3.0383
mv    2.1570
hvc   2.1357
phci  1.8830
vasg  1.8570
tms   1.5705
phcn  1.4475
phct  1.4473
vass  1.3097
tmt   1.2485
phcg  1.1992
mdn   1.1737
tmg   1.0988
abrs  0.9537
plot(rrfImp, top = 20, main='Variable Importance')

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

library(glmnet)
trainData <- read.csv('https://raw.githubusercontent.com/selva86/datasets/master/GlaucomaM.csv')

x <- as.matrix(trainData[,-63]) # all X vars
y <- as.double(as.matrix(ifelse(trainData[, 63]=='normal', 0, 1))) # Only Class

# Fit the LASSO model (Lasso: Alpha = 1)
set.seed(100)
cv.lasso <- cv.glmnet(x, y, family='binomial', alpha=1, parallel=TRUE, standardize=TRUE, type.measure='auc')

# Results
plot(cv.lasso)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# plot(cv.lasso$glmnet.fit, xvar="lambda", label=TRUE)
cat('Min Lambda: ', cv.lasso$lambda.min, '\n 1Sd Lambda: ', cv.lasso$lambda.1se)
df_coef <- round(as.matrix(coef(cv.lasso, s=cv.lasso$lambda.min)), 2)

# See all contributing variables
df_coef[df_coef[, 1] != 0, ]
Min Lambda:  0.01166507 
 1Sd Lambda:  0.2513163

Min Lambda:  0.01166507 
1Sd Lambda:  0.2513163

(Intercept) 3.65
at         -0.17
as         -2.05
eat        -0.53
mhci        6.22
phcs       -0.83
phci        6.03
hvc        -4.15
vass       -23.72
vbrn       -0.26
vars       -25.86
mdt        -2.34
mds         0.5
mdn         0.83
mdi         0.3
tmg         0.01
tms         3.02
tmi         2.65
mv          4.94

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# Load data
trainData <- read.csv("http://rstatistics.net/wp-content/uploads/2015/09/ozone1.csv", stringsAsFactors=F)
print(head(trainData))
  Month Day_of_month Day_of_week ozone_reading pressure_height Wind_speed
1     1            1           4             3            5480          8
2     1            2           5             3            5660          6
3     1            3           6             3            5710          4
4     1            4           7             5            5700          3
5     1            5           1             5            5760          3
6     1            6           2             6            5720          4

  Humidity Temperature_Sandburg Temperature_ElMonte Inversion_base_height
1 20.00000             40.53473            39.77461              5000.000
2 40.96306             38.00000            46.74935              4108.904
3 28.00000             40.00000            49.49278              2693.000
4 37.00000             45.00000            52.29403               590.000
5 51.00000             54.00000            45.32000              1450.000
6 69.00000             35.00000            49.64000              1568.000

  Pressure_gradient Inversion_temperature Visibility
1               -15              30.56000        200
2               -14              48.02557        300
3               -25              47.66000        250
4               -24              55.04000        100
5                25              57.02000         60
6                15              53.78000         60

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# Step 1: 定义基础的模型，只有截距
base.mod <- lm(ozone_reading ~ 1 , data=trainData)  

# Step 2: 包含所有特征的模型
all.mod <- lm(ozone_reading ~ . , data= trainData) 

# Step 3: 进行逐步回归算法，方向设置为双向forward and backward stepwise
stepMod <- step(base.mod, scope = list(lower = base.mod, upper = all.mod), direction = "both", trace = 0, steps = 1000)  

# Step 4: 得到变量列表
shortlistedVars <- names(unlist(stepMod[[1]])) 
shortlistedVars <- shortlistedVars[!shortlistedVars %in% "(Intercept)"] # remove intercept

# Show
print(shortlistedVars)
[1] "Temperature_Sandburg"  "Humidity"              "Temperature_ElMonte"  
[4] "Month"                 "pressure_height"       "Inversion_base_height"

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# install.packages('relaimpo')
library(relaimpo)

# Build linear regression model
model_formula = ozone_reading ~ Temperature_Sandburg + Humidity + Temperature_ElMonte + Month + pressure_height + Inversion_base_height
lmMod <- lm(model_formula, data=trainData)

# calculate relative importance
relImportance <- calc.relimp(lmMod, type = "lmg", rela = F)  

# Sort
cat('Relative Importances: \n')
sort(round(relImportance$lmg, 3), decreasing=TRUE)
Relative Importances: 
Temperature_ElMonte    0.214
Temperature_Sandburg   0.203
pressure_height        0.104
Inversion_base_height  0.096
Humidity               0.086
Month                  0.012

1
2
3
4

bootsub <- boot.relimp(ozone_reading ~ Temperature_Sandburg + Humidity + Temperature_ElMonte + Month + pressure_height + Inversion_base_height, data=trainData,
                       b = 1000, type = 'lmg', rank = TRUE, diff = TRUE)

plot(booteval.relimp(bootsub, level=.95))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

str(trainData)
'data.frame':    366 obs. of  13 variables:
 $ Month                : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Day_of_month         : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Day_of_week          : int  4 5 6 7 1 2 3 4 5 6 ...
 $ ozone_reading        : num  3 3 3 5 5 6 4 4 6 7 ...
 $ pressure_height      : num  5480 5660 5710 5700 5760 5720 5790 5790 5700 5700 ...
 $ Wind_speed           : int  8 6 4 3 3 4 6 3 3 3 ...
 $ Humidity             : num  20 41 28 37 51 ...
 $ Temperature_Sandburg : num  40.5 38 40 45 54 ...
 $ Temperature_ElMonte  : num  39.8 46.7 49.5 52.3 45.3 ...
 $ Inversion_base_height: num  5000 4109 2693 590 1450 ...
 $ Pressure_gradient    : num  -15 -14 -25 -24 25 15 -33 -28 23 -2 ...
 $ Inversion_temperature: num  30.6 48 47.7 55 57 ...
 $ Visibility           : int  200 300 250 100 60 60 100 250 120 120 ...
set.seed(100)
options(warn=-1)

subsets <- c(1:5, 10, 15, 18)

ctrl <- rfeControl(functions = rfFuncs,
                   method = "repeatedcv",
                   repeats = 5,
                   verbose = FALSE)

lmProfile <- rfe(x=trainData[, c(1:3, 5:13)], y=trainData$ozone_reading,
                 sizes = subsets,
                 rfeControl = ctrl)

lmProfile
Recursive feature selection

Outer resampling method: Cross-Validated (10 fold, repeated 5 times) 

Resampling performance over subset size:

 Variables  RMSE Rsquared   MAE RMSESD RsquaredSD  MAESD Selected
         1 5.222   0.5794 4.008 0.9757    0.15034 0.7879         
         2 3.971   0.7518 3.067 0.4614    0.07149 0.3276         
         3 3.944   0.7553 3.054 0.4675    0.06523 0.3708         
         4 3.924   0.7583 3.026 0.5132    0.06640 0.4163         
         5 3.880   0.7633 2.950 0.5525    0.07021 0.4334         
        10 3.751   0.7796 2.853 0.5550    0.06791 0.4457        *
        12 3.767   0.7779 2.869 0.5511    0.06664 0.4424         

The top 5 variables (out of 10):
   Temperature_ElMonte, Pressure_gradient, Temperature_Sandburg, Inversion_temperature, Humidity

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

# 定义control function
ga_ctrl <- gafsControl(functions = rfGA,  # another option is `caretGA`.
                        method = "cv",
                        repeats = 3)

# 遗传算法特征选择
set.seed(100)
ga_obj <- gafs(x=trainData[, c(1:3, 5:13)], 
               y=trainData[, 4], 
               iters = 3,   # normally much higher (100+)
               gafsControl = ga_ctrl)

ga_obj
Genetic Algorithm Feature Selection

366 samples
12 predictors

Maximum generations: 3 
Population per generation: 50 
Crossover probability: 0.8 
Mutation probability: 0.1 
Elitism: 0 

Internal performance values: RMSE, Rsquared
Subset selection driven to minimize internal RMSE 

External performance values: RMSE, Rsquared, MAE
Best iteration chose by minimizing external RMSE 
External resampling method: Cross-Validated (10 fold) 

During resampling:
  * the top 5 selected variables (out of a possible 12):
    Month (100%), Pressure_gradient (100%), Temperature_ElMonte (100%), Visibility (90%), Inversion_temperature (80%)
  * on average, 7.5 variables were selected (min = 5, max = 10)

In the final search using the entire training set:
   * 6 features selected at iteration 3 including:
     Month, Day_of_month, Wind_speed, Temperature_ElMonte, Pressure_gradient ... 
   * external performance at this iteration is

       RMSE    Rsquared         MAE 
     3.6605      0.7901      2.8010 
# Optimal variables
ga_obj$optVariables

1. ‘Month’
2. ‘Day_of_month’
3. ‘Wind_speed’
4. ‘Temperature_ElMonte’
5. ‘Pressure_gradient’
6. ‘Visibility’

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

# Define control function
sa_ctrl <- safsControl(functions = rfSA,
                        method = "repeatedcv",
                        repeats = 3,
                        improve = 5) # n iterations without improvement before a reset

# Genetic Algorithm feature selection
set.seed(100)
sa_obj <- safs(x=trainData[, c(1:3, 5:13)], 
               y=trainData[, 4],
               safsControl = sa_ctrl)

sa_obj
Simulated Annealing Feature Selection

366 samples
12 predictors

Maximum search iterations: 10 
Restart after 5 iterations without improvement (0.2 restarts on average)

Internal performance values: RMSE, Rsquared
Subset selection driven to minimize internal RMSE 

External performance values: RMSE, Rsquared, MAE
Best iteration chose by minimizing external RMSE 
External resampling method: Cross-Validated (10 fold, repeated 3 times) 

During resampling:
  * the top 5 selected variables (out of a possible 12):
    Temperature_ElMonte (73.3%), Inversion_temperature (63.3%), Month (60%), Day_of_week (50%), Inversion_base_height (50%)
  * on average, 6 variables were selected (min = 3, max = 8)

In the final search using the entire training set:
   * 6 features selected at iteration 10 including:
     Month, Day_of_month, Day_of_week, Wind_speed, Temperature_ElMonte ... 
   * external performance at this iteration is

       RMSE    Rsquared         MAE 
     4.0574      0.7382      3.0727 
# Optimal variables
print(sa_obj$optVariables)
[1] "Month"               "Day_of_month"        "Day_of_week"        
[4] "Wind_speed"          "Temperature_ElMonte" "Visibility"         

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

library(InformationValue)
inputData <- read.csv("http://rstatistics.net/wp-content/uploads/2015/09/adult.csv")
print(head(inputData))
  AGE         WORKCLASS FNLWGT  EDUCATION EDUCATIONNUM       MARITALSTATUS
1  39         State-gov  77516  Bachelors           13       Never-married
2  50  Self-emp-not-inc  83311  Bachelors           13  Married-civ-spouse
3  38           Private 215646    HS-grad            9            Divorced
4  53           Private 234721       11th            7  Married-civ-spouse
5  28           Private 338409  Bachelors           13  Married-civ-spouse
6  37           Private 284582    Masters           14  Married-civ-spouse

          OCCUPATION   RELATIONSHIP   RACE     SEX CAPITALGAIN CAPITALLOSS
1       Adm-clerical  Not-in-family  White    Male        2174           0
2    Exec-managerial        Husband  White    Male           0           0
3  Handlers-cleaners  Not-in-family  White    Male           0           0
4  Handlers-cleaners        Husband  Black    Male           0           0
5     Prof-specialty           Wife  Black  Female           0           0
6    Exec-managerial           Wife  White  Female           0           0

  HOURSPERWEEK  NATIVECOUNTRY ABOVE50K
1           40  United-States        0
2           13  United-States        0
3           40  United-States        0
4           40  United-States        0
5           40           Cuba        0
6           40  United-States        0

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

# 选择分类变量，用来计算信息值Info Value.
cat_vars <- c ("WORKCLASS", "EDUCATION", "MARITALSTATUS", "OCCUPATION", "RELATIONSHIP", "RACE", "SEX", "NATIVECOUNTRY")  # get all categorical variables

# Init Output
df_iv <- data.frame(VARS=cat_vars, IV=numeric(length(cat_vars)), STRENGTH=character(length(cat_vars)), stringsAsFactors = F)  # init output dataframe

# 得到每个变量的信息值
for (factor_var in factor_vars){
  df_iv[df_iv$VARS == factor_var, "IV"] <- InformationValue::IV(X=inputData[, factor_var], Y=inputData$ABOVE50K)
  df_iv[df_iv$VARS == factor_var, "STRENGTH"] <- attr(InformationValue::IV(X=inputData[, factor_var], Y=inputData$ABOVE50K), "howgood")
}

# Sort
df_iv <- df_iv[order(-df_iv$IV), ]

df_iv

1

WOETable(X=inputData[, 'WORKCLASS'], Y=inputData$ABOVE50K)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

library(randomForest)
library(DALEX)

# Load data
inputData <- read.csv("http://rstatistics.net/wp-content/uploads/2015/09/adult.csv")

# Train random forest model
rf_mod <- randomForest(factor(ABOVE50K) ~ ., data=inputData, ntree=100)
rf_mod

# Variable importance with DALEX
explained_rf <- explain(rf_mod, data=inputData, y=inputData$ABOVE50K)

# Get the variable importances
varimps = variable_dropout(explained_rf, type='raw')

print(varimps)
Call:
 randomForest(formula = factor(ABOVE50K) ~ ., data = inputData,      ntree = 100) 
               Type of random forest: classification
                     Number of trees: 100
No. of variables tried at each split: 3

        OOB estimate of  error rate: 17.4%
Confusion matrix:
      0    1 class.error
0 24600  120 0.004854369
1  5547 2294 0.707435276


        variable dropout_loss        label
1   _full_model_          852 randomForest
2   EDUCATIONNUM          842 randomForest
3      EDUCATION          843 randomForest
4  MARITALSTATUS          844 randomForest
5         FNLWGT          845 randomForest
6     OCCUPATION          847 randomForest
7            SEX          847 randomForest
8    CAPITALLOSS          847 randomForest
9   HOURSPERWEEK          847 randomForest
10           AGE          848 randomForest
11          RACE          848 randomForest
12     WORKCLASS          849 randomForest
13  RELATIONSHIP          850 randomForest
14 NATIVECOUNTRY          853 randomForest
15      ABOVE50K          853 randomForest
16   CAPITALGAIN          893 randomForest
17    _baseline_          975 randomForest
plot(varimps)