Sum Scores in Twin Growth Curve Models: Practicality Versus Bias

Abstract

To study behavioral or psychiatric phenotypes, multiple indices of the behavior or disorder are often collected that are thought to best reflect the phenotype. Combining these items into a single score (e.g. a sum score) is a simple and practical approach for modeling such data, but this simplicity can come at a cost in longitudinal studies, where the relevance of individual items often changes as a function of age. Such changes violate the assumptions of longitudinal measurement invariance (MI), and this violation has the potential to obfuscate the interpretation of the results of latent growth models fit to sum scores. The objectives of this study are (1) to investigate the extent to which violations of longitudinal MI lead to bias in parameter estimates of the average growth curve trajectory, and (2) whether absence of MI affects estimates of the heritability of these growth curve parameters. To this end, we analytically derive the bias in the estimated means and variances of the latent growth factors fit to sum scores when the assumption of longitudinal MI is violated. This bias is further quantified via Monte Carlo simulation, and is illustrated in an empirical analysis of aggression in children aged 3–12 years. These analyses show that measurement non-invariance across age can indeed bias growth curve mean and variance estimates, and our quantification of this bias permits researchers to weigh the costs of using a simple sum score in longitudinal studies. Simulation results indicate that the genetic variance decomposition of growth factors is, however, not biased due to measurement non-invariance across age, provided the phenotype is measurement invariant across birth-order and zygosity in twins.

Appendices

Appendix 1

Under the assumption that the SLGM is the true data-generating model, summing the individual items at each time point results in a score that is a function of underlying factor loadings as well as the proposed growth in the factor. As an example, consider sum scores of 4 time points. A more explicit formulation of Eq. (9) is then,

$$\left[ {\begin{array}{*{20}{c}} {\mathop \sum \limits_{j=1}^p E[{y_{ij1}}]} \\ {\mathop \sum \limits_{j=1}^p E[{y_{ij2}}]} \\ {\mathop \sum \limits_{j=1}^p E[{y_{ij3}}]} \\ {\mathop \sum \limits_{j=1}^p E[{y_{ij4}}]} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {\mathop \sum \limits_{j=1}^{\text{p}} {\lambda _{j1}}}&0&0&0 \\ 0&{\mathop \sum \limits_{j=1}^{\text{p}} {\lambda _{j2}}}&0&0 \\ 0&0&{\mathop \sum \limits_{j=1}^{\text{p}} {\lambda _{j3}}}&0 \\ 0&0&0&{\mathop \sum \limits_{j=1}^{\text{p}} {\lambda _{j4}}} \end{array}} \right]{\mathbf{\Gamma \upmu }}$$

(38)

In practice, fitting a univariate LGM to the sum score results in,

$${\mathbf{S}}{{\mathbf{S}}_i}={\mathbf{\Gamma }}{{\mathbf{\upmu }}_{\mathbf{s}}}+{{\mathbf{\upvarepsilon }}_{{\mathbf{s}}i}}$$

(14, repeated)

where ${{\mathbf{\upmu }}_{\mathbf{s}}}$ is inflated relative to ${\mathbf{\upmu }}$ as a function of the summed factor loadings and ${{\mathbf{\upvarepsilon }}_{{\mathbf{s}}i}}$ is the sum of $p$ item residuals at each measurement occasion.

Appendix 2

Bias derivations

Bias is first evaluated in the means of the growth factors. Although the focus is on linear growth, these derivations could be extended to curvilinear growth. Recall that Eq. (6) gives $~E\left[ {{{\mathbf{y}}_i}} \right]={\mathbf{\Lambda \Gamma \upmu }}$ because the item intercepts ${\mathbf{\upnu }}$ are set to 0 without loss of generality and $E\left[ {{{\mathbf{\upxi }}_i}} \right]=~{\mathbf{\upmu }}={\text{~}}\left[ {\begin{array}{*{20}{c}} {{\mu _\alpha }}&{{\mu _\beta }} \end{array}} \right]'$. After pre-multiplying both sides of Eq. (6) with transposed loading matrices, the growth factor means can be isolated by pre-multiplying with $\left( \mathbf{\Gamma ^{'} \Lambda ^{'} \Lambda \Gamma } \right)^{{ - 1}}$,

$$\begin{gathered} E\left[ {{\mathbf{y}}_{i} } \right] = {\mathbf{\Lambda \Gamma \upmu }} \hfill \\ {\mathbf{\Gamma }}^{'} {\mathbf{\Lambda }}^{'} E\left[ {{\mathbf{y}}_{i} } \right] = {\mathbf{\Gamma }}^{'} {\mathbf{\Lambda }}^{'} {\mathbf{\Lambda \Gamma \upmu }} \hfill \\ \left( {{\mathbf{\Gamma }}^{'} {\mathbf{\Lambda }}^{'} {\mathbf{\Lambda \Gamma }}} \right)^{{ - 1}} {\mathbf{\Gamma }}^{'} {\mathbf{\Lambda }}^{'} E\left[ {{\mathbf{y}}_{i} } \right] = {\mathbf{\upmu }} \hfill \\ \end{gathered}$$

(39)

As described in the text, a fitted model may include some constrained factor loading matrix ${\mathbf{\tilde{\Lambda }}}$, but $E\left[ {{{\mathbf{y}}_i}} \right]~$is invariant to the constraints of ${\mathbf{\tilde{\Lambda }}}$, giving Eq. (11). This result generalizes to any form of ${\mathbf{\tilde{\Lambda }}}$, and the bias can be computed for the comparison of any misspecified measurement model to a known population value of ${\mathbf{\Lambda }}$.

Similarly, an expression that calculates the variance of the parameter estimates from constrained measurement models can be derived. Equation (7) establishes that the variance in the SLGM model is ${\mathbf{\Sigma }}_{y} = {\mathbf{\Lambda }}({\mathbf{\Gamma \Phi \Gamma ^{\prime}}} + {\mathbf{\Psi }}){\mathbf{\Lambda ^{\prime}}} + {\mathbf{\Theta }}.$ To solve for ${\mathbf{\Phi }}$, pre- and post-multiply both sides by the factor and growth loading matrices as follows,

$$\begin{gathered} {\mathbf{\Sigma }}_{y} = {\mathbf{\Lambda }}\left( {{\mathbf{\Gamma \Phi \Gamma ^{\prime}}} + {\mathbf{\Psi }}} \right){\mathbf{\Lambda ^{\prime}}} + {\mathbf{\Theta }} \hfill \\ {\mathbf{\Sigma }}_{y} = {\mathbf{\Lambda \Gamma \Phi \Gamma ^{\prime}\Lambda ^{\prime}}} + {\mathbf{\Lambda \Psi \Lambda ^{\prime}}} + {\mathbf{\Theta }} \hfill \\ ({\mathbf{\Sigma }}_{y} - {\mathbf{\Lambda \Psi \Lambda ^{\prime}}} - {\mathbf{\Theta }}) = {\mathbf{\Lambda \Gamma \Phi \Gamma ^{\prime}\Lambda ^{\prime}}} \hfill \\ {\mathbf{\Gamma ^{\prime}\Lambda ^{\prime}}}({\mathbf{\Sigma }}_{y} - {\mathbf{\Lambda \Psi \Lambda }}^{\prime} - {\mathbf{\Theta }}){\mathbf{\Lambda \Gamma }} = {\mathbf{\Gamma ^{\prime}\Lambda ^{\prime}\Lambda \Gamma \Phi \Gamma ^{\prime}\Lambda ^{\prime}\Lambda \Gamma }} \hfill \\ \left( {{\mathbf{\Gamma ^{\prime}\Lambda ^{\prime}\Lambda \Gamma }}} \right)^{{ - 1}} {\mathbf{\Gamma ^{\prime}\Lambda ^{\prime}}}({\mathbf{\Sigma }}_{y} - {\mathbf{\Lambda \Psi \Lambda ^{\prime}}} - {\mathbf{\Theta }}){\mathbf{\Lambda \Gamma }}\left( {{\mathbf{\Gamma ^{\prime}\Lambda ^{\prime}\Lambda \Gamma }}} \right)^{{ - 1}} = {\mathbf{\Phi }} \hfill \\ \end{gathered}$$

(40)

By replacing ${\mathbf{\Lambda }}$ with the constrained loading matrix ${\mathbf{\tilde{\Lambda }}}$, ${\mathbf{\hat{\Phi }}}$ can be obtained,

$$\left( {{\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}\tilde{\Lambda }\Gamma }}} \right)^{{ - 1}} {\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}}}({\mathbf{\Sigma }}_{y} - {\mathbf{\Lambda \Psi \Lambda ^{\prime}}} - {\mathbf{\Theta }}){\mathbf{\tilde{\Lambda }\Gamma }}\left( {{\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}\tilde{\Lambda }\Gamma }}} \right)^{{ - 1}} = {\mathbf{\hat{\Phi }}}$$

(41)

$$\left( {{\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}\tilde{\Lambda }\Gamma }}} \right)^{{ - 1}} {\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}}}\left( {{\mathbf{\Lambda \Gamma \Phi \Gamma \Lambda ^{\prime}}} + {\mathbf{\Lambda \Psi \Lambda ^{\prime}}} + {\mathbf{\Theta }} - {\mathbf{\Lambda \Psi \Lambda ^{\prime}}} - {\mathbf{\Theta }}} \right){\mathbf{\tilde{\Lambda }\Gamma }}\left( {{\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}\tilde{\Lambda }\Gamma }}} \right)^{{ - 1}} = {\mathbf{\hat{\Phi }}}$$

(42)

$${\mathbf{\hat{\Phi }}} = \left( {{\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}\tilde{\Lambda }\Gamma }}} \right)^{{ - 1}} {\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}}}\left( {{\mathbf{\Lambda \Gamma \Phi \Gamma ^{\prime}\Lambda ^{\prime}}}} \right){\mathbf{\tilde{\Lambda }\Gamma }}\left( {{\mathbf{\Gamma ^{\prime}\tilde{\Lambda }^{\prime}\tilde{\Lambda }\Gamma }}} \right)^{{ - 1}}$$

(43)

Again, it is important to note that the ${\mathbf{\Lambda }}$’s in the middle parentheses are not the constrained ${\mathbf{\tilde{\Lambda }}}$’s. These represent the population-level covariance structure, which needs to be assumed or fixed to evaluate the bias of a misspecified measurement model.

Now consider the form of ${\mathbf{\tilde{\Lambda }}}$ when fitting a univariate LGM to unweighted sums of items at each time point. The univariate model simply applies the linear growth function to one value at each time point, and the underlying individual item loadings ${\lambda _{jt}}$ are no longer considered once the items are summed to form ${{\mathbf{y}}_{{{\mathbf{s}}_i}}}$. Therefore, the univariate model treats the constrained loading matrix ${\mathbf{\tilde{\Lambda }}}$ essentially as an identity matrix. This results in fitting the univariate growth model of Eq. (1) to the aggregated sums of multiple items at each time point,

$$E[{\mathbf{SS}}_{i} ] = {\mathbf{\Gamma \hat{\upmu }}}_{{\mathbf{s}}}$$

(14, repeated)

where ${\mathbf{S}}{{\mathbf{S}}_i}$ is the aforementioned vector of sum scores for individual $i$ and ${\mathbf{\hat{\upmu }}}_{{\mathbf{s}}}$ is a vector of intercept and slope means estimated when the univariate model is fit to the sum scores. The derivation for the estimated intercept and slope means for the sum score then has the form,

$${\mathbf{\hat{\upmu }}}_{{\mathbf{s}}} = ({\mathbf{\Gamma ^{\prime}\Gamma }})^{{ - 1}} {\mathbf{\Gamma ^{\prime}}}E[{\mathbf{SS}}_{i} ]$$

(15, repeated)

with the true form of $E[{\mathbf{S}}{{\mathbf{S}}_i}]$ given by Eq. (9) under the assumption that the true data-generating process is actually the SLGM for the item-level data. The estimated intercept and slope covariance matrix is similarly derived for sum scores, where ${\mathbf{\tilde{\Lambda }}}$ is treated as identity from Eq. (13) and ${{\mathbf{\Lambda }}_{\mathbf{s}}}$ from Eq. (9) is the population-level ${\mathbf{\Lambda }}$ for the summed items, giving,

$${\mathbf{\hat{\Phi }}}_{{\mathbf{s}}} = \left( {{\mathbf{\Gamma ^{\prime}\Gamma }}} \right)^{{ - 1}} {\mathbf{\Gamma ^{\prime}}}\left( {{\mathbf{\Lambda }}_{{\mathbf{s}}} {\mathbf{\Gamma \Phi \Gamma ^{\prime}\Lambda }}_{{\mathbf{s}}} ^{\prime } } \right){\mathbf{\Gamma }}\left( {{\mathbf{\Gamma ^{\prime}\Gamma }}} \right)^{{ - 1}}$$

(18, repeated)

Appendix 3

The measurement loading matrices that were used in the analytic demonstrations are given below:

$${{\mathbf{\Lambda }}_1}=\left[ {\begin{array}{*{20}{c}} 1&0&0&0 \\ {.7}&0&0&0 \\ {.5}&0&0&0 \\ 1&0&0&0 \\ {.7}&0&0&0 \\ {.5}&0&0&0 \\ 0&1&0&0 \\ 0&{.7}&0&0 \\ 0&{.5}&0&0 \\ 0&1&0&0 \\ 0&{.7}&0&0 \\ 0&{.5}&0&0 \\ 0&0&1&0 \\ 0&0&{.7}&0 \\ 0&0&{.5}&0 \\ 0&0&1&0 \\ 0&0&{.7}&0 \\ 0&0&{.5}&0 \\ 0&0&0&1 \\ 0&0&0&{.7} \\ 0&0&0&{.5} \\ 0&0&0&1 \\ 0&0&0&{.7} \\ 0&0&0&{.5} \end{array}} \right]\quad {{\mathbf{\Lambda }}_2}=\left[ {0\begin{array}{*{20}{c}} 1&0&0&0 \\ {.7}&0&0&0 \\ {.7}&0&0&0 \\ {.7}&0&0&0 \\ {.7}&0&0&0 \\ {.7}&0&0&0 \\ 0&1&0&0 \\ 0&{.7}&0&0 \\ 0&{.7}&0&0 \\ 0&{.7}&0&0 \\ 0&{.6}&0&0 \\ 0&{.6}&0&0 \\ 0&0&1&0 \\ 0&0&{.7}&0 \\ 0&0&{.7}&0 \\ 0&0&{.7}&0 \\ 0&0&{.5}&0 \\ 0&0&{.3}&0 \\ 0&0&0&1 \\ 0&0&0&{.7} \\ 0&0&0&{.7} \\ 0&0&0&{.7} \\ 0&0&0&{.5} \\ 0&0&0&{.3} \end{array}} \right]\quad {{\mathbf{\Lambda }}_3}=\left[ {\begin{array}{*{20}{c}} 1&0&0&0 \\ {.7}&0&0&0 \\ {.7}&0&0&0 \\ {.3}&0&0&0 \\ {.3}&0&0&0 \\ {.5}&0&0&0 \\ 0&1&0&0 \\ 0&{.5}&0&0 \\ 0&{.5}&0&0 \\ 0&{.5}&0&0 \\ 0&{.5}&0&0 \\ 0&{.5}&0&0 \\ 0&0&1&0 \\ 0&0&{.5}&0 \\ 0&0&{.5}&0 \\ 0&0&{.5}&0 \\ 0&0&{.5}&0 \\ 0&0&{.5}&0 \\ 0&0&0&1 \\ 0&0&0&{.3} \\ 0&0&0&{.3} \\ 0&0&0&{.7} \\ 0&0&0&{.7} \\ 0&0&0&{.5} \end{array}} \right]$$

These correspond to the description of the three scenarios described in the text.

Appendix 4