# edps 589
# Fall 2018
# C.J.Anderson
#
# 4-way table & log-linear models
#
##############################################################
# 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") )
count <- 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) )

##############################################################
# Loglinear Model (GP,MEP)                                   #
##############################################################	 

############################################################


logit1b<- lrm(mstaus~ems+pms+gender,data=gepm,family=binomial)

summary(mod1 <- glm (count ~ mstatus + ems + pms + gender + gender*pms 
                + mstatus*ems + mstatus*pms + mstatus*ems*pms,
             data=gepm, family=poisson) )
			 
1-pchisq(mod1$deviance,mod1$df.residual)

# 95% ci for regression coefficients
exp(coef(mod1))
exp(confint(mod1)) 

gepm$fit1 <- mod1$fitted

# This give the covariance matrix of estimated coefficients
covb <- vcov(mod1)
sub_covb <- c(covb[7,7],covb[7,10],covb[10,7],covb[10,10])
sub.covb <- matrix(sub_covb,nrow=2,ncol=2)

# --- you have to figure out linear combination of parameter
#     estimates to give the odds ratio that your interested in 

# 95% ci for regression coefficients
exp(coef(mod1))
exp(confint(mod1)) 

###############################################################
# Logit models                                                #
#   This isn't completely general                             #
###############################################################
# deal with duplicity of p's from using glm for logit model for
#   data in data frame with counts or weights for each line of
#   data.

# Run this but change glm output to logit.mod
   # run 1st
	rm(gepm.o)
	
   # run the rest down to "###"
 	
    # - get numbers to use to sum pairs
		every2 <- seq(from=1,to=length(gepm$counts),by=2)		
    # - sum pairs
		n_i <- gepm$counts[every2] + gepm$counts[every2+1]
    # - put in order
		gepm.o <- gepm[order(gepm$mstatus),]
	# - put into data file
		gepm.o$n_i <- rep(n_i,2)
    # - need pi and (1-pi) in correct rows	
	    gepm$p.fit <- logit.mod$fitted 		
		gepm.o$p.fit <- ifelse(gepm.o$mstatus=="divorce",(1-gepm.o$p.fit),gepm.o$p.fit )
    # Fitted counts
        gepm.o$fit.count <- gepm.o$p.fit*n_i
	# Compute test statistics
	
     (X2 <- sum((gepm.o$counts-gepm.o$fit.count)**2/gepm.o$fit.count))
     (p_x2 <- 1-pchisq(X2,logit.mod$df.residual))

     (G2 <- 2*sum(gepm.o$counts*log(gepm.o$counts/gepm.o$fit.count)))
     (p_g2 <- 1-pchisq(G2,logit.mod$df.residual))
###   


### E,P,G 
summary( logit1 <- glm(mstatus ~ ems + pms + gender,weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit1
# Run code above for test statistics

summary(loglinear1 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems
                        + mstatus*ems + mstatus*pms + mstatus*gender,
                         data=gepm, family=poisson) ) 

### E,GP 
summary( logit2 <- glm(mstatus ~ ems + pms + gender + gender*pms,weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit2
# Run code above

summary(loglinear2 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+mstatus*gender*pms,
                         data=gepm, family=poisson) ) 

### EG,P 
summary( logit3 <- glm(mstatus ~ ems + pms + gender + gender*ems,weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit3
# Run code above

						 
summary(loglinear3 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+mstatus*gender*ems,
                         data=gepm, family=poisson) ) 
### EG,P 
summary( logit4<- glm(mstatus ~ ems + pms + gender + pms*ems,weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit4
# Run code above
						 
summary(loglinear4 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+mstatus*pms*ems,
                         data=gepm, family=poisson) ) 
						 
### EG,PG 
summary( logit5<- glm(mstatus ~ ems + pms + gender + gender*ems + gender*pms,
                  weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit5
# Run code above
						 
summary(loglinear5 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+ mstatus*gender*ems + mstatus*gender*pms,
                         data=gepm, family=poisson) ) 
						 
						 						 
### EP,PG 
summary( logit6<- glm(mstatus ~ ems + pms + gender + gender*pms + ems*pms,
                  weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit6
# Run code above

				 
summary(loglinear6 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+ mstatus*gender*pms + mstatus*gender*pms,
                         data=gepm, family=poisson) ) 
						 
### EP,EG 
summary( logit7 <- glm(mstatus ~ ems + pms + gender + gender*ems + ems*pms,
                  weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit7
# Run code above
		 
summary(loglinear7 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+ mstatus*gender*ems + mstatus*ems*pms,
                         data=gepm, family=poisson) ) 
						 
### EP,EG,PG
summary( logit8 <- glm(mstatus ~ ems + pms + gender + gender*ems + ems*pms + pms*gender,
                  weights=counts,data=gepm,family=binomial(link="logit")) )
logit.mod<-logit8
# Run code above
		 
summary(loglinear8 <-  glm (count ~ mstatus + ems + pms + gender + gender*pms +gender*ems 
                        + pms*ems + gender*pms*ems + mstatus*ems + mstatus*pms + mstatus*gender
						+ mstatus*gender*ems + mstatus*ems*pms,
                         data=gepm, family=poisson) )