Masahiro KURODA
6.1
Introduction
This chapter presents analysis of log-linear models for contingency tables. In analysis of contingency tables, we are interested in associations and interactions among the categorical variables or interpretations of the parameters describing the model structure. Then our goal is to find a best model such that the model structure is simple and the model has few parameters. The log-linear models which is used in analysis of contingency tables are a generalized linear model for counted data and can easily describe the variety of associations and interactions among the variables. To search a best model, we assess the effects on interaction terms in log-linear models by goodness of fit tests. The methodology for analysis of contingency tables is described in many books, for example, Bishop et al. (1975), Everitt (1977) and Agresti (2002). This chapter is organized as follows: Section 6.2 introduces log-linear models and generalized linear models. Moreover we provide the procedures to find the best model by using goodness of fit tests for independence and comparing two log-linear models. Section 6.3 presents the XploRe functions to make inferences for log-linear models. Contingency table analysis using XploRe are illustrated in Section 6.4.
6.2
Log-linear models
Let Y = ( Y 1 , Y 2 , . . . , Y D ) be categorical variables. Then a rectangular ( N × D ) data matrix consisting of N observations on Y can be rearranged as a D -way contingency table with cells defined by joint levels of the variables. Let n ij ...t denote the frequency for a cell Y = (i,j,...,t) and n = { nij...t} . Suppose that Y has a multinomial distribution with an unknown parameter θ = {θ ij...t}, where θ ij ... t ≥ 0 and ∑ θ ij ... t = 1. The log-linear model is expressed in the form
where X is a D × r design matrix and λ is an r × 1 parameter vector. When Y has a Poisson distribution, the log-linear model is re-written by
where m = { m ij...t = Nθij...t } is the vector of expected frequencies.
6.2.1
Log-linear models for two-way contingency tables
Consider an I × J contingency table. The log-linear model is represented by
for all i and j , under the constraints of the λ terms to sum to zero over any subscript such as
The log-linear model given by (6.3) is called the saturated model or the full model that there is statistically dependence between Y 1 and Y 2 . By analogy with analysis of variance models, we define the overall mean by
the main effects of Y 1 and Y 2 by
and the two-factor effect between Y 1 and Y 2 by
Then the main and two-factor effects are determined by the odds and odds ratios, and can be written by
and
For the independence model that Y 1 is statistically independent of Y 2 , the cell probability θ ij can be factorized into the product of marginal probabilities θ i+ and θ+j, that is,
q=1 θ ij . Then the two-factor effect is θ ij and θ +j
so that the log-linear model for the independence model is expressed by
for all i and j .
6.2.2
Log-linear models for three-way contingency tables
For an I × J × K contingency table, the saturated log-linear model for the contingency table is
for all i, j and k. The λ terms are also satisfied the constraints that
We define the λ terms as follows: The overall mean is given by
The main effects of Y 1 , Y 2 and Y 3 are
Each interaction effect is given by
and
Like log-linear models for two-ways contingency tables, the λ terms in the log- linear models for three-way contingency tables directly relate to the odds and odds ratios. Here we introduce an important class of independence models that are called hierarchical log-linear models. In the hierarchical models, if high-order λ terms with certain variables are contained in the model, all lower-order λ terms with these same variables are included. For instance, when a log-linear model con- tains {λ ij }, the log-linear model also contains {λ i} and {λ ij}. We assume that all log-linear models are the hierarchical model. Table 6.1 is the list of the hierarchical log-linear models for three-way contingency tables. Then in- terpretations of parameters in the log-linear models refer to the highest-order terms. In log-linear models for conditional independence models, the two-factor ef- fects correspond to partial associations.  For instance, the log-linear model [Y 1 Y 2 ][Y 2 Y 3 ] permits two-factor terms for associations between Y 1 and Y 2 , and Y 2 and Y 3 , but does not contain a two-factor term for an association between Y 1 and Y 3 . Then the log-linear model [Y 1 Y 2 ][Y 2 Y 3 ] specifies conditional indepen- dence between Y 1 and Y 3 given Y 2 . In the log-linear model [Y 1 Y 2 ][Y 1 Y 3 ][Y 2 Y 3 ] called the no three-factor interaction model , there exists conditional depen- dence for all three pairs. Then the no three-factor interaction model has equal conditional odds ratios between any two variables at each level of the third variable. For example, the conditional odds ratio of Y 1 to Y 2 in the k-th level of Y 3 does not depend on k , and is given by
for i = 1, . . . , I − 1, j = 1, . . . , J − 1 and all k . With multi-way contingency tables, the independence models are more compli- cated than the models for three-way contingency tables. The log-linear models can also describe easily several models for multi-way contingency tables. The basic principles of log-linear models for three-way contingency tables can be extended readily to multi-way contingency tables.
6.2.3
Generalized linear models
The log-linear model is a special case of generalized linear models (McCullagh and Nelder, 1989). For cell frequencies in contingency tables, the generalized
No three-factor interaction [Y 1 Y 2] [Y 1 Y 3] [Y 2 Y 3]
linear models assume a Poisson distribution as the link function.  Thus the log-linear models are given by equation (6.2). Consider a 2 × 3 contingency table. From the constraints
the parameter vector is identified by Thus the log-linear (6.2) can be written as
Table 6.1: Independence models for three-way contingency tables When the maximum likelihood estimates (MLEs) of λ can not be found directly, iterative algorithms such as the Newton-Raphson and Fisher-scoring algorithms or the iterative proportional fitting procedure are applied. To compute the MLEs λ in log-linear models, XploRe uses the Newton-Raphson and Fisher- scoring algorithms.
6.2.4
Fitting to log-linear models
Chi-squared goodness of fit tests
To test a log-linear model against the saturated model, we estimate the ex- pected frequencies of the log-linear model and evaluate the adequacy by the Pearson chi-squared statistic. When the MLEs Â in a log-linear model are obtained, the expected frequencies are estimated from
To assess a log-linear model fitting to the data by comparing n to ˆ m , the Pearson chi-squared statistic
is computed. As another measure of goodness of fit, the likelihood ratio test statistic is used. This test statistic is computed from If the sample size is sufficiently large, χ 2 and G 2 have an asymptotic chi-squared distribution with degrees of freedom ( df ) equal to the difference of the number of free parameters in the saturated model and a log-linear model. Then the chi-squared goodness of fit test can be conducted by the value of χ 2 or G 2 . Moreover the likelihood ratio test statistic G 2 can be used to compare two log- linear models M 1 and M 2 . Then M 2 is nested in M 1 , such that every nonzero λ terms in M 2 is contained in M 1 . For example, the log-linear model [Y 1 Y 2 ][Y 3 is the nested model in the log-linear model [ Y 1 Y 2 ][ Y 2 Y 3 ] and these models are expressed by
Thus M 2 is simpler than M 1 and M 1 must hold when M 2 holds. Assuming that M 1 holds, we test whether M 2 fits the data as well as M 1 . To compare two models, the following test statistic is used:
where G 2 and G 2 are the likelihood ratio test statistics for M 1 and M 2 . Then △G 2 also has an asymptotic chi-squared distribution with df equal to df for M 1 – df for M 2 . When the value of △ G 2 is in a critical region, we conclude that M 2 provides a better description of the data. Furthermore △ G 2 is computed to compare a nested model in M 2 with M 2 . If the value of △G 2 is out a critical region, we re-compare another nested model in M 1 with M 1 . Repeating goodness of fit tests to compare nested models, we find a best model. As another criteria to compare nested models, the Akaike information criteria (AIC) and Bayesian information criteria (BIC) are well known.
Model residuals
The goodness of fit statistic gives the summary of how a log-linear model fits to the data. We examine lack of fit by comparing observed data to the fitted data individually. For cell ( i, j ) in a two-way contingency table, the Pearson standardized residual is defined by The Pearson residual is also related to the Pearson chi-squared test statistics by When a log-linear model holds, the residuals { e ij } have an approximate normal distribution with mean 0. Then, by checking whether the Pearson residuals are larger than about ±2 that is standard normal percentage points, we detect the presence of the data that are influential on the fit of a log-linear model.
6.3
Inference for log-linear models using XploRe
Estimation of the parameter vector λ
6.3.1
To make inferences for log-linear models, we use the functions in the glm library. The library is available by library("glm") The parameter vector λ can be estimated by using the glmest function glmest("polog", x, n) where polog is a Poisson distribution with logarithm link function, x is the design matrix and n is the cell frequencies for contingency tables. Executing lambda = glmest("polog", x, n) the estimates of λ are assigned to the variable lambda. Then lambda contains the following output: b : the estimated parameter vector λ bv : the estimated covariance of b stat : several statistics The expected frequencies m are also computed from m = exp(x*lambda.b) Moreover the glmest function and other functions in glm library can be also specified several options by defining opt with glmopt. For the optional param- eters, refer to Härdle et al. (1999) or Help function in XploRe.
6.3.2
Computing statistics for the log-linear models
A number of statistical characteristics can be computed using the glmstat function. Then statistical characteristics can be obtained from stat = glmstat("polog", x, m, lambda.b, lambda.bv) and are stored in the output stat: df : degrees of freedom deviance : the deviance of the estimated model pearson : the Pearson statistic loglik : the log-likelihood of the estimated model, using the estimated dis- persion parameter aic, bic : Akaike's AIC and Schwarz' BIC criterion, respectively r2, adr2 : the (pseudo) coefficient of determination and its adjusted version, respectively it : the number of iterations needed ret : the return code nr : the number of replicated observation in x, if they were searched for.
6.3.3
Model comparison and selection
The computation of the likelihood ratio test statistic for comparing nested models can be performed by the glmlrtest function: {lr, pvalue} = glmlrtest(loglik2, df2, loglik1, df1) where loglik1 and loglik2 are the log-likelihoods for the log-linear models M 1 and M 2 . Note that M 2 must be the nested model in M 1 . These values are obtained from the glmstat function. The augments df1 and df2 are dfs for each model. Executing the above call, the test statistic 1r and the p -value pvalue are yielded. Moreover, to select the best model automatically, XploRe has the glmselect, glmforward and glmbackward functions. The glmselect function performs a complete search model selection, the glmforward and glmbackward functions do the forward and backward search model selections, respectively. The syntax of these functions is the same as glmest. Note that best models found by these functions are not always hierarchical log-linear models. Then we repeat to compute the likelihood ratio statistics for comparing the nested models and finally find the best model. When the parameters λ 0 , {λ ; A} and {λ β} are contained in all models, the optional parameter fix that specifies them is described as follows: opt = glmopt("fix", 11213) To search the best model by using the backward search model selection, we type select = glmbackward("polog", x, m, opt) Then the output list select consists of five components: best : the five best models bestcrit : a list containing bestcrit.aic and bestcrit.bic , the Akaike and Schwarz criteria for the five best models bestord : the best models of each order beatordcrit : like bestcrit , but for the best model for each order bestfit : containing bestfit.b, bestfit.bv and bestfit.stat , the estima- tion results for the best model
6.4
Numerical analysis of contingency tables
Testing independence
6.4.1
Chi-squared test The data in Table 6.2 are a cross-sectional study of malignant melanoma taken from Roberts et al. (1981) and treated in Dobson (2001). For the two-way
table, we are interested in whether there exists a association between Tumor type and Site. Table 6.3: Expected frequencies for the independence model
Let m = {mij} be the expected frequencies for the contingency table. The log-linear model that Tumor type is independent of Site is expressed by log m i for all i and j . From the constrain entified by the parameter vector To find the expected frequencies, we estimate the MLEs Â using the following statements: library("glm") n=#(22,2,10,16,54,115,19,33,73,11,17,28) x=read("design.dat") lambda = glmest("polog", x, n) where design.dat is specified by 1 -1 -1 -1 -1 1 -1 -1 -1 1 1 -1 -1 -1 0 0 0 -1 -1 1 1 1 1 1 0 0 1 1 0 0 0 0 1 1 0 -1 1 0 1 0 1 1 0 1 0 0 -1 1 0 0 1 1 0 0 1 1 1 0 0 1 0 1 0 1 -1 0 1 0 1 The expected frequencies shown in Table 6.3 can be obtained by m = exp(x*lambda.b) and are compared to the data in Table 6.2 by using χ 2 . The value of χ 2 is computed from (6.7) 0 lambda.stat stat = glmstat("polog", x, m, lambda.b, lambda.bv) stat.pearson Then χ 2 of 65.8 is very significant compared to the chi-square distribution with 6 df and indicates that the independence model does not fit to the data. We can conclude that there exists the association between Tumor type and Site. Note that the function crosstable provides the chi-squared statistic for testing independence for two-way contingency tables.
Model residuals

Table 6.4: Pearson residuals for the independence model Table 6.4 shows the Pearson standardized residuals for the fit of the indepen- dence model. The values are easily computed from e = (n-m)/sqrt(m) or We can see that the residual for Hutchinson's melanotic freckle and Head & neck reflects the overall poor fit, because the value of 6.75 2 = 45.56 is related to χ² = 65.8.
6.4.2
Model comparison
Chi-squared test
The data in Table 6.5 summarize to a survey the Wright State University school of Medicine and the United Health Services in Dayton, Ohio. The analysis for the contingency table is given in Agresti (2002). For the three-way table, we search the best model by using the likelihood ratio tests.
Table 6.6 shows the expected frequencies for log-linear models of no three-factor interaction and conditional independence models. The expected frequencies for each model are computed by using glmest. The expected frequencies for the log-linear model [ AC ][ AM ][ CM ] are found using the following statements: library("glm") n=#(911, 538, 44, 456, 3, 43, 2, 279) x=read("design.dat") lambda = glmest("polog", x, n) m = exp(x*lambda.b)
Then, under the constraints with the λ terms, the parameter vector λ is iden- tified by
and the design matrix x is specified by 1 1 1 1 1 1 1
To compute the expected frequencies of the nested models in the log-linear model [ AC ][ AM ][ CM ], we delete the columns of x corresponding to λ setting to zero in these models and then execute the above statements. For example, deleting the seventh column of x, we can obtain the expected frequencies of the log-linear model [ AC ][ AM ]. The command
x[,1|2|3|4|5|6]
produces the design matrix for the log-linear model [ AC ][ AM ]. Table 6.7 shows results of the likelihood ratio and Pearson chi-squared tests for log-linear models.  The statements to compute the values of G 2 for the saturated model M 1 and a log-linear model M 2 are stat1 = glmstat("polog", x1, n, lambda1.b, lambda1.bv) df1 = rows(lambda1.b) stat2 = glmstat("polog", x2, n, lambda2.b, lambda2.bv) df2 = rows(lambda2.b) {lr,pvalue}=glmlrtest(stat2.loglik, df2, stat1.loglik, df1) lr pvalue where the design matrix x1 for the saturated model is specified by
The value of χ 2 is also computed from lambda = glmest("polog",x,n) lambda.stat Then the values of G 2 and χ 2 or p-value indicate that the model [ AC ][ AM ][ CM ] fits well to the data.
Model residuals
To examine lack of fit to the data, we analyze the residuals for each log-linear model. Table 6.8 shows the Pearson standardized residuals for the log-linear models. All residuals for the log-linear model [ AC ][ AM ][ CM ] are very small and demonstrate that the model well fits to the data. On the other hand, the residuals for conditional independence models indicate poorly fit to the data. In particular, the extremely large residuals of -34.604 for the model [ AC ][ AM ] and of 34.935 for the model [ AC ][ CM ] cause the lack of fit to the data.
Table 6.8: The Pearson standardized residuals for log-linear models

Test for partial associations
Moreover, it is possible to compare nested log-linear models by testing par- tial associations. We test to compare the model [ AC ][ AM ] with the model [AC][AM][CM]. Then the test examines whether there exists a partial associ- ation between Alcohol use and Cigarette use, that is, Each log-linear model is expressed by
From Table 6.7, △G² = 497.4 – 0.4 = 497.0 and df = 2−1 = 1, so that △G 2 provides strong evidence of a partial association between Cigarette use and Marijuana use. We can also test for partial asso- ciations by comparing the models [ AM ][ CM ] and [ AC ][ CM ] with the model [AC][AM][CM].
Search for the best model
Next we illustrate the best model search using the glmbackward function.  Using opt = glmopt("fix", 1121314) we specify to contain λ 0 , { λ ậ } , { λ f } and { λ f } in the model. To choose a best model, we execute select=glmbackward("polog", x, n, opt) select.best where x is the design matrix for the saturated model. Then select.best displays the five best models for the data:
1 2 3 4 5 0 7 8 In the above output, each row corresponds to the parameter vector λ in the saturated log-linear model as follows:
Those components that are not contained in a log-linear model are indicated by zero. The first column shows the no three-interaction model, since the row 8 is zero. The second column represents the saturated model. The last three columns are not the hierarchical models. Therefore the model [ AC ][ AM ][ CM ] is also selected as the best model. The output select.bestfit includes all estimation results with the best model.
Bibliography
Agresti, A. (2002). Categorical Data Analysis 2nd edition, John Wiley & Sons, New York. Bishop, Y.M., Fienberg, S. and Holland, P. (1975). Discrete Multivariate Anal- ysis, M.I.T. Press, Cambridge, MA. Dobson, A.J. (2001). An Introduction to Generalized Linear Models 2nd edi- tion, Chapman & Hall, London. Everitt, B.S. (1977). The analysis of Contingency Tables, Chapman & Hall, London. Härdle, W., Klinke, S. and Müller, M. (1999). XploRe—Learning Guide , Springer- Verlag, Berlin. McCullagh, P. and Nelder, J.A. (1989).  Generalized Linear Models 2nd edi- tion, Chapman & Hall, London. Roberts, G., Martyn, A.L., Dobson, A.J. and McCarthy, W.H. (1981). Tumour thickness and histological type in malignant melanoma in New South Wales, Australia, 1970-76, Pathology 13: 763-770.
Agresti, A. (2002). Categorical Data Analysis 2nd edition, John Wiley & Sons, New York. Bishop, Y.M., Fienberg, S. and Holland, P. (1975). Discrete Multivariate Anal- ysis, M.I.T. Press, Cambridge, MA. Dobson, A.J. (2001). An Introduction to Generalized Linear Models 2nd edition, Chapman & Hall, London. Everitt, B.S. (1977). The analysis of Contingency Tables , Chapman & Hall, London. Härdle, W., Klinke, S. and Müller, M. (1999). XploRe—Learning Guide, Springer-Verlag, Berlin. McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models 2nd edition, Chapman & Hall, London. Roberts, G., Martyn, A.L., Dobson, A.J. and McCarthy, W.H. (1981). Tu- mour thickness and histological type in malignant melanoma in New South Wales, Australia, 1970-76, Pathology 13: 763-770.