#
# Epage 137
#
# 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.
# 
#-----

# version$language == "R" for R version$language == NULL for SPlus

if(is.null(version$language) == FALSE) require(alr3) else library(alr3)

# Oehlert's cake baking data
#cakes <- read.table("data/cakes.txt", header=TRUE)
if(is.null(version$language) == FALSE) data(cakes)
attach(cakes)

# the full model with an interaction term X1:X2:
#
m1 <- lm(Y ~ X1 + X2 + I(X1^2) + I(X2^2) + X1:X2, data=cakes)
s1 <- summary( m1 )
s1 

#--
# the delta method for the estimates of tilde{X}_1 and tilde{X}_2
#-- 

# Extract the estimated coefficients:
#
b0  <- coef(m1)[1]
b1  <- coef(m1)[2]
b2  <- coef(m1)[3]
b11 <- coef(m1)[4]
b22 <- coef(m1)[5]
b12 <- coef(m1)[6]

# Extract the covariance matrix for the beta_hat estimates.  This is \sigma^2 (X' X)^{-1}
#
Var <- vcov(m1)

# Compute the needed analytic derivatives (the gradient) of \tidle{X}_1:
# 
x1m <- "-(2*b1*b22+b2*b12)/(4*b11*b22 - b12^2)"
x1m.expr <- parse(text=x1m)

# this will evaluate x1m at the estimated above beta's:
eval(x1m.expr) 

# this will compute the analytic derivatives:
derivs <- c(D(x1m.expr,"b0"),D(x1m.expr,"b1"),D(x1m.expr,"b2"),D(x1m.expr,"b11"),D(x1m.expr,"b22"),D(x1m.expr,"b12"))
derivs

# this will compute the numeric gradient at our estimated beta_hat: 
eval.derivs<-c(eval(D(x1m.expr,"b0")),eval(D(x1m.expr,"b1")),eval(D(x1m.expr,"b2")),eval(D(x1m.expr,"b11")),eval(D(x1m.expr,"b22")),eval(D(x1m.expr,"b12")))

# now we can evalaute the standard error of \tilde{X}_1:  
sqrt(t(eval.derivs) %*% Var %*% eval.derivs)

# or we can use the ALR3 package "delta.method" (but we need to specify the parameter names correctly):
delta.method(m1,"-(2*b1*b4+b2*b5)/(4*b3*b4 - b5^2)")


# Now lets test a more general model where each term can have a different value:
#

# first create a factor for the block:
cakes$X1X2 <- X1 * X2                     # put in a product factor 
cakes$B    <- factor( cakes$block, ordered=FALSE )
mg <- lm( Y ~ -1 + B + B:X1 + B:X2 + I(X1^2) + I(X2^2) + X1X2, data=cakes, na.action=na.omit)

# the more restrictive model is m1 from above

anova(mg,m1)