Showing posts with label machine learning. Show all posts
Showing posts with label machine learning. Show all posts

Sunday, February 28, 2021

Dealing with skewed predictor variables in prediction problems

Some datasets may have (heavily) skewed predictors. This may increase bias as well as variance of models assuming linear associations between predictor and response. Predictors can then be scaled using non-linear transformations in order to have more linear associations between predictors and response. Here Andrew Gelman discusses some potential advantages of log-transforming predictors when associations between predictor and response are non-linear. 

Many decision-tree methods, however, are insensitive to the scaling of the predictors. They thus do not require you to pre-process your predictors. Some suggestions for methods that are invariant with respect to monotonous transformations (e.g., log transform, squaring) of the predictors: 

  •  Single decision tree created using the CART or MOB algorithms, which are e.g. implemented in the rpart and partykit R packages. The advantage of the latter is that it does not have a preference towards splitting on variables with a higher number of possible split points. Package partykit also implements the conditional inference tree algorithm, which is a good method, but it employs linear association tests for variable selection, which are sensitive to the scaling of the predictors. 
  • Decision-tree ensemblse based on the previously mentioned tree algorithms, as implemented in the randomForest and mobForest packages. Or a CART-based gradient boosting decision-tree ensemble based on the previously mentioned tree akgorithms, as implemented in package gbm. 
  •  Prediction rule ensemble based on MOB or CART. This will strike a balance between the ease of interpretability of a single tree, and the higher predictive accuracy of a decision-tree ensemble. This method is implemented in package pre. Specify use.grad=FALSE in order to derive the rules using MOB, specify tree.unbiased=FALSE to derive the rules using CART (otherwise, default is to use conditional inference trees).

Sunday, September 25, 2016

Effect of step size in forward stagewise regression

Forward stagewise regression is a linear model selection algorithm. It is a modification of least angle regression (LARS; Efron, Hastie, Johnstone & Tibshirani, 2002). It aims to select only a subset of X (predictor) variables, for efficient prediction of a Y (response) variable.

The forward stagewise regression  (FSR) algorithm roughly works as follows:
  1. Initilize by setting i = 0, coefficient vector b_i = 0, and the residual vector rY_i = Y - X * b_i.
  2. Increase (or decrease*) the coefficient of the predictor variable most highly correlated with rY_i by a small step, yielding the current coefficient vector b_{i+1}.
  3. Eat potato chips (optional)
  4. Calculate rY_{i+1} = Y - X * b_{i+1}.
  5. Set i = i + 1, repeat steps 2 - 4 until convergence (i.e., consecutive b_i no longer change).
In step 2, the coefficient is increased or decreased, depending on the sign of the correlation between the predictor variable and rY_i. The step size, that is, the size of the increase (or decrease) in step 2, can be provided by the user, or the optimal value can be determined by cross validation. Step 3 is optional, and should only be performed on a small subset of the iterations.

Although FSR can be performed with the R package lars, I  wanted to program it myself, to check out the effect of different step sizes, and to see how the results compare to ordinary least squares regression (OLS).  I created a function for performing FSR (called swReg), a function for plotting the coefficient paths (called plot.swReg), and applied it on the Boston housing dataset. Here are the code and the results:


R functions for performing FSR, plotting coefficient paths and computing cross-validated prediction error
swReg <- function(X, Y, stepsize = .1, threshold = .1) {
  sX <- scale(X)
  sY <- scale(Y)
  b <- vector(length = ncol(X))
  rY <- sY
  sgn <- 0
  lastsgn <- 0
  index <- 0
  lastindex <- 0
  iteration <- 0
  path <- list()
  while(abs(max(t(rY) %*% sX)) > threshold & (index != lastindex | sgn == lastsgn)) {
    iteration <- iteration + 1
    lastindex <- index
    lastsgn <- sgn
    index <- which(abs(t(rY) %*% sX) == max(abs(t(rY) %*% sX)))
    sgn <- sign((t(rY) %*% sX)[index])
    b[index] <- b[index] + sgn*stepsize
    rY <- sY - sX %*% b
    path[[iteration]] <- b
  }
  bUnstand <- (b/attr(sX, "scaled:scale"))*attr(sY, "scaled:scale") 
  intercept <- attr(sY, "scaled:center") - bUnstand %*% attr(sX, "scaled:center")
  bUnstand <- c(intercept, bUnstand)
  names(bUnstand) <- c("(Intercept)", colnames(X))
  
  return(list(unstandardized.coef = bUnstand, standardized.coef = b,
            iteration = iteration, stepsize = stepsize, threshold = threshold, 
            coef.path = data.frame(matrix(
              unlist(path), ncol = ncol(X), nrow = iteration, byrow = TRUE, 
              dimnames = list(1:iteration, colnames(X)))),
            data = list(X = X, Y = Y))
         )  
}

plot.swReg <- function(object, legend = TRUE) {
  plot(object$coef.path[,1], type = "l", 
  ylim = c(min(object$coef.path), max(object$coef.path)), 
  ylab = "coefficient", xlab = "iteration", 
  main = "Standardized coefficient paths", 
  sub = paste("stepsize =", object$stepsize))
  for (i in 2:ncol(X)) {lines(object$coef.path[,i], col = i)}
  if(legend) {
    legend("topleft", legend = colnames(object$data$X), cex = .5, 
    lty = 1, col = 1:ncol(object$coef.path), y.intersp = .25, 
    x.intersp = .25, bty = "n", seg.len = .5)
  }
}

predict.swReg <- function(object, newdata = object$data$X) {
  cbind(rep(1, times = nrow(newdata)), as.matrix(newdata)) %*% object$unstandardized.coef
}

swReg.xval <- function(object, k = 10, seed = 42) {
  set.seed(seed)
  ids <- peperr::resample.indices(n = nrow(object$data$X), sample.n = 10, method = "cv")
  models <- list()
  xvaldatasets <- list()
  for (i in 1:k){
    xvaldatasets[[i]] <- list(train = list(X = object$data$X[ids$sample.index[[i]],],
                                           Y = object$data$Y[ids$sample.index[[i]]]),
                              test = list(X = object$data$X[ids$not.in.sample[[i]],],
                                          Y = object$data$Y[ids$not.in.sample[[i]]]))
    models[[i]] <- swReg(xvaldatasets[[i]]$train$X, xvaldatasets[[i]]$train$Y, 
                         stepsize = object$stepsize, threshold = object$threshold)
    xvaldatasets[[i]]$test$xvalpredY <- predict.swReg(models[[i]], newdata = xvaldatasets[[i]]$test$X)
  }
  dataset <- cbind(object$data$X, fold = NA, Y = object$data$Y, Ycvpred = NA)
  coefficients <- list()
  for(i in 1:k){
    dataset[,"fold"][ids$not.in.sample[[i]]] <- i
    dataset[,"Ycvpred"][ids$not.in.sample[[i]]] <- xvaldatasets[[i]]$test$xvalpredY 
    coefficients[[i]] <- list(unstandardized.coef = models[[i]]$unstandardized.coef, 
                              standardized.coef = models[[i]]$standardized.coef)
  }
  return(list(dataset = dataset, cv.coefficients = coefficients))
} 

Example: Boston housing data
## Data preparation:
library(MASS)
X <- as.matrix(Boston[,-14])
Y <- Boston$medv

## Run forward stagewise regression:
tmp1 <- swReg(X, Y, st = .2, thr = .2)
tmp2 <- swReg(X, Y, st = .1, thr = .1)
tmp3 <- swReg(X, Y, st = .01, thr = .01)
tmp4 <- swReg(X, Y, st = .001, thr = .001)

## Plot coefficient paths:
par(mfrow=c(2,2))
plot.swReg(tmp1)
plot.swReg(tmp2)
plot.swReg(tmp3)
plot.swReg(tmp4)


Coefficient paths


What we see is that the paths actually look very similar, but less iterations are required before convergence with larger step sizes. Also, larger steps sizes provide sparser solutions (i.e., less non-zero coefficients in the final solution). Note that the algorithm is applied on standardized X and Y values, so the step size could be interpreted as the minimal correlation between a predictor and the response variable that the user finds 'relevant'. The FSR function above provides both standardized and unstandardized regression coefficients.


Comparison with OLS


What we see is: the larger the step size, the sparser the final solution. With step size equal to .2, only 5 predictor variables are selected. With step size equal to .001, all variables are selected. Also, the smaller the step size, the more the final solution resembles the OLS solution. The absolute distances between the OLS and FSR solutions decrease with step size:







References 

Efron, B., Hastie, T., Johnstone, I., & Tibshirani, R. (2004). Least angle regression. The Annals of Statistics, 32(2), 407-499. link to pdf

Tuesday, April 21, 2015

Subsampling for cross validation in R

If you want to obtain cross validated estimates of performance in R, you have to split your dataset into a number of more or less equally sized parts. For leave-on-out cross validation (LOO CV), this is easy (see example below). However, LOO CV generally provides overly optimistic estimates of performance, and a number of folds < no. of observations (e.g., 5 or 10) is preferred. This may be require some custom code, if the number of observations is not a multiple of 5 or 10. Luckily, the resample.indices() function of the R-package peperr provides a generic solution to the problem. Below, examples are provided for a linear regression model, but the point is that lm() can be replaced by any model building function.


R code example for k-fold cross validation:
library(peperr)
# simulate original dataset with 10 variables and 500 observations:
dataset <- data.frame(matrix(rnorm(5000), ncol=10, nrow=500))
# set number of folds:
k <- 10
# get observation ids for in folds:
ids <- resample.indices(n=nrow(dataset), sample.n=k, method="cv")
# create a list for gathering results:
xvaldatasets <- list()
# for each fold k:
for (i in 1:k){
  # make datasets of training and test observations:
  xvaldatasets[[i]] <- list(train=dataset[ids$sample.index[[i]],],
                            test=dataset[ids$not.in.sample[[i]],])
  # fit model on training observations, and use it to make predictions for test observations:
  xvaldatasets[[i]]$test$xvalpred <- predict(lm(X10 ~ ., data=xvaldatasets[[i]]$train),
                                          newdata=xvaldatasets[[i]]$test)
}

# evaluate predictive accuracy in each fold k:
xvalcors <- list()
for (i in 1:k){
  # for example, by calculating the correlation between the predicted and observed values of X10:
  xvalcors[[i]] <- cor(xvaldatasets[[i]]$test$X10, xvaldatasets[[i]]$test$xvalpred)
}
# print correlations between true and predicted values for each fold:
unlist(xvalcors)
# calculate mean:
mean(unlist(xvalcors))
# calculate SD:
sd(unlist(xvalcors))




R code example for leave-one-out cross validation:
# simulate original dataset with 10 variables and 500 observations:
dataset <- data.frame(matrix(rnorm(5000), ncol=10, nrow=500))
# create objects for gathering results:
xvaldatasets <- list()
xvalpreds <- rep(NA,  times=ncol(dataset))
# for each fold k:
for (i in 1:nrow(dataset)){
# make datasets of training and test observations:
xvaldatasets[[i]] <- list(train=dataset[-i,],                           
                          test=dataset[i,])
# fit model on training observations, and use it to make predictions for test observations:
xvalpreds[i] <- predict(lm(X10 ~ ., data=xvaldatasets[[i]]$train),
                        newdata=xvaldatasets[[i]]$test)
}
# evaluate predictive accuracy:
cor(xvalpreds, dataset$X10)

Thursday, December 19, 2013

Tikz: simple decision tree (or CART tree) example



I created the above diagram of a decision tree using the tikz package for LaTeX. It took me a while to find a simple, clear way to produce a decision tree. This is the code I used:

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{positioning}
\newdimen\nodeDist
\nodeDist=35mm
\begin{document}

\begin{tikzpicture}[
    node/.style={%
      draw,
      rectangle,
    },
  ]

    \node [node] (A) {IDS $>$ 13.5?};
    \path (A) ++(-135:\nodeDist) node [node] (B) {exit 1};
    \path (A) ++(-45:\nodeDist) node [node] (C) {LCIanx $>$ 0.3604?};
    \path (C) ++(-135:\nodeDist) node [node] (D) {exit 2};
    \path (C) ++(-45:\nodeDist) node [node] (E) {exit 3};

    \draw (A) -- (B) node [left,pos=0.25] {no}(A);
    \draw (A) -- (C) node [right,pos=0.25] {yes}(A);
    \draw (C) -- (D) node [left,pos=0.25] {no}(A);
    \draw (C) -- (E) node [right,pos=0.25] {yes}(A);
\end{tikzpicture}

\end{document}

Tuesday, July 9, 2013

Error due to partykit

I'm a frequent user of the R package party. Recently, the R package partykit was created by the authors of party, for 'representing, summarizing, and visualizing tree-structured regression and classification models'. And all of a sudden, the where function didn't function anymore, and the following error message gets printed:

Error in function (classes, fdef, mtable)  :
  unable to find an inherited method for function "where", for signature "constparty"


This is due to changes in the structure of the BinaryTree Class in partykit package. See following code for fixing this problem.

library(party)
library(partykit)
airq <- subset(airquality, !is.na(Ozone))
airct <- ctree(Ozone ~ ., data = airq)
where(airct)
# Error!
# Fix error:
detach("package:partykit", unload=TRUE)
airct <- ctree(Ozone ~ ., data = airq) # you have to rebuild the tree
where(airct) # now it works!

Tuesday, June 12, 2012

Notes on RuleFit

Edit (dec 24 2017): I have developed an R package for deriving prediction rule ensembles: pre (available from CRAN). It provides most of the functionality of the Rulefit package, with some improvements/adjustments. The main differences with Rulefit are: 1) pre derives prediction rules using unbiased tree growing algorithms from package partykit, instead of CART trees, which have biased variable selection; 2) pre is completely R-based, so easier to use the rules, coefficients, importance etc. for further computations (but also a bit slower); 3) Rulefit allowed for a bagging and/or boosting approach to generating the initial ensemble of prediction rules, pre allows for bagging and/or boosting and/or random forest approaches.

 

Original post:

RuleFit provides an ensemble of regression functions and CART tree derived prediction rules. It's one of the very few rule prediction rule ensemble methods that can be used for both classification and regression (c.f., C4.5). The paper describing the methodology and program can be found here (Friedman, J. H., & Popescu, B. E. (2008). Predictive learning via rule ensembles. The Annals of Applied Statistics, 916-954). The TMVA (Toolkit for Multivariate Analysis) package in R provides an implementation of rulefit as well, but only for classification.

  

Arguments of the rulefit function

sparse relates to \lambda in Friedman & Popescu (2005)
memory.par is \nu in Friedman & Popescu (2005)
samp.fract is \eta in Friedman & Popescu (2005)
tree.size is \overline{L} in Friedman & Popescu (2005)
inter.supp is \kappa in Friedman & Popescu (2005)

The penalty function in the Friedman & Popescu (2005) paper shows only the lasso penalty (equation 4), but elastic net penalty is implemented in RuleFit, as well. If you set the sparse argument of the rulefit function to low values, (e.g., sparse=.0000001), the model will be built using the ridge penalty, resulting in a much larger number of terms in the model. 

samp.fract=min(1,(11*sqrt(neff)+1)/neff), by default. When all observation weights are equal, neff = n for regression; so samp.fract=min(1,(11*sqrt(n)+1)/n). By default, samp.fract becomes smaller when n increases. It becomes less than 1 about when neff>122.
 
test.frac should be a value between 0.1 and 0.5, and is 0.2 by default. Test sample is used for what some would refer to as validation sample: determining the optimal value of \lambda, the penalty parameter.

If test.reps is set to a value >0, then [value]-fold cross validation is used for determining the optimal value of \lambda. By default, test.reps=round(min(20,max(0.0,5200/neff-2))): it is at most 20 and at least 0. By default, the value of test.resps decreases when neff (n) increases. By default, test.reps < 20 when neff > 237, and test.reps < 4 when neff < 866.  

Using the max.trms argument is very useful in controlling the size of the final ensemble. However, the maximum number specified is an approximate maximum number of terms: in most cases, the ensemble will be somewhat bigger, especially when test.reps > 1.

The sparse argument selects the sparse regression method used determining the weights of prediction functions in the final ensemble (thus, it also selects the prediction functions: prediction functions with zero weight are not included in the final ensemble). For ensembles with human-interpretable sizes, use sparse=3: it uses forward stepwise (regression) or stagewise (classification) regression (see Hastie & Taylor, 2007).

Hastie, T., Taylor, J., Tibshirani, R., & Walther, G. (2007). Forward stagewise regression and the monotone lasso. Electronic Journal of Statistics, 1, 1-29.


Other functions

The rfxeval function uses all data, to estimate the expected extra-sample error Err, the average generalization error when the method f^{hat}(X) is applied to an independent test sample (see paragraph 7.10 from Hastie, Tibshirani and Friedman, 2008).



Troubleshooting

Do not use rfxeval() before rfpred() or rules(): it seems to change the current model in the RuleFit home directory. If rfxval() is used, rebuild the rulefit model using rulefit(x,y), before using rules() and rfpred(). I'm not sure this only involves the rules and rfpred functions.

You can obtain a list (instead of a plot) of the variable importances by using the following code:
imp <- varimp(plot=F); imp

The rules() function provides a list of prediction functions, printed in a command prompt window. This is not very convenient for editing, copying and pasting. However, when the rules() function is used, a file named "rulesout.hlp" is created in the working directory, which can be opened using any text editor. Be aware that R and Rulefit will hang, when the final ensemble consists of < 10 prediction functions. I think this is due to the defaults for the rules function (rules(begin=1, end=begin+9). There's three ways to solve this:
1) open windows task manager, and end the process called 'rf_go.exe'; find the file 'rules.out' in the working directory and open it in a text editor.
2) open windows task manager, and end the process called 'rf_go.exe'; make sure the R working directory is the dierctory from which rulefit is run; type 'readLines("rulesout.hlp")' in R.
3) type "rules(1, [some number])", where [some number] is the number of terms in the current RuleFit model.


Lambda value

According to Hastie, Tibshirani and Friedman (2009), cross validation should be used to estimate the shrinkage, or smoothing, parameters. To use cross-validation for determining the lambda parameter (i.o.w., the size of the final model) in the rulefit function, use the test.reps argument. It provides the shrinkage parameter used for model selection, but this is done implicitly, and currently there's no way to obtain the actual lambda value used.

Hastie, T., Tibshirani, R., Friedman, J., Hastie, T., Friedman, J., & Tibshirani, R. (2009). The elements of statistical learning (Vol. 2, No. 1). New York: Springer.


Range of prediction rules: missings

According to the rulefit help file, missing predictor variable values (NAs) are coded 9.0e30 (observations with missing output variable values are of no use to the algorithm, and are excluded from the analysis). In the description of the prediction rules, you may come across something like:

                                                                        
"Rule   1:     2  variables"                                               
"     support =  0.3496      coeff =   1.538      importance =   100.00          "     V1:  range = -0.9900E+36   1.500"                             
"     V2:  range =   1.500      0.9000E+31"                        
"     V3:  range = 0.9000E+31   0.9900+36" 
 

This rule would not apply for observations with missing values for V2 (range does not include values equal to or higher than 0.9e31), but does apply for observations with missing values for V3 (range includes values equal to or higher than 0.9e31).