The gepR package

The gepR package is an experimental data mining tool. It uses the Gene Expression Programming (GEP) technique to find functional relationships between the output variable and input variables. The computing code is implemented in C with parallel support enabled by OpenMP.

library(gepR)

Dataset

Some example datasets are loaded with the package. For now, using the heating dataset, we will fit a function between the Heating Load of a house and various characteristics of the house. There are 768 rows and 8 independent variables in th data set. The first column Heating_Load is the responses variable. Here is how the heating data looks like:

head(heating)
  Heating_Load Relative_Compactness Surface_Area Wall_Area Roof_Area
1        15.55                 0.98        514.5     294.0    110.25
2        15.55                 0.98        514.5     294.0    110.25
3        15.55                 0.98        514.5     294.0    110.25
4        15.55                 0.98        514.5     294.0    110.25
5        20.84                 0.90        563.5     318.5    122.50
6        21.46                 0.90        563.5     318.5    122.50
  Overall_Height Orientation Glazing_Area Glazing_Area_Distribution
1              7           2            0                         0
2              7           3            0                         0
3              7           4            0                         0
4              7           5            0                         0
5              7           2            0                         0
6              7           3            0                         0

Split training and validation

First split the datasets into two parts for training and validation.

train_idx <- sample(1:nrow(heating), round(0.5*nrow(heating)))
valid_idx <- setdiff(1:nrow(heating), train_idx)

train_y = heating[train_idx, 1]
train_x = heating[train_idx, 2:ncol(heating)]
valid_y = heating[valid_idx, 1]
valid_x = heating[valid_idx, 2:ncol(heating)]

Fit a function using GEP

We will first load the gepR.dll using gep_load_dll(). When no parameter is supplied (as in this case), it looks for the gepR.dll in the default location. Alternatively, a path and file name can be supplied as a charater string in the parameter of the function call. Next, we run the gep_train() function by specifying the data for dependent (y) and independent (x) variables. We then tell the algorithm to run for 1000 iterations per round (it’ll run 3 rounds by default as specified by the maxpass parameter), the goal is to achieve a 0.95 R-square on the training data, and let it run on 10 threads (by default nthreads is 4). There are many more parameters to customize the algorithm’s behavior but we will accept the default for now. For more information, here are the docs for gep_train() and gep_score(). 

gep_load_dll()  # Load the gepR.dll from the default location
gepmod <- gep_train(y = train_y, x = train_x, maxiter = 1000, goal=0.95, nthreads = 10)

    Starting round 1 of 3
    Iter 1: e	 R-square = 0.787
    Iter 41: (log(cbrt(c))*cbrt(log(e)))	 R-square = 0.831
    Iter 44: (log(cbrt(c))*cbrt(round(e)))	 R-square = 0.834
    Iter 172: (cbrt(cos(round(d)))-g)	 R-square = 0.845
    Iter 174: (cbrt(g)-cos(round(d)))	 R-square = 0.857
    Iter 494: (e*((c+h)*exp(g)))	 R-square = 0.916
    Iter 500: (e*((c+e)*exp(g)))	 R-square = 0.916
    Iter 506: (e*((c*e)*exp(g)))	 R-square = 0.918
    Cannot meet the R-square target of 0.950 in 1000 iterations.
    New variable i = 0.001151 * *e**Eceg + (7.686344) added.
    Starting round 2 of 3
    Iter 1: i	 R-square = 0.918
    Iter 46: (i+abs(log(abs(d))))	 R-square = 0.919
    Iter 48: (i+abs(cbrt(d)))	 R-square = 0.920
    Iter 52: (i*abs(cbrt(d)))	 R-square = 0.925
    Iter 64: (sin(abs(g))+sqrt(abs(i)))	 R-square = 0.929
    Iter 728: (sqrt(i)+sqrt(g))	 R-square = 0.933
    Iter 735: (sqrt(i)+cbrt(g))	 R-square = 0.935
    Iter 797: (sqrt(i)+sin(cbrt(g)))	 R-square = 0.936
    Cannot meet the R-square target of 0.950 in 1000 iterations.
    New variable j = 9.228363 * +QSiDg + (-25.186019) added.
    Starting round 3 of 3
    Iter 1: round(j)	 R-square = 0.939
    Iter 96: round(round((g+(j+a))))	 R-square = 0.939
    Iter 164: round((log(cos(a))+round(j)))	 R-square = 0.941
    Iter 218: (round(round(j))+cos(sin(c)))	 R-square = 0.942
    Iter 226: (round(round(j))+cos(sqrt(c)))	 R-square = 0.945
    Iter 777: (round(j)+cos((c*d)))	 R-square = 0.946
    Iter 786: (round(j)+cos((c*a)))	 R-square = 0.951
    Stopping criteria met at iteration 786, R-square = 0.951
    (round(j)+cos((c*a)))
    R-square target reached at round 3. Done.
    GEP model saved to file gep_sol.dat

Predict the new data

We can supply the model and newdata to the gep_score() function for prediction.

predict_y <- gep_score(valid_x, gepmod)
(RMSE <- sqrt(sum((predict_y - valid_y)^2)/length(valid_idx)))
[1] 2.382726
plot(predict_y, valid_y)

Multiple Linear Regression (as a comparison)

GEP fits a composite linear-nonlinear regression model for the data. For comparison, let’s fit a linear regression model using the lm() function.

linmod <- lm(Heating_Load ~ ., data = heating, subset = train_idx)
predict_y_lm <- predict(linmod, newdata = valid_x)
Warning in predict.lm(linmod, newdata = valid_x): prediction from a rank-
deficient fit may be misleading
(RMSE <- sqrt(sum((predict_y_lm - valid_y)^2)/length(valid_idx)))
[1] 2.876745
plot(predict_y_lm, valid_y)

It is good practice to unload the DLL by the end of the session, though it is OK not to do so.

gep_unload_dll()

End Note

This is an experimental package and it is not publicly available yet. The main reason is that I have not figured out a general, cross-platform installation method for an average user. It needs some customized compiling and linking operations to make the DLL. I believe there is a way to make it easier, but have not had time to dig it out. Contact me for more information: yanchaoliu@wayne.edu