The Tiger Mom and A Clockwork Orange

True disciple is doing what you want.

A wise friend once told me that. Amy Chua, better known as the Tiger Mother, wrote about discipline (from a different point of view) in Why Chinese Mothers Are Superior.

What Chinese parents understand is that nothing is fun until you're good at it. To get good at anything you have to work, and children on their own never want to work, which is why it is crucial to override their preferences. [...] Tenacious practice, practice, practice is crucial for excellence; rote repetition is underrated in America. Once a child starts to excel at something -- whether it's math, piano, pitching or ballet -- he or she gets praise, admiration and satisfaction. This builds confidence and makes the once not-fun activity fun. This in turn makes it easier for the parent to get the child to work even more.

For what it's worth, Chua's book is apparently less strident and more nuanced than the WSJ article. Anyway, like her methods or not, I have a lot of sympathy for a parent trying to teach her kids about delayed gratification, that you can do difficult things if you try, and that hard work pays off.

If it's true that mastering a complex skill takes 10,000 hours of practice, then the persistence to push through those hours is a fairly important lesson to learn early. Recent research has caused a reappraisal in how much talent arrises from innate genius versus how much is the product of effort, practice and persistence.

What brought to mind my old friend's remark about disciple was Paul Buchheit's take: motivation can be either intrinsic or extrinsic. Amy Chua is teaching her kids to be extrinsically motivated, to respond to the praise and admiration of others. You do it because you are told to. You put your energy into chasing the approval of external authorities. In contrast, he describes intrinsic motivation like this:

To the greatest extent possible, do whatever is most fun, interesting, and personally rewarding (and not evil).

Follow your heart, as hippy moms tell their children. Buccheit says, "I'm kind of lazy, or maybe I lack will power or discipline or something. Either way, it's very difficult for me to do anything that I don't feel like doing." Sounds familiar. "The intrinsic path to success is to focus on being the person that you are, and put all of your energy and drive into being the best possible version of yourself."

The difference between intrinsic and extrinsic motivation is easily recognized in the moral dimension. In A Clockwork Orange, Anthony Burgess imagines the transfer of aesthetic sense from creation to violence. Deprived of outlet, creativity turns destructive. The main thrust of the story is an examination of attempts to impose an external morality by force versus growing an internal morality.

Paul Buccheit was the software developer that originated Google's gmail. For myself, and I'm sure lots of others, a key attraction to programming was the ability to create in a powerful medium without asking anyone's permission. The creative freedom, the feeling that the authorities hadn't (yet) figured out how to lock things down was incredibly inspiring.

It's easy to see how that attraction to technology meshes with Daniel Pink's elements of motivation — autonomy, mastery, and purpose. If you've got root and a compiler, you've got autonomy. And it's all about a pissing contest of mastery. (This might partially explain the gender ratio in the field.) And technology is rife with appeals to higher purpose, from the open-source movement to the digital media that helped fuel the revolutions in Tunisia and Egypt.

The Tiger Mom demands mastery before autonomy, leaving purpose firmly in the hands of the parent. Buccheit puts autonomy first, trusting in a natural sense of direction to lead to mastery and purpose. In terms of Maslow's hierarchy of needs, Amy Chua has the "esteem" level covered, but stops short of the top level - intrinsic self-directed creativity.

Some suggest that American society erects a border fence at the entry to the highest level. Well, nobody ever mentions why it's always drawn as a pyramid, but there's probably a reason. Not everyone gets to be at the top. Usually, that's the realm only of the elite.

David Cain of Raptitude writes under the title How to Make Trillions of Dollars, that the fundamentals of being a self-directed person are these:

Creativity. Curiosity. Resilience to distraction. Patience with others.

Cain defines self-reliance as “an unswerving willingness to take responsibility for your life, regardless of who had a hand in making it the way it is”.

Again, we're basically talking about the top levels of the pyramid. But, I like the addition of resilience to distraction. Discipline is a hard sell in a culture the promotes immediate gratification and nonstop indulgence. It's hard to hear your own voice over the clamor of consumer culture and expectations from family, boss and everyone else. Hearing it is nearly impossible while drowning in distractions like twitter and facebook. And here's something else to remember about online amusements:

If you're not paying, you're not the customer; you're the product.

At a recent data mining conference I saw rooms full of marketers ready to slice, dice and mash up your personal data to more precisely target advertising. Resisting this attack is an essential skill of modern life. A healthy cynicism is a necessary defense mechanism. Hearing yourself think is only going to get harder.

The flawed idea implicit in consumer culture is what you consume is what you are. But, valuing consuming over creating or doing is inevitably a dead end. The reason I'm not so hot on the iPad is that it's a device for consuming. The old macs were (marketed as) tools for programmers, graphic artists, musicians and film-makers -- in other works doers, builders and creators.

Purpose has to come from values. The recent travails of the financial sector show what happens when motivations or at least incentives become disconnected from morals and values.

Of course, lots of technology has a purpose no higher than selling golf clubs on the internet. And technology, itself, can be a distraction. It's easy to get caught up in a rat race of the latest whizzy buzzword laden language, tool or application framework dujour. For years, I've had a half-joking theory that the true purpose of the internet is to absorb the excess productivity of mankind.

Creative, conceptual work driven by autonomy, mastery and purpose pursued with uninterrupted concentration. That begins to answer the question, how do we get some motivation, apply it to something good and inspire the same in those around us, especially our ungrateful screaming offspring.

You can argue one way or another about whether a 13 year old has the foresight to be intrinsically motivated. I certainly didn't have the wherewithal at that age to set a long term goal. And pushing intrinsic motivation on your kids sounds to me something like imposing democracy by force. I doesn't make that much sense.

It takes determination to undertake exhausting frustrating efforts whose payoff is distant and uncertain. Most things that are worth doing are hard. The drive and courage to try anyway is what makes real progress possible. That doesn't come easily, and neither does the judgement necessary to gauge what is worth while against the scale of your own values.

From my current position, I don't particularly want to lecture anyone on how to succeed in life. I respond negatively to coercion and can be a champion slacker, both of which were much to the detriment of my academic career. Still, here it is:

• Value creation over consumption.
• Surround yourself with creators.
• Do what you love, do it a lot, and do it hard.

Using R for Introductory Statistics, Chapter 5, Probability Distributions

In Chapter 5 of Using R for Introductory Statistics we get a brief introduction to probability and, as part of that, a few common probability distributions. Specifically, the normal, binomial, exponential and lognormal distributions make an appearance.

For each distribution, R provides four functions whose names start with the letters d, p, q or r followed by the family name of the distribution. For example, rnorm produces random numbers drawn from a normal distribution. The letters stand for:

d	density/mass function
p	probability (cumulative distribution function) P(X <= x)
q	quantiles, given q, the smallest x such that P(X <= x) > q
r	random number generation

Normal

The Gaussian or normal distribution has a prominent place due to the central limit theorem. It is widely used to model natural phenomena like variations in height or weight as well as noise and error. The 68-95-99.7 rule says:

• 68% of the data falls within 1 standard deviation of the mean
• 95% of the data falls within 2 standard deviations of the mean
• 99.7% of the data falls within 3 standard deviations of the mean

To plot a normal distribution, define some points x, and use dnorm to generate the density at those points.

x <- seq(-3,3,0.1)
plot(x=x, y=dnorm(x, mean=0, sd=1), type='l')

Binomial

A Bernoulli trial is an experiment which can have one of two possible outcomes. Independent repeated Bernoulli trials give rise to the Binomial distribution, which is the probability distribution of the number of successes in n independent Bernoulli trials. Although the binomial distribution is discrete, in the limit as n gets larger, it approaches the normal distribution.

x <- seq(0,20,1)
plot(x=x, y=dbinom(x,20,0.5))

Uniform

A Uniform distribution just says that all allowable values are equally likely, which comes up in dice or cards. Uniform distributions come in either continuous or discrete flavors.

Log-normal

The log-normal distribution is a probability distribution of a random variable whose logarithm is normally distributed. If X is a random variable with a normal distribution, then Y = exp(X) has a log-normal distribution. Analogously to the central limit theorem, the product of many independent random variables multiplied together tends toward a lognormal distribution. It can be used to model continuous random quantities whose distribution is skewed and non-negative, for example income or survival.

samples <- rlnorm(100, meanlog=0, sdlog=1)
par(fig=c(0,1,0,0.35))
boxplot(samples, horizontal=T, bty="n", xlab="log-normal distribution")
par(fig=c(0,1,0.25,1), new=T)
s <- seq(0,max(samples),0.1)
d <- dlnorm(s, meanlog=0, sdlog=1)
hist(samples, prob=T, main="", col=gray(0.9), ylim=c(0,max(d)))
lines(density(samples), lty=2)
curve(dlnorm(x, meanlog=0, sdlog=1), lwd=2, add=T)
rug(samples)

Exponential

The exponential distribution is used to model the time interval between successive random events such as time between failures arising from constant failure rates. The following plot is generated by essentially the same code as above.

More on probability distributions

• R's d, p, q and r methods for families of distributions
• An Introduction to R, Chapter 8: Probability distributions
• CRAN Task View: Probability Distributions
• John D Cook's map of the relationships between probability distributions.

More Using R for Introductory Statistics

• Using R for Introductory Statistics, Chapters 1 and 2
• Using R for Introductory Statistics 3.1
• Using R for Introductory Statistics 3.2
• Using R for Introductory Statistics 3.3
• Using R for Introductory Statistics, Chapter 3.4
• Using R for Introductory Statistics, Chapter 4
• Using R for Introductory Statistics, Chapter 4, Model Formulae
• Using R for Introductory Statistics, Chapter 5

Annotated source code

We programmers are told that reading code is a good idea. It may be good for you, but it's hard work. Jeremy Ashkenas has come up with a simple tool that makes it easier: docco. Ashkenas is also behind underscore.js and coffeescript, a dialect of javascript in which docco is written.

Interesting ways to mix prose and code have appealed to me ever since I first discovered Mathematica's live notebook, which lets you author documents that combine executable source code, typeset text and interactive graphics. For those who remember the early 90's chiefly for their potty training, running Mathematica on the Next pizza boxes was like a trip to the future. Combining the quick cycles of a Read-evaluate-print-loop with complete word processing and mathematical typesetting encourages you to keep lovely notes on your thinking and trials and errors.

Along the same lines, there's Sweave for R and sage for Python.

Likewise, one of the great innovations of Java was Javadoc. Javadoc doesn't get nearly enough credit for the success of Java as a language. It made powerful API's like the collections classes a snap and even helped navigate the byzantine complexities of Swing and AWT.

These days, automated documentation is expected for any language. Nice examples are: RubyDoc, scaladoc, Haddock (for Haskell). Doxygen works with a number of languages. Python has pydoc, but in practice seems to rely more on the library reference. Anyway, there are a bunch, and if your favorite language doesn't have one, start coding now.

The grand-daddy of these ideas is Donald Knuth's literate programming.

I believe that the time is ripe for significantly better documentation of programs, and that we can best achieve this by considering programs to be works of literature. Hence, my title: "Literate Programming."

Let us change our traditional attitude to the construction of programs: Instead of imagining that our main task is to instruct a computer what to do, let us concentrate rather on explaining to human beings what we want a computer to do.

The practitioner of literate programming can be regarded as an essayist, whose main concern is with exposition and excellence of style. Such an author, with thesaurus in hand, chooses the names of variables carefully and explains what each variable means. He or she strives for a program that is comprehensible because its concepts have been introduced in an order that is best for human understanding, using a mixture of formal and informal methods that reinforce each other.

Indeed, Ashkenas references Knuth, calling docco "quick-and-dirty, hundred-line-long, literate-programming".

This goodness needs to come to more language. There's a ruby port called rocco by Ryan Tomayko. And for Clojure there's marginalia.

I love the quick-and-dirty aspect and that will be the key to encouraging programmers to do more documentation that looks like this. I hope they build docco, or something like it, into github. Maybe one day there will be a Norton's anthology of annotated source code.

Vaguely related

• literate-programming.rb
• The Semicolon Wars Every programmer knows there is one true programming language. A new one every week
• The Factor language's Visual REPL
• It would be awesome if similar work enables rich interactions in the Clojure REPL
• IPython, a rich Python shell

Using R for Introductory Statistics, Chapter 5

Any good stats book has to cover a bit of basic probability. That's the purpose of Chapter 5 of Using R for Introductory Statistics, starting with a few definitions:

Random variable

A random number drawn from a population. A random variable is a variable for which we define a range of possible values and a probability distribution. The probability distribution specifies the probability that the variable assumes any value in its range.

The range of a discrete random variable is a finite list of values. For any value k in the range, 0 ≤ P(X=k) ≤ 1. The sum over all values k in the range is 1.

R sample function implements these definitions. For example, let's ask R to create a six-sided die for us.

p.die <- rep(1/6,6)
sum(p.die)

Now, let's roll it 10 times.

die <- 1:6
sample(die, size=10, prob=p.die, replace=T)

Now, let's roll 1000 dice and plot the results.

s <- table(sample(die, size=1000, prob=p.die, replace=T))
lbls = sprintf("%0.1f%%", s/sum(s)*100)
barX <- barplot(s, ylim=c(0,200))
text(x=barX, y=s+10, label=lbls)

Expected value

Expected value (or population mean) of a discrete random variable X is the weighted average of the values in the range of X.

In the case of our six-sided die, the expected value is 3.5, computed like so:

sum(die*p.die)

Things change a bit when we move from discrete to continuous random variables. A continuous random variable is described by a probability density function. If f(x) is the probability density of a random variable X, P(X≤b) is the area under f(x) and to the left of b. The total area under f(x) = 1. f(x) ≥ 0 for all possible values of X. P(X>b) = 1 - P(X≤b). The expected value is the balance point where exactly half of the total area under f(x) is to the right and half is to the left.

A random sample is a sequence of independent identically distributed random variables. A value derived from a random sample, such as sample mean, sample standard deviation, etc. is called a statistic. When we compute statistics of samples, our hope is that the sample statistic is not too far off from the equivalent measurement of the whole population.

The interesting thing is that derived statistics are also random variables. If we role our die several times, we have taken a random sample of size n. That sample can be summarized by computing its mean, denoted by X bar.

The sample mean is, itself, a random variable with its own distribution. So, let's take a look at that distribution. Because we'll use it later, let's define a function that generates a bunch of samples of size n and computes their means. It returns a vector of sample means.

generate.sample.means <- function(n) {
  sample.means <- numeric()
  for (i in 1:1000) { 
    sample.means <- append(sample.means, sum(sample(die, size=n, prob=p.die, replace=T))/n)
  }
  return (sample.means)
}

sample.means <- generate.sample.means(100)
plot(density(sample.means), main="Distribution of sample means",xlab="sample mean", col="orange")

Not coincidentally, it fits pretty closely to a normal distribution (dashed blue line). It's mean is about the same as the parent population, namely right around 3.5. The standard deviation of the sample means can be derived by dividing the standard deviation of the parent population by the square root of the sample size: σ / √ n.

Let's compare the mean and sd of our sample with predicted values. We compute the standard deviation of our parent population of die rolls by squaring the deviation from the mean for each possible value and averaging that. Dividing that by the square root of our sample size gives the predicted sd of our sample means, about 0.17, which is about spot on with the actual sd.

> mean(sample.means)
[1] 3.49002
> sd(sample.means)
[1] 0.1704918
> sqrt(sum( (1:6-3.5)^2 ) / 6) / sqrt(100)
[1] 0.1707825

To overlay the normal distribution on the plot above, we used R's dnorm function like this:

x = seq(3,4,0.01)
lines(x=x,y=dnorm(x,mean=3.5,sd=0.1707825), col=rgb(0x33,0x66,0xAA,0x90,maxColorValue=255), type="l", lty=2)

Inspection of the formula for the standard deviation of sample means supports our common sense intuition that a bigger sample will more likely reflect the whole population. In particular, as the size of our sample goes up, our estimated mean is more likely to be closer to the parent population mean. This idea is known as the law of large numbers. We can show that it works by creating similar plots with increasing n.

sample.means <- generate.sample.means(100)
plot(density(sample.means), main="Distribution of sample means", xlab="sample mean", col="yellow", xlim=c(3.2,3.8), ylim=c(0,8))
sample.means <- generate.sample.means(500)
lines(density(sample.means), col="orange")
sample.means <- generate.sample.means(1000)
lines(density(sample.means), col="red")
legend(3.6,7,c("n=100","n=500","n=1000"), fill=c("yellow", "orange", "red"))

What we've just discovered is the central limit theorem, which states that for any parent population with mean μ and standard deviation σ, the sampling distribution for large n is a normal distribution with mean μ and standard deviation σ / √ n. Putting that in terms of the standard normal distribution Z gives:

It's a little surprising that the normal distribution arises out of a means of samples from a discrete uniform distribution. More surprisingly, samples from any parent distribution give rise to the normal distribution in exactly the same way. Next, we'll look at several widely used families of distributions and R's functions for working with them.

Notes

The graph above shaded in orange serving as an example probability density function is produced with the following R code:

plot(x=c(seq(0,5,0.01)), y=c(dlnorm(x=seq(0,5,0.01),meanlog=0,sdlog=1)), type="l", xlab='', ylab='', yaxt="n", xaxt="n", bty="n")
polygon(x=c(seq(0,2,0.01),2), y=c(dlnorm(x=seq(0,2,0.01),meanlog=0,sdlog=1),0), col="orange")
mtext('b', at=c(2), side=1)
text(0.6,0.2,"P(X≤b)")
abline(c(0,0),c(0,5))

More Using R for Introductory Statistics

• Using R for Introductory Statistics, Chapters 1 and 2
• Using R for Introductory Statistics 3.1
• Using R for Introductory Statistics 3.2
• Using R for Introductory Statistics 3.3
• Using R for Introductory Statistics, Chapter 3.4
• Using R for Introductory Statistics, Chapter 4
• Using R for Introductory Statistics, Chapter 4, Model Formulae

Using R for Introductory Statistics, Chapter 4, Model Formulae

Several R functions take model formulae as parameters. Model formulae are symbolic expressions. They define a relationship between variables rather than an arithmetic expression to be evaluated immediately. Model formulae are defined with the tilde operator. A simple model formula looks like this:

response ~ predictor

Functions that accept formulae typically also take a data argument to specify a data frame in which to look up model variables and a subset argument to select certain rows in the data frame.

We've already seen model formula used for simple linear regression and with plot and boxplot, to show that American cars are heavy gas guzzlers. Two common uses of formula are:

• y ~ x where x and y are numeric
• x ~ f where x is numeric and f is a factor

The Lattice graphics package can accept more complicated model formulas of this form:

response ~ predictor | condition

We'll try this out with a dataset called kid.weights from the UsingR package. We get age, weight, height and gender for 250 kids ranging from 3 month to 12 years old.

library(UsingR)
library(lattice)
dim(kid.weights)
[1] 250   4

We expect weight and height to be related, but we're wondering if this relationship changes over time as kids grow. Often, when we want to condition on a quantitative variable (like age), we turn it into a categorical variable by binning. Here, we'll create 4 bins by taking age in 3 year intervals.

age.classes = cut(kid.weights$age/12, 3*(0:4))
unique(age.classes)
[1] (3,6]  (6,9]  (9,12] (0,3] 
Levels: (0,3] (3,6] (6,9] (9,12]

With age as a factor, we can express our question as the model formula:

height ~ weight | age.classes

The lattice graphics function xyplot accepts this kind of formula and draws a panel for each level of the conditioning variable. The panels contain scatterplots of the response and predictor, in this case height and weight, divided into subsets by the conditioning variable. The book shows a little trick that let's us customize xyplot, adding a regression line to each scatterplots.

plot.regression = function(x,y) {
  panel.xyplot(x,y)
  panel.abline(lm(y~x))
}

We pass the helper function plot.regression as a custom panel function in xyplot.

xyplot( height ~ weight | age.classes, data=kid.weights, panel=plot.regression)

There's quite a bit more to model formulae, but that's all I've figured out so far.

More on formulae

• Chapter 11 Statistical models in R from An Introduction to R
• R's model formula
• I() can be used insulate arithmetic expressions within formulae.

Cloud bioinformatics

Personally, the thing I love about cloud computing is never having to ask permission. There's no ops guy or pointy-haired boss between me and the launch-instance button. As lovely as that is, the cloud is also a powerful tool for scientific computing, particularly bioinformatics.

Next-gen sequencing, which can produce gigabytes per day, is one factor pushing bioinformatics into the cloud. Data analysis is now the major bottleneck for sequencing-based experiments. Labs are finding out that generating sequencing data is getting to be cheaper than processing it. According to Dave O’Connor Lab at the University of Wisconsin's Department of Pathology and Laboratory Medicine, "There is a real disconnect between the ability to collect next-generation sequence data (easy) and the ability to analyze it meaningfully (hard)."

O'Connor's group works with LabKey Software, a Seattle-based bioinformatics software company founded by the Fred Hutchinson Cancer Research Center. LabKey develops open-source data management software for proteomics, flow cytometry, plate-based assay, and HIV vaccine study data, described in a presentation by Lead Developer Adam Rauch. Their technology stack seems to include: Java, Spring, GWT, Lucene and Gauva (aka Google Collections). LabKey integrates with the impressive Galaxy genomics workflow system and the Trans-Proteomic Pipeline (TPP).

A good part of modern biology boils down to mining biological data, with the goal of correlating sequence, transcription or peptides to outputs like function, phenotype or disease. Machine learning and statistical modeling tend toward long-running CPU-intensive jobs that get run intermittently as new data arrives, making them ideal candidates for the cloud.

Amazon's EC2 seems to be better positioned than either Microsoft's Azure or Google's AppEngine for scientific computing. Amazon has been ahead of the curve in seeing the opportunity in genomic data overload. Microsoft has made some welcome efforts to attract scientific computing, including the Microsoft Biology Foundation and grants for scientific computing in Azure. But they're fighting a headwind arising from proprietary licensing and a closed ecosystem. Oddly, considering Google's reputation for openness, AppEngine looks surprisingly restrictive. Research computing typically involves building and installing binaries, programming in an odd patchwork of languages and long running CPU intensive tasks, none of which is particularly welcome on AppEngine. Maybe Google has a better offering in the works?

It's worth noting that open-source works without friction in cloud environments while many proprietary vendors have been slow to adapt their licensing models to on-demand scaling. For example, lots of folks are using R for machine learning in the cloud, while MatLab is still bogged down in licensing issues. The not-having-to-ask-permission aspect is lost.

According to Xconomy, Seattle has a growing advantage in the cloud. There are several Seattle companies operating in the bioinformatics and cloud spaces. Sage Bionetworks, also linked to FHCRC, was founded by Eric Schadt, also of Pacific Biosciences, and Stephen Friend former founder of Rosetta Inpharmatics. Revolution Analytics sells a scalable variant of R for all kinds of applications including life sciences. Seattle hosts a lot of activity in analytics, cloud computing and biotechnology, which will keep Seattle on the technology map for some time to come.

Using R for Introductory Statistics, Chapter 4

Chapter 4 of Using R for Introductory Statistics gets us started working with multivariate data. The question is: what are the relationships among the variables? One way to go about answering it is by pairwise comparison of variables. Another technique is to divide the data into categories by the values of some variables and analyze the remaining variables within each category. Different facets of the data can be encoded with color, shape and position to create visualizations that show graphically the relationships between several variables.

Taking variables one or two at a time, we can rely on our previous experience and apply our toolbox of univariate and bivariate techniques, such as histograms, correlation and linear regression.

We can also hold some variables constant and analyze the remaining variables in that context. Often, this involves conditioning on a categorical variable, as we did in Chapter 3 by splitting marathon finishing time into gender and age classes. As another example, the distribution of top speeds of italian sports cars describes a dependent variable, top speed, conditioned on two categorical variables, country of origin (Italy) and category (sports car). Because they're so familiar, cars make a great example.

R comes with a dataset called mtcars based on Motor Trend road tests for 32 cars in the 1973-74 model year. They recorded 11 statistics about each model of car. We can get a quick initial look using the pairs function which plots a thumbnail scatterplot for every pair of variables. Pairs is designed to work on numbers, but they've coded categorical values as integers in this data, so it works.

> names(mtcars)
 [1] "mpg"  "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear" "carb"
> pairs(mtcars)

Question 4.7 asks us to describe any trends relating weight, fuel efficiency, and number of cylinders. They also make the distinction between American made cars and imports, another categorical value along with cylinders. Let's make a two-panel plot. First, we'll make a boxplot comparing the distribution of mileage for imports and domestics. In the second panel, we'll combine all four variables.

# make a copy of mtcars 'cause we're going to add a column
cars <- mtcars

# add origin column as a factor to cars data frame
imports <- c(1:3, 8:14, 18:21, 26:28, 30:32)
origin <- rep("domestic", nrow(mtcars))
origin[imports] <- "import"
cars$origin <- factor(origin, levels=c('import', 'domestic'))

# make a vector of colors to color the data points
us.col <- rgb(0,0,255,192, maxColorValue=255)
im.col <- rgb(0,255,0,192, maxColorValue=255)
col <- rep(us.col, nrow(cars))
col[imports] <- im.col

# set up a two panel plot
par(mfrow=c(1,2))
par(mar=c(5,4,5,1)+0.1)

# draw boxplot in first panel
boxplot(mpg ~ origin, data=cars, col=c(im.col, us.col), outpch=19, outcol=c(us.col, im.col), ylab="mpg")
grid(nx=NA, ny=NULL)

# draw scatterplot in second panel
par(mar=c(5,0.5,5,2)+0.1)
plot(mpg~wt, data=cars, col=col, yaxt='n', pch=as.character(cars$cyl), xlab="weight (thousands of lbs)")
grid(nx=NA, ny=NULL)

# fit a line describing mpg as a function of weight
res <- lm(mpg ~ wt, data=cars)
abline(res, col=rgb(255,0,0,64, maxColorValue=255), lty=2)

# return parameters to defaults
par(mar=c(5,4,4,2)+0.1)
par(mfrow=c(1,1))

Domestics fair worse, but why? For one thing, the most fuel efficient cars are all light imports with 4 cylinder engines. Domestic cars are heavier with bigger engines and get worse milage. Other factors are certainly involved, but weight does a pretty good job of explaining fuel consumption, with a correlation of almost 87%.

> cor(cars$wt, cars$mpg)
[1] -0.8676594

Of course, fuel economy is not all there is to life on the open road. What about speed? We have quarter mile times, which probably measure acceleration better than top speed. We might hunt for variables that explain quarter mile performance by doing scatterplots and looking for correlation. The single factor that best correlates with qsec is horsepower. But a combination of factors does noticeably better - the power to weight ratio.

rm(origin)
attach(cars)

# transmission type encoded by color
palette <- c("#23809C","#7A1305")
col <- sapply(cars$am, function(x) palette[as.integer(x)+1])

# 5 panels
par(mfrow=c(1,5), omi=c(0,0.6,0.5,0.2))
par(mar=c(5,0.5,4,0.5)+0.1)

# weight
plot(qsec ~ wt, ylab="", xlab="weight", col= col, pch=as.character(cars$cyl))
mtext(sprintf("%.3f",cor(qsec,wt)), side=3)
grid(nx=NA, ny=NULL)

# displacement
plot(qsec ~ disp, ylab="", yaxt='n', xlab="displacement", col= col, pch=as.character(cars$cyl))
mtext(sprintf("%.3f",cor(qsec,disp)), side=3)
axis(2,labels=F, tick=T)
grid(nx=NA, ny=NULL)

# displacement / weight
plot(qsec ~ I(disp/wt), ylab="", yaxt='n', xlab="disp/wt", col= col, pch=as.character(cars$cyl))
mtext(sprintf("%.3f",cor(qsec,disp/wt)), side=3)
axis(2,labels=F, tick=T)
grid(nx=NA, ny=NULL)

# power
plot(qsec ~ hp, ylab="", yaxt='n', xlab="hp", col= col, pch=as.character(cars$cyl))
mtext(sprintf("%.3f",cor(qsec,hp)), side=3)
axis(2,labels=F, tick=T)
grid(nx=NA, ny=NULL)

# power / weight
plot(qsec ~ I(hp/wt), ylab="", yaxt='n', xlab="hp/wt", col= col, pch=as.character(cars$cyl))
mtext(sprintf("%.3f",cor(qsec,hp/wt)), side=3)
axis(2,labels=F, tick=T)
grid(nx=NA, ny=NULL)

# restore defaults
par(mar=c(5,4,4,2)+0.1)
par(mfrow=c(1,1), omi=c(0,0,0,0))

# add titles and legend
title("What factors influence acceleration?")
mtext("quarter mile time in seconds", side=2, padj=-4)
legend(x=88,y=21.5,c('automatic','manual'), fill=palette[c(1,2)], cex=0.6)

detach(cars)

The numbers above the scatterplots are correlation with qsec. Using position, character, color, multi-part graphs and ratios we've managed to visualize 5 variables, which show increasingly good correlation with the variable we're trying to explain.

Here's one theory that might emerge from staring at this data: 4-bangers with automatic transmission are slow. Here's another theory: there's an error in the data. Look at the slow 4 cylinder way at the top. It's quarter mile is nearly three seconds longer than the next slowest car. An outlier like that seems to need an explanation. That car, according to mtcars, is the Mercedes 230. But, the 230 is listed right next to the 240D - D for diesel. The 240D is a solid car. Many are still running. But they're famously slow. What are the odds that the times for these two cars got transposed?

We can check if the data supports our theories about how cars work. For example, we might guess that car makers are likely to put bigger engines in heavier cars to maintain adequate performance at the expense of gas mileage. Comparing weight with numbers of cylinders in the scatterplot above supports this idea. Displacement measures the total volume of the cylinders and that must be closely related to the number of cylinders. Try plot(disp ~ as.factor(cyl), data=mtcars). Displacement and carburetors are big determinants of the horsepower of an engine. Statisticians might be horrified, but try this: plot(hp ~ I(disp * carb), data=mtcars). Multiplying displacement by carburetion is a quick and dirty hack, but in this case, it seems to work out well.

Chapter 4 introduces parts of R's type system, specifically, lists and data.frames, along with subsetting operations and the apply family of functions. I don't go into it here because that was the first thing I learned about R and if you're a programmer, you'll probably want to do the same. One thing the book doesn't cover at all, so far as I can tell, is clustering.

Take a look at the plots above and see if you can't see the cars clustering into familiar categories: the two super-fast sports cars, the three land-yacht luxury cars weighing in at over 5000 pounds a piece, the 8-cylinder muscle cars, and the 4-cylinder econo-beaters. We don't have to do this by eye, because R has several clustering algorithms built in.

Hierarchical clustering (hclust) works by repeatedly merging the two most similar items or existing clusters together. Determining 'most similar' requires a measure of distance. In general, coming up with a good distance metric takes careful thought specific to the problem at hand, but let's live dangerously.

> plot(hclust(dist(scale(mtcars))), main='Hierarchical clustering of cars', sub='1973-74 model year', xlab='Cars', ylab="")

Not bad at all. The scale function helps out here by putting the columns on a common center and scale so the dist function ends up giving equal weight to each variable.

It's easy to analyze the bejeezes out of something you already you already know, like cars. The trick is to get a similar level of insight out of something entirely new.

Links to more...

• Basic R skills: splitting and plotting