# Example 1: Table 3.2

religion.counts<-c(178,570,138,138,648,252,108,442,252)
table.3.2<-cbind(expand.grid(list(Highest.Degree=c("<HS","HS or JH", "Bachelor or Grad"),
Religious.Beliefs=c("Fund", "Mod", "Lib"))),count=religion.counts) 

(table.3.2.array<-tapply(table.3.2$count,table.3.2[,1:2], sum) )

(row.sum  <- apply(table.3.2.array,1,sum))
(col.sum  <- apply(table.3.2.array,2,sum))
(table.3.2.exp <- outer(row.sum,col.sum) / sum(table.3.2.array))

sum( (table.3.2.exp - table.3.2.array)^2  / table.3.2.exp )  #  X^2 stat

( LR.null <- sum(table.3.2.array*log(table.3.2.exp))   )
( LR.alt  <- sum(table.3.2.array*log(table.3.2.array)) )

 -2 * (LR.null - LR.alt)  #  G stat

# R also knows how to test independence ...

(X.sq<-chisq.test(table.3.2.array)) 
X.sq$exp
X.sq$res
sum(X.sq$res^2)

glm(count~Highest.Degree + Religious.Beliefs, family="poisson", data=table.3.2)$deviance 


#  Example 2: Table 3.3 and 3.4

schizo.counts<-c(90,12,78,13,1,6,19,13,50) 
table.3.3<- cbind(expand.grid(list(Origin=c("Biogenic", "Environmental","Combination"), 
School=c("Eclectic","Medical","Psychoanalytic"))), count=schizo.counts) 

chisq.test(tapply(table.3.3$count,table.3.3[,1:2], sum))

(fit.3.3 <-glm(count~Origin+School, family="poisson", data=table.3.3))
 fit.3.3$deviance

#	Partitioning the Deviance  --  A

table.3.3.1  <- cbind(expand.grid(list(Origin=c("Bio", "Env"),School=c("Ecl","Med"))), count=c(90,12,13,1))
table.3.3.2  <- cbind(expand.grid(list(Origin=c("Bio+Env","Combination"), School=c("Eclectic","Medical"))), count=c(102,78,14,6))
table.3.3.3  <- cbind(expand.grid(list(Origin=c("Bio","Env"), School=c("Ecl+Med","Psy"))), count=c(103,13,19,13))
table.3.3.4  <- cbind(expand.grid(list(Origin=c("Bio+Env","Com"), School=c("Ecl+Med","Psy"))), count=c(116,84,32,50))
 

tapply(table.3.3$count,table.3.3[,1:2], sum) 

tapply(table.3.3.1$count,table.3.3.1[,1:2], sum) 
tapply(table.3.3.2$count,table.3.3.2[,1:2], sum) 
tapply(table.3.3.3$count,table.3.3.3[,1:2], sum) 
tapply(table.3.3.4$count,table.3.3.4[,1:2], sum) 

(dev.1 <- glm(count~Origin+School, family="poisson", data=table.3.3.1)$deviance )
(dev.2 <- glm(count~Origin+School, family="poisson", data=table.3.3.2)$deviance )
(dev.3 <- glm(count~Origin+School, family="poisson", data=table.3.3.3)$deviance )
(dev.4 <- glm(count~Origin+School, family="poisson", data=table.3.3.4)$deviance ) 

dev.1 + dev.2 + dev.3 + dev.4
fit.3.3$deviance


#	Partitioning the Deviance  --  B

table.3.3.5  <- cbind(expand.grid(list(Origin=c("Bio","Env","Com"),School=c("Ecl","Med"))),count=c(90,12,78,13,1,6))
table.3.3.6  <- cbind(expand.grid(list(Origin=c("Bio","Env","Com"),School=c("Ecl+Med","Psy"))),count=c(103,13,84,19,13,50))

tapply(table.3.3.5$count,table.3.3.5[,1:2], sum) 
tapply(table.3.3.6$count,table.3.3.6[,1:2], sum) 


(dev.5 <- glm(count~Origin+School, family="poisson", data=table.3.3.5)$deviance )
(dev.6 <- glm(count~Origin+School, family="poisson", data=table.3.3.6)$deviance )

dev.5 + dev.6
fit.3.3$deviance

#	Res table 

tapply(chisq.test(tapply(table.3.3.6$count,table.3.3.6[,1:2], sum) )$res,table.3.3.6[,1:2], sum) 



#  Example 3: Why consider independence


(p.mat <- cbind(c(.05,.25),c(.15,.55)))	# almost independence delta = 1/6

pp.hat <- rep(NA,10000)
pi.hat <- rep(NA,10000)

n <- 500
for(i in 1:10000)
{
	count.mat <- array(table(sample(letters[1:4],replace=T,size=n,prob=p.mat)),dim=c(2,2))
	pp.hat[i] <- count.mat[1,1] / n
	pi.hat[i] <- sum(count.mat[1,])*sum(count.mat[,1]) / n^2
}

boxplot(list(pp.hat,pi.hat),names=c("Cell count","Marginal mult"),main="p.11 estimation")
t.test(pp.hat)

boxplot(list(log(pp.hat),log(pi.hat)),names=c("Cell count","Marginal mult"),main="Log(p.11) estimation")



log(p.mat[1,1])
(1-p.mat[1,1]) / (p.mat[1,1] * n)   			#  var of cell count estimator
var(log(pp.hat))

( pi.11 <- sum(p.mat[1,])*sum(p.mat[,1]))   		#  marginal mult estimator expectation

(var.pi <- (1-sum(p.mat[,1]))/(sum(p.mat[,1])*n) + (1-sum(p.mat[1,]))/(sum(p.mat[1,])*n) - 2*(1/6)/n)
var(log(pi.hat))



(mean(log(pi.hat)) - log(p.mat[1,1]))^2			# Bias
(-1/6)^2								# Bias approx.

mean(log(pi.hat/p.mat[1,1])^2)   				# MSE est. 2
mean(log(pp.hat/p.mat[1,1])^2)   				# MSE est. 1
var(log(pi.hat)) +(mean(log(pi.hat)) - log(p.mat[1,1]))^2		# MSE  = Var + Bias





delta.mat <- function(mat,delta=0)
{
	d <- delta * mat[1,1]
	mat[1,1] <- mat[1,1] + d
	mat[2,2] <- mat[2,2] + d
	mat[1,2] <- mat[1,2] - d
	mat[2,1] <- mat[2,1] - d

	return(mat)
}


#	simulation 1:

mse.ratio <- array(dim=c(5,5))
sig.05    <- array(0,dim=c(5,5))

pp.hat <- rep(NA,1000)
pi.hat <- rep(NA,1000)

(mat <- outer(c(.2,.8),c(.3,.7)))
delta.mat(mat,delta=-1/6)

delta.vec <- c(0,0.05,0.1,0.2,0.3)
n.vec     <- c(100,200,400,700,1100)

for(i1 in 1:5)
for(i2 in 1:5)
{
sim.mat <- delta.mat(mat,delta=delta.vec[i1])
for(i3 in 1:1000)
{
	count.mat <- array(table(sample(letters[1:4],replace=T,
	size=n.vec[i2],prob=sim.mat)),dim=c(2,2))
	sig.05[i1,i2] <- sig.05[i1,i2] + (chisq.test(count.mat)$p.value < 0.05)/1000
	pp.hat[i3] <- count.mat[1,1] / n.vec[i2]
	pi.hat[i3] <- sum(count.mat[1,])*sum(count.mat[,1]) / n.vec[i2]^2
}
	mse.ratio[i1,i2] <- mean((pi.hat - sim.mat[1,1])^2)/mean((pp.hat - sim.mat[1,1])^2)	
}
par(mfcol=c(1,2))
contour(delta.vec,n.vec,mse.ratio,xlab="Delta",ylab="Sample size - n",main="MSE Ratio: Marg/Cell")
contour(delta.vec,n.vec,sig.05,xlab="Delta",ylab="Sample size - n",main="Power:   Pr(p.value<0.05)")


#	simulation 2:

mse.ratio.2 <- array(dim=c(5,5))
sig.05.2    <- array(0,dim=c(5,5))

pp.hat <- rep(NA,1000)
pi.hat <- rep(NA,1000)

(mat <- outer(c(.2,.8),c(.3,.7)))
delta.mat(mat,delta=-.01)

delta.vec <- c(0,0.05,0.1,0.2,0.3)
k.vec     <- 2:6
n <- 800

for(i1 in 1:5)
for(i2 in 1:5)
{
mat <- outer(rep(1,k.vec[i2])/k.vec[i2],rep(1,k.vec[i2])/k.vec[i2])
sim.mat <- delta.mat(mat,delta=delta.vec[i1])
for(i3 in 1:1000)
{
	count.mat <- array(table(sample(paste(1:k.vec[i2]^2),replace=T,
	size=n.vec[i2],prob=sim.mat)),dim=c(k.vec[i2],k.vec[i2]))
	sig.05.2[i1,i2] <- sig.05.2[i1,i2] + (chisq.test(count.mat)$p.value < 0.05)/1000
	pp.hat[i3] <- count.mat[1,1] / n.vec[i2]
	pi.hat[i3] <- sum(count.mat[1,])*sum(count.mat[,1]) / n.vec[i2]^2
}
	mse.ratio.2[i1,i2] <- mean((pi.hat - sim.mat[1,1])^2)/mean((pp.hat - sim.mat[1,1])^2)	
}

par(mfcol=c(1,2))
contour(delta.vec,k.vec,mse.ratio.2,xlab="Delta",ylab="Number of rows/cols  - k",main="MSE Ratio: Marg/Cell")
contour(delta.vec,k.vec,sig.05.2,xlab="Delta",ylab="Number of rows/cols  - k",main="Power: P(p.value<0.05)")