#
#  MATH3330 Section B    October 2001
#  Week7b_Coffee_Part2.SSC
#
#

# Read data.frame in Week5b_Coffee_Part1.SSC

summary(coffee)
plot( coffee , ask = F)


panel.xyell <- function(x, y, ...) { # better name
	panel.xyplot(x,y,...)
	data.ell(x,y,...)
}

pairs( coffee[,-1] , panel = panel.xyell )
cor(( coffee[,-1] ) )


# Coffee alone

fit.c <- lm ( Heart ~ Coffee, data = coffee )
summary( fit.c )

# 95% confidence interval for the effect of Coffee

coef(fit.c)[2]
coef(fit.c)["Coffee"]

t.val <- qt( .975, df = 18 ) 

summary( fit.c)$coef
names( summary( fit.c))
sd <- summary( fit.c )$coef [ 'Coffee','Std']

ci.width <- t.val * sd

lower.limit <- coef(fit.c)[2] - ci.width
upper.limit <- coef(fit.c)[2] + ci.width

list( lower.limit, upper.limit )

# writing a function to do all this

names(fit.c)
fit.c$df.residual

names(summary(fit.c))
summary(fit.c)$coef


ci <- function( fit, term , conf = .95) {
	df <- fit$df.residual
	co <- coef(fit)[term]
	nam <- names(coef(fit)[term])
	sd <- summary(fit)$coef[term,'Std']
	t.val <- qt( 1 - (1-conf)/2, df = df )
	data.frame( name = nam, conf = conf, value = co, SE = sd,
		half.width = sd*t.val, lower.bound = co - sd*t.val,
		upper.bound = co + sd*t.val)
}

names(coef(fit.c)[2])

ci(fit.c, 'Coffee')			# 95% is default
ci(fit.c, 'Coffee', .99)

#
#  Multiple regression
#

fit.cs <- lm ( Heart ~ Coffee + Stress, coffee )

summary( fit.cs )

ci( fit.cs, 'Coffee' )
ci( fit.cs, 'Stress' )

rbind( ci( fit.c, 'Coffee' ), ci(fit.cs , 'Coffee'))