#
# Epage 92
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

library(alr3)

###################################################################
# BGS for BGSgirls only
if(is.null(version$language) == FALSE) data(BGSgirls)
use <- c(2,4,8,12)
# Fig. 4.1 "../EPS-figs/bgs31.eps"
if(is.null(version$language) == FALSE) pairs(BGSgirls[,use],gap=0.4,cex.labels=1.5)
if(is.null(version$language) == TRUE) pairs(BGSgirls[,use],labels=1.5)

attach(BGSgirls)

# the direct model "Model 1" in the text:
m_direct = lm( Soma ~ WT2 + WT9 + WT18 )
s_direct = summary(m_direct)

# print the correlation matrix
print(cor(BGSgirls[,use]),digits=4)
# create some linear combinations
BGS       <- BGSgirls # Required for SPlus, not needed for R
BGS$DW9   <- BGS$WT9-BGS$WT2
BGS$DW18  <- BGS$WT18-BGS$WT9
BGS$DW218 <- BGS$WT18-BGS$WT2
m1 <- if(is.null(version$language) == FALSE)
 lm(Soma~WT2+WT9+WT18+DW9+DW18,BGS) else
 lm(Soma~WT2+WT9+WT18+DW9+DW18,BGS,singular.ok=TRUE)
s1 <- summary(m1)
print(s1,signif.stars=FALSE)
m2 <- lm(Soma~WT2+DW9+DW18+WT9+WT18,BGS,singular.ok=TRUE)
m3 <- lm(Soma~WT2+WT9+WT18+DW9+DW18,BGS,singular.ok=TRUE)
if(is.null(version$language) == FALSE) {
  ans <- data.frame(m1=coef(m1),m2=coef(m2)[c(1,2,5,6,3,4)],m3=coef(m1))
  print(ans,digits=3) } else {
  print(coef(m1)); print(coef(m2)); print(coef(m3)) }
s2 <- summary(m2)
print(s2,signif.stars=FALSE)

s3 <- summary(m3)

# create the linear combinations desired for this problem:
#
AVE   <- ( WT2 + WT9 + WT18 ) / 3
LIN   <- WT18-WT9
QUAD  <- WT2 - 2 * WT9 + WT18

m4 <- lm(Soma~AVE+LIN+QUAD) 
s4 <- summary(m4)
print(s4,signif.stars=FALSE)

# print R^2 and sigma^2 the residual sigmas
#
c(s_direct$r.squared,s1$r.squared,s2$r.squared,s3$r.squared)

c(s_direct$sigma,s1$sigma,s2$sigma,s3$sigma)