Mapping Geographic Trends in Early Childhood Social, Emotional, and Behavioural Difficulties in Glasgow: 2010–2017

Measuring variation in childhood mental health supports the development of local early intervention strategies. The methodological approach used to investigate mental health trends (often determined by the availability of individual level data) can affect decision making. We apply two approaches to identify geographic trends in childhood social, emotional, and behavioural difficulties using the Strengths and Difficulties Questionnaire (SDQ). SDQ forms were analysed for 35,171 children aged 4–6 years old across 180 preschools in Glasgow, UK, between 2010 and 2017 as part of routine monitoring. The number of children in each electoral ward and year with a high SDQ total difficulties score (≥15), indicating a high risk of psychopathology, was modelled using a disease mapping model. The total difficulties score for an individual child nested in their preschool and electoral ward was modelled using a multilevel model. For each approach, linear time trends and unstructured spatial random effects were estimated. The disease mapping model estimated a yearly rise in the relative rate (RR) of high scores of 1.5–5.0%. The multilevel model estimated an RR increase of 0.3–1.2% in average total scores across the years, with higher variation between preschools than between electoral wards. Rising temporal trends may indicate worsening social, emotional, and behavioural difficulties over time, with a faster rate for the proportion with high scores than for the average total scores. Preschool and ward variation, although minimal, highlight potential priority areas for local service provision. Both methodological approaches have utility in estimating and predicting children’s difficulties and local areas requiring greater intervention.


Introduction
It is estimated that one in every five children under 7 years old meets the criteria for a mental health diagnosis [1]. Identifying mental health problems when they first become evident in younger children [2] and delivering early interventions can reduce the individual and societal impact [3]. We can assess the presence of social, emotional, and behavioural difficulties (e.g., rating a child's willingness to help others, their relationships with peers, emotional regulation, and inattention), which indicate if a child is at risk of mental health problems [4,5]. Childhood social, emotional, and behavioural problems can predict school exclusion [6] and later mental health conditions [7] at an individual level, and lead to higher healthcare and social welfare costs at a societal level [8][9][10]. Identifying patterns and trends in early childhood mental health can support the development of targeted interventions.
This paper focuses on geographic patterns in mental health according to where a child lives. Geographic variation in outcomes can be expected as the residential environment influences the day-to-day social and physical experiences of the child and their caregiver.

1.
Map the overall risk of preschool children's social and emotional difficulties in an entire city from 2010-2017; 2.
Analyse the relationship between preschool children's social, emotional, and behavioural difficulties in neighbouring areas across a whole city from 2010-2017; 3.
Investigate the extent to which the relationship between geography and preschool children's social, emotional, and behavioural difficulties is explained by demographics;

4.
Explore the temporal and spatiotemporal trends in preschool children's social, emotional, and behavioural difficulties across an entire city from 2010-2017.

Setting
The ChiME (Child Mental Health in Education) study took place in Glasgow City, Scotland, from 2010-2017, the most populated city in Scotland, United Kingdom. It has been shown that the levels of social, emotional, and behavioural difficulties in Glasgow are similar to the rest of the UK at preschool age [28].

Study Design
The ChiME data were routinely collected by the Glasgow City Council, providing a whole population sample from 2010-2017 (derived from the Triple P parenting intervention [29]). Previous spatial analysis using a subset of ChiME data is published elsewhere [21].
Preschool staff were asked to complete the teacher version of the Strengths and Difficulties Questionnaire (SDQ) [30]. Forms were completed on paper (2010), through a mix of paper and electronic submissions (2011) or solely electronically (2012-2017) using the education services' information management system. For 2010 and 2011, the SDQ version used was for 4-16-year-olds. This was changed to the version for 2-4-year-olds in subsequent years due to staff perception that question wording for conduct problems in that version was more developmentally appropriate.
All children living in Glasgow, who were attending a local authority preschool or private preschool with funded places in 2010-2017, and who were due to start school the following year, were eligible. During the study period, there were 99,346 children aged 4-5 years in Glasgow [31]. In total, data were collected from 41,128 children from 182 preschools in the study. Children in private nurseries without funded places or those not in preschool were not part of the study. Approximately 90% of 4-year-old children in Scotland are registered in some form of early learning and childcare. There is no relationship between being registered in a preschool and area deprivation in Scotland [32].

Variables
Biological sex at birth and date of birth were obtained from educational services' administrative databases and were linked to the ChiME database. Age was calculated in years from planned start in primary school (1 August following data collection). Each child was assigned a Scottish Index of Multiple Deprivation (SIMD) [33] score based on their residential postcode and cohort. As a result of the local authority's considerable levels of socioeconomic deprivation, scores were categorised into quintiles within Glasgow rather than for the whole of Scotland to form a local index. In the absence of household data, deprivation quintiles were used as a proxy for household socioeconomic status.

Outcome
In the SDQ, there are 25 questionnaire items scored on a three-point Likert scale (0, 1, 2) across five domains: emotional symptoms, conduct problems, hyperactivity/inattention, peer problems, and prosocial behaviour. The first four of these domains contribute to a "total difficulties" score out of 40. Scoring was based on the code provided on the SDQ website [34]. New banding classifies scores from 0-10 as close to average, 11-14 as slightly raised, 15-17 as high, and 18 and above as very high for 2-4-year-olds in the UK [35]. High scores indicate a high likelihood of psychiatric diagnosis.

Geographical Data
Residential postcodes were situated within the 21 electoral wards of Glasgow at that time. Electoral wards were aligned with local planning and policy making and were updated to reflect population change. Two new ward boundaries were introduced after the study period in 2017, creating 23 boundaries, with a new ward created in both the east and west of the city. The ward boundaries were obtained from SASPAC (Small Area Statistics PACkage) at https://saspac.org/ (accessed on 2 July 2020), using the boundaries at 2011.

Statistical Methods
Descriptive analyses compared the median total difficulties scores and the proportion of children with high scores by demographic characteristics, wards, and years.
Models were estimated using Bayesian inference [36], which is increasingly applied in child psychology research to support the analysis of complex data structures [37]. Analysis was conducted in R version 4.0.1 [38] and using package R-INLA (www.r-inla.org, accessed on 2 February 2020) [39] with weakly informative priors [40].
Both models followed the spatiotemporal formulation by Bernardinelli et al. [41]. According to this formulation, we included a linear time trend common to all wards in the city, an overall spatial random intercept (describing how wards deviate from the city-wide average), and a differential time trend (random slope) estimating how much a ward deviates from the linear trend. The intercept and slope may be positively correlated, meaning wards with high difficulties have a higher differential trend, or negatively correlated, where wards with high difficulties have a lower differential trend.
Spatial correlation of residuals aggregated by ward can be measured through Moran's I statistic [42]. For the multilevel model, we generated residuals by simulating values from the distribution and scaled them using R package DHARMa [43]. The value of I estimated the correlation between adjacent wards-using queen adjacency (wards sharing a border or point were considered neighbouring). The value of I indicates where neighbouring wards were similar (I ≈ 1), dissimilar (I ≈ −1), or had no association (I ≈ 0). Moran's I test assessed the null hypothesis that the value of I that was observed was equal to zero and dictated the structure of the spatial random intercept in the model.

Disease Mapping Model
The first outcome of interest was the number of children with or at high risk of psychopathology in an area. The outcome was defined as the number of children in each ward and year with a high total difficulties score (≥15). We assumed a Poisson distribution for the outcome, as high scores are discrete and relatively rare.
The model was specified as follows: where Y j was the number of high scoring children in ward j = 1, . . . , 21, E j was the population of the ward (used as an offset), and θ j was the probability of having a high score. In the model, β 0 was the intercept, β 1 was the main linear trend across the eight yearly cohorts, and u 0j the ward-level spatial intercept. There was no evidence of a spatial correlation in the model; therefore, the spatial effect was unstructured to model spatial heterogeneity. The random slope u 1j allowed for separating linear trends for each ward.
Although shown here to be independent of random intercepts for simplicity, correlated intercepts and slopes were assessed during model building. Coefficient β k was associated with ward characteristic k. Ward characteristics were the proportion of children outside the expected age for school start (under 4.5 years or over 5.5 years), the proportion of children in the most deprived quintile, and the proportion of boys. The parameters were added to the model in the order described and were retained based on the Deviance Information Criterion (DIC) [44], with lower values indicating better fit.

Multilevel Model
The outcome of interest was the total difficulties score for an individual child, which was modelled using a multilevel model. An advantage of using this outcome is that we can examine how scores change on average, rather than solely focusing on whether a child has reached a cut-off. This still has clinical relevance as population average scores relate to rising rates of psychiatric diagnosis [45]. Furthermore, by retaining the individual level data, a multilevel approach that considers multiple contextual effects can be applied.
Normative SDQ data of children in the UK show that individual scores are skewed, with many zeros and large variance requiring a distribution that can model zero inflation and overdispersion, respectively [46]. Therefore, we assumed a zero-inflated negative binomial distribution [47,48]. The mean of the model is λ, p is the zero-inflation parameter, and r is the overdispersion parameter.
Here, child i = 1, . . . , N was nested within ward j = 1, . . . , 21 and preschools k = 1, . . . , 180. In the model, β 0 was the intercept, β 1 the coefficient for the overall time trend, v 0j the ward-level spatial intercept, and α k a preschool effect. There was no evidence of a residual spatial correlation of ward level residuals in the model; therefore, the spatial effect was unstructured. The preschool effect and ward were cross-classified (i.e., they were not hierarchically nested). Random slope v 1j allowed for separating linear trends for each ward. As in the disease mapping model, correlated random intercepts and slopes were considered. Coefficient β ϕ was associated with individual-level variables ϕ-age (centred on mean age of 59 months and squared), sex, and deprivation quintile. The parameters were added to the model in the order described and were retained based on the Deviance Information Criterion (DIC) [44].

Inference
For both modelling approaches, relative rates (RR) were derived from exponentiated parameters. For example, RR = 1.10 would mean there was a 10% higher rate or score, while RR = 0.90 would mean it was 10% lower. In the disease mapping approach, RR > 1 represented an increase in the number of high scoring children relative to the population of the ward. In the multilevel approach, RR > 1 indicated an increase in the scores of individual children, on average.
The exponentiated intercepts for the ward effects (exp u 0j and exp v 0j in Equations (1) and (2), respectively, and school effects (exp(a k ) in Equation (2)) represent the RR for each ward and school compared with the average. The RR of the random slopes (exp u 1j and exp v 1j ) in Equations (1) and (2), respectively, described how much steeper the linear time trend was in each ward compared with the overall time trend.
Once the covariates were added to the models, the random effects described the unexplained effects after adjustment. Uncertainty in the ward level RR was represented using the exceedance probability, i.e., the probability that the corresponding RR was greater than 1. Exceedance probability values of 0.8 or above were considered to indicate high certainty in the elevated RR.

Descriptive Data
Of the 41,128 in the initial dataset, children were excluded if they lived outside Glasgow, were missing date of birth, were under 4 or over 6 years old at primary school entry, had a missing or invalid postcode, or were missing a total difficulties score (n = 5957), resulting in 35,171 included children in 180 preschools. Of the included preschools, 60 returned data every year, while the others were involved for a subset of the study years. The characteristics of the preschools by the number of years they were involved in the study is shown in Supplementary Table S1. Table 1 shows the distribution of demographic characteristics in the included sample and among those with high scores (SDQ ≥ 15; representing 9.0% of the total sample). Children who were outside the expected age range (4.5-5.5), those living in more deprived areas, and boys were observed to have higher levels of difficulties. The proportion of children with high scores increased over time. The distribution of demographics against total scores and for each year are displayed in Supplementary Figure S1 and Supplementary Table S2. In wards 3 and 6, over half of the children were in the most deprived quintile (Supplementary Table S3). While in wards 11 and 17, just under half of children were in the least deprived quintile (Supplementary Table S3). The median score per ward had a narrow range (3 to 5) and there was some geographical clustering of similar median scores (Figure 1, left). The proportion of children with high scores ranged between wards from 6.37% to 11.9% (Figure 1, right). The three areas with the highest rate of high scores were found in the east of the city. Wards 3, 19, and 21 had rates of 11.4%, 11.2%, and 11.9%, respectively. There were several areas in the south of the city (wards 13, 18, and 20) where under 7% of the children had high scores. For outcomes by ward and cohort, see Supplementary Figures S2 and S3.  Table 2 shows the results from the unadjusted and adjusted final disease mapping models. Random slopes were removed from the model as they did not improve fit (DIC = 1023 vs. DIC = 1022) compared with the model without random slopes. For each model, an increase in the relative rate (RR) of high scores per cohort was identified, with an estimated yearly rise of between 1.5% and 5.0% after adjustment. Of the demographics, only the proportion of boys in a ward was associated with the RR of high scoring children in each ward. Although the credible intervals for the other covariates include the null value of one, their inclusion influences the ward variance and removing them did not improve model fit, so they were retained in the model.

Disease Mapping Model Results
With adjustment for covariates, the mean ward variance fell from 0.026 to 0.018, suggesting that some of the heterogeneity at ward level was due to the demographic composition of the wards. Spatial heterogeneity is demonstrated by the spatial plots in Figure 2. This is most pronounced in ward 3, where RR fell from 21.6% to 9.2%. After adjustment, there was still high certainty (>80%) of increased RR in the south-east (wards 5, 9, and 10) and in the east (wards 2, 3, 19, and 21). With the exception of ward 5, these were mostly wards associated with moderate to high levels of deprivation.  Table 2 shows the results from the unadjusted and adjusted final disease mapping models. Random slopes were removed from the model as they did not improve fit (DIC = 1023 vs. DIC = 1022) compared with the model without random slopes. For each model, an increase in the relative rate (RR) of high scores per cohort was identified, with an estimated yearly rise of between 1.5% and 5.0% after adjustment. Of the demographics, only the proportion of boys in a ward was associated with the RR of high scoring children in each ward. Although the credible intervals for the other covariates include the null value of one, their inclusion influences the ward variance and removing them did not improve model fit, so they were retained in the model. With adjustment for covariates, the mean ward variance fell from 0.026 to 0.018, suggesting that some of the heterogeneity at ward level was due to the demographic composition of the wards. Spatial heterogeneity is demonstrated by the spatial plots in Figure 2. This is most pronounced in ward 3, where RR fell from 21.6% to 9.2%. After adjustment, there was still high certainty (>80%) of increased RR in the south-east (wards 5, 9, and 10) and in the east (wards 2, 3, 19, and 21). With the exception of ward 5, these were mostly wards associated with moderate to high levels of deprivation. CrI-credible interval; DIC-deviance information criterion.

Figure 2.
Relative rate increase in number of children with high total SDQ scores before (A) and after (C) and adjustment and exceedance probability before (B) and after (D) and adjustment for sex, age, cohort, and deprivation.  Table 3 shows the results from the unadjusted and adjusted final multilevel model. The preschool effect improved model fit (DIC = 197,533) compared with the model with only a ward effect (DIC = 198,974). There was no evidence of an interaction between wards and preschools. The random slope did not further improve model fit (DIC = 197,527) and thus was excluded. The average total difficulties score was considerably higher for boys compared with girls by 34.4-39.6%. For every squared monthly deviation in age from the Figure 2. Relative rate increase in number of children with high total SDQ scores before (A) and after (C) and adjustment and exceedance probability before (B) and after (D) and adjustment for sex, age, cohort, and deprivation.  Table 3 shows the results from the unadjusted and adjusted final multilevel model. The preschool effect improved model fit (DIC = 197,533) compared with the model with only a ward effect (DIC = 198,974). There was no evidence of an interaction between wards and preschools. The random slope did not further improve model fit (DIC = 197,527) and thus was excluded. The average total difficulties score was considerably higher for boys compared with girls by 34.4-39.6%. For every squared monthly deviation in age from the mean of 59 months, RR increased by 0.3-0.4%. Compared with those in the least deprived group, the average total difficulties score consistently increased with rising deprivation. The largest increase was found for the most deprived group, with an average score of 19.3-29.3% higher than the least deprived. There were no interactions between covariates. The credible intervals for the temporal trend showed a 0.3-1.2% increase in total difficulties on average across the years. There was more variation between preschools than between wards. Adding covariates to the model had a marginal effect on variance estimates. Figure 3 shows the RR associated with living in each ward before and after adjustment for covariates. Before adjustment for covariates, only two wards had a RR greater than 1 with a high certainty (>80%) in the north-east and south-west (wards 5 and 16). After adjustment, there was at least 80% probability that the RR was greater than one for four wards in the centre, south, south-west, and north-east of the city (wards 1, 5, 16, and 18). Relative rate increase in average total difficulties scores before (A) and after (C) and adjustment and exceedance probability before (B) and after (D) and adjustment for sex, age, cohort, and deprivation.

Sensitivity Analyses
Because of the changing characteristics of included preschools over the years of the study (shown in Supplementary Table S1), we conducted a sensitivity analysis to investigate whether the temporal trend held for the nurseries that took part in the study every year. Among these preschools, there was still an increasing linear trend among in the disease mapping (RR = 1.054 (1.032-1.076)) and multilevel model (RR = 1.0014 (1.020-1.009)).
Furthermore, we investigated the impact of the choice of disease mapping outcome on the model. The disease mapping approach was repeated using the median difficulties

Sensitivity Analyses
Because of the changing characteristics of included preschools over the years of the study (shown in Supplementary Table S1), we conducted a sensitivity analysis to investigate whether the temporal trend held for the nurseries that took part in the study every year. Among these preschools, there was still an increasing linear trend among in the disease mapping (RR = 1.054 (1.032-1.076)) and multilevel model (RR = 1.0014 (1.020-1.009)).
Furthermore, we investigated the impact of the choice of disease mapping outcome on the model. The disease mapping approach was repeated using the median difficulties per ward as the outcome. The variation between wards was estimated at zero (0.000 (0.000-0.074)), meaning there was no heterogeneity in the median difficulties per ward.

Discussion
Our study aimed to identify, map, and model geographic trends in social, emotional, and behavioural difficulty scores of the SDQ in early childhood using differing methodological approaches. A key finding was an overall increase in the number of high scoring children (i.e., children with likely mental health difficulties) in each ward over time from the disease mapping model. We also identified a citywide increase over time in the child score, on average, through the multilevel modelling approach. The rate of increase in the average score was lower than that of the proportion of high scores.
Existing evidence of decreasing temporal trends in teacher-rated SDQ scores for UK children is largely outside this age group, limiting comparability [49,50]. Decreasing trends in parent-rated SDQ in 4-12-year-old girls were found in Scotland from 1995-2014 [51]. The contrasting results in this study may point to differences between teacher-raters and parentraters, our focus on sub-national outcomes for a whole population, or other factors outside the scope of this study. Other factors that might be relevant include increasing numbers of children entering early years' provision in Glasgow between 2010 and 2014 [32], changes in approaches to the management of council-run nurseries, and implementation of the Scottish Government's Early Years Framework, which began around 2010 [52]. Internationally, there have been varying trends in emotional and conduct problems in young children that may reflect cross-cultural or methodological differences [53].
The geographic pattern differed between the approaches. For example, ward 18 (which has low levels of deprivation) was associated with increased average SDQ score in the multilevel model and a low rate of high scoring children in the disease mapping model. The opposite effect was found in ward 9, where there were high levels of deprivation. A single ward was identified with certainty in both modelling approaches as having a level of difficulties in excess of the city as a whole (ward 5, where most children are in the middle quintile for deprivation). There may be unmeasured covariates that would help explain this pattern [54].
Both approaches found minimal variation between wards, as shown by the small variance terms. It could be that wards were not granular enough to capture neighbourhood variability. Despite this, a small number of wards were found to be up to 30% above average for high scoring children and up to 10% higher for average scores. There was no evidence of spatial correlation between neighbouring wards in either model. This implies that an approach focusing on clustered regions would not be appropriate in this case, as discussed by Peel et al. for Canberra, Australia [19,55]. We did not find evidence to support the use of random slopes, meaning the linear temporal trends were consistent across all wards in the city. However, in the multilevel model, wards where children had increased scores on average differed from the geographic pattern of the 2010-2012 subset of the data [21]. It may be that covering a longer time period provided a more robust estimate in this paper.
The multilevel model approach meant the role of preschool and ward could be included. Preschool variance was larger than ward variance, reflecting variation in the measurement of SDQ between the preschools and the need for whole-preschool approaches to support the children in greatest need ahead of school start. This further supports the importance of preschool context in understanding the variation in child mental health outcomes and developing place-based interventions that are informed by preschool-level outcomes, or delivery interventions via preschools [56].
A limitation of the disease mapping model is the fact that characteristics were aggregated by ward. Defining the variables in this way is likely to dilute any existing variations, giving less power to find an effect. In line with prior knowledge, biological sex was identified as important in both approaches. This may reflect differences in how teachers perceive social, emotional, and behavioural problems in young girls and boys [57]. The multilevel model found an association with age and socioeconomic deprivation. This demographic trend is consistent with other settings [21,22,58]. Of the residual variation observed, there were some wards with high deprivation (e.g., wards 20) that were better than expected. This may be partly due to an effect described by Kershaw et al. [59] in Canada as "off-diagonal" places, where there is a buffering process at the neighbourhood level that reduces the effect of deprivation on development. Conversely, wards 1 and 5 had higher RRs than expected considering their low levels of deprivation (under 10% of children in these wards are in the most deprived quintile). Off-diagonal effects point to the potential presence of other neighbourhood level factors that contribute to spatial variation.
Although the outcomes in each approach are not directly comparable, they both relate to clinical disorder rates at a population [45]. There may be different mechanisms at play for average scoring children compared with those in the high score category, which likely contributed to the difference between trends in our approaches [60].
While there are several population-level child health surveillance programmes in the UK, the outcomes collected [61] and coverage can vary-with children living in the most deprived areas being the least likely to take part [62]. The ChiME study is a unique resource in its ability to address questions about long-term geographic variation in population mental health in young children, in a UK city. This allows the results of this study to be placed into context with international population studies, contributing to our understanding of spatial patterns in child mental health and the role of demographics.

Conclusions
Our interpretation of the results is that, on average, children's scores have increased marginally over time, while, at the same time and partly as a result, more children are reaching the high score cut-off. Both effects are consistent across electoral wards, supporting a citywide intervention. The variation between wards can nevertheless be used to identify a small number of priority areas.
The current approach to policy, such as the Scottish National Performance Framework, focuses on excessive rates of high scores [63]. According to the framework, if parentrated scores increase by more than 1% for 3 years, performance is worsening at a local level. An alternative explanation may be that there is increased reporting of mental health problems [53]. Using the Scottish Government's guideline, the current study found a worsening performance in the number of high scoring children across 8 years, prompting the need for further action, especially (but not exclusively) in neighbourhoods with high socioeconomic deprivation. Off-diagonal effects show, however, that focusing entirely on deprivation may not always lead to effective place-based interventions [14], which should be directed towards children with the greatest difficulties and should be delivered at preschool-level.
We argue that multilevel modelling of individual scores may be beneficial in identifying priority preschools and areas, warranting further investigation as rising scores transition over the high score cut-off. Therefore, using both approaches together gives deeper insight into the trends in early childhood mental health, which are potentially useful in supporting decision-making.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10. 3390/ijerph191811520/s1. Figure S1: SDQ outcomes by covariates. Figure S2: Median SDQ scores by Electoral Ward and Cohort. Figure S3: Proportion of children with high SDQ scores by Electoral Ward and Cohort. Author Contributions: L.T. and P.W. were involved in the ChiME conception, design, and delivery. L.M., G.C. and M.H. made contributions to the design of the ChiME study. S.J.E.B. provided input on the methodology and analysis of the data. S.O. conducted the analysis and prepared the first draft of the manuscript. All of the authors were involved in the interpretation of results, editing, and reviewing the final manuscript. All authors have read and agreed to the published version of the manuscript.