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

##
## a: Goodness-of-Fit (GoF) test
##

## First of two examples
## In this context, a GoF test is basically a multinomial test

## dice rolls; is the die fair?
o <- c(119, 90, 86, 116, 103, 86)

## summary information
c <- length(o)
n <- sum(o)

## H0: the die is fair
p <- rep(1/c, c)
e <- n*p

## test stat
x2 <- sum((o - e)^2/e)

## when simulating, each die role can be done as
sample(1:c, 1, prob=p)

## n die rolls as
table(sample(1:c, n, replace=TRUE))

## or equivalently
drop(rmultinom(1, size=n, prob=rep(1/c, c)))

## I'll use the first one in this example and the
## second one in the next one

## simulation loop, slower than usual because n is big
N <- 100000
X2s <- rep(NA, N)
for(i in 1:N) {
  Os <- table(sample(1:c, size=n, replace=TRUE))
  X2s[i] <- sum((Os - e)^2/e)
}

## visual
hist(X2s, main="") 
abline(v=x2, col=2, lwd=2)
legend("topright", "obs", lwd=2, col=2, bty="n")

## p-value
## Chi-squared tests (like GoF) are always one-sided
mean(X2s > x2)

## math way (via CLT-like approximation)
pchisq(x2, c - 1, lower.tail=FALSE)

##
## Second of two examples: Is it Gaussian?
##

## scores on a national math test

nels <- read.csv("nels_math.csv")   ## download from book web page
y <- nels[nels[,1] == 1, 2]         ## first school only
n <- length(y)

## H0 it's reasonable to model these data with a Gaussian
## What parameters should we use? The best ones (MLE) of course.
ybar <- mean(y)
s2 <- var(y)

## Now the question is more specific.
## H0: this Gaussian is a good model for these data

## turn into c classes and probabilities
c <- 10
bins <- seq(min(y), max(y), length=c + 1)  ## c + 1 endpoints -> c bins
bins[1] <- -Inf
bins[c + 1] <- Inf
o <- p <- rep(NA, c)
for(i in 1:c) {
  p[i] <- pnorm(bins[i + 1], ybar, sqrt(s2)) 
  p[i] <- p[i] - pnorm(bins[i], ybar, sqrt(s2))
  o[i] <- sum(y < bins[i + 1] & y >= bins[i])
}

## Now we've basically got a multinomial test
## H0: do these probabilities jive with those observations
rbind(o, p)

## better visual provided in the book
mids <- bins[-11] + diff(bins)/2
matplot(mids, cbind(o/sum(o), p), pch=20:21, ylab="probs")
legend("bottom", c("obs-p", "H0-p"), pch=20:21, col=1:2)

## test stat
e <- p*n
x2 <- sum((o - e)^2/e)

## simulation, faster because n is small
for(i in 1:N) {
  Os <- rmultinom(1, size=n, prob=p)
  X2s[i] <- sum((Os - e)^2/e)
}

## p-valuec
mean(X2s > x2)

## math way
pchisq(x2, c - 1, lower.tail=FALSE)

##
## B: Homogeneity and Independence
##

## Two examples

##
## Example 1: 
##

## company job satisfaction survey
o <- rbind(
  c(2, 6, 20, 12, 6),
  c(4, 6, 20, 14, 7),
  c(2,  12,  28, 30, 11),
  c(1, 2, 18, 22, 12),
  c(2, 3, 4, 6, 8))
colnames(o) <- c("miser", "discon", "content", "happy", "smitten")
rownames(o) <- c("<10y", "10-20y", "20-30y", "30-40y", ">40y")

## this makes sense for a Homogeneity test because tenure
## stratification (rows) isn't random, but reported satsifaction 
## (cols) is

## row and col totals
n <- rowSums(o)
c <- colSums(o)
ntot <- sum(n)

## a pretty looking table isn't strictly necessary, but doesn't hurt
tab <- rbind(o, c)
tab <- cbind(tab, c(n, ntot))
colnames(tab) <- c(colnames(o), "n") 
rownames(tab) <- c(rownames(o), "c")
tab

## following the formula to calculate expectaitons
e <- matrix(NA, nrow(o), ncol(o))
for(i in 1:nrow(o)) {
  for(j in 1:ncol(o)) {
    e[i,j] <- n[i]*c[j]/ntot
  }
}

## alternatively via outer product
efast <- outer(n, c/ntot)
  
## test statistic is vectorized
x2 <- sum((o - e)^2/e)
x2

## now thinking about simulation ...

## if every row has the same probabilities, then we can
## sample with categories 1 ... c with prob p=c/ntot
## as follows (check this is the MLE)
Cs <- sample(1:length(c), ntot, prob=c/ntot, replace=TRUE)

## we must arrange those into rows, so it helps to have
## a vector of similar length
ns <- rep(1:length(n), n)

## then we can build a table as follows
Os <- table(ns, factor(Cs, levels=1:length(c)))

## the "factor" step is a precaution in case not
## all categories are represented

## now put that in a loop
for(i in 1:N) {
  Cs <- sample(1:length(c), ntot, prob=c/ntot, replace=TRUE)
  Os <- table(ns, factor(Cs, levels=1:length(c)))
  Es <- outer(rowSums(Os), colSums(Os)/ntot)
  X2s[i] <- sum((Os - Es)^2/Es)
}

## visualize
hist(X2s)
abline(v=x2, col=2, lwd=2)

## p-value (fail to reject)
mean(X2s > x2)

## compare to the math way
pchisq(x2, (nrow(o) - 1)*(ncol(o) - 1), lower.tail=FALSE)

##
## Example 2: 
##

## data on crimes assault crimes
assault <- read.csv("assault.csv")  ## download from book web page
  
## look at just two of the variables
o <- table(assault$locale, assault$income)
o <- o[,c(1,3,4,2)]   ## change column order for increasing income
o

## this is good for an Independence test because both location
## and income are random for any particular crime

## form full table
r <- rowSums(o)
c <- colSums(o)
ntot <- sum(r)
tab <- rbind(o, c)
tab <- cbind(tab, c(r, ntot))
colnames(tab) <- c(colnames(o), "r") 
rownames(tab) <- c(rownames(o), "c")
tab

## calculate test statistic
e <- outer(r, c/ntot)
x2 <- sum((o - e)^2/e)
x2

## to simulate, we must separatly sample from rows and cols
Rs <- sample(1:length(r), ntot, prob=r/ntot, replace=TRUE)
Cs <- sample(1:length(c), ntot, prob=c/ntot, replace=TRUE)
table(factor(Rs, levels=1:length(r)), factor(Cs, levels=1:length(c)))

## simulation  (slower than usual because ntot is large)
for(i in 1:N) {
  Rs <- sample(1:length(r), ntot, prob=r/ntot, replace=TRUE)
  Cs <- sample(1:length(c), ntot, prob=c/ntot, replace=TRUE)
  Os <- table(factor(Rs, levels=1:length(r)), 
              factor(Cs, levels=1:length(c)))
  Es <- outer(rowSums(Os), colSums(Os)/ntot)
  X2s[i] <- sum((Os - Es)^2/Es)
}

## p-value (easy reject)
mean(X2s > x2)

## compare to the math way
pchisq(x2, (nrow(o) - 1)*(ncol(o) - 1), lower.tail=FALSE)