# Example: comparing two means 

Anor <- read.table("http://stat4ds.rwth-aachen.de/data/Anorexia.dat",
                   header=TRUE)
# Get difference post-pre treatment for the group cb and c
cogbehav <- Anor$after[Anor$therapy == "cb"] -
  Anor$before[Anor$therapy == "cb"]
control <- Anor$after[Anor$therapy == "c"] -
  Anor$before[Anor$therapy == "c"]

# Get the 95% CI via t.test function
res <- t.test(cogbehav, control, var.equal = TRUE, conf.level = 0.95)
res$conf.int


n1 <- length(cogbehav)
n2 <- length(control)
v1 <- var(cogbehav)
v2 <- var(control)

s2 <- ((n1 - 1) * v1 + (n2 - 1) * v2)/(n1 + n2 - 2)

CI <- mean(cogbehav) - mean(control) +
  c(-1,1) * qt(0.975, df = n1 + n2 - 2) *  sqrt(s2 * (1/n1 + 1/n2))
CI


n1 <- length(cogbehav)
n2 <- length(control)
summary(lm1)$coefficients[2,1] + 
  qt(0.975, df = n1 + n2 -2) *
  summary(lm1)$coefficients[2,2]
summary(lm1)$coefficients[2,1] - 
  qt(0.975, df = n1 + n2 -2) *
  summary(lm1)$coefficients[2,2]


res.two <- t.test(cogbehav, control, var.equal = TRUE)
res.two


res.one <- t.test(cogbehav, control, var.equal = TRUE, alternative = "greater")
res.one


s <- sqrt(s2 * (1/n1 +1/n2))
testStat <- (mean(cogbehav) - mean(control))/s
## two-sided
2 * pt(testStat, df = n1+n2-2, lower = FALSE )
summary(lm1)$coefficients[2,4] 

## one-sided
pt(testStat, df = n1+n2-2, lower = FALSE )


library(RColorBrewer)
plotclr <- brewer.pal(6, "YlOrRd")

curve(dt(x, n1 + n2 - 2), xlim = c(-5, 5), ylim = c(0, 0.4),
      main = "p-values and rejection region", col = "blue",
      lwd = 2, xlab = "x-y",  ylab = expression(t[13]),  yaxs="i")
cord.x <- c(qt(0.95, n1 + n2 - 2),seq(qt(0.95, n1 + n2 - 2), 5, 0.01), 5)
cord.y <- c(0, dt(seq(qt(0.95, n1 + n2 - 2), 5, 0.01), 13), 0)
polygon(cord.x, cord.y, col = plotclr[3], border = NA )

abline(v = res.one$statistic, lty = 2, lwd = 2, col = "red")
text(0, 0.2, paste("Accept", expression(H0)))
text(2.7, 0.08, paste("Reject", expression(H0)))
text(as.double(res.one$statistic) - 0.15, 0.02, "t", 
     col = "red", cex = 1.2)