# 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")