The twins were treated as individuals by adjusting standard errors for the clustering of observations in twin pairs using the type=complex specification and the cluster option implemented in Mplus 5.2 (Muthén and Muthén, 2007). This option uses the pseudomaximum likelihood method in which the asymptotic covariance matrix of the estimates is obtained by a robust sandwich estimator (Asparouhov, 2005). Bivariate cross-lagged path models were fitted separately for each substance use variable (any drinking, drinking to intoxication, smoking, and daily smoking) with educational achievement across the four study waves, adjusting for covariates. Besides autoregressive and cross-lagged paths between successive time points, the models included additional autoregressive paths on variables at waves 3 and 4 from the equivalent measures two waves earlier (e.g., a path from drinking frequency at age 14 to drinking frequency in young adulthood). Such paths may be used to improve model fit in models with three or more measurements (Tucker et al., 2013). Residual correlations at each wave were also modeled in order to account for the remaining cross-sectional correlation between educational achievement and substance use, unexplained by the cross-lagged associations. The full estimated model is presented as a path diagram in Figure 1.
