#
# Epage 
#
# 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)
}

# simulation example
set.seed(1013185)
case1 <- data.frame(x1=rnorm(100),x2=rnorm(100),
                    x3=rnorm(100),x4=rnorm(100))
e <- rnorm(100)
case1$y <- 1 + case1$x1 + case1$x2 + e
m1 <- lm(y~x1+x2+x3+x4,data=case1)


X <- as.matrix(case1[,-5]) # change from data.frame to a matrix, drop y
Var2 <- matrix(c(1,   0, .95,   0,
                 0,   1,   0,-.95,
               .95,   0,   1,   0,
                 0,-.95,   0,   1), ncol=4)
s1 <- chol(Var2)  # cholesky factor of Var2
X <- X %*% s1
dimnames(X)[[2]] <- paste("x",1:4,sep="")
case2 <- data.frame(X)
case2$y <- 1 + case2$x1 + case2$x2 + e
m2 <- lm(y~x1+x2+x3+x4,data=case2)

sum1 <- summary(m1)
sum2 <- summary(m2)
varb1 <- vcov(m1)
varb2 <- vcov(m2)
sum1
100*varb1
sum2
100*varb2

# Case 2, but with (much) larger sample size
m <- 1100
X2 <- matrix(rnorm(4*m),ncol=4) %*% s1
if(is.null(version$language) == FALSE) dimnames(X2)[[2]] <- paste("x",1:4,sep="")
beta <- c(1,1,0,0)
e <- rnorm(m)
Y2 <- 1 + X2 %*% beta + e
#  The following doesn't work in Splus.  I don't know why:
require(xtable)
if(is.null(version$language) == FALSE) {
  m3 <- lm(Y2 ~ X2) 
  sum3 <- summary(m3)
  varb3 <- sum3$sigma^2 * sum3$cov.unscaled
  if(is.null(version$language) == FALSE) xtable(sum3)  else sum3
  if(is.null(version$language) == FALSE) xtable(m*varb3) else m*varb3
}
# Case 2, but fitting only with the variables x1 and x4

# y ~ x1 + x4 => everything is significant
m2a <- update(m2,~.-x2-x3)
sum2a <- summary(m2a)
if(is.null(version$language) == FALSE) xtable(sum2a) else sum2a

# y ~ x3 + x4 => everything is significant
m2b <- update(m2,~.-x1-x2)
sum2b <- summary(m2b)
if(is.null(version$language) == FALSE) xtable(sum2b) else sum2b
sum2b