#
# Note one could also use the t.test function in the stats package 
#
# 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(DAAG)
library(lattice)

# Part (a):
#
# trt = 0 => control (not treated)
# trt = 1 => treated
#
m_0 = mean( nswdemo[ nswdemo$trt == 0, ]$re75 )
s_0 = sd( nswdemo[ nswdemo$trt == 0, ]$re75 )
n_0 = length( nswdemo[ nswdemo$trt == 0, ]$re75 )

m_1 = mean( nswdemo[ nswdemo$trt == 1, ]$re75 )
s_1 = sd( nswdemo[ nswdemo$trt == 1, ]$re75 )
n_1 = length( nswdemo[ nswdemo$trt == 1, ]$re75 )

print( "Part (a): Control group confidence interval:" )
print( m_0 + qt( c(0.025,0.975), n_0-1 ) * s_0 / sqrt(n_0) )

print( "Part (a): Treated group confidence interval:" )
print( m_1 + qt( c(0.025,0.975), n_1-1 ) * s_1 / sqrt(n_1) )

# Part (b):
#
# trt = 0 => control (not treated)
# trt = 1 => treated
#
m_0 = mean( nswdemo[ nswdemo$trt == 0, ]$re78 )
s_0 = sd( nswdemo[ nswdemo$trt == 0, ]$re78 )
n_0 = length( nswdemo[ nswdemo$trt == 0, ]$re78 )

m_1 = mean( nswdemo[ nswdemo$trt == 1, ]$re78 )
s_1 = sd( nswdemo[ nswdemo$trt == 1, ]$re78 )
n_1 = length( nswdemo[ nswdemo$trt == 1, ]$re78 )

print( "Part (b): Control group confidence interval:" )
print( m_0 + qt( c(0.025,0.975), n_0-1 ) * s_0 / sqrt(n_0) )

print( "Part (b): Treated group confidence interval:" )
print( m_1 + qt( c(0.025,0.975), n_1-1 ) * s_1 / sqrt(n_1) )

# The difference between in mean income between the treated and controls in 1978:
#
d_bar = m_0 - m_1 # the mean difference (if negative the treated group has a larger mean)
sed = sqrt( (s_0/sqrt(n_0))^2 + (s_1/sqrt(n_1))^2 ) # standard error of the difference

print( "Part (b): Confidence interval of the difference" )
print( d_bar + qt( c(0.025,0.975), n_0 + n_1 - 2 ) * sed )

# Plot a box and whiskers plot:
#

#postscript("../../WriteUp/Graphics/Chapter4/prob_1_bwplot.eps", onefile=FALSE, horizontal=FALSE)
nswdemo$trt = as.factor( nswdemo$trt )
bwplot( trt ~ re78, data=nswdemo, xlab="re78 value", ylab="Control (0); Treated (1)" )
#dev.off()