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

##
## a: tests of variance
##

## 
## one-sample test
##

## ice-cream scoops
y <- c(4.0, 5.5, 4.9, 3.8, 4.3, 4.4, 5.1, 4.1, 4.0, 3.0, 3.7, 3.4, 3.6,
       3.0, 3.5, 4.0)
n <- length(y)

## calculate observed statistic
s2 <- var(y)

## test H0: sigma2 = 0.5^2, vs. H1: sigma2 != 0.5^2
sigma2 <- 0.5^2

## need a mu for simulation, but it doesn't matter what you choose
mu <- 4

## same simulation as Chapter 2, except different stat
N <- 100000
S2s <- rep(NA, N) 
for(i in 1:N) {
  Ys <- rnorm(n, mu, sqrt(sigma2))
  S2s[i] <- var(Ys)                 ## different stat here
}

## visualize, looks like a reject to me
hist(S2s, main="", xlim=range(c(s2, S2s)))
abline(v=c(s2, quantile(S2s, mean(S2s > s2))), lty=1:2, col=2, lwd=2)
legend("topright", c("s2", "reflect"), col=2, lty=1:2, lwd=2, bty="n")

## p-value calculation
2*mean(S2s > s2)

## math version: easier to use transformation
v <- s2*(n - 1)/sigma2

## visual based on chi-squared
x <- seq(0, 3*n, length(1000))
plot(x, dchisq(x, n-1), type="l", lwd=2, main="")
abline(v=v, lty=1, col=2, lwd=2)

## p-value based on chi-squared
2*pchisq(v, n - 1, lower.tail=FALSE) 


##
## confidence intervals
##

## the math way
q <- qchisq(c(0.025, 0.975), n - 1)
s2*(n - 1)/rev(q)   ## note this reversal

##
## now two samples
##

## suppose some were actually sorbet
x <- y[10:16]       ## sorbet
s2x <- var(x)
nx <- length(x)
y <- y[1:9]         ## ice cream
ny <- length(y)
s2y <- var(y)

## test statistic
f <- s2y/s2x

## doesn't matter what sigma2 = sigma2x = sigma2y is
sigma2 <- 1
Fs <- rep(NA, N)
for(i in 1:N) {
  Ys <- rnorm(ny, mu, sqrt(sigma2))
  Xs <- rnorm(nx, mu, sqrt(sigma2))
  s2ys <- var(Ys)
  s2xs <- var(Xs)
  Fs[i] <- s2ys/s2xs
}

## visual
hist(Fs[Fs < 8], main="", xlim=range(c(0,8)))  ## focus histogram near 0
abline(v=c(f, quantile(Fs, mean(Fs > f))), lty=1:2, col=2, lwd=2)

## p-value
2*mean(Fs > f)

## math way

## visual
fgrid <- seq(0, 8, length=1000)
plot(fgrid, df(fgrid, ny - 1, nx - 1), type="l", lwd=2, main="")
abline(v=f, lty=1, col=2, lwd=2)
legend("topright", c("f", "reflect"), lty=1:2, col=2, lwd=2, bty="n")

## p-value
2*pf(f, ny - 1, nx - 1, lower.tail=FALSE)

##
## b: ANOVA
## 

## list representation
w <- c(3.5, 3.4, 3.7, 3.4, 3.9, 4.2, 4.0, 4.0, 3.8, 3.8)
ylist <- list(cream=y, sorbet=x, fatfree=w)
ylist

## benefits of list representation
nj <- sapply(ylist, length)
m <- length(nj)
n <- sum(nj)

##
## summary statistics
##

ybarj <- sapply(ylist, mean)
ybarj

s2j <- sapply(ylist, var)
s2.pool <- sum((nj - 1)*s2j)/(n - m)
s2.pool

ybar <- mean(unlist(ylist))
s2 <- var(unlist(ylist))
c(ybar=ybar, s2=s2)

##
## ANOVA decomposition
##

sse <- sum((nj - 1)*s2j)
ssr <- sum(nj*(ybarj - ybar)^2)
sst <- ssr + sse

## coefficient of determination
r2 <- ssr/sst
r2

## f-stat
msr <- ssr/(m - 1)
mse <- sse/(n - m)
f <-  msr/mse

##
## testing via simulation
##

N <- 1000000
mu <- sigma2 <- 1               ## any mu and sigma2 values work
R2s <- Fs <- rep(NA, N)         ## storage for both stats
for(i in 1:N) {
  
  ## draw under null hypothesis
  ## essentially a "for" loop {
  Ys <- lapply(nj, rnorm, mean=ybar, sd=sqrt(s2))
  
  ## within-group calculations
  Ybarjs <- sapply(Ys, mean)
  S2js <- sapply(Ys, var)
  ## } where the "for" loop would end
  
  ## grand average
  Ybars <- sum(nj*Ybarjs)/sum(nj)
  
  ## sums of squares
  SSEs <- sum((nj - 1)*S2js)
  SSRs <- sum(nj*(Ybarjs - Ybars)^2)
  
  ## sampled R^2 statistic
  SSTs <- SSEs + SSRs
  R2s[i] <- SSRs/SSTs
  
  ## sampled F statistic
  MSRs <- SSRs/(m - 1)
  MSEs <- SSEs/(n - m)
  
  ## f statistic
  Fs[i] <- MSRs/MSEs
}

## visual
par(mfrow=c(1, 2))
hist(R2s, main="")
abline(v=r2, col=2, lwd=2)
hist(Fs, main="")
abline(v=f, col=2, lwd=2)
legend("topright", "observed stat", lty=1, lwd=2, col=2,  bty="n")

## sim p-values
pval.mc <- c(r2.mc=mean(R2s > r2), f.mc=mean(Fs > f))
pval.mc

## math way
pf(f, m-1, n-m, lower.tail=FALSE)