# Hald data
 
hald <- scan("")
  7.0000000      26.000000      6.0000000      60.000000      78.500000    
  1.0000000      29.000000      15.000000      52.000000      74.300000    
  11.000000      56.000000      8.0000000      20.000000      104.30000    
  11.000000      31.000000      8.0000000      47.000000      87.600000    
  7.0000000      52.000000      6.0000000      33.000000      95.900000    
  11.000000      55.000000      9.0000000      22.000000      109.20000    
  3.0000000      71.000000      17.000000      6.0000000      102.70000    
  1.0000000      31.000000      22.000000      44.000000      72.500000    
  2.0000000      54.000000      18.000000      22.000000      93.100000    
  21.000000      47.000000      4.0000000      26.000000      115.90000    
  1.0000000      40.000000      23.000000      34.000000      83.800000    
  11.000000      66.000000      9.0000000      12.000000      113.30000    
  10.000000      68.000000      8.0000000      12.000000      109.40000    


hald <- matrix(hald,ncol=5,byrow=T)
hald <- data.frame(hald)
names(hald) <-  c("ta","ts","taf","ds","heat") # Assign names to the columns

# Hald data
# ta        = '% tricalcium aluminate'
# ts        = '% tricalcium silicate'        
# taf       = '% tetracalcium alumino ferrite'
# ds        = '%  dicalcium silicate'     
# heat      = 'heat (calories /g cement)';

# Let's fit a regression predicting heat from ta and ts.
# Here's the R function that does the regression

lmhald <- lm(heat ~ ta + ts,data=hald)

lmhald  # Prints the result of the regression

summary(lmhald) # Gives more results from the regression

attributes(lmhald)
lmhald$coef

names(summary(lmhald))

summary(lmhald)$sigma

summary(lmhald)$fstatistic

summary(lmhald)$cov.unscaled  # Gives estimated covariance matrix of beta-hat

# Same analysis with matrices
# Test H_0: beta_1 = beta_2 = 0

xmat <- cbind(rep(1,13),hald$ta,hald$ts)
xpx <- t(xmat) %*% xmat
xpx
xpxinv <- solve(xpx)
xpxinv
cmat <- cbind(c(0,1,0),c(0,0,1))
t(cmat)

ctxpxinvc <- t(cmat) %*% xpxinv %*% cmat
ctxpxinvc
middle <- solve(ctxpxinvc)
middle

fnum <- t( t(cmat)%*%lmhald$coef) %*% middle %*% ( t(cmat)%*%lmhald$coef) /2
fnum
fstat <- fnum / (summary(lmhald)$sigma)^2
fstat

1-pf(fstat,2,10) # gives p-value

# Creates the contours of a 95% confidence region for C'beta
# CI is all (beta1,beta2) such that Fstat < critpt

critpt <- qf(.95,2,10)  # gives (1-alpha) critical point
critpt

# We want to draw the ellipse 
# t(C beta - C \hat{\beta}) middle  (C'beta - C'\hat{\beta})/(2 MSE) = critpt

eigmid <- eigen(middle)  # Find eigenvalues, vectors of middle
eigmid

# Standard form for ellipse: x^2/a^2 + y^2/b^2 =1
# Major axis: a= 1/sqrt(390.1705/ (critpt*5.79*2) ) = 0.3489541
# Minor axis: b= 1/sqrt(2930.7526/ (critpt*5.79*2)) = 0.1273227
# Eccentricity = c/a = sqrt(a^2 – b^2)/a = .931
# Direction of major axis is [-.995, .0993]’
# Direction of minor axis is [-.0993, -.995]’

# Axes length, rotation

axlength <- sqrt(critpt* (summary(lmhald)$sigma)^2 * 2/eigmid$val)

plot.ellipse <- function(mu,amat,const)  {
#
# plot ellipse given by (x - mu)' amat (x-mu) = const
# 
# create a circle: y' I y = 1, where y = (y1[i], y2[i])'
#
 ix <- 0:100
 y1 <- cos((ix * 2 * pi)/100)
 y2 <- sin((ix * 2 * pi)/100)
#
# now get the lengths and rotation for the axes 
#
 eiga <- eigen(amat)
# if (min(eiga$val) <= 0) stop("not all eigenvalues are positive")
 axlength <- sqrt(const/eiga$val)
 wmat <- cbind(y1 * axlength[1], y2 * axlength[2]) 
#
# each row of wmat is w = D^{-1/2} sqrt(const) y;
# we now have w' D w = const
# now rotate the adjusted ellipse
# Let x = row of xmat; want x' P D P'x = const, so let x' = w' P'
#
 xmat <- wmat %*% t(eiga$vec )
#
# now shift the points by mu
#
 xmat <- xmat + matrix(mu,ncol=2,nrow=length(y1),byrow=T)
 plot(xmat[, 1], xmat[, 2], type = "n",xlab="beta1",ylab="beta2",
    main="95% confidence region for beta")
 lines(xmat[, 1], xmat[, 2]) 
 points(mu[1],mu[2],pch=3) # plot mean
}  # End of function plot.ellipse

plot.ellipse(c(1.4683,.6623),middle,2*critpt* (summary(lmhald)$sigma)^2)


# Alternatively: use the contour function to draw the joint confidence region
# Note: very sloppy code and inefficient

x1 <- 1+(0:100)/100
x2 <- x1-1

z <- matrix(nrow=length(x1),ncol=length(x2))
for (i in 1:length(x1)) {
  for (j in 1:length(x2)) {
z[i,j] <- t(c(x1[i]-1.4683, x2[j]-.6623)) %*% 
   middle %*% c(x1[i]-1.4683, x2[j]-.6623) /(2*5.79)
}}
contour(x1,x2,z,levels=critpt,xlab="beta1",ylab="beta2")
points(1.4683,.6623)

# Condition indices in R

xmat <- cbind(rep(1,13),hald$ta,hald$ts)
xpx <- t(xmat) %*% xmat # X'X
xpxdiag <- diag(xpx) # Extract diagonal entries of X'X
xscale <- xmat %*% diag(1/sqrt(xpxdiag))
xscale
xpxs <- t(xscale) %*% xscale

eigxpxs <- eigen(xpxs)
eigxpxs
sqrt(eigxpxs$values[1]/eigxpxs$values) # Condition indices

# Variance decomposition proportions 
# V(betahat_k) = sigma^2 \sum_{j=1}^p V_{kj}^2/(jth eigenvalue)
# V(betahat) = sigma^2 V D^{-2} V', V is matrix of eigenvalues

# To get variance decomposition proportions in SAS format, want the entries
# in column k to be vector of V_{kj}^2/(jth eigenvalue) for j=1,2,3
# So look at transpose of eigenvector matrix

vardecomp <- t(eigxpxs$vectors)^2 / matrix(eigxpxs$values,nrow=nrow(xpxs),ncol=ncol(xpxs),byrow=F)
vardecsum <- apply(vardecomp,2,sum)
vardecsum <- matrix(vardecsum,nrow=nrow(vardecomp),ncol=ncol(vardecomp),byrow=T)
vardecompprop <- vardecomp/vardecsum
vardecompprop

# To get collinoint, repeat with non-intercept columns after centering

xmat <- xmat[,2:3]
xmean <- apply(xmat,2,mean)
xmat <- sweep(xmat,2,xmean)

xpx <- t(xmat) %*% xmat
xpxdiag <- diag(xpx)
xscale <- xmat %*% diag(1/sqrt(xpxdiag))
xscale
xpxs <- t(xscale) %*% xscale

eigxpxs <- eigen(xpxs)
eigxpxs
sqrt(eigxpxs$values[1]/eigxpxs$values) # Condition indices

# Variance decomposition proportions 


vardecomp <- t(eigxpxs$vectors)^2 / matrix(eigxpxs$values,nrow=nrow(xpxs),ncol=ncol(xpxs),byrow=F)
vardecsum <- apply(vardecomp,2,sum)
vardecsum <- matrix(vardecsum,nrow=nrow(vardecomp),ncol=ncol(vardecomp),byrow=T)
vardecompprop <- vardecomp/vardecsum
vardecompprop