loglin is used to fit log-linear models to multidimensional
  contingency tables by Iterative Proportional Fitting.
loglin(table, margin, start = rep(1, length(table)), fit = FALSE,
       eps = 0.1, iter = 20, param = FALSE, print = TRUE)a contingency table to be fit, typically the output from
    table.
a list of vectors with the marginal totals to be fit.
(Hierarchical) log-linear models can be specified in terms of these
    marginal totals which give the ‘maximal’ factor subsets contained
    in the model.  For example, in a three-factor model,
    list(c(1, 2), c(1, 3)) specifies a model which contains
    parameters for the grand mean, each factor, and the 1-2 and 1-3
    interactions, respectively (but no 2-3 or 1-2-3 interaction), i.e.,
    a model where factors 2 and 3 are independent conditional on factor
    1 (sometimes represented as ‘[12][13]’).
The names of factors (i.e., names(dimnames(table))) may be
    used rather than numeric indices.
a starting estimate for the fitted table.  This optional
    argument is important for incomplete tables with structural zeros
    in table which should be preserved in the fit.  In this
    case, the corresponding entries in start should be zero and
    the others can be taken as one.
a logical indicating whether the fitted values should be returned.
maximum deviation allowed between observed and fitted margins.
maximum number of iterations.
a logical indicating whether the parameter values should be returned.
a logical.  If TRUE, the number of iterations and
    the final deviation are printed.
A list with the following components.
the Likelihood Ratio Test statistic.
the Pearson test statistic (X-squared).
the degrees of freedom for the fitted model. There is no adjustment for structural zeros.
list of the margins that were fit.  Basically the same
    as the input margin, but with numbers replaced by names
    where possible.
An array like table containing the fitted values.
    Only returned if fit is TRUE.
A list containing the estimated parameters of the
    model.  The ‘standard’ constraints of zero marginal sums
    (e.g., zero row and column sums for a two factor parameter) are
    employed.  Only returned if param is TRUE.
The Iterative Proportional Fitting algorithm as presented in
  Haberman (1972) is used for fitting the model.  At most iter
  iterations are performed, convergence is taken to occur when the
  maximum deviation between observed and fitted margins is less than
  eps.  All internal computations are done in double precision;
  there is no limit on the number of factors (the dimension of the
  table) in the model.
Assuming that there are no structural zeros, both the Likelihood
  Ratio Test and Pearson test statistics have an asymptotic chi-squared
  distribution with df degrees of freedom.
Note that the IPF steps are applied to the factors in the order given
  in margin.  Hence if the model is decomposable and the order
  given in margin is a running intersection property ordering
  then IPF will converge in one iteration.
Package MASS contains loglm, a front-end to
  loglin which allows the log-linear model to be specified and
  fitted in a formula-based manner similar to that of other fitting
  functions such as lm or glm.
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988). The New S Language. Wadsworth & Brooks/Cole.
Haberman, S. J. (1972). Algorithm AS 51: Log-linear fit for contingency tables. Applied Statistics, 21, 218--225. 10.2307/2346506.
Agresti, A. (1990). Categorical data analysis. New York: Wiley.
loglm in package MASS for a
  user-friendly wrapper.
glm for another way to fit log-linear models.
# NOT RUN {
## Model of joint independence of sex from hair and eye color.
fm <- loglin(HairEyeColor, list(c(1, 2), c(1, 3), c(2, 3)))
fm
1 - pchisq(fm$lrt, fm$df)
## Model with no three-factor interactions fits well.
# }
Run the code above in your browser using DataLab