Statistical Modeling of Interlocus Interactions in a Complex Disease: Rejection of the Multiplicative Model of Epistasis in Type 1 Diabetes

Corresponding author: Heather J. Cordell, Department of Medical Genetics, University of Cambridge, Wellcome Trust Centre for Molecular Mechanisms in Disease, Wellcome Trust/MRC Bldg., Addenbrooke's Hospital, Hills Rd., Cambridge CB2 2XY, United Kingdom., heather.cordell{at}cimr.cam.ac.uk (E-mail)

Communicating editor: Z-B. ZENG

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

In general, common diseases do not follow a Mendelian inheritance pattern. To identify disease mechanisms and etiology, their genetic dissection may be assisted by evaluation of linkage in mouse models of human disease. Statistical modeling of multiple-locus linkage data from the nonobese diabetic (NOD) mouse model of type 1 diabetes has previously provided evidence for epistasis between alleles of several Idd (insulin-dependent diabetes) loci. The construction of NOD congenic strains containing selected segments of the diabetes-resistant strain genome allows analysis of the joint effects of alleles of different loci in isolation, without the complication of other segregating Idd loci. In this article, we analyze data from congenic strains carrying two chromosome intervals (a double congenic strain) for two pairs of loci: Idd3 and Idd10 and Idd3 and Idd5. The joint action of both pairs is consistent with models of additivity on either the log odds of the penetrance, or the liability scale, rather than with the previously proposed multiplicative model of epistasis. For Idd3 and Idd5 we would also not reject a model of additivity on the penetrance scale, which might indicate a disease model mediated by more than one pathway leading to ß-cell destruction and development of diabetes. However, there has been confusion between different definitions of interaction or epistasis as used in the biological, statistical, epidemiological, and quantitative and human genetics fields. The degree to which statistical analyses can elucidate underlying biologic mechanisms may be limited and may require prior knowledge of the underlying etiology.

MUCH effort has been invested in the mapping and identification of loci that predispose to common multifactorial diseases, such as type 1 diabetes. However, little information is available as to the nature of interaction between genes. It is hoped that the identification of the mode of gene interaction will facilitate understanding of the pathological mechanisms involved in complex diseases, as well as the further identification of disease susceptibility loci. One major challenge in complex disease genetics is to be able to detect, at a reasonable statistical level, genes that alone have small effects on the disease phenotype. The chances of detecting such small effects may be increased when the interaction of one such gene with another is taken into consideration. For example, in type 1 diabetes (CORDELL et al. 1995 , CORDELL et al. 2000 ), type 2 diabetes (COX et al. 1999 ), and inflammatory bowel disease (CHO et al. 1998 ), identification of the most suitable model of interaction between loci provided increased evidence for linkage at one locus when the interaction at another locus was taken into consideration.

The modeling of interlocus effects in human complex disease is still in its infancy (COX et al. 1999 ), but there is some doubt that it offers significant advantages except under special conditions (LEAL and OTT 2000 ). It is challenging, owing to the size of the data set required to provide power to distinguish between different genetic models and the fact that other disease loci are segregating in the study population and environmental factors may be operating. Nevertheless, two-locus genetic modeling has been reported, for example, between IDDM1 and IDDM2 and between IDDM1 and IDDM4, using affected sib-pair data (CORDELL et al. 1995 ). In that data set, the IDDM1/IDDM2 joint effect was well described by the multiplicative model, while the IDDM1/IDDM4 effect followed a heterogeneity model.

The study of interlocus interactions in complex disease has been confused by differences in definition and terminology between biologists, epidemiologists, and statisticians, and between quantitative and human geneticists. The term "epistatic" was originally introduced by BATESON 1909 to describe a "masking" effect in which a factor at one Mendelian locus prevents another from manifesting its effect. This definition is perhaps closest to the phenomenon of interest to a biologist or biochemist investigating interaction between proteins in a mechanistic sense. In fact, the concept of biological interaction is often not precisely defined (GREENLAND and ROTHMAN 1998 ) but usually corresponds to such a situation in which the qualitative nature of the mechanism of action of each of two factors is affected by the presence or absence of the other (SIEMIATYCKI and THOMAS 1981 ).

In quantitative genetics the term epistatic has classically been used to refer to a deviation from additivity in the effects of alleles at different loci with respect to prediction of a quantitative phenotype. This definition is due to FISHER 1918 , but note that this is not equivalent to the BATESON 1909 definition, as pointed out by R. C. Punnett in his review of the FISHER 1918 article (NORTON and PEARSON 1976 ). Epistasis in the FISHER 1918 sense is similar to the usual concept of statistical interaction: departure from a (specific) linear model describing the relationship between predictive factors (here assumed to be alleles at different genetic loci) and an outcome or phenotype of interest. Note that with this definition, the choice of scale becomes important as factors that are additive with respect to an outcome measured on one scale may exhibit interaction when a different, transformed scale is used (FRANKEL and SCHORK 1996 ; GREENLAND and ROTHMAN 1998 ).

In human genetics, three main models of gene interaction for the penetrance (the probability of developing disease given genotype) are commonly considered (RISCH 1990 ): a heterogeneity model (RISCH 1990 ; NEUMAN and RICE 1992 ), which is a generalization of a model in which an individual becomes affected through possessing a predisposing genotype at either of two loci; an additive model, which has been shown to approximate the heterogeneity model when modeling familial relative risks (RISCH 1990 ; CORDELL et al. 1995 ); and a multiplicative model (HODGE 1981 ). The additive and heterogeneity models are usually assumed to represent nonepistatic models and to correspond to a situation in which the biological pathways involved in disease are at some level separate or "independent." The multiplicative model is usually considered to be an epistatic model in which the loci and pathways involved are not independent; note, however, that a multiplicative model can be considered to be an additive model when transformed to the logarithmic scale. In a statistical sense, therefore, the multiplicative model signifies independent additive effects on a logarithmic scale. A fourth model is one of additivity on a liability or probit scale, where loci contribute to an underlying, unobserved, continuous trait in an additive fashion and development of disease occurs if this trait exceeds a certain threshold (PEARSON 1900 ; FISHER 1930 ; WRIGHT 1934A , WRIGHT 1934B ). The question that inevitably arises is whether there is some "natural" scale (such as the probit scale, or alternatively the raw penetrance scale) on which models have a specific biological or causal interpretation. We postpone further discussion of this question until later.

One approach to the dissection of complex disease genetics (and the modeling of gene interactions in complex diseases) is to exploit a rodent model of a human complex disease, such as the model of human type 1 diabetes, the nonobese diabetic (NOD) mouse. Spontaneous diabetes in this inbred mouse strain has a similar etiology to human type 1 diabetes, both in terms of the major physiological features of the disease and also some shared genetic determinants, notably at the MHC (MAKINO et al. 1980 ; LYONS and WICKER 1999 ). The advantage that NOD mouse genetic analysis holds is that large backcross and intercross pedigrees can be bred and phenotyped under stable environmental conditions, allowing the identification of disease linkages with convincing statistical support (LYONS and WICKER 1999 ). It is notable that our previous modeling of multiple locus linkage data from such an NOD cross provided evidence for the multiplicative model of epistasis on the penetrance scale (RISCH et al. 1993 ). This model gave very similar results to the classic polygenic threshold liability model, which is often assumed in the inheritance of complex traits. However, in this previous study, all of the complex gene-gene interactions were assessed simultaneously and the nature of interactions between specific pairs of loci was not examined in detail.

To define interactions between specific pairs of loci, double congenic mouse strains may be developed. For standard single congenic strains, specific chromosome intervals from one inbred strain [the donor, in this case the diabetes-resistant strains C57BL/6 (B6) or C57BL/10 (B10)] are introgressed onto the background of the recipient strain (the diabetes-susceptible NOD strain) by backcrossing (WICKER et al. 1994 ). Conversely, NOD chromosome segments can be introgressed into the diabetes-resistant strain to try to convert a resistant strain into a sensitive one (YUI et al. 1996 ). Allelically variable markers are used to guide the strain construction and to make sure that after several generations of backcrossing only the desired chromosome segment is of B6 or B10 origin. The influence of this introgressed ("congenic") segment on disease frequency can then be evaluated by creating the homozygous strain and monitoring disease in a cohort of NOD congenic mice compared with a cohort of NOD parental mice. If an introgressed interval contains an Idd locus then the NOD congenic strain will be protected from diabetes (WICKER et al. 1994 ; LYONS and WICKER 1999 ). As an extension of this strategy, double congenic strains have now been developed, where a single strain possesses two well-defined congenic regions derived from two separate single congenic strains. These double congenic strains allow a more specific assessment of the mode of interaction between two chromosome regions (compared to backcross data), as the genetic background is not a variable and sufficient numbers of mice can be bred to achieve enough statistical power to distinguish between models. In the work presented here, we report on the analysis of previously published data from an NOD double congenic strain for Idd3 and Idd10 on chromosome 3 (WICKER et al. 1994 ), and on data from an NOD double congenic strain for Idd5 and Idd3 on chromosomes 1 and 3, respectively (HILL et al. 2000 ). The diabetes frequencies in both single and double congenic strains were modeled mathematically on linear, log odds, and liability scales using several statistical models. Given our previous multiple-locus results, we were surprised to find strong evidence against multiplicativity for both pairs of loci, even though models that were additive on a log odds or liability scale fitted the data well.

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

The development of the congenic strains analyzed here has been previously described by WICKER et al. 1994, PODOLIN et al. 1997 , and HILL et al. 2000 , but no specific modeling of epistasis was undertaken by these authors. Here we model the joint action of two pairs of type 1 diabetes susceptibility loci, Idd3 and Idd10 and Idd3 and Idd5, in causing disease. For Idd3 and Idd10, F₁ data were used so that animals could be classified into nine categories: homozygous for NOD at both loci, homozygous for NOD at Idd3 and heterozygous for NOD/B6 at Idd10, homozygous for NOD at Idd3 and homozygous for B6 at Idd10, heterozygous at Idd3 and homozygous for NOD at Idd10, heterozygous at both loci, heterozygous at Idd3 and homozygous for B6 at Idd10, homozygous for B6 at Idd3 and homozygous for NOD at Idd10, homozygous for B6 at Idd3 and heterozygous at Idd10, and homozygous for B6 at both loci. Two different strains, NOD.B6Idd10^R1 and NOD. B6Idd10^R2, were used to generate the Idd10 congenics (PODOLIN et al. 1997 ). These strains differ in their development of diabetes due to differences in the length of the congenic region around the Idd10 locus, with the NOD.B6Idd10^R1 congenic region in fact corresponding to two Idd loci, Idd10 and Idd17 (PODOLIN et al. 1997 ).

Counts of the number of mice developing diabetes in the nine categories for Idd3 and Idd10 are given in Table 1. Counts for comparable NOD double congenic strains for Idd3 and Idd5 (HILL et al. 2000 ) are given in Table 2. For Idd3 and Idd5, data were available only in the homozygous categories, and the non-NOD-derived allele for Idd5 was derived from the B10 rather than the B6 strain. Given data in the form of Table 1 or Table 2, mathematical models for the joint effect of alleles at the two loci were fitted.

Models for penetrance:
Consider modeling the data for the nine genotype categories in Table 1. Let p_ij be the probability that an animal develops disease given that it has genotype i at locus 1 and j at locus 2, where i and j take values from 0 to 2 corresponding to the number of B (B6 or B10) alleles in the genotype. The "raw" estimate of p₀₀ for strain NOD. B6Idd10^R1 in Table 1 is therefore 63/81 or 0.778, for example. Since p_ij is a probability, any realistic model requires that any estimate of p_ij be constrained to lie in the interval [0, 1] for all i, j.

Additive model: The additive model as used in human genetics is usually parameterized as p_ij = x_i + y_j, where x_i and y_j are parameters to be estimated representing the contributions of the different genotypes at loci 1 and 2, respectively (RISCH 1990 ). This model is equivalent to the standard quantitative genetics model

(COCKERHAM and ZENG 1996 ), where now µ, a₁, d₁, a₂, d₂ represent genetic parameters to be estimated corresponding to the mean effect and additive and dominance effects at loci 1 and 2, and x_i and z_i are dummy variables taking values x_i = 1 for NN, 0 for NB, and -1 for BB, and z_i = -0.5 for NN, 0.5 for NB, and -0.5 for BB at locus i. Here y is a quantitative phenotype that we assume for the moment corresponds to the genotype-specific penetrance p_ij. Note that the quantities denoted x have different meanings in the two parameterizations, which may be a little confusing but have been chosen to correspond to the formulations previously used in the literature. Some elementary algebra verifies that the two formulations are completely equivalent, illustrating that there are in fact only five free parameters corresponding to µ, a₁, d₁, a₂, d₂ in formulation 2 (i.e., there is a linear constraint between the six parameters x₀, x₁, x₂, y₀, y₁, y₂ in formulation 1).

Heterogeneity model: In human genetics, the additive model is often considered to be a good approximation to a model of genetic heterogeneity, in which loci 1 and 2 are considered to be independent causes of disease (RISCH 1990 ). In the heterogeneity model, p_ij is written as p_ij = x_i + y_j - x_iy_j (NEUMAN and RICE 1992 ). When x_i and y_j take values x₀ = y₀ = 1, x₁ = y₁ = x₂ = y₂ = 0, or alternatively x₀ = y₀ = x₁ = y₁ = 1, x₂ = y₂ = 0, this leads to penetrance matrices of the form given in Table 3, which clearly correspond to a "classical" heterogeneity model in which an individual can be affected through possessing a predisposing genotype at either locus. When the penetrances are not constrained to be 0 or 1, however, the biological interpretation of the generalized heterogeneity model p_ij = x_i + y_j - x_iy_j is less clear. RISCH 1990 and CORDELL et al. 1995 have shown that the heterogeneity and additive models for the penetrance give very similar results when used to model familial relative risks of disease. This may not be the case, however, when modeling the penetrances directly; e.g., note that the penetrance matrices in Table 3 cannot be achieved using an additive model. SCHORK et al. 1993 have also considered heterogeneity models that have a parameterization equivalent to the RISCH 1990 parameterization if the parameter in the SCHORK et al. 1993 formulation is set equal to 1.

The heterogeneity model used in human genetics does not have a direct equivalent in the quantitative genetics literature. Note, however, that the model may be written as log (1 - p_ij) = log(1 - x_i) + log(1 - y_j); i.e., it may be considered to be an additive or nonepistatic model on the log(1 - p_ij) scale. This illustrates that, like the additive model, it must have only five free parameters.

Multiplicative model: The multiplicative model (HODGE 1981 ) may be written as p_ij = x_iy_j. This model is usually considered to represent a form of epistasis between loci 1 and 2, but note that on a logarithmic scale it is equivalent to an additive model and, therefore, again has only five free parameters.

General (epistatic) model: Restricted models for the penetrance may be compared to a general epistatic model that corresponds to a saturated model in which we estimate nine parameters: p₁₁, p₁₂, p₁₃, p₂₁, p₂₂, p₂₃, p₃₁, p₃₂, and p₃₃. In the quantitative genetics literature, this model is usually written

(1)

(COCKERHAM and ZENG 1996 ), where x_i and z_i are dummy variables as defined previously. The parameters µ, a₁, d₁, a₂, d₂ are as described previously, while i_aa, i_ad, i_da, and i_dd correspond to interaction effects. The advantage of this formulation is that, provided the parameter estimates are approximately independent, we can estimate and test each effect in the full model directly; i.e., we can examine the interactions between specific alleles at the two loci rather than merely testing the overall existence of any interaction between two loci. Again some elementary algebra verifies that the two formulations are equivalent.

Application to data in Table 2:
The data in Table 2 can be modeled in a similar way to that of Table 1. Note that the data in Table 2 have less degrees of freedom (d.f.), with the general model having four free parameters and the restricted models having three free parameters. The standard quantitative genetics model (WEBER et al. 1999 ; ZENG et al. 2000 ) for these data would be

where x_i takes values 0.5 for NN and -0.5 for BB at locus i, and µ, a₁, a₂, and i_aa are the genetic parameters to be estimated, with a₁ and a₂ corresponding to the main effects of loci 1 and 2 and i_aa corresponding to the interaction effect. Since Table 2 contains data only for animals homozygous at both loci, we are unable to model the differences in effect between having a single copy or two copies of a particular allele at a locus (i.e., dominance effects).

Models for the log odds:
In epidemiological studies, rather than modeling the penetrance directly, a more common measure is the natural logarithm of the odds, log(p/(1 - p)), which has the advantage of not being constrained to the interval [0, 1]. The standard epidemiological procedure is to fit an additive model to the log odds so that we write log(p_ij/(1 - p_ij)) = x_i + y_j or, equivalently, p_ij = . This model allows us to model the effects at loci 1 and 2 as independent (in a statistical sense) additive effects on the log odds scale; note, however, that it leads to interactive effects (epistasis) on the penetrance scale. We may also consider fitting heterogeneity and multiplicative models to the log odds in the same way as to the penetrances, although it is unclear what biological meaning should be attached to such models. Note that if all the penetrances are small, an additive model for the log odds should be equivalent to a multiplicative model for the penetrance since x_i + y_j = log() log(p_ij) and so p_ij e^xie^yj = X_iY_j, say.

Liability models:
Another model we consider is a liability or probit model similar to that described by RISCH et al. 1993 . This model corresponds to the classical polygenic threshold model (PEARSON 1900 ; FISHER 1930 ; WRIGHT 1934A , WRIGHT 1934B ) in which each disease susceptibility locus contributes to W (an underlying unobserved continuous variable) in an additive fashion. Around each genotype contribution is a residual variance, which here, as in RISCH et al. 1993 , is assumed to equal 1.0. A threshold V is defined so that the area beyond the threshold corresponds to the overall frequency of affected individuals. The probabilities p_ij can be written p_ij = 1 - (V, W, 1) = 1 - (V, x_i + y_j, 1), where (V, µ, ²) is defined to be the cumulative distribution function of a random variable that is normally distributed with mean µ and variance ², i.e., the probability that such a variable takes a value less than V.

The liability model used here is slightly more general than the one described by RISCH et al. 1993 , which was parameterized by a single parameter for each locus, leading to two rather than three free parameters for their data, which were in a similar form to Table 2. [Note also that RISCH et al. 1993 model heterozygote effects rather than B10 homozygous effects since they consider data generated from a backcross rather than from a congenic strain.] The RISCH et al. 1993 parameterization essentially forces the value of W in category four of Table 2 to equal -1 times the value of W in category one. Our parameterization here is easier to extend to the data for all nine categories in Table 1. In addition it has the following useful property: The probabilities (V, x_i + y_j, 1) can be written as (-x_i -y_j, -V, 1) and it can be shown that fitting the model p_ij = 1 - (-x - y_j, -V, 1) is invariant to the choice of V; i.e., we can use the standard cumulative normal distribution function (-x_i -y_j, 0, 1) and estimation of V is not required.

A natural extension to the additive liability model just described would be to consider heterogeneity or multiplicative effects for the liability; e.g., p_ij = 1 - (V, x_i + y_j - x_iy_j, 1) or p_ij = 1 - (V, x_iy_j, 1). Unfortunately these formulations do not have the invariant property of the additive liability model and, moreover, do not have an obvious genetic interpretation. We therefore do not present results from these analyses. Note that although the additive liability model is additive on the liability scale, it leads to interactive effects (epistasis) on the penetrance scale.

Fitting the likelihood:
Given a model for p_ij and data in all relevant genotype categories, the likelihood for the data may be written as

where a_ij and u_ij are the numbers of affected and unaffected animals in genotype category ij, respectively. E.g., for fitting an additive model to the data for each strain in Table 2 the likelihoods may be written

• NOD: (x₀ + y₀)⁵⁵(1 - x₀ - y₀)¹⁸

• Idd5 congenic strain: (x₀ + y₂)⁴²(1 - x₀ - y₂)⁴⁸

• Idd3 congenic strain: (x₂ + y₀)¹²(1 - x₂ - y₀)⁴⁷

• Idd3/5 double congenic strain: (x₂ + y₂)²(1 - x₂ - y₂)⁸⁹.

The overall likelihood may be calculated as the product of the likelihoods for each of the four strains. The likelihood may then be maximized with respect to the parameters to be estimated. Standard statistical theory predicts that twice the difference between the natural logarithms of the maximized likelihoods for nested models should be distributed as a ² with degrees of freedom equal to the difference in the number of estimated parameters. The maximized log-likelihood for a restricted model (additive, heterogeneity, multiplicative) can therefore be compared to that for the general unrestricted (saturated) model, allowing a test for the goodness of fit of the restricted model to be performed.

Generalized linear models:
It is worth noting that the models described here can all be considered as generalized linear models (NELDER and WEDDERBURN 1972 ; MCCULLAGH and NELDER 1989 ). In the generalized linear model framework, the probability p_ij is related to a linear predictor via a link function: g(p) = µ + x_i + y_j. Here the link functions g(p) correspond to the previously described models as follows

• Additive model for penetrance: g(p) = p

• Multiplicative model for penetrance: g(p) = log(p)

• Heterogeneity model for penetrance: g(p) = -log(1 - p)

• Additive model for liability: g(p) = probit(p)

• Additive model for log odds: logistic regression: g(p) = logit(p).

These models can be fitted using standard routines in most statistics packages (e.g., PROC GENMOD in SAS; glm in S-Plus; GLM in Stata, GLIM, etc.). In the epidemiological literature a number of more general parametric families of link function have been proposed (THOMAS 1981 ; GUERRERO and JOHNSON 1982 ; BRESLOW and STORER 1985 ).

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Modeling the joint effects of Idd3 and Idd10:
Table 4 shows the results for fitting the models described above to the data from Table 1 using the strain NOD. B6Idd10^R1. Table 5 shows the results using the strain NOD.B6Idd10^R2. Results are given in terms of fitted values for the penetrances, and differences in -2 ln likelihoods and P values for rejection of the models are compared to the general model. The asymptotic P value assumes the likelihood-ratio statistic has a standard ² distribution that may not be true if cell counts are too small or if maximization is carried out subject to nonstandard constraints, e.g., constraining p_ij to be between 0 and 1 for the penetrance models. We therefore also present empirical P values for rejecting each specific restricted model, which were calculated by simulating data under the null penetrances for that restricted model, and we note how often the simulated -2 ln likelihood difference exceeded the observed difference. The fitted values for the penetrances for some of these models are displayed graphically in Fig 1, allowing comparisons to be made between the shapes of the graphs for the restricted models and the general model. In the discussion below, we assume a 5% significance level for acceptance/rejection of models; i.e., models with a P value of <0.05 are rejected, although note that we include exact P values where possible.

View larger version (28K): In this window In a new window Download PPT slide   Figure 1. Fitted penetrances under different genetic models as a function of Idd3/Idd10 genotype. Data were generated using the strain NOD.B6Idd10R1.

From Table 4, using the NOD.B6Idd10^R1 strain we find no evidence to reject an additive model for either the log odds or the liability (P values are 0.12 and 0.20, respectively). There is evidence against all three models (additive, heterogeneity, and multiplicative)for the penetrance. In addition there is evidence against a heterogeneity model for the log odds and strong evidence (P = 7 x 10^-16) against a multiplicative model for the log odds. The empirical P values are very close to the asymptotic P values. It is interesting that a multiplicative model for the penetrance is rejected (P = 3 x 10^-5) even though an additive model for the liability is accepted. This contrasts with the results of RISCH et al. 1993 , who found that the additive liability model gave predictions reasonably close to a multiplicative model for the penetrances when applied to backcross data in four genotype categories at a variety of Idd loci. These authors also found that a multiplicative model of epistasis (on the penetrance scale) provided a good fit to the backcross data they analyzed. The differences between our results and those of RISCH et al. 1993 may result partly from consideration of different susceptibility loci, partly from the fact that we have congenic data as opposed to backcross data, and partly from the fact that we are modeling data from all nine genotype categories as opposed to four categories. These last two considerations in particular might be expected to provide a more powerful procedure for discriminating between and rejecting specific genetic models, since congenic and double congenic data allow us to estimate penetrances directly (as opposed to backcross data, which allows estimation of relative penetrances only), while consideration of all possible genotype categories allows us to model differences in effect between a single and a double dose of the predisposing alleles. Note that, in practice, it would not be possible to analyze data at these two loci using the methods of RISCH et al. 1993 with a backcross experiment, since Idd3 and Idd10 are linked and so their segregation is not independent: The methods proposed by RISCH et al. 1993 would not be applicable without some modification to take into account the linkage between the loci. This illustrates one advantage of using a congenic approach.

Using the NOD.B6Idd10^R2 strain (Table 5) we again find no evidence to reject an additive model for either the log odds or the liability. In this strain an additive model for the penetrance is also accepted (P = 0.17), while heterogeneity and multiplicative models for either the penetrance or the log odds are again rejected. The rejection of a heterogeneity model for the penetrance while an additive model is accepted illustrates a difference between these two models when considered on the penetrance scale. The acceptance of the additive model for the penetrance in this strain, which was rejected when using the strain NOD.B6Idd10^R1, may result from the fact that for the strain NOD.B6Idd10^R2 there are two genotype categories missing. From Table 4, these genotype categories contribute to causing the additive model to be rejected when using the NOD.B6Idd10^R1 strain, since the estimated penetrances in the additive model are quite different from the observed penetrances (those estimated in the general model). However, if these categories are dropped in the analysis of the NOD.B6Idd10^R1, we still find evidence against an additive model (P = 0.001). The difference in acceptance/rejection of the additive model is therefore more likely to result from the previously demonstrated difference in diabetes development between the two strains (seen in row 3 of Table 1, 33 vs. 49%) owing to the Idd17 effect.

Table 6 shows the parameter estimates from the general model for strain NOD.B6Idd10^R1 using the quantitative genetics parameterization given in Equation 1. It is not possible to estimate the parameters of this model for strain NOD.B6Idd10^R2 since there are nine parameters to estimate but only seven cells containing data. From Table 6 we find that when modeling the penetrance, there are significant additive x additive (i_aa) and dominance x additive (i_da) interaction effects. When modeling the log odds or liability there are no significant interaction effects, as expected from the fact that the additive model fits well in these cases. It is interesting to note that, regardless of whether the penetrance, log odds, or liability are modeled, the dominance effect d₂ at locus 2 (Idd10) is not significantly different from 0.

The interpretation of the parameter estimates in Table 6 as true additive/dominance/interaction effects relies on the model being orthogonal so that the parameter estimates do not change according to which other parameters are currently being estimated. This is true for a designed experiment such as a backcross and is found to be approximately true for the congenic data analyzed here (results not shown), although it may not hold in general for such data owing to the fact that genotype frequencies may occur in unbalanced proportions and in addition we are analyzing a binary outcome (i.e., a proportion). For nondesigned experiments such as analysis of data from human studies, this approximation is unlikely to hold and careful attention to fitting models in an appropriate hierarchy will be required; e.g., dominance effects may not be included without the relevant additive effects, or interaction effects without the relevant main effects, etc.

Modeling the joint effects of Idd3 and Idd5:
Table 7 and Table 8 show the results of fitting models to the data from Table 2 (using the congenic strains for Idd3 and Idd5). We find that the data are well modeled by either an additive model for the penetrance, an additive model for the log odds, or an additive model for the liability: No significant interaction effects are observed on any of the scales examined. This means that we would not reject the hypothesis that these two loci act additively on the penetrance scale (i.e., with no epistasis), but note that a slightly better fit is provided by an additive model on the liability scale, which on the penetrance scale would include epistasis. The additive model on the log odds scale also fits adequately, which is consistent with results presented by GHOSH et al. 1993 , who fitted logistic models to an ordinal insulitis phenotype in backcross data, and found no evidence for an interaction between Idd3 and Idd5 on the log odds scale. Note that again the multiplicative penetrance model is rejected (P = 0.005), indicating a difference between the multiplicative model and the additive liability model and also indicating that the multiplicative model does not here provide a good fit to the action of Idd3 and Idd5.

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

In this analysis we have fitted models for the joint action of two pairs of NOD diabetes loci, Idd3 and Idd10 and Idd3 and Idd5. The locus Idd3, which we believe to be the interleukin-2 (IL2) locus (DENNY et al. 1997 ), is common to both. Even though we had the NOD.Idd3 Idd10 double congenic data in 1994 (WICKER et al. 1994 ), only in this present article have we analyzed the data in a more formal mathematical way. To our knowledge, only four other groups have evaluated interlocus effects in congenic strains, in models of diabetes (FOX et al. 2000 ), lupus (MOHAN et al. 1999 ), cancer (FIJNEMAN et al. 1996 ), and hypertension (RAPP et al. 1998 ). However, none has carried out the type of modeling we describe here. In all of those reports epistasis has either been inferred (MOHAN et al. 1999 ; FOX et al. 2000 ) or described mathematically (FIJNEMAN et al. 1996 ; RAPP et al. 1998 ) but using simpler functions than those described here. Our consideration of the joint effects of Idd loci within the NOD mouse model further validates this strain as a model of human type 1 diabetes. Methods for analysis of interlocus effects in human complex disease are still being developed. Although the genes predisposing to the disease in humans may not be identical to those in the NOD mouse, the hope is that the underlying genetic basis in terms of number of genes involved and their influence on physiological disease processes may present similarities between the two species.

Our results suggest that the joint effects of Idd3 and Idd10 follow a model that is additive on either the log odds or liability scale, but epistatic on the penetrance scale. The action of these two loci does not appear to be multiplicative, in spite of previous suggestions that most loci involved in type 1 diabetes will contribute to disease in a multiplicative manner (RISCH et al. 1993 ). Moreover, a multiplicative model is also rejected for the action of Idd3 and Idd5. Although multiplicative genetic models were strongly rejected in all our analyses, it was not always possible to distinguish between alternative models that fit the data; e.g., for Idd3 and Idd5 the joint action could be well modeled by either an additive model for the penetrance, an additive model for the log odds, or an additive model for the liability.

Having accepted and rejected specific models for the joint action of two loci, the question arises as to the interpretation (biological or otherwise) of these results. The biologist is interested primarily in mechanisms and pathways, but the detection of a statistical interaction does not necessarily imply interaction on the biological or mechanistic level (WITTE 1998 ). Moreover, since the presence of a statistical interaction depends on the scale of measurement (i.e., whether we choose to model the penetrance, the log odds, or the liability), it is unclear what the biological interpretation should be.

Although this issue of interpretation has not previously been discussed in detail in the genetics literature, it has historically received much attention in the epidemiological literature (see GREENLAND and ROTHMAN 1998 ). ROTHMAN et al. 1980 define four types of interaction: statistical, biological, public health, and individual decision making, and note that the scale of interest may be different for each type. A statistical interaction, for instance, refers to departure from additivity of effects on any chosen outcome scale and depends on the statistical model postulated. It can be argued that the scale itself has little biological meaning: The preferred statistical procedure would be to choose a scale that removes interaction to generate the most parsimonious model (ELSTON 1961 ; COX 1984 ; BRESLOW and STORER 1985 ). This procedure makes sense when the goal is simple description of observed phenomena or prediction of outcome (ROTHMAN et al. 1980 ). For public health purposes, an interaction between two factors relates to the public health consequences of joint exposure to the factors, which may be considered to be directly proportional to the resulting number of cases in the population. It is usually argued that, for public health purposes, the natural scale is the incidence rate, with interaction between risk factors equivalent to a departure from additivity of incidence rate differences (BLOTT and DAY 1979 ; ROTHMAN et al. 1980 ). Similar economic considerations usually underlie the choice of scale in quantitative genetics, where traits such as total yield, milk production, etc. are of interest. In individual (personal) decision making, ROTHMAN et al. 1980 argue that similar reasoning leads to the absolute risk of disease being the relevant outcome scale.

The phenomenon of most interest in the present context is biological interaction and the degree to which it can be elucidated by statistical analysis. Unfortunately, this turns out to be the most complex of the interactions considered. The problem is that any given data pattern can usually be obtained from a number of dissimilar mechanisms or models for disease development (SIEMIATYCKI and THOMAS 1981 ; THOMPSON 1991 ). For instance, consider the following five models of disease causation: model 1, a "no-hit" model in which those who become diseased are those who fail to experience one or more occurrences of a beneficial event, and two factors act additively in increasing the rate at which beneficial events occur; model 2, a "single-hit" model in which occurrence of a single adverse event is sufficient for development of disease, and two factors show more-than-additive effects (i.e., one factor augments the biologic effect of another) in terms of increasing the rate at which events occur [such effects are sometimes called "synergistic" although BLOTT and DAY 1979 propose that this term be reserved for the public health concept of interaction]; model 3, each factor affects a different stage of a multistage process in which the pathogenic process involves transition from a normal state to stage 1, followed by subsequent transition from stage 1 to disease; model 4, two factors have additive effects on an intervening variable that bears an exponential relationship to incidence of disease; model 5, two factors have less-than-additive ("antagonistic") effects on an intervening variable that bears an exponential power relationship to incidence of disease. THOMPSON 1991 shows that each of these very different causal models results in a multiplicative statistical model for the effects of the two risk factors. SIEMIATYCKI and THOMAS 1981 describe additional models in which multiplicative and additive statistical models may arise from either interacting or noninteracting biological models.

The problem of interpretability is compounded by the low statistical power for detecting statistical interactions even when they are present and by the "discretizing" of an underlying continuous variable that can influence the presence of interaction (ROTHMAN and KELLER 1972 ). The fact that we are dealing with multifactorial discrete traits with reduced penetrance also adds to the complexity. For instance, if the penetrances p_ij are constrained to be either 0 or 1, the heterogeneity model on the penetrance scale as described earlier has a natural biological interpretation, but for a complex disease this constraint would seem unrealistic. GREENLAND and ROTHMAN 1998 show that for a discrete trait, departures from additivity in average risk correspond to biological interaction in a general sense (via many possible mechanisms), although no departure from additivity does not necessarily imply absence of biological interaction. These results in terms of absolute risk would support the idea of using the penetrance as the outcome of interest in our analyses; however, GREENLAND and ROTHMAN 1998 make assumptions that are essentially equivalent to assuming complete penetrance within each genotype category, which is again unrealistic for complex disease (and moreover patently contradicted by our data). Even if the underlying model in terms of all genetic and environmental effects could be expressed in these terms, analysis of any subset of these factors would integrate out the effects of the other factors, leading to reduced penetrance and a violation of the assumptions of GREENLAND and ROTHMAN 1998 .

The overall conclusion from this is that from the numerical data alone it will usually be impossible to discern how or even whether two risk factors interact in any biologically meaningful way. The assumed biological interpretation of models given in RISCH 1990 is likely to be unrealistically simplistic. Nevertheless, there may be some value in modeling epistasis. If a prior biological model can be postulated in detail, it may be possible to infer a statistical model and determine how well the data fit the model. For instance, with models of multistage carcinogenesis, specific scenarios for the action of etiological factors and in particular whether they act in the same step or at different steps in the process can be shown to correspond to particular mathematical models (ROTHMAN et al. 1980 ; SIEMIATYCKI and THOMAS 1981 ). In reality, however, it may be that many statistical models fit the data equally well. In this case, allowing for different modes of interaction between potential disease loci can lead to increased power to detect any one of the loci. Simulation studies (CORDELL et al. 1995 , CORDELL et al. 2000 ; LEAL and OTT 2000 ) suggest that this increase in power may be relatively modest. Nevertheless, in analysis of real data for type 1 diabetes (CORDELL et al. 1995 , CORDELL et al. 2000 ), type 2 diabetes (COX et al. 1999 ), and inflammatory bowel disease (CHO et al. 1998 ), increased evidence for linkage at one locus was seen when the interaction at another locus was taken into consideration. Finally, identification of the most parsimonious statistical model for the joint effects of alleles at several loci provides a means for improved prediction of phenotype, compared to considering the loci in isolation. Understanding the joint action of the loci may also lead to improved targeting of interventions; e.g., for the loci considered here, the fact that the disease frequency in the double congenic strain is so much lower than that in either of the single congenic strains suggests that intervention via relatively few gene products may be of value. The question of true biological interaction, for instance, between alleles at the HLA and insulin genes in human type 1 diabetes, remains of paramount interest. However, further research concerning the interpretation and significance of results will be required if statistical analyses such as those presented here are to be used to elucidate the underlying biologic mechanisms involved in complex binary traits.

	FOOTNOTES

¹ Present address: Department of Immunology, The Scripps Research Institute, 10550 N. Torrey Pines Rd., La Jolla, CA 92037.
² Present address: The Breakthrough Toby Robins Breast Cancer Research Centre, Institute of Cancer Research, Mary-Jean Mitchell Green Bldg., Chester Beatty Laboratories, Fulham Rd., London SW3 6JB, United Kingdom.

	ACKNOWLEDGMENTS

We thank colleagues at Cambridge and Case Western Reserve Universities for fruitful discussions. This work was supported by grants from the Wellcome Trust, the U.K. Medical Research Council, the Juvenile Diabetes Fund, and Diabetes U.K. Some of the results of this article were obtained by using the program package S.A.G.E., which is supported by a U.S. Public Health Service resource grant (1 P 41 RR03655) from the National Center for Research Resources.

Manuscript received December 18, 2000; Accepted for publication February 5, 2001.

	LITERATURE CITED

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

BATESON, W., 1909 Mendel's Principles of Heredity. Cambridge University Press, Cambridge/London/New York.

BLOTT, N. J. and N. E. DAY, 1979  Synergism and interaction: are they equivalent? Am. J. Epidemiol. 110:99-100[Free Full Text].

BRESLOW, N. E. and B. E. STORER, 1985  General relative risk functions for case-control studies. Am. J. Epidemiol. 122:149-162[Abstract/Free Full Text].

CHO, J. H., D. L. NICOLAE, L. H. GOLD, C. T. FIELDS, and M. C. LABUDA et al., 1998  Identification of novel susceptibility loci for inflammatory bowel disease on chromosomes 1p, 3q, and 4q: evidence for epistasis between 1p and IBD1. Proc. Natl. Acad. Sci. USA 95:7502-7507[Abstract/Free Full Text].

COCKERHAM, C. C. and COCKERHAM, C. C.Z-B. ZENG, 1996  Design III with marker loci. Genetics 143:1437-1456[Abstract].

CORDELL, H. J., J. A. TODD, S. T. BENNETT, Y. KAWAGUSHI, and M. FARRALL, 1995  Two-locus maximum lod score analysis of a multifactorial trait: joint consideration of IDDM2 and IDDM4 with IDDM1 in type 1 diabetes. Am. J. Hum. Genet. 57:920-934[Medline].

CORDELL, H. J., G. C. WEDIG, K. B. JACOBS, and R. C. ELSTON, 2000  Multilocus linkage tests based on affected relative pairs. Am. J. Hum. Genet. 66:1273-1286[Medline].

COX, D. R., 1984  Interaction. Int. Stat. Rev. 52:1-31.

COX, N. J., M. FRIGGE, D. L. NICOLAE, P. CONCANNON, and C. L. HANIS et al., 1999  Loci on chromosomes 2 (NIDDM1) and 15 interact to increase susceptibility to diabetes in Mexican Americans. Nat. Genet. 21:213-215[Medline].

DENNY, P., C. J. LORD, N. J. HILL, J. V. GOY, and E. R. LEVY et al., 1997  Mapping of the insulin-dependent diabetes locus, Idd3, to a 0.35 cM interval containing the Interleukin-2 gene. Diabetes 46:695-700[Abstract].

ELSTON, R. C., 1961  On additivity in the analysis of variance. Biometrics 17:209-219.

FIJNEMAN, R. J., S. S. DE VRIES, R. C. JANSEN, and P. DEMANT, 1996  Complex interactions of new quantitative trait loci, Sluc1, Sluc2, Sluc3, and Sluc4, that influence the susceptibility to lung cancer in the mouse. Nat. Genet. 14:465-467[Medline].

FISHER, R. A., 1918  The correlation between relatives on the supposition of Mendelian inheritance. Trans. R. Soc. Edinb. 52:399-433.

FISHER, R. A., 1930 The Genetical Theory of Natural Selection. Oxford University Press, London/New York/Oxford.

FOX, C. J., A. D. PATERSON, S. M. MORTIN-TOTH, and J. S. DANSKA, 2000  Two genetic loci regulate T cell dependent islet inflammation and drive autoimmune diabetes pathogenesis. Am. J. Hum. Genet. 67:67-81[Medline].

FRANKEL, W. N. and N. J. SCHORK, 1996  Who's afraid of epistasis? Nat. Genet. 14:371-373[Medline].

GHOSH, S., S. M. PALMER, N. R. RODRIGUES, H. J. CORDELL, and C. M. HEARNE et al., 1993  Polygenic control of autoimmune diabetes in nonobese diabetic mice. Nat. Genet. 4:404-409[Medline].

GREENLAND, S., and K. J. ROTHMAN (Editors), 1998 Concepts of interaction, pp. 329–342 in Modern Epidemiology, Ed. 2. Lippincott-Raven Publishers, Philadelphia.

GUERRERO, V. M. and R. A. JOHNSON, 1982  Use of Box-Cox transformation with binary response models. Biometrika 69:309-314[Abstract/Free Full Text].

HILL, N. J., P. A. LYONS, N. ARMITAGE, J. A. TODD, and L. S. WICKER et al., 2000  NOD Idd5 locus controls insulitis and diabetes and overlaps the orthologous CTLA4/IDDM12and NRAMP1 loci in humans. Diabetes 49:1744-1747[Abstract].

HODGE, S., 1981  Some epistatic two-locus models of disease. I. Relative risks and identity by descent distributions in affected sib pairs. Am. J. Hum. Genet. 33:382-395.

LEAL, S. M. and J. OTT, 2000  Effects of stratification in the analysis of affected-sib-pair data: benefits and costs. Am. J. Hum. Genet. 66:567-575[Medline].

LYONS, P. A., and L. S. WICKER, 1999 Localising quantitative trait loci in the NOD mouse model of type 1 diabetes, pp. 208–225 in Genes and Genetics of Autoimmunity, Vol. 1, edited by A. N. THEOFILOPOULOS. Karger, Basel.

MCCULLAGH, P., and J. A. NELDER, 1989 Generalized Linear Models, Ed. 2. Chapman & Hall, London.

MAKINO, S., K. KUNIMOTO, Y. MURAOKA, Y. MIZUSHIMA, and K. KATAGIRI et al., 1980  Breeding of a non-obese, diabetic strain of mice. Exp. Anim. 29:1-13.

MOHAN, C., L. MOREL, P. YANG, H. WATANABE, and B. CROKER et al., 1999  Genetic dissection of lupus pathogenesis: a recipe for nephrophilic autoantibodies. J. Clin. Invest. 103:1685-1695[Medline].

NELDER, J. A. and R. W. M. WEDDERBURN, 1972  Generalized linear models. J. R. Stat. Soc. Ser. A 135:370-384.

NEUMAN, R. J. and J. P. RICE, 1992  Two-locus models of disease. Genet. Epidemiol. 9:347-365[Medline].

NORTON, B. and E. S. PEARSON, 1976  A note on the background to and refereeing of R. A. Fisher's 1918 paper `The correlation between relatives on the supposition of Mendelian inheritance'. Notes Rec. R. Soc. Lond. 31:151-162[Medline].

PEARSON, K., 1900  Mathematical contributions to the theory of evolution VIII: on the inheritance of characters not capable of exact quantitative measurement. Philos. Trans. A 195:79-150.

PODOLIN, P. L., P. DENNY, C. J. LORD, N. J. HILL, and J. A. TODD et al., 1997  Congenic mapping of the insulin dependent diabetes (Idd) gene, Idd10, localises two genes mediating the Idd10 effect and eliminates the candidate Fcgr1. J. Immunol. 159:1835-1843[Abstract].

RAPP, J. P., M. R. GARRETT, and A. Y. DENG, 1998  Construction of a double congenic strain to prove an epistatic interaction on blood pressure between rat chromosomes 2 and 10. J. Clin. Invest. 101:1591-1595[Medline].

RISCH, N., 1990  Linkage strategies for genetically complex traits. I. Multilocus models. Am. J. Hum. Genet. 46:222-228[Medline].

RISCH, N., S. GHOSH, and J. A. TODD, 1993  Statistical evaluation of multiple locus linkage data in experimental species and relevance to human studies: application to murine and human IDDM. Am. J. Hum. Genet. 53:702-714[Medline].

ROTHMAN, K. and A. KELLER, 1972  The effect of joint exposure to alcohol and tobacco on risk of cancer of the mouth and pharynx. J. Chron. Dis. 25:711-716[Medline].

ROTHMAN, K. J., S. GREENLAND, and A. M. WALKER, 1980  Concepts of interaction. Am. J. Epidemiol. 112:467-470[Free Full Text].

SCHORK, N. J., M. BOEHNKE, J. D. TERWILLIGER, and J. OTT, 1993  Two-trait-locus linkage analysis: a powerful strategy for mapping complex genetic traits. Am. J. Hum. Genet. 53:1127-1136[Medline].

SIEMIATYCKI, J. and D. C. THOMAS, 1981  Biological models and statistical interactions: an example from multistage carcinogenesis. Int. J. Epidemiol. 10:383-387[Abstract/Free Full Text].

THOMAS, D. C., 1981  General relative risk models for survival time and matched case-control studies. Biometrics 37:673-686.

THOMPSON, W. D., 1991  Effect modification and the limits of biological inference from epidemiologic data. J. Clin. Epidemiol. 44:221-232[Medline].

WEBER, K., R. EISMAN, S. HIGGINS, L. KUHL, and A. PATTY et al., 1999  An analysis of polygenes affecting wing shape on chromosome three in Drosophila melanogaster. Genetics 153:773-786[Abstract/Free Full Text].

WICKER, L. S., J. A. TODD, J.-B. PRINS, P. L. PODOLIN, and R. J. RENJILIAN et al., 1994  Resistance alleles at two non-MHC-linked insulin dependent diabetes loci on chromosome 3, Idd3 and Idd10, protect NOD mice from diabetes. J. Exp. Med. 180:1705-1713[Abstract/Free Full Text].

WITTE, J. S., 1998 Gene-environment interaction, pp. 1613–1614 in Encyclopedia of Biostatistics, edited by P. ARMITAGE and T. COLTON. John Wiley & Sons, Chichester, UK.

WRIGHT, S., 1934a  An analysis of variability in number of digits in an inbred strain of guinea pigs. Genetics 19:506-536[Free Full Text].

WRIGHT, S., 1934b  The results of crosses between inbred strains of guinea pigs, differing in number of digits. Genetics 19:537-551[Free Full Text].

YUI, M. A., K. MURALIDHARAN, B. MORENO-ALTAMIRANO, G. PERRIN, and K. CHESTNUT et al., 1996  Production of congenic mouse strains carrying NOD-derived diabetogenic genetic intervals: an approach for the genetic dissection of complex traits. Mamm. Genome 7:331-334[Medline].

ZENG, Z-B., J. LIU, L. F. STAM, C-H. KAO, and J. M. MERCER et al., 2000  Genetic architecture of a morphological shape difference between two Drosophila species. Genetics 154:299-310[Abstract/Free Full Text]