# ------------------------------------------------
# Program: Univariate4.S  
#  Author: Steven M. Boker
#    Date: Tue Sep 28 10:03:34 EST 2004
#
# Provides examples of the use of QQ plots for the 
# visualization of univariate distributions 
#
# ------------------------------------------------

# ------------------------------------------------
# Quantile-Quantile plot of sepal.length of Setosa versus Versicolor 

graphsheet(height=6.4,width=6.4)
qq(variety ~ sepal.length, 
         data=iris,
         subset=variety=="Setosa" | variety=="Versicolor",
         aspect=1)

# ------------------------------------------------
#  Create a sample from each of two Normal distributions
#  with different means and compare them using a QQ plot.

x <- c(rnorm(300, mean=0, sd=1), rnorm(300, mean=2, sd=1))

# Create a category vector

xtype <- as.factor(c(rep(0,300), rep(1,300)))

# Make the QQ plot

graphsheet(height=6.4,width=6.4)
qq(xtype ~ x, aspect=1)


# ------------------------------------------------
#  Create a sample from each of two Normal distributions
#  with different variance and compare them using a QQ plot.

x <- c(rnorm(300, mean=0, sd=1), rnorm(300, mean=0, sd=2))

# Create a category vector

xtype <- as.factor(c(rep(0,300), rep(1,300)))

# Make the QQ plot

graphsheet(height=6.4,width=6.4)
qq(xtype ~ x, aspect=1)

# ------------------------------------------------
#  Create a sample from each of two Normal distributions
#  with different mean & variance and compare them using a QQ plot.

x <- c(rnorm(300, mean=0, sd=1), rnorm(300, mean=2, sd=2))

# Create a category vector

xtype <- as.factor(c(rep(0,300), rep(1,300)))

# Make the QQ plot

graphsheet(height=6.4,width=6.4)
qq(xtype ~ x, aspect=1)

# ------------------------------------------------
#  Create an unknown mixture distribution sample 
#  from each of two Normal distributions
#  with different means and compare them 
#  to the normal distribution using a QQ plot.

x <- c(rnorm(300, mean=0, sd=1), rnorm(300, mean=8, sd=1))

# Create a category vector

xtype <- as.factor(c(rep(0,300), rep(1,300)))

# Make the QQ plot

graphsheet(height=6.4,width=6.4)
qqmath(~ x, 
         distribution=qnorm, data=iris,
         prepanel = prepanel.qqmathline, 
	 panel = function(x, y) {
	   panel.qqmathline(y, distribution = qnorm)
	   panel.qqmath(x, y)
	 },
         aspect=1,
         xlab = list(label="Normal Distribution",cex=1.75),
         ylab=list(label="Mixture with Differing Means",cex=1.75),
         scales=list(cex=1.25))

# ------------------------------------------------
#  Create an unknown mixture distribution sample 
#  from each of two Normal distributions
#  with different mean & variance and compare them 
#  to the normal distribution using a QQ plot.

x <- c(rnorm(300, mean=0, sd=1), rnorm(300, mean=0, sd=4))

# Create a category vector

xtype <- as.factor(c(rep(0,300), rep(1,300)))

# Make the QQ plot

graphsheet(height=6.4,width=6.4)
qqmath(~ x, 
         distribution=qnorm, data=iris,
         prepanel = prepanel.qqmathline, 
	 panel = function(x, y) {
	   panel.qqmathline(y, distribution = qnorm)
	   panel.qqmath(x, y)
	 },
         aspect=1,
         xlab = list(label="Normal Distribution",cex=1.75),
         ylab=list(label="Mixture with Differing Variances",cex=1.75),
         scales=list(cex=1.25))

# ------------------------------------------------
#  Create a sample from each of two different types of distributions
#  compare them using a QQ plot.

x <- c(rnorm(300, mean=0, sd=1), rexp(300))

# Create a category vector

xtype <- as.factor(c(rep(0,300), rep(1,300)))

# Make the QQ plot

graphsheet(height=6.4,width=6.4)
qq(xtype ~ x, aspect=1)

# ------------------------------------------------
#  Subtract a grand mean and create a boxplot.

graphsheet(height=6.4,width=6.4)
bwplot(iris$variety ~ iris$sepal.length, 
       aspect=1)

graphsheet(height=6.4,width=6.4)
newSepalLength <- iris$sepal.length - mean(iris$sepal.length)
bwplot(iris$variety ~ newSepalLength, 
       aspect=1)

# ------------------------------------------------
#  Subtract marginal means and create a boxplot.

oneway(sepal.length~variety, spread = 1, location=mean, data=iris)

graphsheet(height=6.4,width=6.4)
bwplot(variety ~ oneway(sepal.length~variety)$residuals, 
       data=iris,
       aspect=1)

# ------------------------------------------------
#  Checking for homogeneity of residuals and pooled residuals.

resSepLen <- oneway(sepal.length~variety, data=iris)$residuals
graphsheet(height=6.4,width=6.4)
qqmath(~ resSepLen | iris$variety, 
	distribution = substitute(function(p) quantile(resSepLen, p)),
	panel=function(x,y){
		panel.grid()
		panel.abline(0, 1)
		panel.qqmath(x, y)
	},
	aspect=1,
	layout=c(2,2),
	xlab = "Pooled Residual Sepal Length (cm)",
	ylab = "Residual Sepal Length (cm)")

# ------------------------------------------------
#  Checking for normality of pooled residuals.

graphsheet(height=6.4,width=6.4)
qqmath(~ oneway(sepal.length~variety)$residuals,
	data = iris,
	prepanel = prepanel.qqmathline, 
	panel = function(x, y) {
	    panel.qqmathline(y, distribution = qnorm)
	    panel.qqmath(x, y)
	},
	aspect=1,
	xlab = "Unit Normal Quantiles",
	ylab = "Residual Sepal Length (cm)")