# NOT RUN {
data(Macdonell)
# display the frequency table
xtabs(frequency ~ finger+round(height,3), data=Macdonell)
## Some examples by james.hanley@mcgill.ca October 16, 2011
## http://www.biostat.mcgill.ca/hanley/
## See: http://www.biostat.mcgill.ca/hanley/Student/
###############################################
## naive contour plots of height and finger ##
###############################################
# make a 22 x 42 table
attach(Macdonell)
ht <- unique(height)
fi <- unique(finger)
fr <- t(matrix(frequency, nrow=42))
detach(Macdonell)
dev.new(width=10, height=5) # make plot double wide
op <- par(mfrow=c(1,2),mar=c(0.5,0.5,0.5,0.5),oma=c(2,2,0,0))
dx <- 0.5/12
dy <- 0.5/12
plot(ht,ht,xlim=c(min(ht)-dx,max(ht)+dx),
ylim=c(min(fi)-dy,max(fi)+dy), xlab="", ylab="", type="n" )
# unpack 3000 heights while looping though the frequencies
heights <- c()
for(i in 1:22) {
for (j in 1:42) {
f <- fr[i,j]
if(f>0) heights <- c(heights,rep(ht[i],f))
if(f>0) text(ht[i], fi[j], toString(f), cex=0.4, col="grey40" )
}
}
text(4.65,13.5, "Finger length (cm)",adj=c(0,1), col="black") ;
text(5.75,9.5, "Height (feet)", adj=c(0,1), col="black") ;
text(6.1,11, "Observed bin\nfrequencies", adj=c(0.5,1), col="grey40",cex=0.85) ;
# crude countour plot
contour(ht, fi, fr, add=TRUE, drawlabels=FALSE, col="grey60")
# smoother contour plot (Galton smoothed 2-D frequencies this way)
# [Galton had experience with plotting isobars for meteorological data]
# it was the smoothed plot that made him remember his 'conic sections'
# and ask a mathematician to work out for him the iso-density
# contours of a bivariate Gaussian distribution...
dx <- 0.5/12; dy <- 0.05 ; # shifts caused by averaging
plot(ht,ht,xlim=c(min(ht),max(ht)),ylim=c(min(fi),max(fi)), xlab="", ylab="", type="n" )
sm.fr <- matrix(rep(0,21*41),nrow <- 21)
for(i in 1:21) {
for (j in 1:41) {
smooth.freq <- (1/4) * sum( fr[i:(i+1), j:(j+1)] )
sm.fr[i,j] <- smooth.freq
if(smooth.freq > 0 )
text(ht[i]+dx, fi[j]+dy, sub("^0.", ".",toString(smooth.freq)), cex=0.4, col="grey40" )
}
}
contour(ht[1:21]+dx, fi[1:41]+dy, sm.fr, add=TRUE, drawlabels=FALSE, col="grey60")
text(6.05,11, "Smoothed bin\nfrequencies", adj=c(0.5,1), col="grey40", cex=0.85) ;
par(op)
dev.new() # new default device
#######################################
## bivariate kernel density estimate
#######################################
if(require(KernSmooth)) {
MDest <- bkde2D(MacdonellDF, bandwidth=c(1/8, 1/8))
contour(x=MDest$x1, y=MDest$x2, z=MDest$fhat,
xlab="Height (feet)", ylab="Finger length (cm)", col="red", lwd=2)
with(MacdonellDF, points(jitter(height), jitter(finger), cex=0.5))
}
#############################################################
## sunflower plot of height and finger with data ellipses ##
#############################################################
with(MacdonellDF,
{
sunflowerplot(height, finger, size=1/12, seg.col="green3",
xlab="Height (feet)", ylab="Finger length (cm)")
reg <- lm(finger ~ height)
abline(reg, lwd=2)
if(require(car)) {
dataEllipse(height, finger, plot.points=FALSE, levels=c(.40, .68, .95))
}
})
############################################################
## Sampling distributions of sample sd (s) and z=(ybar-mu)/s
############################################################
# note that Gosset used a divisor of n (not n-1) to get the sd.
# He also used Sheppard's correction for the 'binning' or grouping.
# with concatenated height measurements...
mu <- mean(heights) ; sigma <- sqrt( 3000 * var(heights)/2999 )
c(mu,sigma)
# 750 samples of size n=4 (as Gosset did)
# see Student's z, t, and s: What if Gosset had R?
# [Hanley J, Julien M, and Moodie E. The American Statistician, February 2008]
# see also the photographs from Student's notebook ('Original small sample data and notes")
# under the link "Gosset' 750 samples of size n=4"
# on website http://www.biostat.mcgill.ca/hanley/Student/
# and while there, look at the cover of the Notebook containing his yeast-cell counts
# http://www.medicine.mcgill.ca/epidemiology/hanley/Student/750samplesOf4/Covers.JPG
# (Biometrika 1907) and decide for yourself why Gosset, when forced to write under a
# pen-name, might have taken the name he did!
# PS: Can you figure out what the 750 pairs of numbers signify?
# hint: look again at the numbers of rows and columns in Macdonell's (frequency) Table III.
n <- 4
Nsamples <- 750
y.bar.values <- s.over.sigma.values <- z.values <- c()
for (samp in 1:Nsamples) {
y <- sample(heights,n)
y.bar <- mean(y)
s <- sqrt( (n/(n-1))*var(y) )
z <- (y.bar-mu)/s
y.bar.values <- c(y.bar.values,y.bar)
s.over.sigma.values <- c(s.over.sigma.values,s/sigma)
z.values <- c(z.values,z)
}
op <- par(mfrow=c(2,2),mar=c(2.5,2.5,2.5,2.5),oma=c(2,2,0,0))
# sampling distributions
hist(heights,breaks=seq(4.5,6.5,1/12), main="Histogram of heights (N=3000)")
hist(y.bar.values, main=paste("Histogram of y.bar (n=",n,")",sep=""))
hist(s.over.sigma.values,breaks=seq(0,4,0.1),
main=paste("Histogram of s/sigma (n=",n,")",sep=""));
z=seq(-5,5,0.25)+0.125
hist(z.values,breaks=z-0.125, main="Histogram of z=(ybar-mu)/s")
# theoretical
lines(z, 750*0.25*sqrt(n-1)*dt(sqrt(n-1)*z,3), col="red", lwd=1)
par(op)
#####################################################
## Chisquare probability plot for bivariate normality
#####################################################
mu <- colMeans(MacdonellDF)
sigma <- var(MacdonellDF)
Dsq <- mahalanobis(MacdonellDF, mu, sigma)
Q <- qchisq(1:3000/3000, 2)
plot(Q, sort(Dsq), xlab="Chisquare (2) quantile", ylab="Squared distance")
abline(a=0, b=1, col="red", lwd=2)
# }
Run the code above in your browser using DataLab