### Class examples from Week 7: Comparing Groups
### Ken Benoit Mar 2010

## Levin and Fox example p146
a1 <- c(82,83,82,80,83)
b1 <- c(78,77,76,78,76)
mean(a1)
mean(b1)
sd(a1)
sd(b1)
# plot densities, mean, and 95% CI
plot(density(a1),col="blue",xlim=c(70,90), main="First set")
lines(density(b1),col="red")
abline(v=mean(a1),col="blue")
abline(v=t.test(a1)$conf.int, col="blue", lty="dashed")
abline(v=mean(b1),col="red")
abline(v=t.test(b1)$conf.int, col="red", lty="dashed")

## Levin and Fox example p147
quartz()  # note: comment this line out for non-Mac users!
a2 <- c(90,98,63,74,85)
b2 <- c(70,90,91,56,78)
mean(a2)
mean(b2)
sd(a2)
sd(b2)
# plot densities, mean, and 95% CI
plot(density(a2),col="blue",xlim=c(30,120), main="Second set")
lines(density(b2),col="red")
abline(v=mean(a2),col="blue")
abline(v=t.test(a2)$conf.int, col="blue", lty="dashed")
abline(v=mean(b2),col="red")
abline(v=t.test(b2)$conf.int, col="red", lty="dashed")

# p146 and p147 examples: plot both groups on a "strip chart"
par(mar=c(3,6,2,2))
stripchart(rev(list(a1,b1,a2,b2)),
           method="jitter", jitter=.05, offset=.2,
           group.names=rev(c("Method A1", "Method B1", "Method A2",
             "Method B2")),
           las=1, col=c("red","blue"), pch="o")
abline(v=mean(a1), col="blue", lty="dashed")
abline(v=mean(b1), col="red", lty="dashed")


## Levin and Fox example of difference between means pp162-165
x1 <- c(8,7,6,10,2,1,4,3,5,6,7,5,9,8,10,6,6,7,4,5,9,10,
        2,4,3,5,4,8,7,4,9,8,9,10,6)
x2 <- c(4,3,1,7,6,9,10,4,3,6,4,3,6,8,6,5,7,3,4,6,9,8,4,
        5,8,2,7,1,5,6,4,8,7,6,4)
# step 1: mean of each sample
(x1.mean <- mean(x1))
(x2.mean <- mean(x2))
# step 2: variance of each sample (using population formula)
varpop <- function(x) { sum(x^2)/length(x) - mean(x)^2 }
varpop(x1)
varpop(x2)
# step 3: find standard error of difference
serrdiff <- function(x1,x2) {
  n1 <- length(x1)
  n2 <- length(x2)
  sqrt((n1*varpop(x1) + n2*varpop(x2))/(n1+n2-2) * ((n1+n2)/(n1*n2)))
}
(s <- serrdiff(x1,x2))
# step 4: compute t-ratio
(t.calc <- (x1.mean - x2.mean)/s)
# step 5: determine critical value for t
df <- length(x1) + length(x2) - 2
(t.crit <- qt(1-.05/2, df))  # "correct" df not available in LF table
qt(1-.05/2, 60)  # same value as L&F p164
# step 6: compare the calculated and table t-values
abs(t.calc) > t.crit
## using the R shortcut of t.test():
t.test(x1, x2, var.equal=TRUE)


## Levin and Fox example of difference between means
## of same sample measured twice pp166-168
xbefore <- c(2,1,3,3,1,4)
xafter <-  c(1,2,1,1,2,1)
(D <- xbefore - xafter)
sum(D^2)
(N <- length(xbefore)) # note: length of before must = length of after
# step 1: mean for each point in time
(xbefore.mean <- mean(xbefore))
(xafter.mean <- mean(xafter))
# step 2: standard deviation for difference b/t times 1 and 2
(sD <- sqrt(sum(D^2)/N - (xbefore.mean - xafter.mean)^2))
# step 3: standard error of mean difference
(sDbar <- sD/sqrt(N-1))
# step 4: translate into units of SE of mean difference
(t.calc <- (xbefore.mean - xafter.mean)/sDbar)
# step 5: compute degrees of freedom
(df <- N-1)
# step 6: compare calculated t-ratio with critical t-ratio 
(t.crit <- qt(1-.05/2, df))
abs(t.calc) > t.crit
## using the R shortcut of t.test():
t.test(xbefore, xafter, var.equal=TRUE, paired=TRUE)


## Levin and Fox example of difference between proportions pp169-171
n1 <- 180
f1 <- 81
n2 <- 150
f2 <- 48
# step 1: compute two sample proportions and combined sample proportion
(p1 <- (f1/n1))
(p2 <- (f2/n2))
(pstar <- (n1*p1 + n2*p2)/(n1+n2))
# step 2: compute std error of difference
(sdiff <- sqrt(pstar*(1-pstar) * (n1+n2)/(n1*n2)))
# step 3: translate into std error of the difference
(z.calc <- (p1 - p2)/(sdiff))
# step 4: compare calculated z with critical z
(z.crit <- qnorm(1-.05/2))
abs(z.calc) > z.crit
## using the R shortcut of prop.test():
prop.test(c(f1,f2), c(n1,n2))