setwd("/Volumes/Ruby2/Class-G4two/Fall04/BockCh4/")


### compute the betas matrix-wise as per Bock p. 173 ###
X <- matrix(c(rep(1,length(Table41Data$Age)), Table41Data$Age, 
				Table41Data$Age^2), nrow = length(Table41Data$Age))
D <- diag(Table41Data$N, length(Table41Data$Age))
XprimeDX <- t(X) %*% D %*% X
XprimeDXinverse <- solve (XprimeDX)

y.Ponzo <- matrix(Table41Data$PonzoM, nrow = length(Table41Data$Age))
XprimeDy.Ponzo <- t(X) %*% D %*% y.Ponzo
beta.Ponzo <- XprimeDXinverse %*% XprimeDy.Ponzo
print(beta.Ponzo)

y.Poggendorff <- matrix(Table41Data$PoggendorffM, 
		nrow = length(Table41Data$Age))
XprimeDy.Poggendorff <- t(X) %*% D %*% y.Poggendorff
beta.Poggendorff <- XprimeDXinverse %*% XprimeDy.Poggendorff
print(beta.Poggendorff)


### fitted values and residuals as per Bock p. 174 ###
yhat.Ponzo <- X %*% beta.Ponzo
e.Ponzo <- y.Ponzo - yhat.Ponzo
Table412.Ponzo <- matrix(c(Table41Data$Age, y.Ponzo, yhat.Ponzo,
			e.Ponzo), nrow = length(Table41Data$Age))
print(Table412.Ponzo)
plot(Table41Data$Age, y.Ponzo, ylim = c(0,0.6))
lines(Table41Data$Age, yhat.Ponzo)

yhat.Poggendorff <- X %*% beta.Poggendorff
e.Poggendorff <- y.Poggendorff - yhat.Poggendorff
Table412.Poggendorff <- matrix(c(Table41Data$Age, y.Poggendorff, 
			yhat.Poggendorff, e.Poggendorff), nrow = length(Table41Data$Age))
print(Table412.Poggendorff)			
plot(Table41Data$Age, y.Poggendorff, ylim = c(0,2))
lines(Table41Data$Age, yhat.Poggendorff)


### error variance estimates, Bock pl 184ff ###
ssw.Ponzo <- sum((Table41Data$N - 1) * Table41Data$PonzoSD^2)
ssr.Ponzo <- t(beta.Ponzo) %*% t(X) %*% D %*% y.Ponzo
ssg.Ponzo <- t(y.Ponzo) %*% D %*% y.Ponzo
sse.Ponzo <- ssg.Ponzo - ssr.Ponzo
### or
sseAlt.Ponzo <- t(e.Ponzo) %*% D %*% e.Ponzo

### goodness of fit ###
F.Ponzo <- (sse.Ponzo/(length(Table41Data$Age) - ncol(X))) /
			(ssw.Ponzo/(sum(Table41Data$N) - length(Table41Data$Age)))
p.Ponzo <- 1 - pf(F.Ponzo, (length(Table41Data$Age) - ncol(X)), 
				(sum(Table41Data$N) - length(Table41Data$Age)))
print(c(F.Ponzo, p.Ponzo))

ssw.Poggendorff <- sum((Table41Data$N - 1) * Table41Data$PoggendorffSD^2)
ssr.Poggendorff <- t(beta.Poggendorff) %*% t(X) %*% D %*% y.Poggendorff
ssg.Poggendorff <- t(y.Poggendorff) %*% D %*% y.Poggendorff
sse.Poggendorff <- ssg.Poggendorff - ssr.Poggendorff
### or  ###
sseAlt.Poggendorff <- t(e.Poggendorff) %*% D %*% e.Poggendorff

### goodness of fit ###
F.Poggendorff <- (sse.Poggendorff/(length(Table41Data$Age) - ncol(X))) /
			(ssw.Poggendorff/(sum(Table41Data$N) - length(Table41Data$Age)))
p.Poggendorff <- 1 - pf(F.Poggendorff, (length(Table41Data$Age) - ncol(X)), 
				(sum(Table41Data$N) - length(Table41Data$Age)))
print(c(F.Poggendorff, p.Poggendorff))

### alternative variance estimates---Poggendorff only ###
sigmahatWithin.Poggendorff <- ssw.Poggendorff / (sum(Table41Data$N) - 
									length(Table41Data$Age))
print(sigmahatWithin.Poggendorff)
sigmahatPooled.Poggendorff <- (ssw.Poggendorff + sse.Poggendorff) /
		  ((sum(Table41Data$N) -  length(Table41Data$Age)) +
		 		(length(Table41Data$Age) - ncol(X)))
print(sigmahatPooled.Poggendorff)


### and now back to standard errors for the coefficients, combining ###
### Bock p. 185 and p. 175 ###
### note that the coefficients (betas) are linear combinations of the means ###
### beta = My., wehre M is ###
### M <- XprimeDXinverse %*% t(X) %*% D ###
### so, using the within veriance estimate,  ###
### the covariance matrix of the coefficients is ###

Vbeta.Poggendorff <- as.real(sigmahatWithin.Poggendorff) *
								XprimeDXinverse
print(Vbeta.Poggendorff)

### (the result near the bottom of p. 185 of Bock's book appears incorrect) ###
### (appears to be hand-multiplied by 0.01239 vs. 0.01029) ###

SEbeta.Poggendorff <- sqrt(diag(Vbeta.Poggendorff))
print(SEbeta.Poggendorff)

### correlation matrix among the coefficients ###

R.Poggendorff <- solve(diag(sqrt(diag(Vbeta.Poggendorff)), 
							nrow=nrow(Vbeta.Poggendorff), 
							ncol=nrow(Vbeta.Poggendorff))) %*%
				Vbeta.Poggendorff %*%
				solve(diag(sqrt(diag(Vbeta.Poggendorff)), 
							nrow=nrow(Vbeta.Poggendorff), 
							ncol=nrow(Vbeta.Poggendorff)))
							
print(R.Poggendorff)