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"))
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]
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]
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.
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.
[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]
[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]
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]
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]
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.
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
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
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.