# edps 589 # Fall 2018 # C.J.Anderson # # 4-way table .model building # Second part is from lecture on model building # ############################################################## # Set up # ############################################################## library(vcd) library(vcdExtra) library(MASS) library(rms) var.values <- expand.grid(mstatus=c("divorce","married"), ems=c("yes","no"), pms=c("yes","no"), gender=c("woman","man") ) counts <- c(17, 4, 54, 25, 36, 4, 214, 322, 28, 11, 60, 42, 17, 4, 68, 130) gepm <- cbind(var.values,counts) (gepm <- as.data.frame(gepm) ) ############################################################## # LG --> P (i.e., independence) # ############################################################## (tab1 <- xtabs(count ~ gender + pms, data=gepm)) # easy -- test inpdependence # G2 and X2 gp.1 <- loglm(count ~ gender + pms,data=tab1) # X2 only chisq.test(tab1,correct=FALSE) # odds odds(tab1) # Using glm gp.I <- aggregate(counts~ gender + pms,data=gepm,FUN=sum) summary( lm.gp <- glm(counts ~ gender + pms, data=gp.I, family=poisson)) (X2 <- sum((gp.I$counts-lm.gp$fitted)**2/lm.gp$fitted)) ############################################################### # Models GP - > E # ############################################################### (gp.e <- aggregate(counts ~ gender + pms + ems,data=gepm,FUN=sum)) # -- as loglinear models summary(gp.e1 <- glm(counts~gender + pms + gender*pms + ems, data=gp.e,family=poisson)) (X2 <- sum((gp.e$counts-gp.e1$fitted)**2/gp.e1$fitted)) summary(gp.e2<- glm(counts~gender + pms + gender*pms + ems + pms*ems, data=gp.e,family=poisson)) (X2 <- sum((gp.e$counts-gp.e2$fitted)**2/gp.e2$fitted)) summary(gp.e3<- glm(counts~gender + pms + gender*pms + ems + pms*ems + gender*ems, data=gp.e,family=poisson)) (X2 <- sum((gp.e$counts-gp.e3$fitted)**2/gp.e3$fitted)) ################################################################ # Models GPE -> M # ################################################################ summary(gepm1 <- glm(counts ~ gender + pms + ems + gender*pms + gender*ems + pms*ems + gender*ems*pms + mstatus*pms + mstatus*ems, data=gepm,family=poisson)) (X2 <- sum((gepm$counts-gepm1$fitted)**2/gepm1$fitted)) summary(gepm2 <- glm(counts ~ gender+pms+ems+gender*pms+gender*ems+pms*ems+gender*ems*pms + mstatus*pms + mstatus*ems + mstatus*pms*ems, data=gepm,family=poisson)) 1-pchisq(gepm2$deviance,gepm2$df.residual) summary(gepm3 <- glm(counts ~ gender+pms+ems+gender*pms+gender*ems+pms*ems+gender*ems*pms + mstatus*pms + mstatus*ems + mstatus*pms*ems + mstatus*gender, data=gepm,family=poisson)) 1-pchisq(gepm3$deviance,gepm3$df.residual)