## Chapter 5
##
## Code to support Video Vignette
##
## Videos and supporting code are not a complete portrayal
## of Chapter content.

##
## a: two-sample Bernoulli/coin flips
##

## buddy's coin flips, from chapter 1
data <- c("H", "T", "H", "H", "H", "T", "T", "H", "T", "T", "H", "H", "H",
          "T", "T", "H", "H", "H", "H", "T", "H", "H", "H", "T", "H", "H", "T", 
          "H", "T", "T")

## convert to binary
y <- as.numeric(data == "H")
n <- length(y)

## summary
sy <- sum(y)        ## number of heads
ny <- length(y)     ## number of coin flips

## summary of second set of coin flips
sx <- 19
nx <- 53

## does are the two coins the same?

##
## via simulation
##

## if they are (H0) then we can combine the samples to estimate theta
that.pool <- (sy + sx)/(ny + nx)

## a natural statistic is the difference in MLEs
that.d <- sy/ny - sx/nx

## two separate copies of Chapter 1 code
N <- 100000                                 ## number of MC "trials"
That.ds <- rep(NA, N)                       ## storage
for(i in 1:N) {                             ## each MC iteration
  Sy <- sum(rbinom(ny, 1, that.pool))       ## could be combined ...
  Sx <- sum(rbinom(nx, 1, that.pool))       ## ... into a single sample
  That.ds[i] <- Sy/ny - Sx/nx               ## save sampled stat
}

## visualizing the sampling distribution
hist(That.ds, main="", xlim=c(-0.65, 0.65)) 
abline(v=c(that.d, -that.d), lty=1:2, col=2, lwd=2)
legend("topright", c("that.d", "reflect"), col=2, lty=1:2, lwd=2, bty="n")

## p-value calculation
pval.mc <- 2*mean(That.ds >= that.d)
pval.mc

## 
## the math way
##

se <- sqrt(that.pool*(1 - that.pool)*(1/ny + 1/nx))
z <- that.d/se
pval.z <- 2*pnorm(-abs(z))
pval.z

##
## b: two-sample Gaussian
##

## from Chapter 1: heights of women from class
y <- c(65, 65, 69, 67, 66, 66, 67, 64, 68, 68, 68, 69, 70, 65, 70, 
       67, 65, 62, 69, 68, 65, 68, 64, 64)

## summarizing
ny <- length(y)
ybar <- mean(y)
s2y <- var(y)

## heights of Blacksburgh high women's BBall team
x <- c(5*12 + 6, 5*12 + 8, 5*12 + 6, 5*12 + 8, 6*12 + 2, 5*12 + 7, 
       5*12 + 10, 6*12 + 0, 5*12 + 7, 5*12 + 7)

## summarizing
nx <- length(x)
xbar <- mean(x)
s2x <- var(x)

## test statistic
muhat.d <- ybar - xbar

##
## by simulation
##

## two copies of Chapter 1: it doesn't matter what mu is
mu <- 7                                           ## RIP Dad and Mickey
Muhat.ds <- rep(NA, N)                   
for(i in 1:N) {  
  sigma2ys <- (ny - 1)*s2y/rchisq(1, ny - 1)      ## Y population (class)
  Ys <- rnorm(ny, mu, sqrt(sigma2ys))
  sigma2xs <- (nx - 1)*s2x/rchisq(1, nx - 1)      ## X population (bball)
  Xs <- rnorm(nx, mu, sqrt(sigma2xs)) 
  Muhat.ds[i] <- mean(Ys) - mean(Xs)              ## diff in averages           
}

## visualizing
hist(Muhat.ds, main="")
abline(v=c(muhat.d, -muhat.d), lty=1:2, col=2, lwd=2)

## p-value
pval.mc <- 2*mean(Muhat.ds < muhat.d)
pval.mc
legend("topright", c("muhat.d", "reflect"), col=2, lty=1:2, lwd=2, bty="n")

##
## by math: Welch test
##

se2 = s2y/ny + s2x/nx                     ## squared standard error
se <- sqrt(se2)                           ## standard error         
nu <- min(ny, nx) - 1                     ## see book for better calculation
t <- (ybar - xbar)/se                     ## t-stat
2*pt(-abs(t), nu)                         ## p-value

##
## paired t-test
##

## new Xs: heights of moms of women in my class
x <- c(64, 63, 68, 66, 67, 66, 66, 63, 66, 67, 69, 68, 70, 64, 68, 
       66, 65, 61, 69, 66, 61, 67, 63, 64)
nx <- length(x)

## view as paired rows
data.frame(y, x)

## create differences "data"
d <- y - x
n <- length(d)

## summary/estimates from differences
dbar <- mean(d)
s2d <- var(d)

## 
## skipping simulation version
##

## t-test on differences
se.d <- sqrt(s2d/n)
t.d <- dbar/se.d
2*pt(-abs(t.d), n - 1)