Study design and participants

Data were collected through the Pathways in ASD study, a large (n = 421) longitudinal inception cohort study of children with ASD and their families that began in 2005. Details of the study methodology are described elsewhere.^{Reference Vaillancourt, Haltigan, Smith, Zwaigenbaum, Szatmari and Fombonne8} Participants were recruited between ages 2 years and 5 years, across five community and academic Canadian sites, within 4 months of receiving an ASD diagnosis (sites in Halifax, Nova Scotia; Montreal, Quebec; Hamilton, Ontario; Edmonton, Alberta; and Vancouver, British Columbia). Participants were then followed longitudinally through a variety of standardised assessment and report measures, using an accelerated design (i.e. staggered entry age, then synchronisation of assessment time points at age 6). The current study focused on the first eight assessment time points: T ₁ (mean age 3.3 years) through T ₈ (mean age 10.7 years).

All procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008. All procedures involving human participants were initially approved at McMaster University by the Hamilton Integrated Research Ethics Board (REB) (project number 04-353) and subsequently by REBs at each participating site. Written informed consent was obtained from all participating families.

Behavioural assessments

Severity of insistence on sameness behaviour was measured using the Autism Diagnostic Interview-Revised (ADI-R).^{Reference Lord, Rutter and Le Couteur23} The ADI-R is a semi-structured validated parent interview; clinicians quantify current and past behaviour on a scale of 0–3, with 3 representing an ASD-associated behaviour that causes severe impairment. The ADI-R has a widely used 12-item subscale capturing ‘stereotyped and repetitive interests and behaviours’. Factor analysis has been used repeatedly to confirm that items on this subscale reliably load onto one of two underlying constructs: insistence on sameness behaviour and repetitive motor behaviour.^{Reference Szatmari, Georgiades, Bryson, Zwaigenbaum, Roberts and Mahoney14,Reference Bishop, Hus, Duncan, Huerta, Gotham and Pickles24,Reference Leekam, Prior and Uljarevic25} For this study, we used the total raw current ADI-R Insistence on sameness (ADI-R-IS) subscale scores (6 items, total possible range 0–18, actual range 0–12) collected at T ₁, T ₄ and T ₈.

Anxiety symptoms were measured using the DSM-oriented Anxiety Problems subscale of the Child Behavior Checklist (CBCL),^{Reference Achenbach and Dumenci26} which was completed by parents at each study time point (T ₁ through T ₈). There are two versions of the CBCL: one for school-aged (6–18 years) and one for preschool-aged (1.5–5 years) children. The CBCL proprietary software provides age- and gender-normed standardised scores (T-scores) across a range of common childhood psychiatric disorders and symptom profiles, with T-scores above 65 considered to indicate a clinical concern because of elevated symptomatology.^{Reference Achenbach and Dumenci26} The CBCL Anxiety subscale appears to have adequate sensitivity (0.83) and moderate specificity (0.52) as a screening tool to detect anxiety disorders in children with ASD (AUC = 0.71).^{Reference Magyar and Pandolfi27}

Other baseline (T ₁) measures used to characterise the sample and the different trajectory groups were: age (in months), gender, family income on an 11-point scale from <Can$5000 to >Can$80 000 per year; parenting stress levels on the Parenting Stress Index-Short Form (PSI-SF) (total stress score);^{Reference Abidin28} Adaptive Behavior Composite score from the Vineland Adaptive Behavior Scales Second Edition (VABS-II);^{Reference Sparrow, Balla and Chicchetti29} IQ on the Developmental Index of the Merrill–Palmer-Revised Scales of Development (M-P-R);^{Reference Roid and Sampers30} and language abilities on the Preschool Language Scale Fourth Edition (PLS-4).^{Reference Zimmerman, Steiner and Pond31}

Analyses

Analyses were conducted in SAS (version 9.5 for Windows). Preliminary analyses using exploratory factor analysis (EFA) were initially used to examine the extent to which symptoms on both measures separated into distinct constructs at each study time point.

Clusters of participants following similar longitudinal trajectories over time were then identified using a well-established finite-mixture modelling procedure (PROC TRAJ).^{Reference Nagin and Odgers32,Reference Jones and Nagin33} This method is useful to examine heterogeneity and describe multiple distinct patterns of change over time in a population.^{Reference Arrandale, Koehoorn, MacNab and Kennedy34}

Multigroup trajectories were estimated for the continuous ADI-R insistence on sameness subscale raw total scores (at T1, T ₄ and T ₈) and the CBCL anxiety subscale T-scores (at T ₁ through T ₈). For the anxiety trajectories, T-scores (as opposed to raw scores) were used to facilitate clinically interpretable modelling of anxiety across both versions of the CBCL (preschool and school-age). Specifically, as the raw scores and number of items from the two versions are incongruous, and since many anxiety symptoms (e.g. separation anxiety) are expected or normative at younger ages, analyses using T-scores were chosen to examine anxiety symptoms in ASD relative to typically developing norms. We also subsequently examined anxiety trajectories using mean item raw scores as a sensitivity analysis. Note that the CBCL proprietary software truncates CBCL T-scores at 50. To preserve maximum heterogeneity in the data, at each study time point we imputed anxiety T-scores for participants who had initially been assigned a T-score of 50, based on their age, gender and raw score, in keeping with the CBCL scoring procedures. Multiple imputation was considered for this step.

All trajectories of insistence on sameness and anxiety were first modelled using quadratic growth curves. The optimal number of trajectory groups was selected on the basis of inspection of the Akaike information criterion (AIC), Bayesian information criterion (BIC) and log-likelihood values, while considering parsimony, group sizes and previous literature. Once the optimal number of groups was identified, the polynomial order was varied (e.g. linear, cubic) and model fit statistics were further compared. Model adequacy was assessed by considering the odds of correct classification (target ≥5) and the average posterior probability of assignment (APPA) (target >0.70) for each group.

Descriptive characteristics at T ₁ were compared between assigned trajectory groups. As an exploratory analysis, we also examined T ₁ measures as predictors of group assignment.^{Reference Roeder, Lynch and Nagin35} Trajectory group assignment was treated as a probabilistic ordinal outcome, and the odds of transitioning to a higher-severity trajectory group relative to the reference low-severity group was estimated for time-stable baseline covariates, while accounting for uncertainty in group assignment.^{Reference Jones, Nagin and Roeder36}

A joint trajectory model of insistence on sameness and anxiety was subsequently estimated to allow examination of these simultaneous developmental paths. We estimated: (a) the probability of belonging to each anxiety trajectory, conditioned on the insistence on sameness trajectory assignments; (b) the probability of belonging to each insistence on sameness trajectory, conditioned on the anxiety trajectory assignments; and (c) the joint probabilities of belonging simultaneously to all combinations of anxiety–insistence on sameness trajectories.^{Reference Nagin and Tremblay37} Three select joint trajectory groups of interest with unique anxiety profiles to which >10% of the study sample were assigned and a low-risk reference group were compared in terms of descriptive characteristics.

With respect to missing data, PROC TRAJ uses the maximum likelihood method to estimate parameters from all available data and therefore permits inclusion of participants with missing information, provided that it is missing at random (MAR).^{Reference Arrandale, Koehoorn, MacNab and Kennedy34} We conducted preliminary analyses to explore patterns and correlates of missingness. We then conducted sensitivity analyses, including restricting group comparisons to participants with APPA ≥0.90, to account for potential misclassification. We also repeated the unidimensional trajectory analyses by (a) excluding participants with high proportions of missing data (i.e. missing >1 of 3 ADI-R time points, or >3 of 8 CBCL time points) and (b) imputing missing outcome data through predictive mean matching using the following auxiliary variables: age at diagnosis, gender, baseline VABS-II and PSI-SF scores. We also estimated the trajectories while accounting for the possibility that participant attrition is group dependent.^{Reference Haviland, Jones and Nagin38}