########################################################################
## Mon Oct  4 12:35:47 PDT 2004
## Demonstration program for lattice graphics
##
## This presents a simplified version of a few of the graphs from
## William Clevelands book Visualizing Data -- thoughtfully made available
## at http://cm.bell-labs.com/cm/ms/departments/sia/wsc
## This merely scratches the surface of these very cool graphic tools
## To snork this into your emacs buffer type
## C-x i ~carlm/213/LatticeGraphics/demonstration.r
########################################################################


## read in the barley data
barley<-read.table("/hdir/0/carlm/213/LatticeGraphics/barley",header=T)


##
## Let's take a look at what's in barley


names(barley)

summary(barley)

tapply(barley$yield, barley$site,summary)

## some notes on what we just did:
## 1) barley is a dataframe -- which is why we can reference its
## elements  using the $ notation
## 2) the output of the suumary() function is a vector. tapply cannot
## put a vector into an array element so what it produces is a list


## we need to load the lattice library in order to use the
## nifty functions therein


library(lattice)

## a 'dotplot'
## trellis.device() == x11()
dotplot( variety ~ yield | year * site, 
        data = barley,
        sub = list("Figure 1.6",cex=.8),
        xlab = "Barley Yield (bushels/acre)")



  ## barchart presents a little less information
barchart(variety  ~ yield | year* site, 
        data = barley, 
        sub = list("Figure 1.6",cex=.8),
        xlab = "Barley Yield (bushels/acre)")

## groups is a useful argument -- note the auto.key arg
dotplot(variety ~ yield | site, 
        data = barley,col=c('red','blue'),
        groups = year,auto.key=T)


### another dataset
singer<-read.table(file="/hdir/0/carlm/213/LatticeGraphics/singer",header=T)
names(singer)
tapply(singer$height,singer$voice.part,summary)

## simple non-lattice histogram 
hist(singer$height)

histogram(~ height | voice.part,
		  data = singer, 
#		  nint = 17, 
#		  endpoints = c(59.5, 76.5),
		  layout = c(2, 4), 
		  sub = list("Figure 1.2",cex=.8),
		  xlab = "Height (inches)")

## compare to normal distribution
## Q: How do we compare a vector of numbers to a distribution?
## A: plot the sorted values against the quantile of the (normal) distribution of the percentile of each (sorted) value.
bulblife<-rexp(n=500,rate=.1)
hist(bulblife)
## sorted values 
sorted.bl<-sort(bulblife)
## percentiles of sorted values
pt.sorted.bl<-rank(sorted.bl)/length(sorted.bl)
## quantiles of normal of pt.sorted.bl vs bl
plot(qnorm(pt.sorted.bl),sorted.bl)

sorted.bl[c(1,50)]
pt.sorted.bl[c(1,50)]
qnorm(pt.sorted.bl)[c(1,50)]
x11()
qqnorm(bulblife)
## a sort of diagonal line (through the 1st and 3rd quartiles)
qqline(bulblife)

qqmath(~ height | voice.part,type='p',
       data=singer,layout=c(2,4))

## specifying panel function -- this will make more sense after we
## have covererd functions, so store it away in the back of your brain
## for later.

qqmath(~ height | voice.part,
       data=singer,layout=c(2,4),
       panel= function(x,y,...) {
         panel.grid()
         panel.qqmathline(x, distribution= qnorm)
         panel.qqmath(x,y,type='smooth',...)
         panel.qqmath(x,y,type='p', col='green',...)}
       )


data(swiss)
names(swiss)


pairs(swiss)
pairs(swiss,panel=panel.smooth)

pairs(swiss,lower.panel=panel.smooth)


### xyplot and coplot --- similar yet different

coplot(Fertility ~ Education | Catholic,
       data=swiss)



coplot(Fertility ~ Education | Catholic,
       panel=panel.smooth, data=swiss)

xyplot(Fertility ~ Education ,
       data=swiss)

xyplot(Fertility ~ Education | Catholic,
       data=swiss)


###Hmmmm xyplot requires that conditioning variables be factors or shingles

##Here is one way of splitting Catholic into quintiles

CathQ5<-quantile(swiss$Catholic,prob=seq(0,1,.2))
xyplot(Fertility ~ Education | cut(Catholic,breaks=CathQ5),
       data=swiss)

## Here is an easier way of doing the same thing

xyplot(Fertility ~ Education | equal.count(Catholic,number=5),
       data=swiss)



## and here is how you would duplicate coplot

xyplot(Fertility ~ Education |
       shingle(Catholic,interval=co.intervals(Catholic)),
       data=swiss)


## a fancy xyplot using some panel functions
xyplot(Education ~ Infant.Mortality |
       equal.count(Catholic,number=4),
       data=swiss, panel=function(x,y){
         panel.grid()
         panel.xyplot(x,y)
         panel.xyplot(x,y,type="r")
         panel.xyplot(x,y,type="smooth",col='papayawhip')})


#postscript(file="bwtest.ps",horzintal=F)
bwplot( Infant.Mortality ~ equal.count(Agriculture,number=5)  |
       equal.count(Education,number=4),  data=swiss,xlab="Agriculture",
       layout=c(2,2))
#dev.off()