UCBAdmissions: Fitting models with loglm() and glm()

##                 Dept   A   B   C   D   E   F
## Admit    Gender                             
## Admitted Male        512 353 120 138  53  22
##          Female       89  17 202 131  94  24
## Rejected Male        313 207 205 279 138 351
##          Female       19   8 391 244 299 317

Independence model

berk.mod0 <- loglm(~ Dept + Gender + Admit, data=UCBAdmissions)
berk.mod0

## Call:
## loglm(formula = ~Dept + Gender + Admit, data = UCBAdmissions)
## 
## Statistics:
##                   X^2 df P(> X^2)
## Likelihood Ratio 2098 16        0
## Pearson          2000 16        0

mosaic(berk.mod0, gp=shading_Friendly)

Conditional independence in UCB admissions data

This model allows associations of Dept with both Gender and Admit, but asserts that Gender is independent of Admit, given Dept.

berk.mod1 <- loglm(~ Dept * (Gender + Admit), data=UCBAdmissions)
berk.mod1

## Call:
## loglm(formula = ~Dept * (Gender + Admit), data = UCBAdmissions)
## 
## Statistics:
##                   X^2 df P(> X^2)
## Likelihood Ratio 21.7  6  0.00135
## Pearson          19.9  6  0.00284

mosaic(berk.mod1, gp=shading_Friendly)

All two-way model

berk.mod2 <-loglm(~(Admit+Dept+Gender)^2, data=UCBAdmissions)
berk.mod2

## Call:
## loglm(formula = ~(Admit + Dept + Gender)^2, data = UCBAdmissions)
## 
## Statistics:
##                   X^2 df P(> X^2)
## Likelihood Ratio 20.2  5  0.00114
## Pearson          18.8  5  0.00207

mosaic(berk.mod2, gp=shading_Friendly)

compare models

anova(berk.mod0, berk.mod1, berk.mod2, test="Chisq")

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Dept + Gender + Admit 
## Model 2:
##  ~Dept * (Gender + Admit) 
## Model 3:
##  ~(Admit + Dept + Gender)^2 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1     2097.7 16                                    
## Model 2       21.7  6    2075.94        10        0.00000
## Model 3       20.2  5       1.53         1        0.21593
## Saturated      0.0  0      20.20         5        0.00114

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

Fit the same, using glm()

Need to transform the data to a frequency data.frame

berkeley <- as.data.frame(UCBAdmissions)
berk.glm1 <- glm(Freq ~ Dept * (Gender+Admit), data=berkeley, family="poisson")
summary(berk.glm1)

## 
## Call:
## glm(formula = Freq ~ Dept * (Gender + Admit), family = "poisson", 
##     data = berkeley)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.478  -0.414   0.010   0.309   2.232  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           6.2756     0.0425  147.74  < 2e-16 ***
## DeptB                -0.4057     0.0677   -5.99  2.1e-09 ***
## DeptC                -1.5394     0.0831  -18.54  < 2e-16 ***
## DeptD                -1.3223     0.0816  -16.21  < 2e-16 ***
## DeptE                -2.4028     0.1101  -21.82  < 2e-16 ***
## DeptF                -3.0962     0.1576  -19.65  < 2e-16 ***
## GenderFemale         -2.0333     0.1023  -19.87  < 2e-16 ***
## AdmitRejected        -0.5935     0.0684   -8.68  < 2e-16 ***
## DeptB:GenderFemale   -1.0758     0.2286   -4.71  2.5e-06 ***
## DeptC:GenderFemale    2.6346     0.1234   21.35  < 2e-16 ***
## DeptD:GenderFemale    1.9271     0.1246   15.46  < 2e-16 ***
## DeptE:GenderFemale    2.7548     0.1351   20.39  < 2e-16 ***
## DeptF:GenderFemale    1.9436     0.1268   15.32  < 2e-16 ***
## DeptB:AdmitRejected   0.0506     0.1097    0.46     0.64    
## DeptC:AdmitRejected   1.2091     0.0973   12.43  < 2e-16 ***
## DeptD:AdmitRejected   1.2583     0.1015   12.40  < 2e-16 ***
## DeptE:AdmitRejected   1.6830     0.1173   14.34  < 2e-16 ***
## DeptF:AdmitRejected   3.2691     0.1671   19.57  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2650.095  on 23  degrees of freedom
## Residual deviance:   21.736  on  6  degrees of freedom
## AIC: 216.8
## 
## Number of Fisher Scoring iterations: 4

Mosaic plot

Note the use of formula to reorder factors in the mosaic

mosaic(berk.glm1, gp=shading_Friendly, 
       labeling=labeling_residuals, 
       formula=~Admit+Dept+Gender)

The same, displaying studentized residuals.

mosaic(berk.glm1, residuals_type="rstandard", 
       labeling=labeling_residuals, shade=TRUE, 
       formula=~Admit+Dept+Gender, 
       main="Model: [DeptGender][DeptAdmit]")

## All two-way model
berk.glm2 <- glm(Freq ~ (Dept + Gender + Admit)^2, data=berkeley, family="poisson")
summary(berk.glm2)

## 
## Call:
## glm(formula = Freq ~ (Dept + Gender + Admit)^2, family = "poisson", 
##     data = berkeley)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8        9       10       11       12  
## -0.7548   0.9947   1.9645  -3.1577  -0.0340   0.0445   0.1571  -0.2203   1.0127  -0.7384  -0.7437   0.5490  
##      13       14       15       16       17       18       19       20       21       22       23       24  
##  0.0676  -0.0474  -0.0691   0.0508   1.0558  -0.6124  -0.7362   0.4268  -0.2012   0.0511   0.1980  -0.0537  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  6.2715     0.0427  146.85  < 2e-16 ***
## DeptB                       -0.4032     0.0678   -5.94  2.8e-09 ***
## DeptC                       -1.5779     0.0895  -17.63  < 2e-16 ***
## DeptD                       -1.3500     0.0853  -15.83  < 2e-16 ***
## DeptE                       -2.4498     0.1176  -20.84  < 2e-16 ***
## DeptF                       -3.1379     0.1617  -19.40  < 2e-16 ***
## GenderFemale                -1.9986     0.1059  -18.87  < 2e-16 ***
## AdmitRejected               -0.5821     0.0690   -8.44  < 2e-16 ***
## DeptB:GenderFemale          -1.0748     0.2286   -4.70  2.6e-06 ***
## DeptC:GenderFemale           2.6651     0.1261   21.14  < 2e-16 ***
## DeptD:GenderFemale           1.9583     0.1273   15.38  < 2e-16 ***
## DeptE:GenderFemale           2.7952     0.1393   20.07  < 2e-16 ***
## DeptF:GenderFemale           2.0023     0.1357   14.75  < 2e-16 ***
## DeptB:AdmitRejected          0.0434     0.1098    0.40     0.69    
## DeptC:AdmitRejected          1.2626     0.1066   11.84  < 2e-16 ***
## DeptD:AdmitRejected          1.2946     0.1058   12.23  < 2e-16 ***
## DeptE:AdmitRejected          1.7393     0.1261   13.79  < 2e-16 ***
## DeptF:AdmitRejected          3.3065     0.1700   19.45  < 2e-16 ***
## GenderFemale:AdmitRejected  -0.0999     0.0808   -1.24     0.22    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2650.095  on 23  degrees of freedom
## Residual deviance:   20.204  on  5  degrees of freedom
## AIC: 217.3
## 
## Number of Fisher Scoring iterations: 4

mosaic.glm(berk.glm2, residuals_type="rstandard", 
           labeling = labeling_residuals, shade=TRUE,
           formula=~Admit+Dept+Gender, 
           main="Model: [DeptGender][DeptAdmit][AdmitGender]")

anova(berk.glm1, berk.glm2, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: Freq ~ Dept * (Gender + Admit)
## Model 2: Freq ~ (Dept + Gender + Admit)^2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         6       21.7                     
## 2         5       20.2  1     1.53     0.22

Add 1 df term for association of [GenderAdmit] only in Dept A

berkeley <- within(berkeley, dept1AG <- (Dept=='A')*(Gender=='Female')*(Admit=='Admitted'))
berkeley[1:6,]

##      Admit Gender Dept Freq dept1AG
## 1 Admitted   Male    A  512       0
## 2 Rejected   Male    A  313       0
## 3 Admitted Female    A   89       1
## 4 Rejected Female    A   19       0
## 5 Admitted   Male    B  353       0
## 6 Rejected   Male    B  207       0

berk.glm3 <- glm(Freq ~ Dept * (Gender+Admit) + dept1AG, 
                 data=berkeley, family="poisson")
summary(berk.glm3)

## 
## Call:
## glm(formula = Freq ~ Dept * (Gender + Admit) + dept1AG, family = "poisson", 
##     data = berkeley)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8        9       10       11       12  
##  0.0000   0.0000   0.0000   0.0000  -0.0632   0.0827   0.2951  -0.4009   0.5573  -0.4152  -0.4182   0.3051  
##      13       14       15       16       17       18       19       20       21       22       23       24  
## -0.3066   0.2184   0.3204  -0.2314   0.6984  -0.4142  -0.4992   0.2863  -0.4203   0.1086   0.4268  -0.1138  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           6.2383     0.0442  141.16  < 2e-16 ***
## DeptB                -0.3685     0.0688   -5.36  8.5e-08 ***
## DeptC                -1.5021     0.0839  -17.89  < 2e-16 ***
## DeptD                -1.2851     0.0825  -15.58  < 2e-16 ***
## DeptE                -2.3655     0.1108  -21.35  < 2e-16 ***
## DeptF                -3.0590     0.1580  -19.36  < 2e-16 ***
## GenderFemale         -2.8018     0.2363  -11.86  < 2e-16 ***
## AdmitRejected        -0.4921     0.0717   -6.86  6.9e-12 ***
## dept1AG               1.0521     0.2627    4.00  6.2e-05 ***
## DeptB:GenderFemale   -0.3073     0.3124   -0.98     0.33    
## DeptC:GenderFemale    3.4031     0.2461   13.83  < 2e-16 ***
## DeptD:GenderFemale    2.6956     0.2468   10.92  < 2e-16 ***
## DeptE:GenderFemale    3.5233     0.2522   13.97  < 2e-16 ***
## DeptF:GenderFemale    2.7121     0.2479   10.94  < 2e-16 ***
## DeptB:AdmitRejected  -0.0507     0.1118   -0.45     0.65    
## DeptC:AdmitRejected   1.1078     0.0997   11.12  < 2e-16 ***
## DeptD:AdmitRejected   1.1570     0.1038   11.14  < 2e-16 ***
## DeptE:AdmitRejected   1.5816     0.1193   13.25  < 2e-16 ***
## DeptF:AdmitRejected   3.1678     0.1685   18.80  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2650.0952  on 23  degrees of freedom
## Residual deviance:    2.6815  on  5  degrees of freedom
## AIC: 199.7
## 
## Number of Fisher Scoring iterations: 3

mosaic.glm(berk.glm3, residuals_type="rstandard", 
           labeling = labeling_residuals, shade=TRUE,
           formula=~Admit+Dept+Gender, 
           main="Model: [DeptGender][DeptAdmit] + DeptA*[GA]")

anova(berk.glm1, berk.glm3, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: Freq ~ Dept * (Gender + Admit)
## Model 2: Freq ~ Dept * (Gender + Admit) + dept1AG
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)    
## 1         6      21.74                         
## 2         5       2.68  1     19.1  1.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1