10 R packages every data scientist should know about

The yhat blog lists 10 R packages they wish they'd known about earlier. Drew Conway calls them "10 reasons to always start your analysis in R". They're all very useful R packages that every data scientist should be aware of. They are:

1. sqldf (for selecting from data frames using SQL)
2. forecast (for easy forecasting of time series)
3. plyr (data aggregation)
4. stringr (string manipulation)
5. Database connection packages RPostgreSQL, RMYSQL, RMongo, RODBC, RSQLite
6. lubridate (time and date manipulation)
7. ggplot2 (data visulization)
8. qcc (statistical quality control and QC charts)
9. reshape2 (data restructuring)
10. randomForest (random forest predictive models)

You can find links to all of these packages and tips on how to use them at link below.

yhat blog: 10 R packages I wish I knew about earlier

Don’t use correlation to track prediction performance

From: http://www.r-bloggers.com/dont-use-correlation-to-track-prediction-performance/?utm_source=feedburner&utm_medium=email&utm_campaign=Feed%3A+RBloggers+%28R+bloggers%29

Using correlation to track model performance is “a mistake that nobody would ever make” combined with a vague “what would be wrong if I did do that” feeling. I hope after reading this feel a least a small urge to double check your work and presentations to make sure you have not reported correlation where R-squared, likelihood or root mean square error (RMSE) would have been more appropriate.

It is tempting (but wrong) to use correlation to track the performance of model predictions. The temptation arises because we often (correctly) use correlation to evaluate possible model inputs. And the correlation function is often a convenient built-in function.

In fact on the data used to train a model one can prove correlation squared is equal R-squared under very mild assumptions. See Nina Zumel’s “Correlation and R-squared” article for a very good explanation of why this is the case. In fact correlation squared is nearly equal to R-squared for any data set that is exchangeable with the training data. So we should expect correlation squared nearly equal to R-squared on properly prepared test data. And this is the core of the paradox: correlation is perfectly good measure on training and test data, but possibly not on data later encountered in a production environment. We don’t know that data later seen in production is in fact exchangeable with training and test data, that is something we hope for and want to track. We don’t want symmetries in the correlation function hiding a divergence such as the an unexpected changing of units or scale in production data.

The correlation function (which we will call cor(,)) has a huge number of obscuring symmetries: it is unchanged under positive scaling, shifts and the swap of its two arguments. This means it is in fact scoring if some ideal shift plus re-scaling of your model predictions is performing well instead of scoring the predictions you are using. And this is not what you want, models in a production environment are supposed to make actual good predictions. Measurements in production are supposed to tell you if the model or data have drifted (not to merely assume they have not).

Here is some R-code showing symmetries in cor(,):

> y = runif(10)
> x = y + 0.5*runif(10)
> cor(x,y)
[1] 0.8893743
> cor(y,x)
[1] 0.8893743
> cor(10*x,y)
[1] 0.8893743
> cor(x+10,y)
[1] 0.8893743

R-squared (written as a function as rsq(,) has none of these symmetries, it changes under simple alterations of its arguments and can become arbitrarily negative.

Here is some R-code showing the lack of symmetries in rsq(,):

> rsq = function(y,f) { 1 - sum((y-f)^2)/sum((y-mean(y))^2) }
> rsq(x,y)
[1] -0.4966555
> rsq(y,x)
[1] 0.09424879
> rsq(10*x,y)
[1] -9.197255
> rsq(x+10,y)
[1] -2250.407

And here is some R-code to remind you that correlation squared and R-squared do agree on training data:

> model = lm(y~x)
> rsq(y,predict(model))
[1] 0.7909866
> cor(y,predict(model))^2
[1] 0.7909866
> model

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
    -0.3309       1.1432  

If you look at this with an open or learning mind it should seem very strange that a function like cor(,) with a huge number of symmetries is closely associated with a function like rsq(,) with many fewer symmetries. At this point we re-recommend Nina Zumel’s “Correlation and R-squared” article to remind ourselves why correlation squared and R-squared are the same on training data. But they point we want to leave with is that the correlation function is using its many symmetries to evaluate if some simple function of a value is a good prediction (hence correlation is a great way to vet possible model inputs), and correlation is not scoring if the unaltered predictions at hand actually are in fact good.

Related posts:

1. Correlation and R-Squared
2. Level fit summaries can be tricky in R
3. Modeling Trick: Masked Variables

Building a package in RStudio is actually very easy

source: http://www.r-bloggers.com/building-a-package-in-rstudio-is-actually-very-easy/?utm_source=feedburner&utm_medium=email&utm_campaign=Feed%3A+RBloggers+%28R+bloggers%29

So, you’ve written some code and you use it routinely. Now you’ve written some code and you’d like to use version control to ensure that development continues in a robust fashion. You do that and you use Github or something so that not only are changes tracked, but the general public receives the benefit of your knowledge. At the same time, you receive the benefit of their editing and ideas, but in a secure way. Now what? All that useful, well-functioning code needs a home, some documentation and a little dignity. It needs to become a package.

This is actually very easy to do. I had a few false starts, but now that I know what’s going on, I can create an R package faster than you can say Hadley Wickham.

1. Start with a Github repository. This will have a set of R sourcecode files, which collectively do something interesting.
2. Within RStudio, create a new project, using your local Github repository as the project’s source directory.
3. Within the “Build Tools” section of the “Project Options” make sure that you’ve selected “Package” from the “Project Build Tools” dropdown box.
4. In that same dialog box, enter the “–no-examples” flag of the “Check Package” build options. (At this stage, you probably don’t have any examples composed. The default assumes that you do. I don’t know why this is.)
5. If RStudio, didn’t create one, create a subdirectory in your project folder called “R”. Move all your code here. Github will reflect the deletion and new file creation.
6. Create a DESCRIPTION file. There are numerous sources which explain how to do this.
7. Make sure your code doesn’t have any errors in it. Easiest way to do that is to source all of them.
8. Build the package.

And that’s it. If you press CTRL+SHIFT+B (in Windows, at any rate) RStudio will build the package and load it for you. All of your functions are now in memory, but won’t appear in the Workspace pane.

RStudio has fine- if brief- documentation on their site, with links to more detailed guidance. The most useful of these to me, so far, was the wiki for Wickham’s devtools package.

Of course, that’s only about half the story. You still have to write documentation for all your functions, provide a demo file, compose a vignette, etc. Then, you may want to submit to CRAN. For now, I’m just using Github, which will be fine for the forseeable future.

The MRMR project may now be loaded directly into R via Github. The documentation is non-existent for now. It’s a very bare bones application, so could be useful as a template for others. Meanwhile, Github has loads of good examples of R packages.

Want to load MRMR? Instructions below will make that happen.

library(devtools)
install_github(repo = "MRMR", username = "PirateGrunt")
library(MRMR)
df = GetNAICData()
head(df)

Armadillo eigenvalues

source: http://gallery.rcpp.org/articles/armadillo-eigenvalues/

Today a (slightly confused) question on StackOverflow wondered how to access R’s facilities for eigenvalues calculations from C code.

For this, we need to step back and consider how this is done. In fact, R farms the calculation out to the BLAS. On could possibly access R’s functions—but would then have to wrestle with the data input/output issues which make Rcpp shine in comparison. Also, Rcpp gets us access to Armadillo (via the RcppArmadillo) package and Armadillo’s main focus are exactly the linear algebra calculations and decompositions.

And with facilities that were added to Rcpp in the 0.10.* release series, this effectively becomes a one-liner of code! (Nitpickers will note that there are also one include statement, two attributes declarations and the function name itself.)

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::vec getEigenValues(arma::mat M) {
    return arma::eig_sym(M);
}

We can illustrate this easily via a random sample matrix.

set.seed(42)
X <- matrix(rnorm(4*4), 4, 4)
Z <- X %*% t(X)

getEigenValues(Z)

        [,1]
[1,]  0.3319
[2,]  1.6856
[3,]  2.4099
[4,] 14.2100

In comparison, R gets the same results (in reverse order) and also returns the eigenvectors.

eigen(Z)

$values
[1] 14.2100  2.4099  1.6856  0.3319

$vectors
         [,1]     [,2]    [,3]     [,4]
[1,]  0.69988 -0.55799  0.4458 -0.00627
[2,] -0.06833 -0.08433  0.0157  0.99397
[3,]  0.44100 -0.15334 -0.8838  0.03127
[4,]  0.55769  0.81118  0.1413  0.10493

Armadillo has other eigenvector computations too, see its documentation.

Using Eigen for eigenvalues

source: http://gallery.rcpp.org/articles/eigen-eigenvalues/

A previous post showed how to compute eigenvalues using the Armadillo library via RcppArmadillo.

Here, we do the same using Eigen and the RcppEigen package.

#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

using Eigen::Map;                // 'maps' rather than copies 
using Eigen::MatrixXd;                  // variable size matrix, double precision
using Eigen::VectorXd;                  // variable size vector, double precision
using Eigen::SelfAdjointEigenSolver;    // one of the eigenvalue solvers

// [[Rcpp::export]]
VectorXd getEigenValues(Map<MatrixXd> M) {
    SelfAdjointEigenSolver<MatrixXd> es(M);
    return es.eigenvalues();
}

We can illustrate this easily via a random sample matrix.

set.seed(42)
X <- matrix(rnorm(4*4), 4, 4)
Z <- X %*% t(X)

getEigenValues(Z)

[1]  0.3319  1.6856  2.4099 14.2100

In comparison, R gets the same results (in reverse order) and also returns the eigenvectors.

eigen(Z)

$values
[1] 14.2100  2.4099  1.6856  0.3319

$vectors
         [,1]     [,2]    [,3]     [,4]
[1,]  0.69988 -0.55799  0.4458 -0.00627
[2,] -0.06833 -0.08433  0.0157  0.99397
[3,]  0.44100 -0.15334 -0.8838  0.03127
[4,]  0.55769  0.81118  0.1413  0.10493

Eigen has other a lot of other decompositions, see its documentation for more details.

Using the GSL to compute eigenvalues

source: http://www.r-bloggers.com/using-the-gsl-to-compute-eigenvalues/?utm_source=feedburner&utm_medium=email&utm_campaign=Feed%3A+RBloggers+%28R+bloggers%29

Two posts showed how to compute eigenvalues using Armadillo andusing Eigen. As we also looked at using the
GNU GSL, this post will show how to conpute eigenvalues using GSL.

As mentioned in the previous GSL post, we instantiate C language pointers suitable for GSL (here the matrix M). Those must be freed manually, as shown before the return statement.

// [[Rcpp::depends(RcppGSL)]]

#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>

// [[Rcpp::export]]
Rcpp::NumericVector getEigenValues(Rcpp::NumericMatrix sM) {

    RcppGSL::matrix<double> M(sM);  // create gsl data structures from SEXP
    int k = M.ncol();
    Rcpp::NumericVector N(k);   // to store results 

    gsl_vector *eigval = gsl_vector_alloc(k);
    gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(k);
    gsl_eigen_symm (M, eigval, w);
    gsl_eigen_symm_free (w);

    for (int j = 0; j < k; j++) {
        N[j] = gsl_vector_get(eigval, j);
    }
    M.free();                          // important: GSL wrappers use C structure
    return N;    // return vector  
}

We can illustrate this easily via a random sample matrix.

set.seed(42)
X <- matrix(rnorm(4*4), 4, 4)
Z <- X %*% t(X)

getEigenValues(Z)

[1] 14.2100  2.4099  1.6856  0.3319

In comparison, R gets the same results (in reverse order) and also returns the eigenvectors.

eigen(Z)

$values
[1] 14.2100  2.4099  1.6856  0.3319

$vectors
         [,1]     [,2]    [,3]     [,4]
[1,]  0.69988 -0.55799  0.4458 -0.00627
[2,] -0.06833 -0.08433  0.0157  0.99397
[3,]  0.44100 -0.15334 -0.8838  0.03127
[4,]  0.55769  0.81118  0.1413  0.10493

Safely Loading Packages in R

source: http://www.r-bloggers.com/safely-loading-packages-in-r/?utm_source=feedburner&utm_medium=email&utm_campaign=Feed%3A+RBloggers+%28R+bloggers%29

Using R snippets written by other developers can be unendingly maddening.  There are a variety of reasons for this, most of which boil down to a simple issue: most code is written such that a system must be configured in precisely the same way as the code’s author’s machine.  Anyone who’s ever seen a line like this:
read.xls("C:/Users/MCaine/code/R/projecteuler/someotherdirectory/data.xls")
knows what I am talking about. To use this without modification, you must:

1. Use Windows.
2. Have exactly the directory structure specified by the address (which is highly unlikely, unless you were the one who wrote it).
3. Have the gdata package installed and included in the project (which is both unlikely, and difficult to know without already being a regular user of the package).

You can see how it would already be easier to just change the address to whatever works on your machine.  For this problem, I’m afraid I have no easy solution.  My preferred approach is to provide URLs, so that the directory structure doesn’t depend on the user’s machine, but this obviously provides its own host of problems.  However, I can take aim at this third issue.

Package management in R is a silly thing, both because it is so easy and it is so easy to screw up.  Most people who write R write it like analysts: they write exactly enough code the get the desired output on their machine, and leave it at that.  When such R Code is made public, it tends to be very difficult to use it in actual replication.  But there are some simple ways around that.

getPackage <- function(pkg){
  if(!require(pkg, character.only=TRUE)){
    install.packages(pkg)
    library(pkg, character.only=TRUE)
  }
  return(TRUE)
}

Consider this gist.  Instead of calling library(package), which fails if the library is not installed on the user’s machine, getPackage(package) invokes that safer require function, which doesn’t fail when the local machine doesn’t have the requisite package.  Instead, it returns false, which triggers the R script to download the package from the user’s default CRAN mirror and then bring it into the user’s working session.  If anything goes wrong, then it throws an error, but it won’t do that for anything silly like the user not having a package that is a mere two commands away.

One note: because of the esoteric manner in which R treats package names, you must pass this function a string and not a package name.  If you’re not on top of your type-know-how, this means that getPackage(plyr) will fail.  You should instead writegetPackage("plyr").

Now, I’m sure there’s a very good reason hidden deep in R’s core that this is a bad way to do things, but it has saved me time and headaches.  R is recognized as a difficult language both within the programming world (for its strange inconsistencies and its non-grown-up hacker culture), and outside of it (because programming is hard).  I wish that more R functions were written in the defensive way to decrease the cognitive barriers to using R for the latter group.