#####
#####
###     STATISTICAL GRAPHICS FOR VISUALIZING DATA:
###     AN INTRODUCTION TO LATTICE GRAPHICS IN R
###
###     ICPSR Summer Program, Blalock Lecture Series
###     August 10-11, 2009
###
###     William G. Jacoby
###     ICPSR and Michigan State University
#####
#####

###
###   Some simple R objects
###

   x <- 3

   x

   y <- 5

   y

   x + y

###
###   Creating vectors and matrices
###

   x <- c(1, 2, 3, 4)

   x

   y <- c(10, 12, 14, 16)

   y

   xdata1 <- cbind(x, y)

   xdata2 <- rbind(y, x)

   xdata1

   xdata2

###
###     Read data from file "states, with names.txt".
###     Note header record.
###

   states <- read.table(file.choose(), header = T)

   states

###
###   Reading data from a file that does not contain a header record
###

   states2 <- read.table(file.choose())
   
   states2

###
###   Use colnames function to assign names to variables
###

   colnames(states2) <- c("state", "polprior", "party",
      "ideol", "igstr", "govsize", "region", "ptygov")

   states2[1:10, ]

###
###   Using the "foreign" package to
###   read data files that have been saved
###   in a different format-- STATA, in this case
###

   library(foreign)

   states3 <- read.dta(file.choose())

   states3[1:10, ]

###
###   Giving R instructions to find particular variables--
###   The first statement will cause an error!
###

   polprior

   states$polprior

   summary(states$polprior)

   with(states, mean(polprior))

###
###   Attaching data frames can be convenient, but it can 
###   also cause serious confusion and trouble
###

   attach(states)

   polprior

   polprior <- c(1, 2, 3, 4)

   polprior

   states$polprior

   detach(states)

   objects()

   remove(polprior)

###
###   Data structures
###

   is.vector(x)

   is.matrix(xdata1)

   is.matrix(states)

   is.data.frame(states)

   is.vector(states$polprior)

   is.vector(states$region)

   is.factor(states$region)

   xdata3 <- as.data.frame(xdata1)

   is.data.frame(xdata3)

###
###   Which objects current exist in the Workspace?
###

   objects()

###
###   Clean up, and remove unnecessary objects
###

   remove("states2", "states3", "x", "y", "xdata1", "xdata2", "xdata3")
   
###
###   Some very basic data manipulations
###

   str(states)

   fix(states)
 
   states$ptyidl <- (states$party + states$ideol) / 2
   
   states[3,2]
   
   states$polprior[3]
   
   states[states$region == "south",]
   
   states[states$region == "south", 1:4]


   
###
###   Getting help
###

   help(sd)
   
   ?sd
   
   sd

###
###   Help for functions in packages is not available until
###   the package is loaded in the current R session
###

   help(xyplot)
   
   library(lattice)
   
   help(xyplot)
   
###
###   Show some of the general display 
###   functions in the lattice package
###

   stripplot(~states$polprior)

   histogram(~states$polprior)

   densityplot(~states$polprior)

   xyplot(states$polprior ~ states$ideol)
   
###
###   Other ways to do the same things
###

   plot(states$polprior, rep(1, length(states$polprior)))
   
   hist(states$polprior)
   
   plot(density(states$polprior))
   
   plot(states$ideol, states$polprior)

###   
###   Show example of a multipanel display--
###   histograms of state policy priorities,
###   separately for states within each
###   of four regions

   histogram(~states$polprior | states$region)
   
###
###   Use the data= argument
###

   histogram(~polprior, data = states)

   xyplot(polprior ~ ideol, data = states)

###
###   Set aspect ratio explicitly, and
###   add axis labels
###

   histogram(~polprior, data = states,  
      aspect = .75, 
      xlab = "State policy priorities, 1992")

   xyplot(polprior ~ ideol, data = states,
      aspect = 1,
      xlab = "State electorate ideology",
      ylab = "state policy priorities, 1992")
   
###
###   Create a dot plot and modify data
###   to change the plot
###

   dotplot(state ~ polprior, data = states,
      aspect = 1.5,
      xlab = "State policy priorities, 1992")

###
###   Decrease size of scale labels
###

   dotplot(state ~ polprior, data = states,
      aspect = 1.5,
      xlab = "State policy priorities, 1992",
      scales = list(cex = .65)
   )

###
###   reorder factor "state" by priority scores
###

   states$state <- reorder(states$state, states$polprior)

###
###   Create dot plot with reordered factor
###

   dotplot(state ~ polprior, data = states,
      aspect = 1.5,
      xlab = "State policy priorities, 1992",
      scales = list(cex = .65)
   )
   
###
###   The "scales=" parameter can be used to 
###   change the details of the axes. For example,
###   replace the quantitative scale in a histogram
###   with textual labels
###

   histogram(~polprior, data = states,
     aspect = .75,
     xlab = "State policy priorities, 1992"
   )

   histogram(~polprior, data = states,
     aspect = .75,
     xlab = "State policy priorities, 1992",
     scales = list(x = list(tick.number = 3,
       at = c(10, 50, 90),
       labels = c("Particularized Benefits", "Mixture", "Collective Goods")))
   )

###
###   Create two-line textual labels
###

   histogram(~polprior, data = states,
     aspect = .75,
     xlab = "State policy priorities, 1992",
     scales = list(x = list(tick.number = 3,
       at = c(10, 50, 90),
       labels = c("Particularized\nBenefits", 
       "Mixed\nPriorities", 
       "Collective\nGoods")))
   )

###
###   Demonstrate use of logarithmic axes
###

###
###   Read data on campaign spending by
###   congressional challengers in 2006
###

   cong06 <- read.table(file.choose(), header = T)

###
###   Construct scatterplot of incumbent victory margin
###   versus challenger campaign spending, using raw data
###

   xyplot(margin ~ spend06, data = cong06)

###
###   Use "col=" parameter to change color of
###   plotting symbols, making them easier to see
###

   xyplot(margin ~ spend06, data = cong06,
      col = "black")
   
###
###   Use logarithm (base 10) to transform
###   campaign spending values and clarify
###   structure within resulting scatterplot
###

   xyplot(margin ~ log10(spend06), data = cong06,
      col = "black")
   
###
###   Modify horizontal axis
###

   xyplot(margin ~ log10(spend06), data = cong06,
      col = "black",
      scales = list(x = list(tick.number = 3, at = c(1, 2, 3),
         labels = c("$10,000", "$100,000", "$1,000,000"))),
         xlab = "Total campaign spending by incumbent\n(Logarithmic scale)",
      ylab = "Incumbent margin of victory")

###
###   Slightly fancier version, with intermediate tick
###   marks to emphasize logarithmic scale on X axis
###

   xyplot(margin ~ log10(spend06), data = cong06,
      col = "black",
      scales = list(x = list(at = c(log10(seq(from = 10, to =100, by = 10)),
        log10(seq(from = 200, to =1000, by = 100))),
         labels = c("$10,000", rep("", 8),"$100,000", 
         rep("", 8), "$1,000,000"))),
         xlab = "Total campaign spending by incumbent\n(Logarithmic scale)",
      ylab = "Incumbent margin of victory")
  
###
###   Even though we have been using the congressional
###   spending data, the 1992 state policy data still
###   exist as another object within the Workspace.
###   Therefore, that dataset is still accessible. 

   states[1:10,]

###
###   Show that characteristics of plotting region
###   elements are controlled by panel function
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities"
     )
     
   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     col = "black",
     pch = 4
     )
     
###
###   Preceding changes were actually made within the panel function.
###   The following function creates an identical graph
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=4, col = "black")}
     )

###
###   The following "empty" panel function shows how 
###   the panel function controls the plotting region
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) { }
     )

###
###   Use panel function to change line style 
###   and plotting symbol in dot plot
###

   dotplot(state ~ polprior, data = states,
      aspect = 1.5,
      panel = function (x, y) {
      panel.dotplot(x, y, lty = 2, pch = 1, col = "black")
      },
      xlab = "State policy priorities, 1992",
      scales = list(cex = .65)
   )
   
###
###   More complex use of panel functions to
###   change both point and line color
###   in a dot plot
###

###
###   Create dot plot of state government size
###

   dotplot(state ~ govsize, data = states)

###
###   In preceding dot plot, the factor, "state"
###   has levels ordered by the values of the variable
###   "polprior". Reorder this factor, and also change
###   some details of the dot plot. Note explicit use 
###   of FUN= argument in the "reorder" function.
###

   states$state <- reorder(states$state, states$govsize, FUN = mean)

   dotplot(state ~ govsize, data = states,
      aspect = 1.5,
      scales = list(cex = .65),
      xlab = "Size of state government (employees per capita)")

###
###   Use two panel functions, "panel.segments" and "panel.xyplot"
###   to create a dot plot in which lines run out to plotted points.
###   Note that this approach treats the dot plot as a sort of
###   "disguised scatterplot". 
###
###   The general format for the "panel.segments" function is:
###      panel.segments(x1, y1, x2, y2, ...)
###   Where: x1 and y1 are coordinates for one endpoint of
###   the line segment, x2 and y2 are the coordinates of the other
###   endpoint. This function is vectorized (which is why it works
###   in the dot plot).
###

   n <- length(states$state)

   dotplot(state ~ govsize, data = states,
     aspect = 1.5,
     scales = list(cex = .65),
     xlab = "Size of state government (employees per capita)",
     panel = function (x, y) {
     panel.segments(rep(0, n), as.numeric(states$state),
       x, as.numeric(states$state), lty = 2, col = "gray")
     panel.xyplot(states$govsize, as.numeric(states$state),
       pch = 16, col = "black")
     }
   )

###
###   The line segments encourage visual comparisons
###   of magnitude. In order to do this accurately,
###   the line segments should begin at the zero
###   point. The following function modifies the
###   dot plot to do this.

   dotplot(state ~ govsize, data = states,
      aspect = 1.5,
      scales = list(cex = .65),
      xlab = "Size of state government (employees per capita)",
      panel = function (x, y) {
      panel.segments(rep(0, n), as.numeric(states$state),
        x, as.numeric(states$state), lty = 2, col = "gray")
      panel.xyplot(states$govsize, as.numeric(states$state),
        pch = 16, col = "black")
      },
      xlim = c(-50, 900)
   )

###
###   Add a vertical line at the zero point.
###

   dotplot(state ~ govsize, data = states,
      aspect = 1.5,
      scales = list(cex = .65),
      xlab = "Size of state government (employees per capita)",
      panel = function (x, y) {
      panel.segments(rep(0, n), as.numeric(states$state),
        x, as.numeric(states$state), lty = 2, col = "gray")
      panel.xyplot(states$govsize, as.numeric(states$state),
        pch = 16, col = "black")
      panel.abline(v = 0, lty = 1, col = "gray")
      },
      xlim = c(-50, 900)
   )

###
###   Panel function can be used to plot two variables simultaneously
###

   xyplot(polprior ~ party, data = states,
     aspect = 1,
     xlab = "Orientations of state electorate",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "red")
     panel.xyplot(states$ideol, y, pch=4, col="green")}
     )

###
###   A graph which illustrates twenty possible plotting symbols 
###

   symbols <- c(1:20)
   
   symbols

   xyplot(symbols ~ symbols,
      aspect = 1,
      panel = function(x, y) {
      panel.xyplot(x, y, pch = symbols, col = "black")
      }
   )

###
###   Make symbols larger with "cex" parameter
###

   xyplot(symbols ~ symbols,
      aspect = 1,
      panel = function(x, y) {
      panel.xyplot(x, y, pch = symbols, col = "black",
        cex = 1.25)
      }
   )
      
###
###   Add the symbol numbers to the
###   graph of plotting symbols
###

   xyplot(symbols ~ symbols,
      aspect = 1,
      panel = function(x, y) {
      panel.xyplot(x, y, pch = symbols, col = "black",
        cex = 1.25)
      panel.text(x, y + .6, labels = symbols, cex = .75)
      }
   )

###
###   Panel functions can insert additional elements into the plotting region
###

###
###   A reference line at Y = 50 
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.abline(h = 50, lty = 2, col = "blue")}
     )

###
###   A line showing the OLS fit
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.lmline(x, y, lty = 1, col = "red")}
     )

###
###   A loess curve
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.loess(x, y, span = .75, family="symmetric",
       degree = 1, col = "green")}
     )

###
###   Place labels just to right of each plotted point
###   DISCLAIMER: For illustration only. The results look bad!
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.text(x+15, y, labels = states$state, cex = .75)}
     )

###
###   Use logical subscript operators to label points selectively.
###   Results look much better!
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.text(x[x>400 | y<20]-15, y[x>400 | y<20], 
       labels = states$state[x>400 | y<20], cex = .75)}
     )

###
###   Use panel functions and multiple subscripts
###   to modify plotting symbols selectively
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x[states$party < mean(states$party)], 
       y[states$party < mean(states$party)], 
       pch=4, col = "red")
     panel.xyplot(x[states$party >= mean(states$party)], 
       y[states$party >= mean(states$party)], 
       pch=1, col = "blue")}
     )

###
###   Return to scatterplot with multiple symbols, create a key 
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x[states$party < mean(states$party)], 
       y[states$party < mean(states$party)], 
       pch=4, col = "red")
     panel.xyplot(x[states$party >= mean(states$party)], 
       y[states$party >= mean(states$party)], 
       pch=1, col = "blue")},
     key = list(text = list(text = c("Republican electorates", "Democratic
electorates"), cex = .75),
       points = list(pch = c(4, 1), col = c("red", "blue"), cex = .75),
       space = "top", border = T)
     )

#####
#####
###     Displays involving subsets and multiple, superimposed, groups
#####
#####

###
###   Read data on state policy priorities, 1982-2005
###

policy.1982.2005 <- read.table(file.choose(), header = T)

policy.1982.2005 [1:5, ]

###
###   Create discrete version of year variable
###

policy.1982.2005$decade <- cut(policy.1982.2005$year, 
                             breaks = c(1981, 1989, 1999, 2006),
                             labels = c("1980s", "1990s", "2000s"))

histogram(~policy.priority, data = policy.1982.2005,
            aspect = .67, 
            subset = decade == "1980s",
            xlab = "State policy priorities",
            main = "Distribution for 1980s")

histogram(~policy.priority, data = policy.1982.2005,
            aspect = .67, 
            subset = (year >= 1993 & year <= 2001),
            xlab = "State policy priorities",
            main = "Distribution for Clinton terms")

###
###   Superimposed density plots, by decade
###

densityplot(~policy.priority, data = policy.1982.2005,
   aspect = 2/3,
   plot.points = F, ref = T,
   group = decade)

###
###   Add a key
###

densityplot(~policy.priority, data = policy.1982.2005,
   aspect = 2/3,
   plot.points = F, ref = T,
   group = decade, auto.key = T)

###
###   Additional changes to grouped density plot; note use
###   of "key" rather than "auto.key"
###

densityplot(~policy.priority, data = policy.1982.2005,
   aspect = 2/3,
   plot.points = F, ref = T,
   groups = decade, 
   col = c("red", "blue", "green"),
   key = list(text = list(labels = levels(policy.1982.2005$decade)),
      lines = list(col = c("red", "blue", "green")) ),
   xlab = "State policy priorities"
   )

###
###   Example of "groups=" parameter in a scatterplot
###

   xyplot(polprior ~ party, data = states,
      aspect = 1, groups = region, pch = 1:4,
      col = c("blue", "red", "green", "magenta"),
      key = list(text = list(labels = levels(states$region)),
         points = list(pch = 1:4, col = c("blue", "red", "green",
"magenta"))),
      xlab = "Partisanship of state electorate",
      ylab = "State policy priorities")

#####
#####
###     Some three-dimensional displays
#####
#####

   cloud(polprior ~ igstr + party, data = states,
      xlab = list(label = "Interest group\nstrength in state", cex = .75),
      ylab = list(label = "Partisanship of\nstate electorate", cex = .75),   
      zlab = list(label = "State policy\npriorities", cex = .75)
   )

###
###   Change viewing angle and plot symbol size
###

   cloud(polprior ~ igstr + party, data = states,
      cex = 1.25,
      xlab = list(label = "Interest group\nstrength in state", cex = .75),
      ylab = list(label = "Partisanship of\nstate electorate", cex = .75),   
      zlab = list(label = "State policy\npriorities", cex = .75),
      screen = list(z = 40, x = -80),
   )

   cloud(polprior ~ igstr + party, data = states,
      cex = 1.25,
      xlab = list(label = "Interest group\nstrength in state", cex = .75),
      ylab = list(label = "Partisanship of\nstate electorate", cex = .75),   
      zlab = list(label = "State policy\npriorities", cex = .75),
      screen = list(z = 20, x = -80),
   )

###
###   Estimate regression model and create 
###   wireframe plot of fitted OLS surface
###

   model1 <- lm(polprior ~ igstr + party, data = states)

   summary(model1)

   pty1 <- seq(from = min(states$party), to = max(states$party), length.out =
10)

   ig1 <- seq(from = min(states$igstr), to = max(states$igstr), length.out =
10)

   pred.data <- expand.grid(ig1, pty1)
   
   names(pred.data) <- c("igstr", "party")
   
   pred.data[1:5,]

   pred.data$yhat <- predict(model1, newdata = pred.data   )

   wireframe(yhat ~ igstr * party, data = pred.data,
      zlab = list(label = "State policy priority", 
         cex = .75, rot = 90),
      xlab = list(label = "Interest group\nstrength in state", cex = .75),
      ylab = list(label = "Partisanship of\nstate electorate", cex = .75)
   )

###
###   Change viewing angle and add color gradient
###

   wireframe(yhat ~ igstr * party, data = pred.data,
      shade = TRUE, 
      alpha = .5,
      col = "black", 
      zlab = list(label = "State policy priority", 
         cex = .75, rot = 90),
      xlab = list(label = "Interest group\nstrength in state", cex = .75),
      ylab = list(label = "Partisanship of\nstate electorate", cex = .75),
   screen = list(z = 20, x = -80)
   )

###
###   Create multipanel display, showing 
###   scatterplot of two variables, while 
###   controlling for a third variable
###

###   First, create the conditioning variable, using
###   a special data form called a "shingle".
###   Basically, this is a categorical variable with
###   overlapping categories

   given.party <- equal.count(states$party, number = 6, overlap = .75)

   plot(given.party)

###   Next, create scatterplot, conditioning on the
###   shingle variable

   xyplot(polprior ~ igstr | given.party, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, col = "black")
     panel.lmline(x, y, col = "red")}
     )
     
###
###   Calculate median policy priority scores
###   within each region, and create dot plot
###   of region medians
###

   medvals <- tapply(states$polprior, states$region, FUN = median)
   
   region <- as.vector(names(medvals))
   
   region
   
   region <- reorder(region, medvals)
   
   region
   
   dotplot(region ~ medvals, aspect = 1.5)

###
###   Perform another multiple regression analysis 
###   of state policy priorities, and 
###   use results to construct a residual plot,
###   with loess curve fitted to residuals

   model2 <- lm(polprior ~ ideol + party + igstr + govsize, data = states)

   summary(model2)

   xyplot(residuals(model2) ~ predict(model2),
     aspect = 1,
     panel = function(x, y) {
     panel.xyplot(x, y, col = "black")
     panel.abline(h = 0, lty = 2)
     panel.loess(x, y, span = .75, 
       degree = 1, family = "gaussian")}
     )

###
###   Use the expression function to employ
###   mathematical expressions in axis labels
###

   xyplot(residuals(model2) ~ predict(model2),
     aspect = 1,
     panel = function(x, y) {
     panel.xyplot(x, y, col = "black")
     panel.abline(h = 0, lty = 2)
     panel.loess(x, y, span = .75, 
       degree = 1, family = "gaussian")
     },
     xlab = expression(hat(Y)),
     ylab = expression(e[i]),
     ylim = c(-40, 40)
   )

###
###   A slightly fancier version of the axis labels
###

   xyplot(residuals(model2) ~ predict(model2),
     aspect = 1,
     panel = function(x, y) {
     panel.xyplot(x, y, col = "black")
     panel.abline(h = 0, lty = 2)
     panel.loess(x, y, span = .75, 
       degree = 1, family = "gaussian")
     },
     xlab = expression("Predicted values,"~~ hat(Y)[i]),
     ylab = expression("Residuals,"~~e[i]),
     ylim = c(-40, 40)
   )

###
###   Positioning multiple graphs on a single page:
###   First step is to create the graphs separately,
###   and save each one to an object.
###

   graph1 <- histogram(~polprior, data = states,
                xlab = "State policy priorities, 1992",
                aspect = .75)
   
   graph2 <- densityplot(~polprior, data = states,
                xlab = "State policy priorities, 1992",
                aspect = .75, col = "black")
   
   graph3 <- qqmath(~polprior, data = states,
                distribution = qunif, col = "black",
                xlab = "Quantiles", 
                ylab = "State policy priorities, 1992")
   
   set.seed(1234)
   graph4 <- stripplot(~polprior, data = states, aspect = .75, 
                xlab = "State policy priorities, 1992",
                col = "black", jitter = T)
                
###
###   Second step is to use "print()" function with position=
###   and more= parameters 
###

   print(graph1, position = c(0, .5, .5, 1), more = T)

   print(graph2, position = c(.5, .5, 1, 1), more = T)

   print(graph3, position = c(0, 0, .5, .5), more = T)

   print(graph4, position = c(.5, 0, 1, .5), more = F)

                                
###
###   Experimental evidence indicates that visual perception
###   of curves is optimized when the mean absolute angle of
###   the curve is 45 degrees. The process of "banking" is
###   used to change the aspect ratio to bring the mean absolute
###   angle as close to 45 degrees as possible
###

###
###   A line showing the OLS fit,
###   aspect ratio set to 1.0
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.lmline(x, y, lty = 1, col = "red")}
     )

###  
###   Same scatterplot, with regression line
###   set to 45 degrees. Note use of "aspect="
###   argument and the "prepanel" function
###

   xyplot(polprior ~ igstr, data = states,
     aspect = "xy",
     prepanel = function (x, y) {
     prepanel.lmline(x, y)}, 
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.lmline(x, y, lty = 1, col = "red")}
     )

###
###   A loess curve, with aspect ratio of 1.0
###

   xyplot(polprior ~ igstr, data = states,
     aspect = 1,
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.loess(x, y, span = .75, family="symmetric",
       degree = 1, col = "green")}
     )

###
###   Same scatterplot, banked to 45 degrees (that is, the mean
###   absolute angle of the line segments that make up the loess
###   curve is 45 degrees)
###

   xyplot(polprior ~ igstr, data = states,
     aspect = "xy",
     prepanel = function (x, y) {
     prepanel.loess(x, y, span = .75, 
       family="symmetric", degree = 1)},
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.loess(x, y, span = .75, family="symmetric",
       degree = 1, col = "green")}
     )

###
###   Print the last graph to a PDF file.
###   

pdf(file = "e://lattice graphics demo//icpsr 2007//banked scatter 1.pdf")

   xyplot(polprior ~ igstr, data = states,
     aspect = "xy",
     prepanel = function (x, y) prepanel.loess(x, y,span = .75,
family="symmetric",
       degree = 1),
     xlab = "Strength of interest groups in state",
     ylab = "State policy priorities",
     panel = function (x, y) {
     panel.xyplot(x, y, pch=1, col = "black")
     panel.loess(x, y, span = .75, family="symmetric",
       degree = 1, col = "green")}
     )

dev.off()

###
###   Setting graphical parameters explicitly, by
###   retrieving and modifying trellis settings
###

###
###   Create box plot of yearly policy priority
###   distributions for years 2000 through 2005.
###   Note the use of the "subset=" parameter.
###

bwplot(year ~ policy.priority, data = policy.1982.2005,
   subset = year >= 2000)

###
###   The "bwplot" function converts the numeric variable
###   "year" to a factor; when this happens, the levels of
###   the newly-created factor are set to consecutive integers.
###   In order to prevent this, coerce "year" into a factor;
###   in that case, the values of the variable become 
###   the levels of the factor
###

bwplot(as.factor(year) ~ policy.priority, data = policy.1982.2005,
   subset = year >= 2000)

###
###   Transpose the axes of the box plot display.
###

bwplot(policy.priority ~ as.factor(year), data = policy.1982.2005,
   subset = year >= 2000)


###
###   Assume we want to change several characteristics of this box plot:
###   The dot indicating the median should be a line segment. 
###   We also want to change the color of the plot to black,
###   and change the whiskers to solid (rather than dashed) lines.
###   These properties cannot be modified through the panel
###   function. So, we must modify the trellis settings.
###

###
###   First, determine which trellis settings are available
###

names(trellis.par.get())

###   
###   The "box.dot", "box.rectangle" and "box.umbrella" settings
###   affect various parts of the box plot. Print out the 
###   settings to see which parameters are associated each one.
###

trellis.par.get("box.dot")

###
###   The "str()" function prints the contents of an
###   object in a more efficient format.
###

str(trellis.par.get("box.dot"))

str(trellis.par.get("box.rectangle"))

str(trellis.par.get("box.umbrella"))

###
###   We can change the median dot to a line segment, using the
###   "pch=" parameter of "box.dot". We can change the dashed
###   lines to solid lines, using the "lty=" parameter of 
###   "dot.umbrella". And, all three settings have a "col="
###   parameter that can be used to change the color.
###

###
###   Use the "par.settings" argument in the "bwplot" function
###   to change the details of this single display.
###

bwplot(as.factor(year) ~ policy.priority, data = policy.1982.2005,
   subset = year >= 2000, col = "black",
   par.settings = list(
     box.dot = list(pch = "|", col = "black"),
     box.rectangle = list(col = "black"),
     box.umbrella =  list(col = "black", lty = 1))
)

###
###   Note that the outliers are still plotted in cyan.
###   We need to use the setting, "plot.symbol", to
###   change the color.
###

bwplot(as.factor(year) ~ policy.priority, data = policy.1982.2005,
   subset = year >= 2000, col = "black",
   par.settings = list(
     box.dot = list(pch = "|", col = "black"),
     box.rectangle = list(col = "black"),
     box.umbrella =  list(col = "black", lty = 1),
     plot.symbol = list(col = "black") )
)

###
###   Add axis labels
###

bwplot(as.factor(year) ~ policy.priority, data = policy.1982.2005,
   subset = year >= 2000, col = "black",
   par.settings = list(
     box.dot = list(pch = "|", col = "black"),
     box.rectangle = list(col = "black"),
     box.umbrella =  list(col = "black", lty = 1),
     plot.symbol = list(col = "black") ),
   xlab = "State policy priority scores",
   ylab = "Year"
)

###
###   Trellis settings can be changed directly, so the
###   modifications last for the remainder of the R session.
###

###
###   Retrieve the relevant settings, and assign each one
###   to a new object.
###

new.dot <- trellis.par.get("box.dot")

new.rectangle <- trellis.par.get("box.rectangle")

new.umbrella <- trellis.par.get("box.umbrella")

new.symbol <- trellis.par.get("plot.symbol")

new.dot$pch <- "|"

new.dot$col <- "black"

new.rectangle$col <- "black"

new.umbrella$col = "black"

new.umbrella$lty <- 1

new.symbol$col <- "black"

trellis.par.set(box.dot = new.dot,
   box.rectangle = new.rectangle,
   box.umbrella = new.umbrella,
   plot.symbol = new.symbol)

###
###   The preceding trellis settings have been
###   modified "globally," so the changes will affect
###   all subsequent graphs that use these settings.
###

bwplot(policy.priority ~ as.factor(year), data = policy.1982.2005,
   subset = year >= 1995 & year <= 1999,
   xlab = "Year",
   ylab = "State policy priorities"
)

bwplot(as.factor(year) ~ policy.priority, data = policy.1982.2005,
   subset = year >= 1990 & year <= 1994,
   xlab = "State policy priorities",
   ylab = "Year"
)

###
###   A subtle point: The "plot.symbol" setting was changed,
###   and this will affect all graphs that use plotting
###   symbols, not just box plots. So, note that the 
###   following function does NOT use the "col="
###   argument to change the plotting symbols from
###   cyan to black.
###

xyplot(polprior ~ govsize, data = states)