# Edps 589
# Fall 2018
# c.j.anderson
#
# Continuation Ratio example  
#
# Employment x father education x race
#
# Can do this using poisson model or binary logit model,
#  which is easier probabably depends on how data are 
#  formated.   I did first ratio both ways, but 2nd one
#  using family=binomial....practice reshaping data.
#
library(vcd)        # for inputting data & xtabs
library(vcdExtra)

var.levels <- expand.grid(
  employ=c("in_school","working","inactive"),
  race=c("White","Black"),
  father=c("<=12",">12")  )
 
counts <- c(204, 195, 131, 100, 53, 67,
             78,  90,  28,  18,  5,  9)
			 
(cr <- cbind(var.levels,counts))

summary(cr$counts)
###########################
## in_school vs working   #
###########################
# Subset
sub1 <- cr[which(cr$employ=="in_school" | cr$employ=="working"),]

############################################
#  In school vs Working as loglinear models#
############################################
# (R,F)
summary(glm(counts ~ race + father + race*father + employ + employ*race + employ*father,
        data=sub1, family=poisson))
# (R)	
summary(glm(counts ~ race + father + race*father + employ + employ*race ,
        data=sub1, family=poisson))
# (F)		
summary(glm(counts ~ race + father + race*father + employ + employ*father,
        data=sub1, family=poisson))
# null			
summary(glm(counts ~ race + father + race*father + employ ,
        data=sub1, family=poisson))

############################################
#  In school vs Working as logit models    #
############################################
# change to wide format
rsub1 <- reshape(sub1, idvar = c("race","father"), 
                       timevar = "employ", 
					   direction = "wide")

# give names that are shorter
names(rsub1) <- c("race","father","in_school","working")

# total is ncases per race x father combination
rsub1$total <- rsub1$in_school + rsub1$working

# Check format
rsub1

# Model fitting
# R + F
summary(r1.rf <- glm(working/total ~ race + father,weights=total,data=rsub1,family=binomial))

# R 
summary(r1.r <-  glm(working/total ~ race ,weights=total,data=rsub1,family=binomial))
anova(r1.r,r1.rf,test="LRT")

# F 
summary(r1.f <- glm(working/total ~ father,weights=total,data=rsub1,family=binomial))
anova(r1.f,r1.rf,test="LRT")

# null
summary(r1.null <- glm(working/total ~ 1,weights=total,data=rsub1,family=binomial))

#####################################################
#  In school or Working vs In active as logit models#
#####################################################

cr$employ2 <- ifelse(cr$employ=="inactive",0,1)

cr.tab <- xtabs(counts ~ race + father + employ2,data=cr)
cr.df <- as.data.frame(cr.tab)
r.cr <- reshape(cr.df, 
                idvar = c("race","father"), 
				timevar = "employ2", 
				direction = "wide")
				
r.cr$total <- r.cr$Freq.0 + r.cr$Freq.1

summary(r2.rf <- glm(Freq.1/total ~ race + father,weights=total,data=r.cr,family=binomial))

summary(r2.r <- glm(Freq.1/total ~ race ,weights=total,data=r.cr,family=binomial))
anova(r2.r,r2.rf,test="LRT")

summary(r2.f <- glm(Freq.1/total ~ father ,weights=total,data=r.cr,family=binomial))
anova(r2.f,r2.rf,test="LRT")