Classical Models for Twin Data: The Case of Categorical Data

ABSTRACT

The classical models ACE and ADE were used in the 1990’s to estimate heredity of a phenotype from data on monozygotic and dizygotic twins. The author extended these models to a model called ACDE with four parameters instead of only three. In that paper, the data were assumed to be continuous. This paper considers the same models in the case where the data is categorical. It is showed how these models can be estimated by maximum likelihood. An example is given based on twin data on BMI from the UK. This is the same data as in the previous paper but in categorized form.

KEYWORDS:

• Categorical data
• twin data
• twin models
• heredity

Twin data models

Models for twin data or other types of family relationships were common in the 1980’s and 1990’s see e.g., Neale and Cardon (1992). One has data on monozygotic (MZ) twins and dizygotic (DZ) twins. It is assumed that MZ twins have all their genes in common and dz twins have half of their genes in common. We observe a phenotype on each twin. A phenotype may be anything that one can measure or observe on each twin, such as symptoms (e.g., allergy), illnesses (e.g.,cancer, diabetes, asthma), physical measures (e.g., weight, height), personality traits (e.g.,nervousness), longevity (how long you live) etc.

Jöreskog (2020) covered the case when the phenotypes were continuous and showed how the classical models ACE, ADE, and ACDE could be estimated. This paper covers the case when the phenotypes are categorical. This requires a different methodology.

Data

Thera are $T_{mz}$ monozygotic twins and $T_{dz}$ dizygotic twins. Each twin in a twin pair responds on a set of $c$ ordered categories. Typically, $c=2,3,4,5$. So the data is represented by a two-way contingency table, one for the MZ group and one for the DZ group:

$\begin{array}{rl}& N_{mz}=\left[\begin{array}{cccc}{n}_{11}^{(mz)}& n_{12}^{(mz)}& \cdots & n_{1c}^{(mz)}\\ n_{21}^{(mz)}& n_{22}^{(mz)}& \cdots & n_{2c}^{(mz)}\\ ⋮& ⋮& ⋮⋮& ⋮\\ n_{c1}^{(mz)}& n_{c2}^{(mz)}& \cdots & n_{cc}^{(mz)}\end{array}\right],\\ & N_{dz}=\left[\begin{array}{cccc}{n}_{11}^{(dz)}& n_{12}^{(dz)}& \cdots & n_{1c}^{(dz)}\\ n_{21}^{(dz)}& n_{22}^{(dz)}& \cdots & n_{2c}^{(dz)}\\ ⋮& ⋮& ⋮⋮& ⋮\\ n_{c1}^{(dz)}& n_{c2}^{(dz)}& \cdots & n_{cc}^{(dz)}\end{array}\right]\end{array}$

where $n_{ij}$ is the number of twins in the sample where twin 1 responds on category $i$ and twin 2 responds on category j. We have

$T_{mz}=\sum _{i=1}^c\sum _{j=1}^c{n}_{ij}^{(mz)},T_{dz}=\sum _{i=1}^c\sum _{j=1}^c{n}_{ij}^{(dz)}$

Model

Let ${\pi }_{ij}$ be the probability of a response in categories $i$ and $j$ as described in the previous section. The model is

(1) $\begin{array}{rl}& {\pi }_{ij}^{(mz)}=\underset{{\tau }_{i-1}}{\overset{{\tau }_i}{\int }}\underset{{\tau }_{j-1}}{\overset{{\tau }_j}{\int }}{\varphi }_2(u,v,{\rho }_{mz})dudv\\ & {\pi }_{ij}^{(dz)}=\underset{{\tau }_{i-1}}{\overset{{\tau }_i}{\int }}\underset{{\tau }_{j-1}}{\overset{{\tau }_j}{\int }}{\varphi }_2(u,v,{\rho }_{dz})dudv,\end{array}$(1)

where

${\varphi }_2(u,v,\rho )=\frac{1}{2\pi \sqrt{1-{\rho }^2}}{e}^{-\frac{1}{2(1-{\rho }^2)}(u^2-2\rho uv+v^2)},$

is the standard bivariate normal density with correlation $\rho$. The thresholds

${\tau }_0=-\mathrm{\infty }<{\tau }_1<{\tau }_2<\dots <{\tau }_{c-1}<{\tau }_c=\mathrm{\infty }$

are assumed to be the same within and between groups. If there are $c$ categories, there are $c-1$ strictly inceasing thresholds. This choice of thresholds ensures that the only difference between the MZ and DZ group is manifested through the difference between ${\rho }_{mz}$ and ${\rho }_{dz}$.

Note that ${\pi }_{ij}={\pi }_{ji},$ for both MZ and DZ. Since the marginal density functions of ${\varphi }_2(u,v;\rho )$ is the standard normal density function independent of $\rho$, it follows from (1) that the univariate row and column marginals

$\sum _{j=1}^c{\pi }_{ij}=\sum _{j=1}^c{\pi }_{ji},i=1,2,\cdots ,c,$

are equal and the same for both MZ and DZ. Denoting this common marginal as

$\alpha =({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_c),$

we have ${\alpha }_i>0$ and ${\sum }_{i=1}^c{\alpha }_i=1$ Then

${\tau }_i={\mathrm{\Phi }}^{-1}({\alpha }_1+{\alpha }_2+\cdots +{\alpha }_i),i=1,2,\dots ,c-1.$

Reversely,

${\alpha }_i=\mathrm{\Phi }({\tau }_i)-\mathrm{\Phi }({\tau }_{i-1}),i=1,2,\dots ,c.$

$\mathrm{\Phi }$ is the standard normal distribution function and ${\mathrm{\Phi }}^{-1}$ is its inverse function. Note that the ${\alpha }_i$, and hence, the ${\tau }_i$ do not depend on ${\rho }_{mz}$ or ${\rho }_{dz}$.

Here we use the same notation and definitions of $h^2$, $c^2$, $d^2$, and $e^2$ as in the previous paper (Jöreskog, 2020). Since there is no scale of the categorical measurement, it is assumed that the total variance $h^2+c^2+d^2+e^2=1$. So the only equations we have are

(2) ${\rho }_{mz}=h^2+c^2+d^2$(2)

(3) ${\rho }_{dz}=\frac{1}{2}{h}^2+c^2+\frac{1}{4}{d}^2$(3)

This is two equations in the three unknowns $h^2,c^2,d^2$. There is no unique solution for $h^2,c^2,d^2$. However, if $d^2=0$, we can solve for $h^2$ and $c^2$. This gives

$h^2=2({\rho }_{mz}-{\rho }_{dz}),c^2=2{\rho }_{dz}-{\rho }_{mz}$

corresponding to the ACE model. Similarly, if $c^2=0$, we can solve for $h^2$ and $d^2$. This gives

$h^2=4{\rho }_{dz}-{\rho }_{mz},d^2=2{\rho }_{mz}-4{\rho }_{dz}$

corresponding to the ADE model.

As argued in the previous paper, neither $c^2=0$ nor $d^2=0$ are reasonable assumptions. Therefore, an alternative model is considered as follows

Equations (2)(2) ${\rho }_{mz}=h^2+c^2+d^2$(2) and (3(3) ${\rho }_{dz}=\frac{1}{2}{h}^2+c^2+\frac{1}{4}{d}^2$(3) ) can be written in matrix form as

(4) $\left(\begin{array}{c}{\rho }_{mz}\\ {\rho }_{dz}\end{array}\right)=\left(\begin{array}{ccc}1& 1& 1\\ \frac{1}{2}& 1& \frac{1}{4}\end{array}\right)\left(\begin{array}{c}{h}^2\\ c^2\\ d^2\end{array}\right)$(4)

Let A be the matrix in (4) and let $\mathbf{B}$ be

$\mathbf{B}={\mathbf{A}}^{\prime }(\mathbf{A}{\mathbf{A}}^{\prime }{)}^{-1}.$

Using paper and pencil algebra, gives

$\mathbf{B}=\left(\begin{array}{cc}\frac{1}{2}& -\frac{2}{7}\\ -\frac{1}{2}& 1\frac{3}{7}\\ 1& -1\frac{1}{7}\end{array}\right).$

The ACDE model is defined by

$\left(\begin{array}{c}{h}^2\\ c^2\\ d^2\end{array}\right)=\mathbf{B}\left(\begin{array}{c}{\rho }_{mz}\\ {\rho }_{dz}\end{array}\right)$

In particalar,

$h^2={\rho }_{mz}-\frac{2}{7}{\rho }_{dz}$

Estimation

Assuming multinomial distributions of the $n_{ij}$ and independent groups, the log likelihood function is

$lnL=\sum _{i=1}^c\sum _{j=1}^c{n}_{ij}^{(mz)}ln{\pi }_{ij}^{(mz)}+\sum _{i=1}^c\sum _{j=1}^c{n}_{ij}^{(dz)}ln{\pi }_{ij}^{(dz)}.$

In the following the shorter notation ${\sum }_{ij}$ is used instead of ${\sum }_{i=1}^c{\sum }_{j=1}^c$. Let

$\tau =({\tau }_1,{\tau }_2,\dots ,{\tau }_{c-1})$

Assuming that the samples are sufficiently large so that all $n_{ij}>0$, it is convenient to minimize the fit function

$F(\theta )=F_{mz}(\tau ,{\rho }_{mz})+F_{dz}(\tau ,{\rho }_{dz})$

where

$F_{mz}(\tau ,{\rho }_{mz})=\sum _{ij}{n}_{ij}^{(mz)}ln{p}_{ij}^{(mz)}-\sum _{ij}{n}_{ij}^{(mz)}ln{\pi }_{ij}^{(mz)}$

and

$F_{dz}(\tau ,{\rho }_{dz})=\sum _{ij}{n}_{ij}^{(dz)}ln{p}_{ij}^{(dz)}-\sum _{ij}{n}_{ij}^{(dz)}ln{\pi }_{ij}^{(dz)}$

and where

$p_{ij}^{(mz)}=n_{ij}^{(mz)}/T_{mz},p_{ij}^{(dz)}=n_{ij}^{(dz)}/T_{dz}.$

The parameter vector $\theta$ is

$\theta =({\tau }_1,{\tau }_2,\dots ,{\tau }_{c-1},{\rho }_{mz},{\rho }_{dz})$

The function $F(\theta )$ which is non-negative, is to be minimized with respect to $\theta$. $F(\theta )$ can be minimized numerically using first and second derivatives. At the minimum, the inverse of the expected Hessian matrix can be used as an estimate of the asymptotic covariance matrix of $\hat{\theta }$, and

$C=2F(\hat{\theta })$

can be used as a chi-square statistic with $2(c^2-1)-(c+1)$ degrees of freedom to test the fit of the model. This is an LR (likelihood ratio) statistic. Alternatively, one can use the GF (goodness-of-fit) statistic

$G=\sum _{ij}\frac{{(n_{ij}^{(mz)}-T_{mz}{\hat{\pi }}_{ij}^{(mz)})}^2}{T_{mz}{\hat{\pi }}_{ij}^{(mz)}}+\sum _{ij}\frac{{(n_{ij}^{(dz)}-T_{dz}{\hat{\pi }}_{ij}^{(dz)})}^2}{T_{dz}{\hat{\pi }}_{ij}^{(dz)}},$

which is also approximately distributed as ${\chi }^2$ with the same degrees of freedom.

Example

In this example we use the same BMI data from the UK as in the previous paper, but the measurement from each twin is classified into three categories: normal, overweight, and obese on the basis of their BMI. So $c=3$, $T_{mz}=794$ and $T_{dz}=758$. The contingency tables are

$N_{mz}=\left(\begin{array}{ccc}325& 78& 8\\ 60& 150& 38\\ 7& 40& 88\end{array}\right)N_{dz}=\left(\begin{array}{ccc}202& 96& 46\\ 97& 112& 56\\ 19& 64& 66\end{array}\right).$

In terms of proportions these are

$\begin{array}{rl}& P_{mz}=\left(\begin{array}{ccc}0.4093& 0.0982& 0.0101\\ 0.0756& 0.1889& 0.0479\\ 0.0088& 0.0504& 0.1108\end{array}\right)\\ & P_{dz}=\left(\begin{array}{ccc}0.2665& 0.1266& 0.0607\\ 0.1280& 0.1478& 0.0739\\ 0.0251& 0.0844& 0.0871\end{array}\right).\end{array}$

Because of the inequality ${\tau }_1<{\tau }_2$ and the restricted range of ${\rho }_{mz}$ and ${\rho }_{dz}$ it is neccessary to have good starting values for the numerical minimization of $F(\theta )$. These are obtained as follows.

Take the average of the row and column margins of $P_{mz}$ and $P_{dz}$. Then form a wheighted average of these averages with weights $w_{mz}=T_{mz}/(T_{mz}+T_{dz})$ and $w_{dz}=T_{dz}/(T_{mz}+T_{dz})$. This gives $\alpha$ as

$\alpha =(0.472,0.339,0.189)$

and estimated ${\tau }_1=-0.070$ and ${\tau }_2=0.882$.

To obtain starting values of ${\rho }_{mz}$ and ${\rho }_{dz}$ minimize $F(\theta )$ with respect to ${\rho }_{mz}$ and ${\rho }_{dz}$ for given ${\tau }_1$ and ${\tau }_2$ separately for MZ and DZ. This is a special case of the polychoric correlation, see Olsson (1979), which gives ${\rho }_{mz}=0.810$ and ${\rho }_{dz}=0.456$. So the starting values for $\theta$ are

${\theta }_0=(-0.070,0.882,0.810,0.456)$

The maximum likelihood estimates are

$\hat{\theta }=(-0.068,0.872,0.811,0.457),$

very close to the starting values. The value of the LR statistic $C$ is 31.50 with 12 degrees of freedom. This has a $P$-value of 0.0017 indicating a not so good fit of the model. The value of the GF statistic $G$ is 31.81, very similar to $C$.

To examine the fit in detail, one can consider the GF residuals

$\frac{n_{ij}-T{\hat{\pi }}_{ij}}{\sqrt{T{\hat{\pi }}_{ij}}}$

separately for MZ and DZ. These are

$\left(\begin{array}{ccc}1.81& 0.49& 1.34\\ -1.60& 0.24& -1.49\\ 0.95& -1.20& -1.09\end{array}\right)$

for MZ and

$\left(\begin{array}{ccc}-1.47& -0.48& 2.70\\ -0.38& 1.15& 0.04\\ -2.15& 1.11& 0.97\end{array}\right)$

for DZ. It is seen that there is only one residual larger than 2 in absolute value. This suggests there are too many twins in the DZ group where twin 1 is in category 1 (normal weight) and twin 2 is in category 3 (obese). There are 46 such twins. According to the model the expected number is 31.

Since the total sample sise $T=1552$ is very large compared to the number of parameters, it is reasonable to consider RMSEA (Root Mean Square Error of Approximation) as a measure of fit, see Browne and Cudeck (1993). This is $R=\sqrt{(C-d)/T}=0.03$ indicating that the model fits at least approximately.

The estimated asymptotic covariance matrix is

$TAcov(\hat{\theta })=\left(\begin{array}{cccc}1.2104& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ 0.7182& 1.4200& \hspace{0.17em}& \hspace{0.17em}\\ 0.0719& -0.1734& 0.5459& \hspace{0.17em}\\ -0.1259& -0.4722& 0.0789& 2.4473\end{array}\right)$

Dividing these numbers by 1552 and taking the square root of the diagonal elements one can obtain the standard errors of the parameter estimates $\hat{\theta }$.

Let

${\mathbf{B}}_{\mathrm{A}\mathrm{C}\mathrm{E}}=\left(\begin{array}{cc}2& -2\\ -1& 2\\ 0& 0\end{array}\right){\mathbf{B}}_{\mathrm{A}\mathrm{D}\mathrm{E}}=\left(\begin{array}{cc}-1& 4\\ 0& 0\\ 2& -4\end{array}\right){\mathbf{B}}_{\mathrm{A}\mathrm{C}\mathrm{D}\mathrm{E}}=\left(\begin{array}{cc}\frac{1}{2}& -\frac{2}{7}\\ -\frac{1}{2}& 1\frac{13}{7}\\ 1& -1\frac{1}{7}\end{array}\right)$

Estimates of $h^2,c^2,d^2$ are obtained from

$\left(\begin{array}{c}{\hat{h}}^2\\ {\hat{c}}^2\\ {\hat{d}}^2\end{array}\right)=\mathbf{B}\left(\begin{array}{c}{\hat{\rho }}_{mz}\\ {\hat{\rho }}_{dz}\end{array}\right)$

and estimates of their asymptotic covariance matrix are obtained from

$TAcov\left(\begin{array}{c}{\hat{h}}^2\\ {\hat{c}}^2\\ {\hat{d}}^2\end{array}\right)=\mathbf{B}\left(\begin{array}{cc}0.5459& \hspace{0.17em}\\ 0.0789& 1.5644\end{array}\right){\mathbf{B}}^{\prime }$

The results are summarized in Table 1.

Table 1. Parameter estimates (0* indicates a fixed value) and standard errors

Parameter	ACE	ADE	ACDE
${\hat{h}}^2$	0.709(0.084)	1.017(0.157)	0.275(0.014)
${\hat{c}}^2$	0.103(0.079)	0*	0.117(0.045)
${\hat{d}}^2$	0*	−0.206(0.159)	0.289(0.047)

The value of ${\hat{h}}^2$ in Table 1 compares well the values of $\hat{H}$ in Table 2 of Jöreskog (2020). The ADE model gives non-admissible results in both papers. The values of ${\hat{h}}^2$ and $\hat{H}$ for the ACE and ACDE models are fairly close. The standard errors are smaller for $\hat{H}$. Intuitively, this is expected since the continuous data has more information than the categorical data. As in the previous paper, the ACDE model gives a more reasonable value of heredity of 28% than 71% as given by the ACE model.