Set Up

Here are libraries that we’ll use.

library(vcd)
## Loading required package: grid
library(vcdExtra)
## Loading required package: gnm
library(MASS)

Data for use in loglinear models

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) )
##    mstatus ems pms gender counts
## 1  divorce yes yes  woman     17
## 2  married yes yes  woman      4
## 3  divorce  no yes  woman     54
## 4  married  no yes  woman     25
## 5  divorce yes  no  woman     36
## 6  married yes  no  woman      4
## 7  divorce  no  no  woman    214
## 8  married  no  no  woman    322
## 9  divorce yes yes    man     28
## 10 married yes yes    man     11
## 11 divorce  no yes    man     60
## 12 married  no yes    man     42
## 13 divorce yes  no    man     17
## 14 married yes  no    man      4
## 15 divorce  no  no    man     68
## 16 married  no  no    man    130

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.

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

summary(mod1 <- glm (counts ~ mstatus + pms + ems + gender + gender*pms
                + mstatus*ems + mstatus*pms + mstatus*pms*ems,
                data=df4, family=poisson) )
## 
## Call:
## glm(formula = counts ~ mstatus + pms + ems + gender + gender * 
##     pms + mstatus * ems + mstatus * pms + mstatus * pms * ems, 
##     family = poisson, data = df4)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.12272  -0.60305   0.00924   0.62703   1.08468  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  2.9270     0.1675  17.470  < 2e-16 ***
## mstatusmarried              -1.0986     0.2981  -3.685 0.000229 ***
## pmsno                        0.7210     0.2178   3.311 0.000929 ***
## emsno                        0.9295     0.1761   5.280 1.29e-07 ***
## genderman                    0.3436     0.1307   2.628 0.008586 ** 
## pmsno:genderman             -1.3106     0.1530  -8.569  < 2e-16 ***
## mstatusmarried:emsno         0.5671     0.3355   1.690 0.091002 .  
## mstatusmarried:pmsno        -0.7922     0.4824  -1.642 0.100564    
## pmsno:emsno                  0.7421     0.2311   3.211 0.001323 ** 
## mstatusmarried:pmsno:emsno   1.7955     0.5121   3.506 0.000454 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.8549  on 15  degrees of freedom
## Residual deviance:    8.1535  on  6  degrees of freedom
## AIC: 112.16
## 
## Number of Fisher Scoring iterations: 4
1-pchisq(mod1$deviance,mod1$df.residual)
## [1] 0.2270702

95% ci for regression coefficients

We can do this in a number of ways.

exp(coef(mod1))
##                (Intercept)             mstatusmarried 
##                 18.6721992                  0.3333333 
##                      pmsno                      emsno 
##                  2.0565333                  2.5333333 
##                  genderman            pmsno:genderman 
##                  1.4100000                  0.2696513 
##       mstatusmarried:emsno       mstatusmarried:pmsno 
##                  1.7631579                  0.4528302 
##                pmsno:emsno mstatusmarried:pmsno:emsno 
##                  2.1002979                  6.0225998
exp(confint(mod1)) 
## Waiting for profiling to be done...
##                                 2.5 %     97.5 %
## (Intercept)                13.2623287 25.6048005
## mstatusmarried              0.1797506  0.5836535
## pmsno                       1.3445827  3.1641322
## emsno                       1.8079908  3.6115977
## genderman                   1.0929295  1.8258637
## pmsno:genderman             0.1993896  0.3633231
## mstatusmarried:emsno        0.9302023  3.4919961
## mstatusmarried:pmsno        0.1685466  1.1412597
## pmsno:emsno                 1.3328079  3.3037901
## mstatusmarried:pmsno:emsno  2.2513211 17.0681176
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

exp(coef(mod1))
##                (Intercept)             mstatusmarried 
##                 18.6721992                  0.3333333 
##                      pmsno                      emsno 
##                  2.0565333                  2.5333333 
##                  genderman            pmsno:genderman 
##                  1.4100000                  0.2696513 
##       mstatusmarried:emsno       mstatusmarried:pmsno 
##                  1.7631579                  0.4528302 
##                pmsno:emsno mstatusmarried:pmsno:emsno 
##                  2.1002979                  6.0225998
exp(confint(mod1)) 
## Waiting for profiling to be done...
##                                 2.5 %     97.5 %
## (Intercept)                13.2623287 25.6048005
## mstatusmarried              0.1797506  0.5836535
## pmsno                       1.3445827  3.1641322
## emsno                       1.8079908  3.6115977
## genderman                   1.0929295  1.8258637
## pmsno:genderman             0.1993896  0.3633231
## mstatusmarried:emsno        0.9302023  3.4919961
## mstatusmarried:pmsno        0.1685466  1.1412597
## pmsno:emsno                 1.3328079  3.3037901
## mstatusmarried:pmsno:emsno  2.2513211 17.0681176

Equivalent logit and loglinear models

Logit (E,P,G) and Loglinear (ME, MP, MG, EPG)

summary(logit1 <- glm(p ~ ems + pms + gender, data=wide.data, weights=ncases, family=binomial))
## 
## Call:
## glm(formula = p ~ ems + pms + gender, family = binomial, data = wide.data, 
##     weights = ncases)
## 
## Deviance Residuals: 
##        2         4         6         8        10        12        14        16  
##  1.10324  -0.96777  -2.12076   0.66225   2.26247  -0.61508  -1.01946  -0.05336  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.1046     0.2746  -7.664 1.80e-14 ***
## emsno         1.5655     0.2466   6.349 2.17e-10 ***
## pmsno         0.8894     0.1666   5.339 9.37e-08 ***
## genderman     0.3057     0.1445   2.116   0.0344 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.956  on 7  degrees of freedom
## Residual deviance:  13.629  on 4  degrees of freedom
## AIC: 56.893
## 
## Number of Fisher Scoring iterations: 4
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) )
## 
## Call:
## glm(formula = counts ~ mstatus + mstatus * ems + mstatus * pms + 
##     ems + pms + gender + ems * pms + ems * gender + pms * gender + 
##     ems * pms * gender + mstatus * gender, family = poisson, 
##     data = df4)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
## -0.40347   1.02681   0.57326  -0.77971   0.90329  -1.91878  -0.50928   0.42333  
##        9        10        11        12        13        14        15        16  
## -0.97188   2.04309   0.40420  -0.46363   0.51324  -0.88084   0.04328  -0.03122  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.92950    0.22025  13.301  < 2e-16 ***
## mstatusmarried           -2.10456    0.27459  -7.664 1.80e-14 ***
## emsno                     0.98045    0.25151   3.898 9.69e-05 ***
## pmsno                     0.49960    0.27207   1.836 0.066313 .  
## genderman                 0.58092    0.27134   2.141 0.032282 *  
## mstatusmarried:emsno      1.56549    0.24658   6.349 2.17e-10 ***
## mstatusmarried:pmsno      0.88936    0.16659   5.339 9.37e-08 ***
## emsno:pmsno               0.99105    0.29996   3.304 0.000954 ***
## emsno:genderman          -0.44916    0.31224  -1.439 0.150288    
## pmsno:genderman          -1.30392    0.38263  -3.408 0.000655 ***
## mstatusmarried:genderman  0.30574    0.14451   2.116 0.034372 *  
## emsno:pmsno:genderman    -0.01418    0.41893  -0.034 0.973006    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.855  on 15  degrees of freedom
## Residual deviance:   13.629  on  4  degrees of freedom
## AIC: 121.63
## 
## Number of Fisher Scoring iterations: 5

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)

summary( logit2 <- glm(p ~ ems + pms + gender + gender*pms,  weights=ncases, data=wide.data, family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + gender * pms, family = binomial, 
##     data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##       2        4        6        8       10       12       14       16  
##  1.2924  -0.4672  -2.1530   0.4614   2.0789  -1.0152  -0.9152   0.2691  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.2277     0.3177  -7.012 2.35e-12 ***
## emsno             1.5699     0.2467   6.365 1.95e-10 ***
## pmsno             1.0258     0.2417   4.243 2.20e-05 ***
## genderman         0.5044     0.2903   1.738   0.0823 .  
## pmsno:genderman  -0.2645     0.3341  -0.792   0.4286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.96  on 7  degrees of freedom
## Residual deviance:  13.00  on 3  degrees of freedom
## AIC: 58.263
## 
## Number of Fisher Scoring iterations: 4
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) ) 
## 
## Call:
## glm(formula = counts ~ mstatus + mstatus * ems + mstatus * pms + 
##     mstatus * gender + mstatus * gender * pms + ems + pms + gender + 
##     gender * pms + gender * ems + pms * ems + gender * pms * 
##     ems, family = poisson, data = df4)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -0.4576   1.2087   0.2696  -0.3816   0.9207  -1.9462  -0.3557   0.2939  
##       9       10       11       12       13       14       15       16  
## -0.9098   1.8693   0.6772  -0.7563   0.4545  -0.7944  -0.2171   0.1591  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     2.94217    0.22040  13.349  < 2e-16 ***
## mstatusmarried                 -2.22772    0.31771  -7.012 2.35e-12 ***
## emsno                           1.00990    0.25336   3.986 6.72e-05 ***
## pmsno                           0.48388    0.27301   1.772  0.07633 .  
## genderman                       0.55718    0.27311   2.040  0.04134 *  
## mstatusmarried:emsno            1.56989    0.24665   6.365 1.96e-10 ***
## mstatusmarried:pmsno            1.02577    0.24174   4.243 2.20e-05 ***
## mstatusmarried:genderman        0.50435    0.29027   1.738  0.08229 .  
## pmsno:genderman                -1.26231    0.38608  -3.270  0.00108 ** 
## emsno:genderman                -0.50363    0.31991  -1.574  0.11542    
## emsno:pmsno                     0.95424    0.30356   3.144  0.00167 ** 
## mstatusmarried:pmsno:genderman -0.26446    0.33406  -0.792  0.42857    
## emsno:pmsno:genderman           0.06429    0.43039   0.149  0.88126    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.9  on 15  degrees of freedom
## Residual deviance:   13.0  on  3  degrees of freedom
## AIC: 123
## 
## Number of Fisher Scoring iterations: 5

Logit( EG,P) and loglinear (MEG, MP)

summary( logit3 <- glm(p ~ ems + pms + gender + gender*ems, weights=ncases, data=wide.data, family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + gender * ems, family = binomial, 
##     data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##       2        4        6        8       10       12       14       16  
##  1.7532  -1.0348  -1.1269   0.3867   1.4109  -0.3428  -1.7345   0.2550  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.5432     0.4044  -6.288 3.21e-10 ***
## emsno             2.0199     0.3913   5.162 2.44e-07 ***
## pmsno             0.8978     0.1666   5.390 7.05e-08 ***
## genderman         1.0856     0.4917   2.208   0.0272 *  
## emsno:genderman  -0.8501     0.5097  -1.668   0.0954 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.956  on 7  degrees of freedom
## Residual deviance:  10.746  on 3  degrees of freedom
## AIC: 56.009
## 
## Number of Fisher Scoring iterations: 4
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) ) 
## 
## Call:
## glm(formula = counts ~ mstatus + mstatus * ems + mstatus * pms + 
##     mstatus * gender + mstatus * gender * ems + ems + pms + gender + 
##     gender * pms + gender * ems + pms * ems + gender * pms * 
##     ems, family = poisson, data = df4)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -0.5722   1.6572   0.6151  -0.8321   0.4214  -1.0452  -0.2984   0.2460  
##       9       10       11       12       13       14       15       16  
## -0.6594   1.2473   0.2229  -0.2604   0.9538  -1.4487  -0.2057   0.1507  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      2.9688     0.2202  13.482  < 2e-16 ***
## mstatusmarried                  -2.5432     0.4045  -6.288 3.22e-10 ***
## emsno                            0.9353     0.2529   3.697 0.000218 ***
## pmsno                            0.5436     0.2720   1.998 0.045682 *  
## genderman                        0.4854     0.2781   1.745 0.080908 .  
## mstatusmarried:emsno             2.0199     0.3913   5.162 2.44e-07 ***
## mstatusmarried:pmsno             0.8978     0.1666   5.390 7.05e-08 ***
## mstatusmarried:genderman         1.0856     0.4917   2.208 0.027252 *  
## emsno:genderman                 -0.3241     0.3217  -1.007 0.313743    
## pmsno:genderman                 -1.4053     0.3885  -3.618 0.000297 ***
## emsno:pmsno                      0.9386     0.3021   3.107 0.001891 ** 
## mstatusmarried:emsno:genderman  -0.8501     0.5097  -1.668 0.095360 .  
## emsno:pmsno:genderman            0.1020     0.4248   0.240 0.810240    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.855  on 15  degrees of freedom
## Residual deviance:   10.746  on  3  degrees of freedom
## AIC: 120.75
## 
## Number of Fisher Scoring iterations: 5

Logit (EG, P) and loglinear (MEG, MP)

summary( logit4<- glm(p ~ ems + pms + gender + pms*ems, weights=ncases, data=wide.data, family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + pms * ems, family = binomial, 
##     data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##       2        4        6        8       10       12       14       16  
## -0.2594  -0.2549  -0.3737   0.2015   0.1725   0.2141   0.4373  -0.3443  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.3049     0.3150  -4.142 3.44e-05 ***
## emsno         0.5962     0.3366   1.771  0.07651 .  
## pmsno        -0.7004     0.4850  -1.444  0.14874    
## genderman     0.3089     0.1458   2.118  0.03415 *  
## emsno:pmsno   1.7999     0.5130   3.509  0.00045 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.95608  on 7  degrees of freedom
## Residual deviance:   0.69784  on 3  degrees of freedom
## AIC: 45.961
## 
## Number of Fisher Scoring iterations: 4
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) ) 
## 
## Call:
## glm(formula = counts ~ mstatus + mstatus * ems + mstatus * pms + 
##     mstatus * gender + mstatus * pms * ems + ems + pms + gender + 
##     gender * pms + gender * ems + pms * ems + gender * pms * 
##     ems, family = poisson, data = df4)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
##  0.11762  -0.23118   0.14538  -0.20932   0.12524  -0.35213  -0.15584   0.12780  
##        9        10        11        12        13        14        15        16  
## -0.09026   0.14699  -0.13621   0.16516  -0.17840   0.39920   0.28061  -0.19949  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  2.8046     0.2283  12.283  < 2e-16 ***
## mstatusmarried              -1.3049     0.3150  -4.142 3.44e-05 ***
## emsno                        1.1646     0.2588   4.500 6.80e-06 ***
## pmsno                        0.7580     0.2809   2.698  0.00697 ** 
## genderman                    0.5447     0.2733   1.993  0.04625 *  
## mstatusmarried:emsno         0.5962     0.3366   1.771  0.07651 .  
## mstatusmarried:pmsno        -0.7004     0.4850  -1.444  0.14874    
## mstatusmarried:genderman     0.3089     0.1458   2.118  0.03415 *  
## emsno:pmsno                  0.6495     0.3130   2.075  0.03798 *  
## pmsno:genderman             -1.2311     0.3828  -3.216  0.00130 ** 
## emsno:genderman             -0.4019     0.3106  -1.294  0.19568    
## mstatusmarried:emsno:pmsno   1.7999     0.5130   3.509  0.00045 ***
## emsno:pmsno:genderman       -0.1030     0.4230  -0.243  0.80763    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.85486  on 15  degrees of freedom
## Residual deviance:    0.69784  on  3  degrees of freedom
## AIC: 110.7
## 
## Number of Fisher Scoring iterations: 4

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

summary(logit5<- glm(p ~ ems + pms + gender + gender*ems + gender*pms,
                  weights=ncases, data=wide.data, family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + gender * ems + gender * 
##     pms, family = binomial, data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##       2        4        6        8       10       12       14       16  
##  1.8791  -0.6166  -1.1765   0.2284   1.2938  -0.6716  -1.6157   0.4993  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.6335     0.4306  -6.116 9.58e-10 ***
## emsno             2.0113     0.3921   5.130 2.90e-07 ***
## pmsno             1.0107     0.2438   4.145 3.39e-05 ***
## genderman         1.2212     0.5358   2.279   0.0227 *  
## emsno:genderman  -0.8210     0.5115  -1.605   0.1085    
## pmsno:genderman  -0.2151     0.3351  -0.642   0.5210    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.956  on 7  degrees of freedom
## Residual deviance:  10.332  on 2  degrees of freedom
## AIC: 57.596
## 
## Number of Fisher Scoring iterations: 4
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) ) 
## 
## Call:
## glm(formula = 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, family = poisson, data = df4)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -0.5994   1.7809   0.3586  -0.5016   0.4430  -1.0900  -0.1765   0.1449  
##       9       10       11       12       13       14       15       16  
## -0.6114   1.1402   0.4423  -0.5054   0.8761  -1.3576  -0.4011   0.2974  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      2.9752     0.2201  13.516  < 2e-16 ***
## mstatusmarried                  -2.6335     0.4306  -6.116 9.59e-10 ***
## emsno                            0.9646     0.2560   3.767 0.000165 ***
## pmsno                            0.5336     0.2729   1.956 0.050518 .  
## genderman                        0.4704     0.2793   1.684 0.092100 .  
## mstatusmarried:emsno             2.0113     0.3921   5.129 2.91e-07 ***
## mstatusmarried:pmsno             1.0107     0.2438   4.145 3.39e-05 ***
## mstatusmarried:genderman         1.2212     0.5358   2.279 0.022657 *  
## emsno:genderman                 -0.3735     0.3308  -1.129 0.258884    
## pmsno:genderman                 -1.3663     0.3926  -3.480 0.000501 ***
## emsno:pmsno                      0.9046     0.3067   2.949 0.003185 ** 
## mstatusmarried:emsno:genderman  -0.8210     0.5115  -1.605 0.108511    
## mstatusmarried:pmsno:genderman  -0.2151     0.3351  -0.642 0.521044    
## emsno:pmsno:genderman            0.1591     0.4333   0.367 0.713504    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.855  on 15  degrees of freedom
## Residual deviance:   10.332  on  2  degrees of freedom
## AIC: 122.34
## 
## Number of Fisher Scoring iterations: 5

Logit (EP, PG) and LogLinear (MEP, MPG)

summary(logit6<- glm(p ~ ems + pms + gender + gender*pms + ems*pms,
                  weights=ncases, data=wide.data, family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + gender * pms + ems * pms, 
##     family = binomial, data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##        2         4         6         8        10        12        14        16  
## -0.10282   0.04534  -0.40883   0.07229   0.06684  -0.03764   0.48337  -0.12296  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.3901     0.3591  -3.871 0.000108 ***
## emsno             0.6090     0.3383   1.800 0.071788 .  
## pmsno            -0.5976     0.5267  -1.134 0.256590    
## genderman         0.4320     0.2837   1.523 0.127778    
## pmsno:genderman  -0.1677     0.3306  -0.507 0.611932    
## emsno:pmsno       1.7809     0.5146   3.460 0.000539 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.95608  on 7  degrees of freedom
## Residual deviance:   0.43965  on 2  degrees of freedom
## AIC: 47.703
## 
## Number of Fisher Scoring iterations: 4
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) ) 
## 
## Call:
## glm(formula = 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, family = poisson, data = df4)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -0.5994   1.7809   0.3586  -0.5016   0.4430  -1.0900  -0.1765   0.1449  
##       9       10       11       12       13       14       15       16  
## -0.6114   1.1402   0.4423  -0.5054   0.8761  -1.3576  -0.4011   0.2974  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      2.9752     0.2201  13.516  < 2e-16 ***
## mstatusmarried                  -2.6335     0.4306  -6.116 9.59e-10 ***
## emsno                            0.9646     0.2560   3.767 0.000165 ***
## pmsno                            0.5336     0.2729   1.956 0.050518 .  
## genderman                        0.4704     0.2793   1.684 0.092100 .  
## mstatusmarried:emsno             2.0113     0.3921   5.129 2.91e-07 ***
## mstatusmarried:pmsno             1.0107     0.2438   4.145 3.39e-05 ***
## mstatusmarried:genderman         1.2212     0.5358   2.279 0.022657 *  
## emsno:genderman                 -0.3735     0.3308  -1.129 0.258884    
## pmsno:genderman                 -1.3663     0.3926  -3.480 0.000501 ***
## emsno:pmsno                      0.9046     0.3067   2.949 0.003185 ** 
## mstatusmarried:emsno:genderman  -0.8210     0.5115  -1.605 0.108511    
## mstatusmarried:pmsno:genderman  -0.2151     0.3351  -0.642 0.521044    
## emsno:pmsno:genderman            0.1591     0.4333   0.367 0.713504    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.855  on 15  degrees of freedom
## Residual deviance:   10.332  on  2  degrees of freedom
## AIC: 122.34
## 
## Number of Fisher Scoring iterations: 5

Logit (EP, EG) and Loglinear (MEP, MEG)

summary( logit7 <- glm(p ~ ems + pms + gender + gender*ems + ems*pms,
                  weights= ncases, data=wide.data, family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + gender * ems + ems * pms, 
##     family = binomial, data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##        2         4         6         8        10        12        14        16  
##  0.13038  -0.32409  -0.11978   0.11917  -0.08189   0.27297   0.13038  -0.20302  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      -1.5199     0.4716  -3.223  0.00127 **
## emsno             0.8279     0.5041   1.642  0.10053   
## pmsno            -0.6147     0.5043  -1.219  0.22290   
## genderman         0.6147     0.5043   1.219  0.22290   
## emsno:genderman  -0.3343     0.5268  -0.635  0.52570   
## emsno:pmsno       1.7048     0.5352   3.186  0.00144 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.95608  on 7  degrees of freedom
## Residual deviance:   0.29001  on 2  degrees of freedom
## AIC: 47.553
## 
## Number of Fisher Scoring iterations: 3
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) ) 
## 
## Call:
## glm(formula = 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, family = poisson, data = df4)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
## -0.05579   0.11784   0.18556  -0.26570   0.03860  -0.11339  -0.09224   0.07546  
##        9        10        11        12        13        14        15        16  
##  0.04379  -0.06919  -0.17327   0.21093  -0.05579   0.11784   0.16507  -0.11818  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      2.8467     0.2341  12.162  < 2e-16 ***
## mstatusmarried                  -1.5199     0.4716  -3.223  0.00127 ** 
## emsno                            1.1169     0.2664   4.193 2.76e-05 ***
## pmsno                            0.7304     0.2808   2.601  0.00930 ** 
## genderman                        0.4772     0.2924   1.632  0.10270    
## pmsno:genderman                 -1.2076     0.3857  -3.131  0.00174 ** 
## emsno:genderman                 -0.3242     0.3333  -0.973  0.33071    
## emsno:pmsno                      0.6783     0.3133   2.165  0.03037 *  
## mstatusmarried:emsno             0.8279     0.5041   1.642  0.10053    
## mstatusmarried:pmsno            -0.6147     0.5043  -1.219  0.22290    
## mstatusmarried:genderman         0.6147     0.5043   1.219  0.22290    
## emsno:pmsno:genderman           -0.1183     0.4244  -0.279  0.78048    
## mstatusmarried:emsno:genderman  -0.3343     0.5268  -0.635  0.52570    
## mstatusmarried:emsno:pmsno       1.7048     0.5352   3.186  0.00144 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.85486  on 15  degrees of freedom
## Residual deviance:    0.29001  on  2  degrees of freedom
## AIC: 112.29
## 
## Number of Fisher Scoring iterations: 3

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

summary( logit8 <- glm(p ~ ems + pms + gender + gender*ems + ems*pms + pms*gender,
                  weights=ncases, data=wide.data ,family=binomial) )
## 
## Call:
## glm(formula = p ~ ems + pms + gender + gender * ems + ems * pms + 
##     pms * gender, family = binomial, data = wide.data, weights = ncases)
## 
## Deviance Residuals: 
##        2         4         6         8        10        12        14        16  
##  0.20163  -0.08550  -0.18209   0.03124  -0.12528   0.07135   0.20163  -0.05307  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      -1.5603     0.4880  -3.197  0.00139 **
## emsno             0.8108     0.5104   1.589  0.11215   
## pmsno            -0.5422     0.5419  -1.000  0.31708   
## genderman         0.6704     0.5282   1.269  0.20437   
## emsno:genderman  -0.2920     0.5411  -0.540  0.58942   
## emsno:pmsno       1.6974     0.5383   3.154  0.00161 **
## pmsno:genderman  -0.1283     0.3387  -0.379  0.70495   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.95608  on 7  degrees of freedom
## Residual deviance:   0.14636  on 1  degrees of freedom
## AIC: 49.41
## 
## Number of Fisher Scoring iterations: 3
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) ) 
## 
## Call:
## glm(formula = 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, family = poisson, 
##     data = df4)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
## -0.08535   0.18268   0.04832  -0.07053   0.05925  -0.17218  -0.02421   0.01975  
##        9        10        11        12        13        14        15        16  
##  0.06724  -0.10570  -0.04565   0.05483  -0.08535   0.18268   0.04304  -0.03105  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     2.85384    0.23409  12.191  < 2e-16 ***
## mstatusmarried                 -1.56030    0.48799  -3.197  0.00139 ** 
## emsno                           1.12856    0.26704   4.226 2.38e-05 ***
## pmsno                           0.71979    0.28183   2.554  0.01065 *  
## genderman                       0.46563    0.29378   1.585  0.11297    
## mstatusmarried:emsno            0.81085    0.51041   1.589  0.11215    
## mstatusmarried:pmsno           -0.54216    0.54191  -1.000  0.31708    
## mstatusmarried:genderman        0.67043    0.52822   1.269  0.20437    
## emsno:pmsno                     0.66544    0.31444   2.116  0.03432 *  
## emsno:genderman                -0.34780    0.33888  -1.026  0.30474    
## pmsno:genderman                -1.18541    0.38994  -3.040  0.00237 ** 
## mstatusmarried:emsno:pmsno      1.69743    0.53825   3.154  0.00161 ** 
## mstatusmarried:emsno:genderman -0.29201    0.54108  -0.540  0.58942    
## mstatusmarried:pmsno:genderman -0.12826    0.33875  -0.379  0.70495    
## emsno:pmsno:genderman          -0.08576    0.43297  -0.198  0.84298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1333.85486  on 15  degrees of freedom
## Residual deviance:    0.14636  on  1  degrees of freedom
## AIC: 114.15
## 
## Number of Fisher Scoring iterations: 3