As noted above, the Whānau Pakari programme evolved through two sequential phases. The initial trial (2012–2014) compared a high-intensity intervention (weekly group sessions during school term time for 12 months plus 6-monthly assessments) with a low-intensity control (6-monthly assessments and advice only) [16]. After the trial, the programme was adopted as the regional childhood obesity service in Taranaki (2014 onwards). The clinical service maintained core elements of the trial intervention while adjusting the weekly session duration to six months due to waning attendance in the second 6-month period. This was supported by evidence suggesting that ≥ 26 hours of contact time (equivalent to approximately 26 weekly one-hour sessions) is required to achieve meaningful BMI SDS reductions [19]. All participants received support from a multidisciplinary team, including a pediatrician, psychologist, dietitian, and physical activity coordinator, with screening and action for any identified weight-related comorbidities [16].

This study analysed combined data from both the trial and the subsequent clinical service to maximise statistical power and capture real-world data. The trial and service data could be integrated because intervention delivery and outcome measurements remained consistent across both phases during the first six months. This methodological consistency allowed for robust data pooling. Furthermore, referrals to the programme were accepted year-round which minimises the risk of seasonal bias. We are also unaware of other factors that may have biased the recruitment of participants into either the trial or the service. Our statistical analyses accounted for possible differential associations between study cohorts (i.e., trial or service).

Eligibility criteria remained constant throughout both programme phases. Children and young people aged approximately 4.0 to < 18.0 years were eligible if their BMI was at or above the 98th percentile for age and sex, or at or above the 91st percentile with concurrent weight-related comorbidities, according to the UK 1990 reference data used for screening in NZ [20]. Participants also had a baseline assessment in January 2012–June 2020 and a follow-up assessment within 6 months ± 5 weeks (i.e., 183 ± 35 days). Exclusion criteria encompassed medical conditions affecting their ability to engage in physical activity, psychological conditions affecting group participation or the absence of committed family support [21].

Holistic clinical assessments occurred at the participant's home or a community venue agreed to by the family. Demographic information collected via caregiver questionnaires included the participant's age, biological sex (male or female) and self-reported ethnicity (based on prioritised ethnicity used by the health board at the time) [22]. Weight was measured using SECA 813 digital scales (SECA, Hamburg, Germany) and height with SECA 213 portable stadiometers. BMI was subsequently calculated. Waist circumference was also measured and the waist-to-height ratio was calculated. Participants were initially screened in clinic using the UK 1990 (UK90 or ‘UK Cole’) growth reference [20]. Weight, height and BMI were converted to age- and sex-specific standard deviation scores (SDS or z-scores) using the age-appropriate WHO growth system (Supplementary Methods; Supplementary Fig. 1).

Dietary patterns (24-h recall) were assessed using the modified Children’s Dietary Questionnaire (CDQ), adapted to the NZ context [23]. Physical activity levels were measured using the Children’s Physical Activity Questionnaire (C-PAQ) [24]. For children aged < 11 years, questionnaires were completed jointly by the caregiver and child. Caregivers reported the duration of sleep during the collection of medical history. Parameters of interest were self-reported (or proxy reported, where age-appropriate): physical activity levels, fruit and vegetable servings per day, time spent on screens, sleep duration and intake of sweet drinks.

The primary outcome was change (Δ) in BMI SDS at 6 months. Secondary outcomes included the likelihood of a reduction in BMI SDS, changes in other anthropometric parameters and changes in diet and lifestyle measures. Participants were stratified by Southern Hemisphere season at baseline, defined by the meteorological criteria [25]: summer (December–February), autumn (March–May), winter (June–August) and spring (September–November). Our analytical approach combined traditional linear models with machine learning techniques.

Demographic and anthropometric characteristics at baseline were summarised as means ± standard deviation (SD), medians [Q1, Q3] or n (%) and compared between seasons using Fisher’s exact tests for categorical variables, one-way ANOVA for continuous variables approximating a normal distribution on visual inspection, or non-parametric Kruskal–Wallis tests for skewed continuous variables.

Univariable (unadjusted) Δ from baseline (overall and within season) were examined using paired t tests, while between-season differences were assessed with one-way ANOVA. Differences were reported as mean Δ with 95% confidence intervals (CI). Linear associations with continuous predictors were analysed using Pearson’s or Spearman’s rank correlation coefficients.

Associations between season and Δ BMI SDS at 6 months were analysed using a generalised linear mixed model (GLMM). The model included family ID as a random intercept, accounting for the clustering of siblings, effectively nesting participant-level data within family clusters. Based on existing evidence of factors influencing childhood obesity intervention outcomes (particularly in NZ), the model included season (the primary exposure of interest), baseline BMI SDS, age at baseline, sex, ethnicity (Māori vs Non-Māori) and programme cohort (service vs trial) as fixed effects. We tested for a season*cohort interaction to assess for differential influences of the controlled trial environment and the routine service delivery on outcome; if non-significant (P ≥ 0.05), the interaction term was removed from the final model. This a priori variable selection approach avoided data-driven model building that could lead to overfitting with our sample size (n = 397). The same analytical approach was applied to other anthropometric outcomes (weight SDS, height SDS and waist-to-height ratio). GLMM results are presented as adjusted means (i.e., least-squares means) with 95% CI; between-group and within-group differences are reported as adjusted mean differences (aMD) with 95% CI.

Seasonal rates of BMI SDS reduction at 6 months were compared with a Fisher's exact test. Further, the likelihood of a BMI SDS reduction was assessed with a multivariable generalised linear model using a modified Poisson procedure with robust error variances [26]. This model included the same covariates as specified for the GLMM, with effect size expressed as the adjusted relative risk (aRR) with 95% CI.

Baseline dietary and lifestyle variables were compared between seasons using Kruskal–Wallis tests, followed by unadjusted pairwise comparisons using Wilcoxon rank-sum tests. Univariable overall and within-season Δ were assessed with Wilcoxon signed-rank tests. Multivariable ranked GLMMs were run to examine possible between-season Δ differences in these parameters compared to baseline and these were constructed as previously described. If there was evidence of between-season differences (P < 0.05 in the ranked GLMM), pairwise effect magnitudes were quantified using Wilcoxon-derived Hodges–Lehmann location shift estimates with 95% CI.

Data were analysed using SAS v9.4 (SAS Institute, Cary, NC, USA). All statistical tests were two-sided with significance set at P < 0.05. Given the exploratory nature of the study, planned pairwise comparisons between seasons were not adjusted for multiple comparisons to avoid inflating the risk of Type II errors, which could obscure clinically relevant patterns [27,28,29].

Given the complex associations and potential non-linear relationships between lifestyle, demographic and clinical factors that may influence changes in BMI SDS, we used random forest modelling to capture non-linear effects and higher-order interactions that may not be readily detected by traditional linear models [30, 31]. Random forests can handle multiple correlated predictors and implicitly model interactions, making them suitable for identifying factors associated with Δ BMI SDS within our pediatric population.

Random forest analyses were conducted in R v4.4.1 [32] using the ranger package for model fitting, with supporting packages for data management and visualisation (readxl, writexl, ggplot2, dplyr, caret, forcats, tibble and stringr). Two random forest models were developed to predict Δ BMI SDS over the 6-month follow-up period. The base model incorporated demographic, clinical and temporal variables (season, programme cohort, sex, ethnicity, age at baseline and baseline BMI SDS). The expanded model included these variables plus dietary and lifestyle factors (i.e., sweet drink intake, physical activity levels, screen time, sleep duration and fruit/vegetable servings per day).

Missing values were imputed within each training fold [mean for continuous variables; mode (most frequent category) for categorical variables] ensuring no information from the corresponding test folds was used during model training. Model performance was evaluated using repeated fivefold cross-validation (10 repeats) to reduce variance and guard against overfitting.

Variable importance was assessed using permutation importance (mean increase in prediction error after permuting each predictor) as implemented in ranger [30]. Importance scores were averaged across cross-validation folds. For interpretability plots, we used out-of-fold (OOF) predictions, whereby each participant’s predicted outcome was generated exclusively from models that did not include that participant in training data.