R talks to Weka about Data Mining

R provides us with excellent resources to mine data, and there are some good overviews out there:

• Yanchang’s website with Examples and a nice reference card
• The rattle-package that introduces a nice GUI for R, and Graham William’s compendium of tools
• The caret-package that offers a unified interface to running a multitude of model builders. And there are other tools out there for data mining, like Weka.

Weka has a GUI and can be directed via the command line with Java as well, and Weka has a large variety of algorithms included. If, for whatever reason, you do not find the algorithm you need being implemented in R, Weka might be the place to go. And the RWeka-package marries R and Weka.

I am not an expert neither in R, nor in Weka, nor in data mining. But I happen to play around with them, and I’d like to share a starter on how to work with them. There is good documentation out there (e.g. Open-Source Machine Learning: R Meets Weka or RWeka Odds and Ends), but sometimes you want to document your own steps and ways of working, and this is what I do.

So, I want to build a classification model for the iris-dataset, based on a tree classifier. Joice is the C4.5 algorithm that I did not find implemented in any standard R package (anybody can help me out?).

We want to predict the class of a flower based on their attributes, namely sepal and petal width and length. The three species we have are “setosa”, “versicolor” and “virginica”. A short summary is given above.



Prediction with J48 (aka C4.5)
We next load the RWeka package.

summary(iris)

## Sepal.Length Sepal.Width Petal.Length Petal.Width 
## Min. :4.30 Min. :2.00 Min. :1.00 Min. :0.1 
## 1st Qu.:5.10 1st Qu.:2.80 1st Qu.:1.60 1st Qu.:0.3 
## Median :5.80 Median :3.00 Median :4.35 Median :1.3 
## Mean :5.84 Mean :3.06 Mean :3.76 Mean :1.2 
## 3rd Qu.:6.40 3rd Qu.:3.30 3rd Qu.:5.10 3rd Qu.:1.8 
## Max. :7.90 Max. :4.40 Max. :6.90 Max. :2.5 
## Species 
## setosa :50 
## versicolor:50 
## virginica :50 
## 
## 
## 

library(RWeka)
We now build the classifier, and this works with the J48(.)-function:



iris_j48 <- J48(Species ~ ., data = iris)
iris_j48

## J48 pruned tree
## ------------------
## 
## Petal.Width <= 0.6: setosa (50.0)
## Petal.Width > 0.6
## | Petal.Width <= 1.7
## | | Petal.Length <= 4.9: versicolor (48.0/1.0)
## | | Petal.Length > 4.9
## | | | Petal.Width <= 1.5: virginica (3.0)
## | | | Petal.Width > 1.5: versicolor (3.0/1.0)
## | Petal.Width > 1.7: virginica (46.0/1.0)
## 
## Number of Leaves : 5
## 
## Size of the tree : 9

summary(iris_j48)

## 
## === Summary ===
## 
## Correctly Classified Instances 147 98 %
## Incorrectly Classified Instances 3 2 %
## Kappa statistic 0.97 
## Mean absolute error 0.0233
## Root mean squared error 0.108 
## Relative absolute error 5.2482 %
## Root relative squared error 22.9089 %
## Coverage of cases (0.95 level) 98.6667 %
## Mean rel. region size (0.95 level) 34 %
## Total Number of Instances 150 
## 
## === Confusion Matrix ===
## 
## a b c <-- classified as
## 50 0 0 | a = setosa
## 0 49 1 | b = versicolor
## 0 2 48 | c = virginica

plot(iris_j48)
plot of chunk build_J48

We can assign the model to an object, and printing the object gives us the tree in “Weka-Output”, summary(.) gives us the Summary of the classification on the training set (again, in Weka-style), and plot(.) allows us to nicely plot it.

Evaluation in Weka
Well, we used the whole dataset now for training, but we actually might want to perform cross-validation. This can be done like this:

eval_j48 <- evaluate_Weka_classifier(iris_j48, numFolds = 10, complexity = FALSE,
seed = 1, class = TRUE)
eval_j48

## === 10 Fold Cross Validation ===
## 
## === Summary ===
## 
## Correctly Classified Instances 144 96 %
## Incorrectly Classified Instances 6 4 %
## Kappa statistic 0.94 
## Mean absolute error 0.035 
## Root mean squared error 0.1586
## Relative absolute error 7.8705 %
## Root relative squared error 33.6353 %
## Coverage of cases (0.95 level) 96.6667 %
## Mean rel. region size (0.95 level) 33.7778 %
## Total Number of Instances 150 
## 
## === Detailed Accuracy By Class ===
## 
## TP Rate FP Rate Precision Recall F-Measure MCC ROC Area PRC Area Class
## 0.980 0.000 1.000 0.980 0.990 0.985 0.990 0.987 setosa
## 0.940 0.030 0.940 0.940 0.940 0.910 0.952 0.880 versicolor
## 0.960 0.030 0.941 0.960 0.950 0.925 0.961 0.905 virginica
## Weighted Avg. 0.960 0.020 0.960 0.960 0.960 0.940 0.968 0.924 
## 
## === Confusion Matrix ===
## 
## a b c <-- classified as
## 49 1 0 | a = setosa
## 0 47 3 | b = versicolor
## 0 2 48 | c = virginica
We see slightly worse results now, as you would suspect.

Using Weka-controls
We used the standard options for th J48 classifier, but Weka allows more. You can acces these with the WOW-function:

WOW("J48")

## -U Use unpruned tree.
## -O Do not collapse tree.
## -C <pruning confidence>
## Set confidence threshold for pruning. (default 0.25)
## Number of arguments: 1.
## -M <minimum number of instances>
## Set minimum number of instances per leaf. (default 2)
## Number of arguments: 1.
## -R Use reduced error pruning.
## -N <number of folds>
## Set number of folds for reduced error pruning. One fold is
## used as pruning set. (default 3)
## Number of arguments: 1.
## -B Use binary splits only.
## -S Don't perform subtree raising.
## -L Do not clean up after the tree has been built.
## -A Laplace smoothing for predicted probabilities.
## -J Do not use MDL correction for info gain on numeric
## attributes.
## -Q <seed>
## Seed for random data shuffling (default 1).
## Number of arguments: 1.
If, for example, we want to use a tree with minimum 10 instances in each leaf, we change the command as follows:

j48_control <- J48(Species ~ ., data = iris, control = Weka_control(M = 10))
j48_control

## J48 pruned tree
## ------------------
## 
## Petal.Width <= 0.6: setosa (50.0)
## Petal.Width > 0.6
## | Petal.Width <= 1.7: versicolor (54.0/5.0)
## | Petal.Width > 1.7: virginica (46.0/1.0)
## 
## Number of Leaves : 3
## 
## Size of the tree : 5
And you see the tree is different (well, it just does not go as deep as the other one..).

Building cost-sensitive classifiers
You might want to include a cost matrix, i.e you want to penalize some wrong classifications, see here. If you think classifying for example a versicolor wrongly is very harmful, you want to penalize such a classification in our example, you can do that easily – you just have to choose a different classifier, namely the “Cost-sensitive classifier” in Weka:

csc <- CostSensitiveClassifier(Species ~ ., data = iris, control = Weka_control(`cost-matrix` = matrix(c(0,
10, 0, 0, 0, 0, 0, 10, 0), ncol = 3), W = "weka.classifiers.trees.J48",
M = TRUE))
But you have to tell the “cost-sensitive-classifier” that you want to use J48 as algorithm, and you have to tell him the cost matrix you want to apply, name ly the matrix of the form

matrix(c(0, 1, 0, 0, 0, 0, 0, 1, 0), ncol = 3)

## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 1 0 1
## [3,] 0 0 0
where you penalize “versicolor” being falsly classified as one of the others by factor 10.

And again we evaluate on 10-fold CV:

eval_csc <- evaluate_Weka_classifier(csc, numFolds = 10, complexity = FALSE,
seed = 1, class = TRUE)
eval_csc

## === 10 Fold Cross Validation ===
## 
## === Summary ===
## 
## Correctly Classified Instances 98 65.3333 %
## Incorrectly Classified Instances 52 34.6667 %
## Kappa statistic 0.48 
## Mean absolute error 0.2311
## Root mean squared error 0.4807
## Relative absolute error 52 %
## Root relative squared error 101.9804 %
## Coverage of cases (0.95 level) 65.3333 %
## Mean rel. region size (0.95 level) 33.3333 %
## Total Number of Instances 150 
## 
## === Detailed Accuracy By Class ===
## 
## TP Rate FP Rate Precision Recall F-Measure MCC ROC Area PRC Area Class
## 0.980 0.070 0.875 0.980 0.925 0.887 0.955 0.864 setosa
## 0.980 0.450 0.521 0.980 0.681 0.517 0.765 0.518 versicolor
## 0.000 0.000 0.000 0.000 0.000 0.000 0.500 0.333 virginica
## Weighted Avg. 0.653 0.173 0.465 0.653 0.535 0.468 0.740 0.572 
## 
## === Confusion Matrix ===
## 
## a b c <-- classified as
## 49 1 0 | a = setosa
## 1 49 0 | b = versicolor
## 6 44 0 | c = virginica

and we see that the “versicolors” are now better predicted (only one wrong, compared to 3 in the normal J48 earlier). But this happened at the expense of more fals classification on “virginica”, where we have now 6 wrongly classified instead of 2.

Creating publication quality graphics using R

As part of a one-day workshop, I have developped an online tutorial on how to create publication quality graphics using R (from an academic point of view).

The tutorial can be found here

http://teachpress.environmentalinformatics-marburg.de/2013/07/creating-publication-quality-graphs-in-r-7/

As mentioned in the tutorial, feel free to send me any feedback, criticism, general comments or bug reports.

Enjoy,

Tim

Btw, the entire tutorial was created using Rmarkdown and knitr. The .Rmd file can be found here

https://github.com/tim-salabim/metvurst/blob/master/markdown/20130617_data_vis_workshop.Rmd

Select operations on R data frames

The R language is weird - particularly for those coming from a typical programmer's background, which likely includes OO languages in the curly-brace family and relational databases using SQL. A key data structure in R, the data.frame, is used something like a table in a relational database. In terms of R's somewhat byzantine type system (which is explained nicely here), a data.frame is a list of vectors of varying types. Each vector is a column in the data.frame making this a column-oriented data structure as opposed to the row-oriented nature of relational databases.

In spite of this difference, we often want to do the same sorts of things to an R data.frame that we would to a SQL table. The R docs confuse the SQL-savvy by using different terminology, so here is a quick crib-sheet for applying SQL concepts to data.frames.

We're going to use a sample data.frame with the following configuration of columns, or schema, if you prefer: (sequence:factor, strand:factor, start:integer, end:integer, common_name:character, value:double) where the type character is a string and a factor is something like an enum. Well, more accurately, value is a vector of type double and so forth. Anyway, our example is motivated by annotation of genome sequences, but the techniques aren't particular to any type of data.

> head(df)
sequence strand start end common_name value
1 chromosome + 1450 2112 yvrO 0.9542516
2 chromosome + 41063 41716 graD6 0.2374012
3 chromosome + 62927 63640 graD3 1.0454790
4 chromosome + 63881 64807 gmd 1.4383845
5 chromosome + 71811 72701 moaE -1.8739953
6 chromosome + 73639 74739 moaA 1.2711058
So, given a data.frame of that schema, how do we do some simple select operations?

Selecting columns by name is easy:

> df[,c('sequence','start','end')]
sequence start end
1 chromosome 1450 2112
2 chromosome 41063 41716
3 chromosome 62927 63640
4 chromosome 63881 64807
5 chromosome 71811 72701
...
As is selecting row names, or both:

> df[566:570,c('sequence','start','end')]
sequence start end
566 chromosome 480999 479860
567 chromosome 481397 480999
568 chromosome 503053 501275
569 chromosome 506476 505712
570 chromosome 515461 514277
Selecting rows that meet certain criteria is a lot like a SQL where clause:

> df[df$value>3.0,]
sequence strand start end common_name value
199 chromosome + 907743 909506 hutU 3.158821
321 chromosome + 1391811 1393337 nadB 3.092771
556 chromosome - 431600 431037 apt 3.043373
572 chromosome - 519043 518186 hbd1 3.077040
For extra bonus points, let's find tRNAs.

> df[grep("trna", df$common_name, ignore.case=T),]
sequence strand start end common_name value
18 chromosome + 115152 115224 Asn tRNA -0.461038128
19 chromosome + 115314 115422 Ile tRNA -0.925268307
31 chromosome + 167315 167388 Tyr tRNA 0.112527023
32 chromosome + 191112 191196 Ser tRNA 0.986357577
...
Duplicate row names

Row names are not necessarily unique in R, which breaks the method shown above for selecting by row name. Take matrix a:

< a = matrix(1:18, nrow=6, ncol=3)
< rownames(a) <- c('a', 'a', 'a', 'b', 'b', 'b')
< colnames(a) <- c('foo', 'bar', 'bat')
< a
foo bar bat
a 1 7 13
a 2 8 14
a 3 9 15
b 4 10 16
b 5 11 17
b 6 12 18
It looks to me like trying to index by the row names just returns the first row of a given name:

< a['a',]
foo bar bat 
1 7 13
< a['b',]
foo bar bat 
4 10 16 
But this works:

< a[rownames(a)=='a',]
foo bar bat
a 1 7 13
a 2 8 14
a 3 9 15

'More Resources:

• R programming for those coming from other languages
• Extract or Replace Parts of a Data Frame
• Subsetting Vectors, Matrices and Data Frames
• R Help search engine
• Abandon hope all ye who enter The R Inferno
• R Programming for Bioinformatics

Help for R, the R language, or the R project is notoriously hard to search for, so I like to stick in a few extra keywords, like R, data frames, data.frames, select subset, subsetting, selecting rows from a data.frame that meet certain criteria, and find.

Linear Mixed Models in R

Linear Mixed Models in R

Introduction:

A brief history of time
At the beginning (1992-1995) I would use SAS (first proc glm, later proc mixed), but things started getting painfully slow and limiting if one wanted to move into animal model BLUP. At that time (1995-1996) I moved to DFREML (by Karen Meyer, now replaced by WOMBAT) and AIREML (by Dave Johnson, now defunct—I mean the program), which were designed for the analysis of animal breeding progeny trials, so it was a hassle to deal with experimental design features. At the end of 1996 (or was it the beginning of 1997) I started playing with ASReml (programed by Arthur Gilmour mostly based on theoretical work by Robin Thompson and Brian Cullis). I was still using SAS for data preparation, but all my analyses went through ASReml (for which I wrote the cookbook), which was orders of magnitude faster than SAS (and could deal with much bigger problems). Around 1999, I started playing with R (prompted by a suggestion from Rod Ball), but I didn’t really use R/S+ often enough until 2003. At the end of 2005 I started using OS X and quickly realized that using a virtual machine or dual booting was not really worth it, so I dropped SAS and totally relied on R in 2009.

Options

As for many other problems, there are several packages in R that let you deal with linear mixed models from a frequentist (REML) point of view. I will only mention nlme (Non-Linear Mixed Effects), lme4 (Linear Mixed Effects) and asreml (average spatial reml). There are also several options for Bayesian approaches, but that will be another post.

nlme is the most mature one and comes by default with any R installation. In addition to fitting hierarchical generalized linear mixed models it also allows fitting non-linear ones. Its main advantages are, in my humble opinion, the ability to fit fairly complex hierarchical models using linear or non-linear approaches, a good variety of variance and correlation structures, and access to several distributions and link functions for generalized models. In my opinion, its main drawbacks are i- fitting cross-classified random factors is a pain, ii- it can be slow and may struggle with lots of data, iii- it does not deal with pedigrees by default and iv- it does not deal with multivariate data.

lme4 is a project led by Douglas Bates (one of the co-authors of nlme), looking at modernizing the code and making room for trying new ideas. On the positive side, it seems to be a bit faster than nlme and it deals a lot better with cross-classified random factors. Drawbacks: similar to nlme’s, but dropping point i- and adding that it doesn’t deal with covariance and correlation structures yet. It is possible to fit pedigrees using the mmpedigree package, but I find the combination a bit flimsy.

ASReml-R is, unsurprisingly, an R package interface to ASReml. On the plus side it i- deals well with cross-classified random effects, ii- copes very well with pedigrees, iii- can work with fairly large datasets, iv-can run multivariate analyses and v- covers a large number of covariance and correlation structures. Main drawbacks are i- limited functionality for non-Gaussian responses, ii- it does not cover non-linear models and iii- it is non-free (as in beer an speech). The last drawback is relative; it is possible to freely use asreml for academic purposes (and there is also a version for developing countries). Besides researchers, the main users of ASReml/asreml-r are breeding companies.

All these three packages are available in Windows, Linux and OS X.

A  simple example:

I will use a traditional dataset to show examples of the notation for the three packages: Yates’ variety and nitrogen split-plot experiment. We can get the dataset from the MASS package, after which it is a good idea to rename the variables using meaningful names. In addition, I will follow Bill Venables’s excellent advice and create additional variables for main plot and subplots, as it is confusing to use the same factor for two purposes (e.g. variety as treatment and main plot). Incidentally, if you haven’t read Bill’s post go and read it; it is one of the best explanations I have ever seen for a split-plot analysis.

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

library(MASS) data(oats) names(oats) = c('block', 'variety', 'nitrogen', 'yield') oats$mainplot = oats$variety oats$subplot = oats$nitrogen summary(oats)  block           variety     nitrogen      yield              mainplot  I  :12   Golden.rain:24   0.0cwt:18   Min.   : 53.0   Golden.rain:24  II :12   Marvellous :24   0.2cwt:18   1st Qu.: 86.0   Marvellous :24  III:12   Victory    :24   0.4cwt:18   Median :102.5   Victory    :24  IV :12                    0.6cwt:18   Mean   :104.0  V  :12                                3rd Qu.:121.2  VI :12                                Max.   :174.0    subplot  0.0cwt:18  0.2cwt:18  0.4cwt:18  0.6cwt:18

The nlme code for this analysis is fairly simple: response on the left-hand side of the tilde, followed by the fixed effects (variety, nitrogen and their interaction). Then there is the specification of the random effects (which also uses a tilde) and the data set containing all the data. Notice that 1|block/mainplot is fitting block and mainplot within block. There is no reference to subplot as there is a single assessment for each subplot, which ends up being used at the residual level.

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

library(nlme) m1.nlme = lme(yield ~ variety*nitrogen,                       random = ~ 1|block/mainplot,                       data = oats) summary(m1.nlme) Linear mixed-effects model fit by REML  Data: oats        AIC      BIC    logLik   559.0285 590.4437 -264.5143 Random effects:  Formula: ~1 | block         (Intercept) StdDev:    14.64496  Formula: ~1 | mainplot %in% block         (Intercept) Residual StdDev:    10.29863 13.30727 Fixed effects: yield ~ variety * nitrogen                                     Value Std.Error DF   t-value p-value (Intercept)                      80.00000  9.106958 45  8.784492  0.0000 varietyMarvellous                 6.66667  9.715028 10  0.686222  0.5082 varietyVictory                   -8.50000  9.715028 10 -0.874933  0.4021 nitrogen0.2cwt                   18.50000  7.682957 45  2.407927  0.0202 nitrogen0.4cwt                   34.66667  7.682957 45  4.512152  0.0000 nitrogen0.6cwt                   44.83333  7.682957 45  5.835427  0.0000 varietyMarvellous:nitrogen0.2cwt  3.33333 10.865342 45  0.306786  0.7604 varietyVictory:nitrogen0.2cwt    -0.33333 10.865342 45 -0.030679  0.9757 varietyMarvellous:nitrogen0.4cwt -4.16667 10.865342 45 -0.383482  0.7032 varietyVictory:nitrogen0.4cwt     4.66667 10.865342 45  0.429500  0.6696 varietyMarvellous:nitrogen0.6cwt -4.66667 10.865342 45 -0.429500  0.6696 varietyVictory:nitrogen0.6cwt     2.16667 10.865342 45  0.199411  0.8428 anova(m1.nlme)                  numDF denDF   F-value p-value (Intercept)          1    45 245.14299  <.0001 variety              2    10   1.48534  0.2724 nitrogen             3    45  37.68562  <.0001 variety:nitrogen     6    45   0.30282  0.9322

The syntax for lme4 is not that dissimilar, with random effects specified using a (1|something here) syntax. One difference between the two packages is that nlme reports standard deviations instead of variances for the random effects.

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

library(lme4) m1.lme4 = lmer(yield ~ variety*nitrogen + (1|block/mainplot),                        data = oats) summary(m1.lme4) Linear mixed model fit by REML Formula: yield ~ variety * nitrogen + (1 | block/mainplot)    Data: oats  AIC   BIC logLik deviance REMLdev  559 593.2 -264.5    595.9     529 Random effects:  Groups         Name        Variance Std.Dev.  mainplot:block (Intercept) 106.06   10.299  block          (Intercept) 214.48   14.645  Residual                   177.08   13.307 Number of obs: 72, groups: mainplot:block, 18; block, 6 Fixed effects:                                  Estimate Std. Error t value (Intercept)                       80.0000     9.1064   8.785 varietyMarvellous                  6.6667     9.7150   0.686 varietyVictory                    -8.5000     9.7150  -0.875 nitrogen0.2cwt                    18.5000     7.6830   2.408 nitrogen0.4cwt                    34.6667     7.6830   4.512 nitrogen0.6cwt                    44.8333     7.6830   5.835 varietyMarvellous:nitrogen0.2cwt   3.3333    10.8653   0.307 varietyVictory:nitrogen0.2cwt     -0.3333    10.8653  -0.031 varietyMarvellous:nitrogen0.4cwt  -4.1667    10.8653  -0.383 varietyVictory:nitrogen0.4cwt      4.6667    10.8653   0.430 varietyMarvellous:nitrogen0.6cwt  -4.6667    10.8653  -0.430 varietyVictory:nitrogen0.6cwt      2.1667    10.8653   0.199 anova(m1.lme4) Analysis of Variance Table                  Df  Sum Sq Mean Sq F value variety           2   526.1   263.0  1.4853 nitrogen          3 20020.5  6673.5 37.6856 variety:nitrogen  6   321.7    53.6  0.3028

For this type of problem, the notation for asreml is also very similar, particularly when compared to nlme.

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(asreml) m1.asreml = asreml(yield ~ variety*nitrogen,                            random = ~ block/mainplot,                            data = oats) summary(m1.asreml)$varcomp                              gamma component std.error  z.ratio constraint block!block.var          1.2111647  214.4771 168.83404 1.270343   Positive block:mainplot!block.var 0.5989373  106.0618  67.87553 1.562593   Positive R!variance               1.0000000  177.0833  37.33244 4.743416   Positive wald(m1.asreml, denDF = "algebraic") $Wald                  Df denDF    F.inc           Pr (Intercept)       1     5 245.1000 1.931825e-05 variety           2    10   1.4850 2.723869e-01 nitrogen          3    45  37.6900 2.457710e-12 variety:nitrogen  6    45   0.3028 9.321988e-01 $stratumVariances                df  Variance block block:mainplot R!variance block           5 3175.0556    12              4          1 block:mainplot 10  601.3306     0              4          1 R!variance     45  177.0833     0              0          1

In this simple example one pretty much gets the same results, independently of the package used (which is certainly comforting). I will soon cover another simple model, but with much larger dataset, to highlight some performance differences between the packages.