---
title: '4-way Table: Logit and Loglinear'
author: "Carolyn J. Anderson"
date: '2022-10-31'
output: 
   html_document:
     toc: yes
---

# Set Up

Here are libraries that we'll use.

```{r libs}
library(vcd)
library(vcdExtra)
library(MASS)
```
## Data for use in loglinear models

```{r loglindata}
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)	
df4 <- cbind(var.values,counts)
(df4 <- as.data.frame(df4) )
```

## Data for use in logit models
Split data into married and divorced and created variables needed in logit fit via GLM.  This set up with make comparing loglinear and logit model results much easier that previous code.

```{r mdata}
m <- df4[which(df4$mstatus=="married"),]
m$freq.married <- m$counts

d <- df4[which(df4$mstatus=="divorce"),]
d$freq.divorced <- d$counts

wide.data <- cbind(m[,c(2,3,4,6)], d[,6])
names(wide.data) <- c("ems", "pms", "gender", "freq.married", "freq.divorced")

wide.data$ncases <- wide.data$freq.married + wide.data$freq.divorced
wide.data$p <- wide.data$freq.married/wide.data$ncases
```

# Loglinear (GP,MEP)

## Fit the model
```{r}
summary(mod1 <- glm (counts ~ mstatus + pms + ems + gender + gender*pms
                + mstatus*ems + mstatus*pms + mstatus*pms*ems,
                data=df4, family=poisson) )
			 
1-pchisq(mod1$deviance,mod1$df.residual)
```

## 95% ci for regression coefficients

We can do this in a number of ways.

```{r}
exp(coef(mod1))
exp(confint(mod1)) 

df4$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 
```

The easy way to 95% ci for odds ratios
```{r}
exp(coef(mod1))
exp(confint(mod1)) 
```

# Equivalent logit and loglinear models
#
## Logit (E,P,G) and Loglinear (ME, MP, MG, EPG)

```{r model1}
summary(logit1 <- glm(p ~ ems + pms + gender, data=wide.data, weights=ncases, family=binomial))

summary(loglin1 <- glm(counts ~  mstatus + mstatus*ems + mstatus*pms + ems + pms + gender + ems*pms + ems*gender + pms*gender + ems*pms*gender + mstatus*gender,  data=df4, family=poisson) )
```

Notice that the residual deviance for both models are equivalent (i.e., $G^2=14.629$, $df=4$).  Also that parameters match up; that is,


logit              | Loglinear             | ASE
-------------------|-----------------------|--------
int =-2.1046       | M (married) = -2.1046 | 0.2746
Ems (no)= 1.5655   | ME = 1.5655           | 0.2466
Pms (no)= 0.8894   | MP = 0.8894           | 0.1666             
Gender (man)=0.3057| MG =0.3057            | 0.1445
-----------------------------------------------------


## Logit (E,GP) and Loglinear (ME, MGP,EGP)

```{r}
summary( logit2 <- glm(p ~ ems + pms + gender + gender*pms,  weights=ncases, data=wide.data, family=binomial) )

summary(loglin2 <-  glm (counts ~ mstatus  + mstatus*ems + mstatus*pms + mstatus*gender	+ mstatus*gender*pms + ems + pms + gender + gender*pms + gender*ems + pms*ems + gender*pms*ems, data=df4, family=poisson) ) 
```

## Logit( EG,P) and loglinear (MEG, MP)
```{r}
summary( logit3 <- glm(p ~ ems + pms + gender + gender*ems, weights=ncases, data=wide.data, family=binomial) )

summary(loglin3 <-glm (counts ~ mstatus + mstatus*ems + mstatus*pms + mstatus*gender+ mstatus*gender*ems 
                       + ems + pms + gender + gender*pms +gender*ems + pms*ems + gender*pms*ems,
                         data=df4, family=poisson) ) 
```

## Logit (EG, P) and loglinear (MEG, MP)
```{r}
summary( logit4<- glm(p ~ ems + pms + gender + pms*ems, weights=ncases, data=wide.data, family=binomial) )

summary(loglin4 <-  glm (counts ~ mstatus + mstatus*ems + mstatus*pms + mstatus*gender+mstatus*pms*ems 
                         + ems + pms + gender + gender*pms +gender*ems + pms*ems + gender*pms*ems ,
                         data=df4, family=poisson) ) 
```

## Logit (EG, PG) and Loglinear (MEG, MPG)

```{r}
summary(logit5<- glm(p ~ ems + pms + gender + gender*ems + gender*pms,
                  weights=ncases, data=wide.data, family=binomial) )

summary(loglin5 <-  glm (counts ~ mstatus  + mstatus*ems + mstatus*pms + mstatus*gender + mstatus*gender*ems + mstatus*gender*pms + ems + pms + gender + gender*pms +gender*ems + pms*ems + gender*pms*ems,
                         data=df4, family=poisson) ) 
```

## Logit (EP, PG) and LogLinear (MEP, MPG)
```{r}
summary(logit6<- glm(p ~ ems + pms + gender + gender*pms + ems*pms,
                  weights=ncases, data=wide.data, family=binomial) )

summary(loglin6 <-glm (counts ~ mstatus + mstatus*ems + mstatus*pms + mstatus*gender+ mstatus*ems*gender + mstatus*gender*pms + ems + pms + gender + gender*pms +gender*ems + pms*ems + gender*pms*ems,
                         data=df4, family=poisson) ) 

```
## Logit (EP, EG) and Loglinear (MEP, MEG)
```{r}
summary( logit7 <- glm(p ~ ems + pms + gender + gender*ems + ems*pms,
                  weights= ncases, data=wide.data, family=binomial) )

summary(loglin7 <-  glm (counts ~ 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=df4, family=poisson) ) 
```

## Logit (EP,EG,PG) and loglinear model (MEP, MEG, MPG)

```{r}
summary( logit8 <- glm(p ~ ems + pms + gender + gender*ems + ems*pms + pms*gender,
                  weights=ncases, data=wide.data ,family=binomial) )

summary(loglin8 <- glm (counts ~ mstatus + mstatus*ems + mstatus*pms + mstatus*gender
						+ mstatus*ems*pms + mstatus*ems*gender + mstatus*pms*gender
						+ ems + pms + gender + gender*pms +gender*ems  + pms*ems + gender*pms*ems ,
                         data=df4, family=poisson) ) 
						 