**everyday analytics**, and kindly contributed to R-bloggers)

#### Introduction

The point is that my line of business requires travel, and sometimes that is a lot of the time, like say almost all of last year. Inevitable comparisons to George Clooney’s character in Up in the Air were made (ironically I started to read that book, then left it on a plane in a seatback pocket), requests about favours involving duty free, and of course many observations and gently probing questions about frequent flier miles (FYI I’ve got more than most people, but a lot less than the entrepreneur I sat next to one time, who claimed to have close to 3 million).

But I digress.

#### Background

In my case this means flying Delta.

So I happened to notice in one of my many visits to Delta’s website that they have data on all of their aircraft in a certain site section. I thought this would be an interesting data set on which to do some analysis, as it has both quantitative and qualitative information and is relatively complex. What can we say about the different aircraft in Delta’s fleet, coming at it with ‘fresh eyes’? Which planes are similar? Which are dissimilar?

Aircraft data card from Delta.com |

The data set comprises 33 variables on 44 aircraft taken from Delta.com, including both quantitative measures on attributes like cruising speed, accommodation and range in miles, as well as categorical data on, say, whether a particular aircraft has Wi-Fi or video. These binary categorical variables were transformed into quantitative variables by assigning them values of either 1 or 0, for yes or no respectively.

#### Analysis

data <- read.csv(file="delta.csv", header=T, sep=",", row.names=1)

# scatterplot matrix of intermediary (size/non-categorical) variables

plot(data[,16:22])

We can see that there are pretty strong positive correlations between all these variables, as all of them are related to the aircraft’s overall size. Remarkably there is an almost perfectly linear relationship between wingspan and tail height, which perhaps is related to some principle of aeronautical engineering of which I am unaware.

The exception here is the variable right in the middle which is the number of engines. There is one lone outlier [Boeing 747-400 (74S)] which has four, while all the other aircraft have two. In this way the engines variable is really more like a categorical variable, but we shall as the analysis progresses that this is not really important, as there are other variables which more strongly discern the aircraft from one another than this.

How do we easier visualize a high-dimensional data set like this one? By using a dimensionality reduction technique like principal components analysis.

*Principal Components Analysis*

Next let’s say I know nothing about dimensionality reduction techniques and just naively apply principle components to the data in R:

# Naively apply principal components analysis to raw data and plot

pc <- princomp(data)

plot(pc)

Taking that approach we can see that the first principal component has a standard deviation of around 2200 and accounts for over 99.8% of the variance in the data. Looking at the first column of loadings, we see that the first principle component is just the range in miles.

# First component dominates greatly. What are the loadings?

summary(pc) # 1 component has > 99% variance

loadings(pc) # Can see all variance is in the range in miles

Importance of components:

Comp.1 Comp.2 Comp.3 Comp.4

Standard deviation 2259.2372556 6.907940e+01 2.871764e+01 2.259929e+01

Proportion of Variance 0.9987016 9.337038e-04 1.613651e-04 9.993131e-05

Cumulative Proportion 0.9987016 9.996353e-01 9.997966e-01 9.998966e-01

Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8

Seat.Width..Club. -0.144 -0.110

Seat.Pitch..Club. -0.327 -0.248 0.189

Seat..Club.

Seat.Width..First.Class. 0.250 -0.160 -0.156 0.136

Seat.Pitch..First.Class. 0.515 -0.110 -0.386 0.112 -0.130 0.183

Seats..First.Class. 0.258 -0.124 -0.307 -0.109 0.160 0.149

Seat.Width..Business. -0.154 0.142 -0.108

Seat.Pitch..Business. -0.514 0.446 -0.298 0.154 -0.172 0.379

Seats..Business. -0.225 0.187

Seat.Width..Eco.Comfort. 0.285 -0.224

Seat.Pitch..Eco.Comfort. 0.159 0.544 -0.442

Seats..Eco.Comfort. 0.200 -0.160

Seat.Width..Economy. 0.125 0.110

Seat.Pitch..Economy. 0.227 0.190 -0.130

Seats..Economy. 0.597 -0.136 0.345 -0.165 0.168

Accommodation 0.697 -0.104 0.233

Cruising.Speed..mph. 0.463 0.809 0.289 -0.144 0.115

Range..miles. 0.999

Engines

Wingspan..ft. 0.215 0.103 -0.316 -0.357 -0.466 -0.665

Tail.Height..ft. -0.100 -0.187

Length..ft. 0.275 0.118 -0.318 0.467 0.582 -0.418

Wifi

Video

Power

Satellite

Flat.bed

Sleeper

Club

First.Class

Business

Eco.Comfort

Economy

This is because the scale of the different variables in the data set is quite variable; we can see this by plotting the variance of the different columns in the data frame (regular scaling on the left, logarithmic on the right):

# verify by plotting variance of columns

mar <- par()$mar

par(mar=mar+c(0,5,0,0))

barplot(sapply(data, var), horiz=T, las=1, cex.names=0.8)

barplot(sapply(data, var), horiz=T, las=1, cex.names=0.8, log='x')

par(mar=mar)

We correct for this by scaling the data using the scale() function. We can then verify that the variances across the different variables are equal so that when we apply principal components one variable does not dominate.

# Scale

data2 <- data.frame(scale(data))

# Verify variance is uniform

plot(sapply(data2, var))

After applying the scale() function the variance is now constant across variables |

Now we can apply principal components to the scaled data. Note that this can also be done automatically in call to the prcomp() function by setting the parameter scale=TRUE. Now we see a result which is more along the lines of something we would expect:

# Proceed with principal components

pc <- princomp(data2)

plot(pc)

plot(pc, type='l')

summary(pc) # 4 components is both 'elbow' and explains >85% variance

Great, so now we’re in business. There are various rules of thumb for selecting the number of principal components to retain in an analysis of this type, two of which I’ve read about are:

- Pick the number of components which explain 85% or greater of the variation
- Use the ‘elbow’ method of the scree plot (on right)

# Get principal component vectors using prcomp instead of princomp

pc <- prcomp(data2)

# First for principal components

comp <- data.frame(pc$x[,1:4])

# Plot

plot(comp, pch=16, col=rgb(0,0,0,0.5))

So what were are looking at here are twelve 2-D projections of data which are in a 4-D space. You can see there’s a clear outlier in all the dimensions, as well as some bunching together in the different projections.

library(rgl)

# Multi 3D plot

plot3d(comp$PC1, comp$PC2, comp$PC3)

plot3d(comp$PC1, comp$PC3, comp$PC4)

*Cluster Analysis*

# Determine number of clusters

wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))

for (i in 2:15) wss[i] <- sum(kmeans(mydata,

centers=i)$withinss)

plot(1:15, wss, type="b", xlab="Number of Clusters",

ylab="Within groups sum of squares")

Clustering without the nstart parameter can lead to variable results for each run |

Clustering with the nstart and iter.max parameters leads to consistent results, allowing proper interpretation of the scree plot |

# From scree plot elbow occurs at k = 4

# Apply k-means with k=4

k <- kmeans(comp, 4, nstart=25, iter.max=1000)

library(RColorBrewer)

library(scales)

palette(alpha(brewer.pal(9,'Set1'), 0.5))

plot(comp, col=k$clust, pch=16)

# 3D plot

plot3d(comp$PC1, comp$PC2, comp$PC3, col=k$clust)

plot3d(comp$PC1, comp$PC3, comp$PC4, col=k$clust)

# Cluster sizes

sort(table(k$clust))

clust <- names(sort(table(k$clust)))

# First cluster

row.names(data[k$clust==clust[1],])

# Second Cluster

row.names(data[k$clust==clust[2],])

# Third Cluster

row.names(data[k$clust==clust[3],])

# Fourth Cluster

row.names(data[k$clust==clust[4],])

[1] “CRJ 100/200 Pinnacle/SkyWest” “CRJ 100/200 ExpressJet”

[3] “E120″ “ERJ-145″

[1] “Airbus A330-200″ “Airbus A330-200 (3L2)”

[3] “Airbus A330-200 (3L3)” “Airbus A330-300″

[5] “Boeing 747-400 (74S)” “Boeing 757-200 (75E)”

[7] “Boeing 757-200 (75X)” “Boeing 767-300 (76G)”

[9] “Boeing 767-300 (76L)” “Boeing 767-300 (76T)”

[11] “Boeing 767-300 (76Z V.1)” “Boeing 767-300 (76Z V.2)”

[13] “Boeing 767-400 (76D)” “Boeing 777-200ER”

[15] “Boeing 777-200LR”

[1] “Airbus A319″ “Airbus A320″ “Airbus A320 32-R”

[4] “Boeing 717″ “Boeing 737-700 (73W)” “Boeing 737-800 (738)”

[7] “Boeing 737-800 (73H)” “Boeing 737-900ER (739)” “Boeing 757-200 (75A)”

[10] “Boeing 757-200 (75M)” “Boeing 757-200 (75N)” “Boeing 757-200 (757)”

[13] “Boeing 757-200 (75V)” “Boeing 757-300″ “Boeing 767-300 (76P)”

[16] “Boeing 767-300 (76Q)” “Boeing 767-300 (76U)” “CRJ 700″

[19] “CRJ 900″ “E170″ “E175″

[22] “MD-88″ “MD-90″ “MD-DC9-50″

Ahhh, that’s the way fly (some day, some day…). This is apparently the plane professional sports teams and the American military often charter to fly – this article in the Sydney Morning Herald has more details.

Top: CRJ100/200. Bottom left: Embraer E120. Bottom right: Embraer ERJ-145. |

I’ve flown many times in the venerable CRJ 100/200 series planes, in which I can assure you there is only economy seating, and which I like to affectionately refer to as “little metal tubes of suffering.”

These are split into two clusters, which seem to again divide the planes approximately by size (both physical and accommodation), though there is crossover in the Boeing craft.

# Compare accommodation by cluster in boxplot

boxplot(data$Accommodation ~ k$cluster,

xlab='Cluster', ylab='Accommodation',

main='Plane Accommodation by Cluster')

# Compare presence of seat classes in largest clusters

data[k$clust==clust[3],30:33]

data[k$clust==clust[4],30:33]

First.Class | Business | Eco.Comfort | Economy | |

Airbus A330-200 | 0 | 1 | 1 | 1 |

Airbus A330-200 (3L2) | 0 | 1 | 1 | 1 |

Airbus A330-200 (3L3) | 0 | 1 | 1 | 1 |

Airbus A330-300 | 0 | 1 | 1 | 1 |

Boeing 747-400 (74S) | 0 | 1 | 1 | 1 |

Boeing 757-200 (75E) | 0 | 1 | 1 | 1 |

Boeing 757-200 (75X) | 0 | 1 | 1 | 1 |

Boeing 767-300 (76G) | 0 | 1 | 1 | 1 |

Boeing 767-300 (76L) | 0 | 1 | 1 | 1 |

Boeing 767-300 (76T) | 0 | 1 | 1 | 1 |

Boeing 767-300 (76Z V.1) | 0 | 1 | 1 | 1 |

Boeing 767-300 (76Z V.2) | 0 | 1 | 1 | 1 |

Boeing 767-400 (76D) | 0 | 1 | 1 | 1 |

Boeing 777-200ER | 0 | 1 | 1 | 1 |

Boeing 777-200LR | 0 | 1 | 1 | 1 |

First.Class | Business | Eco.Comfort | Economy | |

Airbus A319 | 1 | 0 | 1 | 1 |

Airbus A320 | 1 | 0 | 1 | 1 |

Airbus A320 32-R | 1 | 0 | 1 | 1 |

Boeing 717 | 1 | 0 | 1 | 1 |

Boeing 737-700 (73W) | 1 | 0 | 1 | 1 |

Boeing 737-800 (738) | 1 | 0 | 1 | 1 |

Boeing 737-800 (73H) | 1 | 0 | 1 | 1 |

Boeing 737-900ER (739) | 1 | 0 | 1 | 1 |

Boeing 757-200 (75A) | 1 | 0 | 1 | 1 |

Boeing 757-200 (75M) | 1 | 0 | 1 | 1 |

Boeing 757-200 (75N) | 1 | 0 | 1 | 1 |

Boeing 757-200 (757) | 1 | 0 | 1 | 1 |

Boeing 757-200 (75V) | 1 | 0 | 1 | 1 |

Boeing 757-300 | 1 | 0 | 1 | 1 |

Boeing 767-300 (76P) | 1 | 0 | 1 | 1 |

Boeing 767-300 (76Q) | 1 | 0 | 1 | 1 |

Boeing 767-300 (76U) | 0 | 1 | 1 | 1 |

CRJ 700 | 1 | 0 | 1 | 1 |

CRJ 900 | 1 | 0 | 1 | 1 |

E170 | 1 | 0 | 1 | 1 |

E175 | 1 | 0 | 1 | 1 |

MD-88 | 1 | 0 | 1 | 1 |

MD-90 | 1 | 0 | 1 | 1 |

MD-DC9-50 | 1 | 0 | 1 | 1 |

Looking at the raw data, the difference I can ascertain between the largest two clusters is that all the aircraft in the one have first class seating, whereas all the planes in the other have business class instead [the one exception being the Boeing 767-300 (76U)].

#### Conclusions

If I did this again, I would structure the data differently and see what relationships such analysis could draw out using only select parts of the data (*e.g. *aircraft measurements only). The interesting lesson here is that it when using techniques like dimensionality reduction and clustering it is not only important to be mindful of applying them correctly, but also what variables are in your data set and how they are represented.

#### References & Resources

*http://en.wikipedia.org/wiki/Principal_components_analysis*

*http://www.smh.com.au/travel/travel-planning/travel-news/plane-luxury-how-us-sports-stars-fly-20110106-19h6y.html*

Code & Data on github*https://github.com/mylesmharrison/delta_PCA_kmeans*

**leave a comment**for the author, please follow the link and comment on his blog:

**everyday analytics**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...