### plotting the likelihoods for the regression parameters ###

### kind of like Bock's pp. 201-204, except ###
### using (asymptotic) normal instead of exact F, and ###
### (I think Bock may also have lost a sign on the covariance) ###

### the quartz command is unique to Mac OS X ###
### where "quartz" is the name of the windowing system ###

quartz(width=10, height=10)

### this opens a graphics window 720 x 720 pixels (10 * 72 dpi) ###

### below sets it up to make a 2 x 2 array of plots ###

par(mfrow=c(2,2))

fBnormal <- function(x,y) {
   ### what's below is the kernel of the bivariate normal density ###
   ### (no normalization), because we can figure out its ordinate ###
   ### most easily for various included percentages ###
   ### due to the unique factoid about 2df chi-squares that ###
   ### the ordinate of the kernel of the bivariate normal density ###
   ### is equal to 1 - alpha for the chi-square ###
   NormalDensKernel <- exp(-0.5*((((x-xmean)/xsd)^2
                                - 2*r*((x-xmean)/xsd)*((y-ymean)/ysd)
                                + ((y-ymean)/ysd)^2))/(1 - r^2))
 }

### plot for betas 1-2 ###
x <- seq((beta.Poggendorff[2]-3.5*SEbeta.Poggendorff[2]),
		 (beta.Poggendorff[2]+3.5*SEbeta.Poggendorff[2]), length=100)
y <- seq((beta.Poggendorff[1]-3.5*SEbeta.Poggendorff[1]),
		 (beta.Poggendorff[1]+3.5*SEbeta.Poggendorff[1]), length=100)
xmean <- beta.Poggendorff[2]
ymean <- beta.Poggendorff[1]

xsd <- SEbeta.Poggendorff[2]
ysd <- SEbeta.Poggendorff[1]
r <- R.Poggendorff[2,1]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.05, 0.01),
   labels = c("95%", "99%"),
    main = "Probability Contours, Poggendorff beta 1 and 2",
    xlab = "Poggendorff beta 2",
	ylab = "Poggendorff beta 1")

### plot the data upper right ###
plot(Table41Data$Age, y.Poggendorff, ylim = c(0,2))
lines(Table41Data$Age, yhat.Poggendorff)


### plot for betas 1-3 ###
x <- seq((beta.Poggendorff[3]-3.5*SEbeta.Poggendorff[3]),
		 (beta.Poggendorff[3]+3.5*SEbeta.Poggendorff[3]), length=100)
y <- seq((beta.Poggendorff[1]-3.5*SEbeta.Poggendorff[1]),
		 (beta.Poggendorff[1]+3.5*SEbeta.Poggendorff[1]), length=100)
xmean <- beta.Poggendorff[3]
ymean <- beta.Poggendorff[1]

xsd <- SEbeta.Poggendorff[3]
ysd <- SEbeta.Poggendorff[1]
r <- R.Poggendorff[3,1]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.05, 0.01),
   labels = c("95%", "99%"),
    main = "Probability Contours, Poggendorff beta 1 and 3",
    xlab = "Poggendorff beta 3",
	ylab = "Poggendorff beta 1")

### plot for betas 2-3 ###
x <- seq((beta.Poggendorff[3]-3.5*SEbeta.Poggendorff[3]),
		 (beta.Poggendorff[3]+3.5*SEbeta.Poggendorff[3]), length=100)
y <- seq((beta.Poggendorff[2]-3.5*SEbeta.Poggendorff[2]),
		 (beta.Poggendorff[2]+3.5*SEbeta.Poggendorff[2]), length=100)
xmean <- beta.Poggendorff[3]
ymean <- beta.Poggendorff[2]

xsd <- SEbeta.Poggendorff[3]
ysd <- SEbeta.Poggendorff[2]
r <- R.Poggendorff[3,2]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.05, 0.01),
   labels = c("95%", "99%"),
    main = "Probability Contours, Poggendorff beta 2 and 3",
    xlab = "Poggendorff beta 3",
	ylab = "Poggendorff beta 2")


### anticipating material we'll discuss more when we get to ###
### our "confidence envelopes" paper, do a ###
### 2nd paneled set with points and 235 (magic number) lines ###

nlines <- 235

library(MASS)

### mvrnorm() is in the MASS library ###
### generates multivariate normal data ###

SampleOfBetas <- mvrnorm(nlines, beta.Poggendorff, Vbeta.Poggendorff)

quartz(width=10, height=10)
par(mfrow=c(2,2))

### for betas 1-2 ###
x <- seq((beta.Poggendorff[2]-3.5*SEbeta.Poggendorff[2]),
		 (beta.Poggendorff[2]+3.5*SEbeta.Poggendorff[2]), length=100)
y <- seq((beta.Poggendorff[1]-3.5*SEbeta.Poggendorff[1]),
		 (beta.Poggendorff[1]+3.5*SEbeta.Poggendorff[1]), length=100)
xmean <- beta.Poggendorff[2]
ymean <- beta.Poggendorff[1]

xsd <- SEbeta.Poggendorff[2]
ysd <- SEbeta.Poggendorff[1]
r <- R.Poggendorff[2,1]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.05, 0.01),
   labels = c("95%", "99%"),
    main = "Probability Contours, Poggendorff beta 1 and 2",
    xlab = "Poggendorff beta 2",
	ylab = "Poggendorff beta 1")
points(SampleOfBetas[,2], SampleOfBetas[,1], pch=20, col="cyan")

### data upper right ###
plot(Table41Data$Age, y.Poggendorff, ylim = c(0,2))
for (i in 1:nlines) {
	ytemp <- X %*% as.matrix(SampleOfBetas[i,])
	lines(Table41Data$Age, ytemp, col="cyan")
}
lines(Table41Data$Age, yhat.Poggendorff)
points(Table41Data$Age, y.Poggendorff)

### for betas 1-3 ###
x <- seq((beta.Poggendorff[3]-3.5*SEbeta.Poggendorff[3]),
		 (beta.Poggendorff[3]+3.5*SEbeta.Poggendorff[3]), length=100)
y <- seq((beta.Poggendorff[1]-3.5*SEbeta.Poggendorff[1]),
		 (beta.Poggendorff[1]+3.5*SEbeta.Poggendorff[1]), length=100)
xmean <- beta.Poggendorff[3]
ymean <- beta.Poggendorff[1]

xsd <- SEbeta.Poggendorff[3]
ysd <- SEbeta.Poggendorff[1]
r <- R.Poggendorff[3,1]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.05, 0.01),
   labels = c("95%", "99%"),
    main = "Probability Contours, Poggendorff beta 1 and 3",
    xlab = "Poggendorff beta 3",
	ylab = "Poggendorff beta 1")
points(SampleOfBetas[,3], SampleOfBetas[,1], pch=20, col="cyan")

### for betas 2-3 ###
x <- seq((beta.Poggendorff[3]-3.5*SEbeta.Poggendorff[3]),
		 (beta.Poggendorff[3]+3.5*SEbeta.Poggendorff[3]), length=100)
y <- seq((beta.Poggendorff[2]-3.5*SEbeta.Poggendorff[2]),
		 (beta.Poggendorff[2]+3.5*SEbeta.Poggendorff[2]), length=100)
xmean <- beta.Poggendorff[3]
ymean <- beta.Poggendorff[2]

xsd <- SEbeta.Poggendorff[3]
ysd <- SEbeta.Poggendorff[2]
r <- R.Poggendorff[3,2]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.05, 0.01),
   labels = c("95%", "99%"),
    main = "Probability Contours, Poggendorff beta 2 and 3",
    xlab = "Poggendorff beta 3",
	ylab = "Poggendorff beta 2")
points(SampleOfBetas[,3], SampleOfBetas[,2], pch=20, col="cyan")



### a third paneled set has the 95% contours for the 3df? ordinate at 0.02 ###
### (see excel) ###
nlines <- 235
SampleOfBetas <- mvrnorm(nlines, beta.Poggendorff, Vbeta.Poggendorff)

quartz(width=10, height=10)
par(mfrow=c(2,2))

### for betas 1-2 ###
x <- seq((beta.Poggendorff[2]-3.5*SEbeta.Poggendorff[2]),
		 (beta.Poggendorff[2]+3.5*SEbeta.Poggendorff[2]), length=100)
y <- seq((beta.Poggendorff[1]-3.5*SEbeta.Poggendorff[1]),
		 (beta.Poggendorff[1]+3.5*SEbeta.Poggendorff[1]), length=100)
xmean <- beta.Poggendorff[2]
ymean <- beta.Poggendorff[1]

xsd <- SEbeta.Poggendorff[2]
ysd <- SEbeta.Poggendorff[1]
r <- R.Poggendorff[2,1]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.02),
   labels = c("95%"),
    main = "3 df Probability Contours, Poggendorff beta 1 and 2",
    xlab = "Poggendorff beta 2",
	ylab = "Poggendorff beta 1")
points(SampleOfBetas[,2], SampleOfBetas[,1], pch=20, col="cyan")

### data upper right ###
plot(Table41Data$Age, y.Poggendorff, ylim = c(0,2))
for (i in 1:nlines) {
	ytemp <- X %*% as.matrix(SampleOfBetas[i,])
	lines(Table41Data$Age, ytemp, col="cyan")
}
lines(Table41Data$Age, yhat.Poggendorff)
points(Table41Data$Age, y.Poggendorff)
for (i in 1:length(Table41Data$Age)) {
 lines(c(Table41Data$Age[i],Table41Data$Age[i]),
			c((yhat.Poggendorff[i]+1.96*(sqrt(sigmahatPooled.Poggendorff/16))),
			(yhat.Poggendorff[i] - 1.96*(sqrt(sigmahatPooled.Poggendorff/16)))))
}


### for betas 1-3 ###
x <- seq((beta.Poggendorff[3]-3.5*SEbeta.Poggendorff[3]),
		 (beta.Poggendorff[3]+3.5*SEbeta.Poggendorff[3]), length=100)
y <- seq((beta.Poggendorff[1]-3.5*SEbeta.Poggendorff[1]),
		 (beta.Poggendorff[1]+3.5*SEbeta.Poggendorff[1]), length=100)
xmean <- beta.Poggendorff[3]
ymean <- beta.Poggendorff[1]

xsd <- SEbeta.Poggendorff[3]
ysd <- SEbeta.Poggendorff[1]
r <- R.Poggendorff[3,1]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.02),
   labels = c("95%"),
    main = "3 df Probability Contours, Poggendorff beta 1 and 3",
    xlab = "Poggendorff beta 3",
	ylab = "Poggendorff beta 1")
points(SampleOfBetas[,3], SampleOfBetas[,1], pch=20, col="cyan")

### for betas 2-3 ###
x <- seq((beta.Poggendorff[3]-3.5*SEbeta.Poggendorff[3]),
		 (beta.Poggendorff[3]+3.5*SEbeta.Poggendorff[3]), length=100)
y <- seq((beta.Poggendorff[2]-3.5*SEbeta.Poggendorff[2]),
		 (beta.Poggendorff[2]+3.5*SEbeta.Poggendorff[2]), length=100)
xmean <- beta.Poggendorff[3]
ymean <- beta.Poggendorff[2]

xsd <- SEbeta.Poggendorff[3]
ysd <- SEbeta.Poggendorff[2]
r <- R.Poggendorff[3,2]

z <- outer(x, y, fBnormal)

contour(x, y, z, 
   levels = c(0.02),
   labels = c("95%"),
    main = "3 df Probability Contours, Poggendorff beta 2 and 3",
    xlab = "Poggendorff beta 3",
	ylab = "Poggendorff beta 2")
points(SampleOfBetas[,3], SampleOfBetas[,2], pch=20, col="cyan")