## Nov 28, 2007
## AWH3
## this code provides an example of a multivariate regression

## example 7.8 in book
y <- matrix(c(1,-1,4,-1,3,2,8,3,9,2),5,2, byrow=T)
n <- nrow(y)
p <- ncol(y)

z <- cbind(rep(1,n), c(0,1,2,3,4))
r <- ncol(z)-1


## calculate beta.hat
beta.hat <- solve(t(z)%*%z)%*%t(z)%*%y
beta.hat

## calculate y.hat (the fitted values)
y.hat <- z%*%beta.hat
y.hat

## calculate epsilon.hat (the residuals)
eps.hat <- y-y.hat
eps.hat

## calculate the cross products of the residuals
## why is this matrix 2 by 2?  Because p=2!
eps.eps <- t(eps.hat)%*%eps.hat

## Thus Sigma.hat is:
Sigma.hat <- (1/(n-r-1))*eps.eps
Sigma.hat

## The variance for Beta
var.beta.11 <- diag(solve(t(z)%*%z)*Sigma.hat[1,1])
var.beta.22 <- diag(solve(t(z)%*%z)*Sigma.hat[2,2])

## student t-tests for the Betas: Ho: Beta = 0
t.1.obs <- beta.hat[,1]/sqrt(var.beta.11)
t.2.obs <- beta.hat[,2]/sqrt(var.beta.22)

## output
cat("Beta1:", beta.hat[,1], "Std Error:", sqrt(var.beta.11), "Obs t-stat:", t.1.obs, "p-value:",  2*(1-pt(abs(t.1.obs), n-(r+1))))
cat("Beta2:", beta.hat[,2], "Std Error:", sqrt(var.beta.22), "Obs t-stat:", t.2.obs, "p-value:",  2*(1-pt(abs(t.2.obs), n-(r+1))))


#########################################################
#########################################################
## Male and Female turtles revisited with multivariate regression

FT <- read.table("turtlesFemale.txt", header=T)
MT <- read.table("turtlesMale.txt", header=T)

FT
MT

p <- ncol(FT)
n.f <- nrow(FT)
n.m <- nrow(MT)
n <- n.f+n.m

y <- as.matrix(rbind(FT,MT))
y

## ok  specify the z matrix 
## f = 0, m =1
z <- cbind(rep(1, (n.f+n.m)), c(rep(0, n.f), rep(1, n.f))) 
z
r <- ncol(z)-1


## calculate beta.hat
beta.hat <- solve(t(z)%*%z)%*%t(z)%*%y
beta.hat

## calculate y.hat (the fitted values)
y.hat <- z%*%beta.hat
y.hat

## calculate epsilon.hat (the residuals)
eps.hat <- y-y.hat
eps.hat

## calculate the cross products of the residuals
## why is this matrix 2 by 2?  Because p=2!
eps.eps <- t(eps.hat)%*%eps.hat

## Thus Sigma.hat is:
Sigma.hat <- (1/(n-r-1))*eps.eps
Sigma.hat

Sigma.mle <- (1/n)*eps.eps
Sigma.mle

## The covariance for Beta
var.beta.11 <- diag(solve(t(z)%*%z)*Sigma.hat[1,1])
var.beta.22 <- diag(solve(t(z)%*%z)*Sigma.hat[2,2])
var.beta.33 <- diag(solve(t(z)%*%z)*Sigma.hat[3,3])

## student t-tests for the Betas
t.1.obs <- beta.hat[,1]/sqrt(var.beta.11)
t.2.obs <- beta.hat[,2]/sqrt(var.beta.22)
t.3.obs <- beta.hat[,3]/sqrt(var.beta.33)

## output
cat("Beta1:", beta.hat[,1], "Std Error:", sqrt(var.beta.11), "Obs t-stat:", t.1.obs, "p-value:",  2*(1-pt(abs(t.1.obs), n-(r+1))))
cat("Beta2:", beta.hat[,2], "Std Error:", sqrt(var.beta.22), "Obs t-stat:", t.2.obs, "p-value:",  2*(1-pt(abs(t.2.obs), n-(r+1))))
cat("Beta3:", beta.hat[,3], "Std Error:", sqrt(var.beta.33), "Obs t-stat:", t.3.obs, "p-value:",  2*(1-pt(abs(t.3.obs), n-(r+1))))

## from the individual t-tests looks like a difference but lets be a little more precise.

############################################################
############################################################
## Now re-estimate the model, but just for the intercept.
z.con <- matrix(rep(1, (n.f+n.m)), n.f+n.m,1)
z.con
r <- ncol(z.con)-1

## calculate beta.hat
beta.hat.con <- solve(t(z.con)%*%z.con)%*%t(z.con)%*%y
beta.hat.con

## calculate y.hat (the fitted values)
y.hat.con <- z%*%beta.hat.con
y.hat.con

## calculate epsilon.hat (the residuals)
eps.hat.con <- y-y.hat.con
eps.hat.con

## calculate the cross products of the residuals
## why is this matrix 2 by 2?  Because m=2!
eps.eps.con <- t(eps.hat.con)%*%eps.hat.con

## Thus Sigma.hat is:
Sigma.hat.con <- (1/(n-r-1))*eps.eps.con
Sigma.hat.con

Sigma.mle.con <- (1/n)*eps.eps.con
Sigma.mle.con

## calculate the likelihood ratio test using the mles for the Sigmas
alpha <- 0.05
Lambda.obs <- (det(Sigma.mle)/det(Sigma.mle.con))^(n/2)
LRT.obs <- -2*log(Lambda.obs)
LRT.obs

qchisq(1-0.05, 3)

# Thus reject the null hypothesis, thus the covariate for gender is important.