Using the R Lattice package for data analysis

R is this language I feel like I've had to re-learn like 10 times. Every time I stop for more than a month it's like starting from scratch. Very annoying. So I can have a reference and stop repeating myself I'm finally putting some of this into posts. Hopefully if you find yourself in a similar situation this'll help. The below analysis is based on a Coursera class in data modeling using R.

Let's first do a categorical break-out using the base package. We'll separate groups of variables on a single scatter plot:

[code language="splus"] y <- x + rnorm(100) g <- gl(2,50) g <-gl(2, 50, labels = c("male", "female"))

##let's check the structure str(g)

Factor w/ 2 levels "male","female": 1 1 1 1 1 1 1 1 1 1 …

plot(x, y) plot(x, y, type = "n") points(x[g== "male"], y[g=="male"], col = "green") points(x[g== "female"], y[g=="female"], col = "blue", pch=19) title("break out male from female") [/code]

categorical scatter plot, base package
categorical scatter plot, base package

In contrast to the base graphics package, Lattice functions generate plots in one shot rather than building them up piecewise: xyplot, bwplot, histogram, stripplot, dotplot, splom (like pairs in base system), levelplot. Generically,

y ~ x | f * g

  • on the left of the ~ is the y variable, on the right is the x variable
  • after the | are conditioning variables = they are optional; the * indicates an interaction.
  • lattice plots functions DIRECTLY on the graphics device
  • lattice graphics functions return an object of the class trellis
  • the print methods for lattice functions actually do the work of plotting the data on the graphics device
  • lattice functions return "plot objects" that can, in principle, be stored (but it's usually better to just save the code + data)
  • on the command line, trellis objects re auto-printed so it appears the function is plotting the data

The following example plots y vs x conditioned on f:

[code language="splus"] x <- rnorm(100) y <- x + rnorm(100, sd = 0.5) f <- gl(2, 50, labels = c("group 1", "group 2")) xyplot(y ~ x | f) [/code]

xyplot2
xyplot2

Let's now do a little exploration of the environmental dataset. Here's what it looks like:

head(environmental) ozone radiation temperature wind 1 41 190 67 7.4 2 36 118 72 8.0 3 12 149 74 12.6 4 18 313 62 11.5 5 23 299 65 8.6 6 19 99 59 13.8

ozone vs. radiation
ozone vs. radiation
Next, we'll create a shingle variable. We'll see shortly how useful this can be when using lattice for plotting.

[code language="splus"] summary(environmental$temperature) temp.cut <- equal.count(environmental$temperature, 4) # create s shingle variable help(equal.count) [/code]

Below are some details from R help below on shingles:

A shingle is a data structure used in Trellis, and is a generalization of factors to ‘continuous’ variables. It consists of a numeric vector along with some possibly overlapping intervals. These intervals are the ‘levels’ of the shingle. The levels andnlevels functions, usually applicable to factors, also work on shingles. The implementation of shingles is slightly different from S. There are print methods for shingles, as well as for printing the result of levels() applied to a shingle. For use in labelling, theas.character method can be used to convert levels of a shingle to character strings.

Here's what this shingle variable looks like:

[code language="splus"]

temp.cut [/code]

Data: [1] 67 72 74 62 65 59 61 69 66 68 58 64 66 57 68 62 59 73 61 61 67 81 79 76 [25] 82 90 87 82 77 72 65 73 76 84 85 81 83 83 88 92 92 89 73 81 80 81 82 84 [49] 87 85 74 86 85 82 86 88 86 83 81 81 81 82 89 90 90 86 82 80 77 79 76 78 [73] 78 77 72 79 81 86 97 94 96 94 91 92 93 93 87 84 80 78 75 73 81 76 77 71 [97] 71 78 67 76 68 82 64 71 81 69 63 70 75 76 68

Intervals: min max count 1 56.5 76.5 46 2 67.5 81.5 51 3 75.5 86.5 51 4 80.5 97.5 51

Overlap between adjacent intervals: [1] 27 30 31

[code language="splus"] #next let's plot ozone vs. radiation, conditioned on temp.cut xyplot(ozone ~ radiation | temp.cut, data = environmental) [/code]

From the chart below we can see that the relationship between ozone and radiation is dependent on the temperature: the bottom left and right panels (in which temp is lowest and second lowest) show not much of a relationship, whereas the top left - and even more so the top right - shows an increasing relationship between radiation and ozone.

lattice shingle plot
lattice shingle plot

[code language="splus"] #let's modify the layout to make it more intuitive - top-bottom layout, rather than quadrants. xyplot(ozone ~ radiation | temp.cut, data = environmental, layout = c(1, 4)) [/code]

shingle plot using R on environmental dataset
shingle plot using R on environmental dataset

[code language="splus"] xyplot(ozone ~ radiation | temp.cut, data = environmental, layout = c(1, 4), as.table = TRUE) #orders from top to bottom, ascending by temperature - much better [/code]

shingle plot reordered using R and lattice
shingle plot reordered using R and lattice

Finally, we can modify the layout ordering of the first chart to make it more intuitive:

[code language="splus"]

xyplot(ozone ~ radiation | temp.cut, data = environmental, layout = c(2, 2), as.table = TRUE) [/code]

shingle plot using lattice and R on environmental dataset
shingle plot using lattice and R on environmental dataset

Next we will create a custom panel in which we add a regression line to each panel. This is done by creating a custom function.

We can see from the panels below that as the temperature increases, so does the ozone level with respect to radiation….

[code language="splus"] #create a custom panel function to add a regression line to each panel xyplot(ozone ~ radiation | temp.cut, data = environmental, as.table = TRUE, pch=20, panel = function(x, y, ...){ panel.xyplot(x, y, ...) fit panel.abline(fit) }) [/code]

shingle plot with regressions by panel
shingle plot with regressions by panel

But on closer inspection it looks like panel 4 (bottom right) is non-linear. Let's try this again using another function called loess (local polynomial regression fitting).

[code language="splus"]

xyplot(ozone ~ radiation | temp.cut, data = environmental, as.table = TRUE, pch=20, panel = function(x, y, ...){ panel.xyplot(x, y, ...) panel.loess(x, y) }, xlab = "solar radiation", ylab = "ozone (ppb)", main = "Ozone vs. Solar Radiation") [/code]

regressions using polynomial fitting on shingles
regressions using polynomial fitting on shingles

Now it's looking better. In the environmental data set we have another variable, wind. let's put this all together including the wind variable.

[code language="splus"] #first create a wind cut, similar to what we did with temp wind.cut

repeat, but on the wind cut and temp cut

xyplot(ozone ~ radiation | temp.cut * wind.cut, data = environmental, as.table = TRUE, pch=20, panel = function(x, y, …){ panel.xyplot(x, y, …) panel.loess(x, y) }, xlab = "solar radiation", ylab = "ozone (ppb)", main = "Ozone vs. Solar Radiation") [/code]

regression on multiple shingles
regression on multiple shingles

we can see how the relationship changes: high temp and low/med wind is most interesting. the highlighting illustrates where we are on the two conditioning variables, temp and wind.

We can also do other interesting things in Lattice. For example we can use splom function (for scatter plot matrix) to generate a scatter plot matrix of all variables in the data set.

splom(~ environmental)

scatter plot matrix
scatter plot matrix

Another interesting thing we can do is generate a histogram with the conditional variable:

histogram(~ ozone | wind.cut, data = environmental, as.table=TRUE) we can see that as the wind increases the distribution changes as well - ozone is more concentrated in the lower wind ranges.

shingled histogram
shingled histogram

we can of course condition on both wind and temp as with the xyplot function:

histogram(~ ozone | wind.cut * temp.cut, data = environmental, as.table=TRUE)

shingled histogram with interactions
shingled histogram with interactions
We can see in the bottom left panel that when the wind is low and temp high that the values of ozone are spread out, whereas in the top right panel we see that most of the ozone values are zero when the temp is low and the wind is high.

Get notified about new posts