Recall, that the independence can be stated in terms of cell probabilities as a product of marginal probabilities,

${\pi }_{ij}={\pi }_{i+}{\pi }_{+j}i=1,\dots ,I,j=1,\dots ,J$

and in terms of cell frequencies,

${\mu }_{ij}=n{\pi }_{ij}=n{\pi }_{i+}{\pi }_{+j}i=1,\dots ,I,j=1,\dots ,J$

By taking natural logarithms on both sides of "=" sign, we obtain the loglinear model of independence:

$\log ({\mu }_{ij})=\lambda +{\lambda }_i^A+{\lambda }_j^B$

where superscripts A and B are just used to denote the two categorical variables.

This is an ANOVA type-representation where:

• $\lambda$ represents the "overall" effect, or the grand mean of the logarithms of the expected counts, and it ensures that $\sum _i\sum _j{\mu }_{ij}=n$, that is, the expected cell counts under the fitted model add up to the total sample size $n$.
• ${\lambda }_I^A$ represents the "main" effect of variable A, or a deviation from the grand mean, and it ensures that $\sum _j{\mu }_{ij}=n_{i+}$, that is, the marginal totals under the fitted model add up to the observed marginal counts. It represents the effect of classification in row $i$.
• ${\lambda }_J^B$ represents the "main" effect of variable B, or a deviation from the grand mean, and it ensures that $\sum _i{\mu }_{ij}=n+j$. This is the effect of classification in column j.
• ${\lambda }_I^A={\lambda }_J^B=0$, or alternatively, $\sum _i{\lambda }_i^A=\sum _j{\lambda }_j^B=0$, to deal with over-parametrization (see below).

The maximum likelihood (ML) fitted values for the cell counts are the same as the expected (fitted) values under the test of independence in two-way tables, i.e., $E({\mu }_{ij})=n_{i+}{n}_{+j}/n$. Thus, the $X^2$ and $G^2$ for the test of independence are goodness-of-fit statistics for the log-linear model of independence testing that the independence model holds versus that it does not, or more specifically testing that the independence model is true vs. saturated model is true. This model also implies that ALL odds ratios should be equal to 1.

For an $I\times J$ table, and the model is

$\log ({\mu }_{ij})=\lambda +{\lambda }_i^A+{\lambda }_j^B$

There are $I$ terms in the set {${\lambda }_I^A$}, but one of them is redundant, so there are $I-1$ unknown parameters, e.g., {${\lambda }_1^A,\dots ,{\lambda }_{I-1}^A$}, and there are $J$ terms in the set {${\lambda }_J^B$}, and one is redundant, so there are $J-1$ unknown parameters, e.g., {${\lambda }_1^B,\dots ,{\lambda }_{J-1}^B$}. (Why is one of them redundant?) There can be many different parameterizations, but regardless of which set we use, we need to set the constraints to account for redundant parameters. Nonexistence of a unique set of parameters does not mean that the expected cell counts will change depending on which set of parameters is being used. It simply means that the estimates of the effects may be obtained under different sets of constraints, which will lead to different interpretations. But expected cell counts will remain the same.

DUMMY CODING: To avoid over-parametrization, one member in the set of $\lambda$s is fixed to have a constant value, typically 0. This corresponds to using dummy coding for the categorical variables (e.g. A = 1, 0). By default, in SAS PROC GENMOD, the last level is set to 0. So, we have

$\log ({\mu }_{11})=\lambda +{\lambda }_1^A+{\lambda }_1^B$

$\log ({\mu }_{22})=\lambda +0+0=\lambda$

By default, in R glm() the first level of the categorical variable is set to 0. So, we have

$\log ({\mu }_{11})=\lambda +0+0=\lambda$

$\log ({\mu }_{22})=\lambda +{\lambda }_2^A+{\lambda }_2^B$

ANOVA-type CODING: Another way to avoid over-parametrization is to fix the sum of the terms equal to a constant, typically 0. That is the ANOVA-type constraint. This corresponds to using the so-called "effect” coding for categorical variables (e.g. A = 1, 0, −1). By default, SAS PROC CATMOD and R loglin(), use the zero-sum constraint, e.g., the expected cell count in the first cell and the last cell,

$\log ({\mu }_{11})=\lambda +{\lambda }_1^A+{\lambda }_1^B$

$\log ({\mu }_{22})=\lambda -{\lambda }_1^A-{\lambda }_1^B$

We will see more on these with a specific example in the next section.

We can have different parameter estimates (i.e, different values of $\lambda$s) depending on the type of constraints we set. So, what is unique about these parameters that lead to the same inference, regardless of parametrization? The differences, that is the log odds, are unique:

${\lambda }_i^A-{\lambda }_{i^{\prime }}^A$

${\lambda }_j^B-{\lambda }_{j^{\prime }}^B$

where the subscript $i$ denotes one level of categorical variable $A$ and "$i$" denotes another level of the same variable; similarly for $B$.

Thus the odds is also unique!

$$\begin{array}{rl}\log (odds)& =\log \left({\displaystyle \frac{{\mu }_{i1}}{{\mu }_{i2}}}\right)=\log ({\mu }_{i1})-\log ({\mu }_{i2})\\ & =(\lambda +{\lambda }_i^A+{\lambda }_1^B)-(\lambda +{\lambda }_i^A+{\lambda }_2^B)={\lambda }_1^B-{\lambda }_2^B\end{array}$$

$\Rightarrow \text{odds}=\exp ({\lambda }_1^B-{\lambda }_2^B)$

The odds ratio measures the strength of the association and depends only on the interaction terms {${\lambda }_{ij}^{AB}$ }, which clearly does not appear in this model, but we will discuss it when we see the saturated log-linear model.

Do you recall, how many odds ratios do we need to completely characterize associations in $I\times J$ tables?

With the saturated model, the $N=IJ$ counts in the cells are still assumed to be independent observations of a Poisson random variable, but no independence is assumed between the variables $A$ and $B$. The model is expressed as

$\log ({\mu }_{ij})=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_{ij}^{AB}$

Note the additional ${\lambda }_{ij}^{AB}$ interaction term, compared with the independence model. We still have the same ANOVA-like constraints to avoid over-parameterization, and again, this changes with the software used. Ultimately, the inference, in the end, will not depend on the choice of constraints.

Parameter estimates and interpretation

The odds ratio is directly related to the interaction terms. For example, for a $2\times 2$ table,

$$\begin{array}{rl}\log (\theta )& =\log \left({\displaystyle \frac{{\mu }_{11}{\mu }_{22}}{{\mu }_{12}{\mu }_{21}}}\right)\\ & =\\ & ={\lambda }_{11}^{AB}+{\lambda }_{22}^{AB}-{\lambda }_{12}^{AB}-{\lambda }_{21}^{AB}\end{array}$$

How many odds ratios are there? There should be $(I-1)(J-1)$ which should be equal to the unique number of ${\lambda }_{ij}$s in the model.

How many unique parameters are there in this model?

With $IJ=N$ terms, the model is a perfect fit! Moreover, the saturated model includes the independence model as a special case, where the interaction terms are 0.

These models include all lower-order terms that comprise higher-order terms in the model.

(A, B) is a simpler model than (AB)

Interpretation does not depend on how the variables are coded.

We want hierarchical models so that interaction terms represent associations. If there is a significant interaction term, we do NOT look at the lower-order terms, but only interpret the higher-order terms because the value of lower-order terms is coding dependent and can be misleading.

Next, we'll see how these models extend to three-way tables.

This model is the default model in a way that serves for testing of goodness-of-fit of the other models. Recall that the saturated model has the maximum number of parameters and fitting a saturated model is the same as estimating ML parameters of distributions appropriate for each cell of the contingency table.

Main assumption

The $N=IJK$ counts in the cells are assumed to be independent observations of a Poisson random variable.

Model structure

Parameter estimation and interpretation

Model Fit

The saturated model has a perfect fit with $G^2=0$ and df = 0 since the number of cells is equal to the number of unique parameters in the model.

Model Selection

Relevant when comparing to simpler models. The saturated model is the most complex model possible.

Fitting in SAS and R

proc genmod data=berkeley order=data;
class D S A;
model count= D S A D*S D*A S*A D*S*A/dist=poisson link=log;
run;

UCBAdmissions 
berk.data = as.data.frame(UCBAdmissions)
berk.data$Gender = relevel(berk.data$Gender, ref='Female')
berk.data$Dept = relevel(berk.data$Dept, ref='F')
berk.sat = glm(Freq~Admit*Gender*Dept, family=poisson(), data=berk.data)

> summary(berk.sat)

Call:
glm(formula = Freq ~ Admit * Gender * Dept, family = poisson(), 
    data = berk.data)

Deviance Residuals: 
 [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

Coefficients:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     3.17805    0.20412  15.569  < 2e-16 ***
AdmitRejected                   2.58085    0.21171  12.190  < 2e-16 ***
GenderMale                     -0.08701    0.29516  -0.295   0.7682    
DeptA                           1.31058    0.23001   5.698 1.21e-08 ***
DeptB                          -0.34484    0.31700  -1.088   0.2767    
DeptC                           2.13021    0.21591   9.866  < 2e-16 ***
DeptD                           1.69714    0.22204   7.644 2.11e-14 ***
DeptE                           1.36524    0.22870   5.969 2.38e-09 ***
AdmitRejected:GenderMale        0.18890    0.30516   0.619   0.5359    
AdmitRejected:DeptA            -4.12505    0.32968 -12.512  < 2e-16 ***
AdmitRejected:DeptB            -3.33462    0.47817  -6.974 3.09e-12 ***
AdmitRejected:DeptC            -1.92041    0.22876  -8.395  < 2e-16 ***
AdmitRejected:DeptD            -1.95888    0.23781  -8.237  < 2e-16 ***
AdmitRejected:DeptE            -1.42370    0.24250  -5.871 4.33e-09 ***
GenderMale:DeptA                1.83670    0.31672   5.799 6.66e-09 ***
GenderMale:DeptB                3.12027    0.38572   8.090 5.99e-16 ***
GenderMale:DeptC               -0.43376    0.31687  -1.369   0.1710    
GenderMale:DeptD                0.13907    0.31938   0.435   0.6632    
GenderMale:DeptE               -0.48599    0.34151  -1.423   0.1547    
AdmitRejected:GenderMale:DeptA  0.86318    0.40267   2.144   0.0321 *  
AdmitRejected:GenderMale:DeptB  0.03113    0.53349   0.058   0.9535    
AdmitRejected:GenderMale:DeptC -0.31382    0.33741  -0.930   0.3523    
AdmitRejected:GenderMale:DeptD -0.10691    0.34013  -0.314   0.7533    
AdmitRejected:GenderMale:DeptE -0.38908    0.36500  -1.066   0.2864    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2.6501e+03  on 23  degrees of freedom
Residual deviance: 1.4788e-13  on  0  degrees of freedom
AIC: 207.06

The intercept is a normalizing constant and should be ignored. The main effects for D, S, and A are all difficult to interpret and not very meaningful since we have significant two-way and three-way associations; the two- and three-way associations are highly meaningful. For example, the estimated coefficient for the SA association is 0.1889.

Exponentiating this coefficient gives

$\exp (0.1889)=1.208$,

which is the estimated SA odds ratio for Department F since that is the reference department in this analysis. The reference group for S is "female," and the reference group for A is "accept." If we write the $2\times 2$ table for $S\times A$ in Department F, i.e., the "partial" table, with the reference groups in the last row and column, we get

for which the estimated odds ratio is:

$351\times 24/317\times 22=1.208$.

The Wald z-statistic for this coefficient,

$z=0.1889/0.3052$,

which corresponds to Chi-square statistic ${0.62}^2=0.38$ with p-value 0.5359 and indicates that the SA odds ratio for Department F is not significantly different from 1.00 or that the log odds ratio is not significantly different from 0.

The 95% confidence interval for the parameter estimate, that is, for the log odds ratio, is $(-0.4092,0.7870)$. Thus the 95% confidence interval for the odds ratio is

$(\exp (-0.4092),\exp (0.7870))=(0.67,2.20)$.

To get the SA odds ratio for any other department, we have to combine the SA coefficient with one of the DSA coefficients. For example, the SA odds ratio for Department A is

$\exp (0.1889+0.8632)=2.864$.

Based on:

$$\begin{array}{rl}{\theta }_{SA(i="A")}& ={\displaystyle \frac{{\mu }_{ijk}{\mu }_{i{j}^{\prime }{k}^{\prime }}}{{\mu }_{i{j}^{\prime }k}{\mu }_{ij{k}^{\prime }}}}\\ & =\\ & =({\lambda }_{jk}+{\lambda }_{j^{\prime }{k}^{\prime }}-{\lambda }_{j^{\prime }k}-{\lambda }_{j{k}^{\prime }})+({\lambda }_{ijk}+{\lambda }_{i{j}^{\prime }{k}^{\prime }}-{\lambda }_{i{j}^{\prime }k}-{\lambda }_{ij{k}^{\prime }})\\ & =(0.1889+0-0-0)+(0.8632+0-0-0)\end{array}$$

The Wald z-statistic for the first DSA coefficient,

$z=0.8632/0.4027$,

indicates that the SA odds ratio for Department A is significantly different from the SA odds ratio in Department F. To see if the SA odds ratio in Department A is significantly different from 1.00, we would have to compute the standard error for the sum of the two coefficients using the estimated covariance matrix, or refit the model by fixing the level of interest equal to 0.

In many situations, however, we take recourse to the saturated model only as a last resort. As the number of variables grow, saturated models become more and more difficult to interpret. In the following sections, we look at simpler models which are useful in explaining the associations among the discrete variables of interest.

This is the most restrictive model in that all variables are assumed to be jointly independent, regardless of any conditioning. Equivalently, this requires that each two and three-way distribution factors into the product of the marginal distributions involved.

Main assumptions

Model Structure

$\log ({\mu }_{ijk})=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C$

proc genmod data=berkeley order=data;
class D S A;
model count = D S A / dist=poisson link=log;
run;

berk.ind = glm(Freq~Admit+Gender+Dept, family=poisson(), data=berk.data)

> summary(berk.ind)

Call:
glm(formula = Freq ~ Admit + Gender + Dept, family = poisson(), 
    data = berk.data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-18.170   -7.719   -1.008    4.734   17.153  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    4.72072    0.04553 103.673  < 2e-16 ***
AdmitRejected  0.45674    0.03051  14.972  < 2e-16 ***
GenderMale     0.38287    0.03027  12.647  < 2e-16 ***
DeptA          0.26752    0.04972   5.380 7.44e-08 ***
DeptB         -0.19927    0.05577  -3.573 0.000352 ***
DeptC          0.25131    0.04990   5.036 4.74e-07 ***
DeptD          0.10368    0.05161   2.009 0.044533 *  
DeptE         -0.20098    0.05579  -3.602 0.000315 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2650.1  on 23  degrees of freedom
Residual deviance: 2097.7  on 16  degrees of freedom
AIC: 2272.7

Next, let us see an example of the joint independence model.

With three variables, there are three ways to satisfy joint independence. The assumption is that one variable is independent of the other two, but the latter two can have an arbitrary association. For example, if $A$ is jointly independent of $B$ and $C$, which we denote by $(A,BC)$, then $A$ is independent of each of the others, and we can factor the three-way joint distribution into the product of the marginal distribution for $A$ and the two-way joint distribution of $(B,C)$.

Main assumptions

Model Structure

proc genmod data=berkeley order=data;
class D S A;
model count = D S A D*S / dist=poisson link=log lrci type3 obstats;
run;

D*S        DeptA   Male     1    1.9436    0.1268    1.6950    2.1921   234.84

berk.jnt = glm(Freq~Admit+Gender+Dept+Gender*Dept, family=poisson(), data=berk.data)

> summary(berk.jnt)

Call:
glm(formula = Freq ~ Admit + Gender + Dept + Gender * Dept, family = poisson(), 
    data = berk.data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-12.744   -3.208   -0.058    2.495    9.869  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       4.88451    0.05728  85.269  < 2e-16 ***
AdmitRejected     0.45674    0.03051  14.972  < 2e-16 ***
GenderMale        0.08970    0.07492   1.197   0.2312    
DeptA            -1.14975    0.11042 -10.413  < 2e-16 ***
DeptB            -2.61301    0.20720 -12.611  < 2e-16 ***
DeptC             0.55331    0.06796   8.141 3.91e-16 ***
DeptD             0.09504    0.07483   1.270   0.2040    
DeptE             0.14193    0.07401   1.918   0.0551 .  
GenderMale:DeptA  1.94356    0.12683  15.325  < 2e-16 ***
GenderMale:DeptB  3.01937    0.21771  13.869  < 2e-16 ***
GenderMale:DeptC -0.69107    0.10187  -6.784 1.17e-11 ***
GenderMale:DeptD  0.01646    0.10334   0.159   0.8734    
GenderMale:DeptE -0.81123    0.11573  -7.010 2.39e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2650.10  on 23  degrees of freedom
Residual deviance:  877.06  on 11  degrees of freedom
AIC: 1062.1

> fits = fitted(berk.jnt)
> resids = residuals(berk.jnt,type="pearson")
> h = lm.influence(berk.jnt)$hat
> adjresids = resids/sqrt(1-h)
> round(cbind(berk.data$Freq,fits,adjresids),2)
         fits adjresids
1  512 319.90     15.18
2  313 505.10    -15.18
3   89  41.88      9.42
4   19  66.12     -9.42
5  353 217.15     12.59
6  207 342.85    -12.59
...

Next, let us look at what is often the most interesting model of conditional independence.

With three variables, there are three ways to satisfy conditional independence. The assumption is that two variables are independent, given the third. For example, if $A$ and $B$ are conditionally independent, given $C$, which we denote by $(AC,BC)$, then the conditional distribution of $(AB)$, given $C$, can be factored into the product of the two conditional marginals, given $C$.

Main assumptions 

(these are stated for the case $(AC,BC)$ but hold in a similar way for $(AB,BC)$ and $(AB,AC)$ as well):

• The $N=IJK$ counts in the cells are assumed to be independent observations of a Poisson random variable, and
• there is no partial interaction between $A$ and $B$: ${\lambda }_{ij}^{AB}=0$, for all $i,j$, and
• there is no three-way interaction: ${\lambda }_{ijk}^{ABC}=0$ for all $i,j,k$.

Note the constraints above are in addition to the usual set-to-zero or sum-to-zero constraints (present even in the saturated model) imposed to avoid overparameterization.

Model Structure

$\log ({\mu }_{ijk})=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{ik}^{AC}+{\lambda }_{jk}^{BC}$

In the example below, we consider the model where admission status and department are conditionally independent, given sex, which we denote by $(AS,DS)$.

proc genmod data=berkeley order=data;
class D S A;
model count = D S A A*S D*S / dist=poisson link=log lrci type3 obstats;
run;

berk.cnd = glm(Freq~Admit+Gender+Dept+Admit*Gender+Dept*Gender, family=poisson(), data=berk.data)

> fits = fitted(berk.cnd)
> resids = residuals(berk.cnd,type="pearson")
> h = lm.influence(berk.cnd)$hat
> adjresids = resids/sqrt(1-h)
> round(cbind(berk.data$Freq,fits,adjresids),2)
         fits adjresids
1  512 367.28     12.17
2  313 457.72    -12.17
3   89  32.78     12.13
4   19  75.22    -12.13
5  353 249.31      9.91
6  207 310.69     -9.91
...

> anova(berk.cnd, test="LR")
Analysis of Deviance Table

Model: poisson, link: log

Response: Freq

Terms added sequentially (first to last)


             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                            23    2650.10              
Admit         1   230.03        22    2420.07 < 2.2e-16 ***
Gender        1   162.87        21    2257.19 < 2.2e-16 ***
Dept          5   159.52        16    2097.67 < 2.2e-16 ***
Admit:Gender  1    93.45        15    2004.22 < 2.2e-16 ***
Gender:Dept   5  1220.61        10     783.61 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The homogeneous associations model is also known as a model of no-3-way interactions. Denoted by (AB, AC, BC), the only restriction this model imposes over the saturated model is that each pairwise conditional association doesn't depend on the value the third variable is fixed at. For example, the conditional odds ratio between $A$ and $B$, given $C$ is fixed at its first level, must be the same as the conditional odds ratio between $A$ and $B$, given $C$ is fixed at its second level, and so on.

Main assumptions

• The N = IJK counts in the cells are assumed to be independent observations of a Poisson random variable, and
• there is no three-way interaction among the variables: ${\lambda }_{ijk}^{ABC}=0$ for all $i,j,k$.

Model Structure

$\log ({\mu }_{ijk})=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{ij}^{AB}+{\lambda }_{jk}^{BC}+{\lambda }_{ik}^{AC}$

In terms of the Berkeley example, this model implies that the conditional association between department and sex does not depend on the fixed value of admission status, that the conditional association between sex and admission status does not depend on the fixed value of department, and the conditional association between department and admission status does not depend on the fixed value of sex. 

Does this model fit? Even this model doesn't fit well, but it seems to fit better than the previous models. The deviance statistic is $G^2=20.2251$ with df= 5, and the Value/df is 4.0450.

proc genmod data=berkeley order=data;
class D S A;
model count = D S A D*S D*A S*A / dist=poisson link=log lrci type3 obstats;
run;