# Edps 589
# C.J.Anderson
#
# Modeling building for loglinear and logit models:  
#   linear by linear model
#   RC(1) association model
#   Correspondence analysis
#

library(vcd)
library(vcdExtra)
library(logmult)   # easiest way to fit RC(M) model
library(gnm)       # generalized non-linear models (logmult is a wrapper)
library(ca)        # Correspondence analysis


#setwd("C:\\Users\\cja\\Dropbox\\edps 589\\7Model_building")
setwd("D:\\Dropbox\\edps 589\\7Model_building")

(hsb <- read.table("hsb_linear_by_linear_data.txt",header=TRUE))

####################################################################
# Fit various models (as poisson)                                  #
####################################################################

## Independence
summary( independence <- glm(counts ~ ses + hsp, data=hsb, family=poisson) )

## Linear x Linear  
summary( lin.by.lin <- glm(counts ~ ses + hsp  + u*v, data=hsb, family=poisson) )
(a <- anova(independence,lin.by.lin) )
dim(a)
1-pchisq(a[2,4],a[2,3])
exp(lin.by.lin$coefficients[8])

## Uniform
summary( uniform <- glm(counts ~ ses + hsp + row*col , data=hsb, family=poisson) )
(a <- anova(independence,uniform) )
dim(a)
1-pchisq(a[2,4],a[2,3])
exp(uniform$coefficients[8])

hsb.tab <- xtabs(counts ~ ses + hsp, data=hsb)

###
# RC(1) association model via rc
###
summary( rc1 <- rc(hsb.tab, nd = 1, 
   weighting=c("none"),
   rowsup = NULL, colsup = NULL,
   se = c("jackknife"),
   nreplicates = 100, 
   family = poisson, 
   )
)   

plot(rc1, conf.int=0.95,main="RC(1) Association Model + 95% CI on Scale Values")

###
# rc(1) via gnm
###
rc1.gnm <-gnm(counts ~ ses + hsp + Mult(ses,hsp), data=hsb, family=poisson, verbose=TRUE)


################################################################
# Correspondence analysis                                      #
################################################################
par(mfrow=c(1,1))
hsb.tab <- xtabs(counts ~ ses + hsp, data=hsb) # create a 2 way table
prop.table(hsb.tab, 1) # row percentages
prop.table(hsb.tab, 2) # column percentages
fit <- ca(hsb.tab)
print(fit)   # basic results
summary(fit) # extended results
plot(fit,main="Correspondence Analysis")    # symmetric map
plot(fit, mass = TRUE, contrib = "absolute", map = "rowgreen", arrows = c(FALSE, TRUE)) # asymmetric map 

## Really should only "look" at first dimension
chisq.test(hsb.tab)

names(fit)
r.sv <- fit$rowcoord[,1]
r.level <- c(1,1,1)
c.sv <- fit$colcoord[,1]
c.level <- c(1,.95,1.05)

par(mfrow=c(2,1))
plot(r.sv,r.level,
     ylim=c(.5,1.5),
     xlim=c(-1.5,1.5),
     type="n",
     main="Correspondence Analysis (96.84% of X^2)",
     yaxt="n",
     xlab="SES",
	 ylab=" ",
	 frame=FALSE)
text(r.sv,r.level,c("High","Low","Middle"),col="blue")

plot(c.sv,r.level,
     ylim=c(.5,1.5),
     xlim=c(-1.5,1.5),
     type="n",
     yaxt="n",
     xlab="High School Program Types",
     yaxt="n",
     ylab=" ",
     frame=FALSE)
text(c.sv,c.level,c("Academic","General","VoTech"),col="red")