Operating Characteristics of Statistical Methods for Detecting Gene-by-Measured Environment Interaction in the Presence of Gene-Environment Correlation under Violations of Distributional Assumptions.

In the past few decades it has become increasingly clear that any inquiry into the roots of psychopathology such as anxiety or depression, as well as other complex behaviors, requires accounting for possible interactions and correlations between genetic vulnerabilities and environmental factors. Historically, genetically informative models assumed that genetic and environmental influences on a particular trait were static across the population. But growing evidence points both to differential effects of the same environmental exposure across genotypes (gene-by-environment interaction) and to differential environmental exposures across genotypes (gene-by-environment correlation). Several methods have been proposed to model these more complex relationships, particularly within the context of twin and family studies (Eaves & Erkanli, 2003; Price & Jaffee, 2008; Purcell, 2002; Rathouz, Van Hulle, Rodgers, Waldman, & Lahey, 2008). An important methodological consideration with the models proposed by ourselves and others in this recent literature is that these are full probability, non-linear structural equations models (SEM). As such, they are based on distributional assumptions, such as multivariate normality of latent genetic and latent environmental factors. In practice, many phenotypes of interest are not normally distributed. Data may be ordinal (e.g. behavior ratings of impulsivity), skewed (e.g. symptom counts of depression) or censored (e.g. age-of-onset). In contrast to recent work, classical twin and adoption study data analysis methods—which do not posit any interaction effects—rely on linear SEMs. Whereas normality may be a useful working assumption in these models, valid inferences are often based only on assumptions about the first two moments (mean, variance and covariance) of the data. Violations of normality have a negligible effect on parameter estimates in such models, and methods are available to adjust standard errors for bias due to non-normality. In the presence of gene-by-environment interactions, however, because they involve the product or square of latent normal quantities, the manifest variables will be non-normal by construction. Alternatively, when the latent factors are normal and do not interact, but the latent errors or measurement errors are non-normal, the manifest variables will also be non-normal. Therefore, when the scale of measurement of the variable(s) of interest is inherently non-normal, it is questionable as to whether the data can distinguish between these two fundamentally different scenarios. The issue of robustness of current G×M analysis methods to violation of distributional assumptions is therefore critical to behavior genetic designs being used in investigations of psychopathology, in particular where many phenotypes are measured with highly skewed distributions (e.g. symptom counts), and requires thorough exploration before any of these methods can be reliably used in such studies. The goal of the current paper is to explore the existence and severity of consequences of such violations of distributional assumptions on statistical tests and estimators. We consider bivariate behavior genetic designs involving a measured environment M and its potential moderating effects (G×M) on variance components impacting on a phenotype of interest P. The models allow for correlation between M and variance components of P (rGM). The specific question is whether, in the context of the set of models laid out in Purcell (2002) and in Rathouz et al. (2008), the data are able to distinguish between non-normality in manifest variables due to G×M versus that due to measurement properties of the phenotype. In our previous work (Van Hulle, Lahey, & Rathouz, 2013), we evaluated the Type I error rates, power, and performance of the Bayesian Information Criterion (BIC) for testing and comparing a subset of the models proposed in Rathouz et al. (2008), equations for which are shown in Table 1. In that paper, data were simulated under a variety of conditions both with and without G×M interactions. To briefly summarize, we found that: (i) when comparing the Cholesky with G×M model with the various submodels, the false positive rates consistently fell short of the nominal Type I error rates (α =.10, .05, .01); (ii) with larger sample sizes (N=2000),in nearly all cases the correct model had the lowest BIC value across all possible models; (iii) with lower sample sizes (N=500), models specifying non-linear main effects were more difficult to distinguish from models containing interaction effects. In that paper, all simulated latent quantities and error variables were normally distributed, thus matching the distributional assumptions of the models. For the current study, we examined Type I error rates and performance of the Bayesian Information Criterion (BIC) for nested and non-nested models under similar data generating mechanisms as in Van Hulle et al. (2013), but where violations of normality in the measurement process were imposed on either the moderator (M) or the phenotype (P). Table 1 Bivariate variance components models for latent component-by-measured environment (G×M) interaction.