We used two sets of data to explore the relationship between the delivery addresses of seized articles and population characteristics or demographics of the local area.

First, the department’s Mail and Passenger System (MAPS) database provided data about mail articles seized from 1 January 2008 until 31 December 2011, including a variety of characteristics of the article, the most relevant of which were the intended delivery address, which could be used for geocoding and the nature of the seized goods. The second set of data, provided by the Australian Bureau of Statistics (ABS), contained a range of demographic fields of consideration from the 2011 Australian Census. These were provided at the level of statistical unit 2 or SA2 polygons, one of the main structures within the Australian Statistical Geography Standard framework, a framework designed to enable consistent and comparable publication of statistics. SA2s are medium-sized geographical units intended to represent a community that interacts together socially and economically (Australian Bureau of Statistics 2016a). They are the smallest area for the release of many ABS statistics. There are 2,214 SA2 spatial units within Australia with populations ranging from 3,000 to 25,000. Visual comparisons as well as exploratory analysis at the finer grade of Statistical Area 1 (SA1) suggested that the SA2 was a more appropriate unit of analysis due to the increased number of seizures per unit, while not diluting the key patterns and relationships.

There is a large number of Census fields available for the analysis of individuals in the Australian population, therefore the analysis needed to focus on a smaller number of Census fields for the project to be feasible. Profiles were therefore developed of persons or groups associated with the admission of high risk biosecurity goods into Australia. Risk factors for the profiles were based on known intelligence and anecdotal evidence, provided by internal departmental studies, internal Customs studies and broader anecdotal evidence of links between imports and specific demographic groups (Department of Agriculture 2016b, Maller C, Kancans R and Carr A 2007, Kancans R, Stenekes and Benedictos T 2010). In particular, this anecdotal evidence included a number of assumptions linking demographics and importation of particular commodities, in particular, that higher seizure incidences often occur in areas with large numbers of young people attending university.

Anecdotal and empirical research indicated that:

• culturally significant events may be the catalyst to the importation of high risk goods (plant and animal) through the postal system (e.g. cultural and religious festivals or events);
• culturally significant plants may be the catalyst to the importation of high risk plant material through the postal system;
• communication is a barrier to informing and educating persons of culturally and linguistically diverse (CALD) backgrounds about high risk biosecurity behaviour relating to the importation of plant and animal material through the postal system;
• cohorts of CALD persons may distribute biosecurity risk material due to close relational ties and spatial proximity;
• CALD persons studying in Australia are considered a high risk group; and
• CALD persons engaged in agricultural production could be a high biosecurity risk group due to the proximity to commercial agriculture.

The considerations discussed above led to a choice of ABS Census fields regarded as potential risk factors of individuals for the importation of biosecurity risk items which included:

• language spoken at home
• student attending an educational institution (university or other tertiary institutions)
• full-time or part-time attendance status of educational institution
• age in 5 year groupings
• people employed in industries related to horticultural production, including nursery and floriculture, mushroom and vegetable growing and fruit and nut tree growing.

In addition to these fields, the list of census fields was broadened to include those that had not been identified in the anecdotal evidence, but may still be important. These other fields are related to socio-economic circumstances including household characteristics, income levels and employment; those commonly used to describe a population. Land use within the SA2 was also a desirable field, due to the differing impact of biosecurity risk material for certain land uses, for example, primary production or parkland. This information was available by mesh block, the smallest geographic unit available from the ABS (Australian Bureau of Statistics 2016a). These units were aggregated to provide proportions of each land use type for each SA2.

The full set of census fields considered are included in Table 1. For the purpose of analysis, counts of each level of the census field were considered separately as dummy variables, so there were technically 2854 variables used in the model to represent these 32 fields.

In order to be able to relate the seizure counts and these fields, we counted the number of seized articles that were addressed to locations within each SA2 polygon. This count required a match between the delivery addresses and the SA2 locations, which was performed using the G-NAF (Geographic National Address File) provided by Navigate1. This operation took consignee details addresses from items seized at the border and assigned geospatial parameters to it, namely the latitude and longitude and the ABS SA2 and SA1 areas. Overall, the geocoder worked as expected. Although failure rates were initially high, this was, however, more likely due to poor address data (e.g. special or foreign characters not readable by the G-NAF). This poor data quality was fixed through a thorough cleanse of the data. Once the geocoding was completed, this provided a key point of data that was unique across the MAPS data and the census information, namely the SA1 and SA2 locations. Using this standard data field, we were able to join and analyse the data with more accuracy.

Of the 404,873 seizure records in the period of consideration, 307,627 (76%) could be matched to a unique SA2 area. Failures to match were for a range of reasons, including that the international sender address was incorrectly included as the recipient address, that the address was poorly spelled or that the address was fake. Of the 307,627 matched records, 161,390 (52%) contained risk articles that either lacked appropriate Customs declaration or were misdeclared. These records were used for analysis – correctly declared items that were seized were not included.

These seizure counts per SA2 polygon were then matched with the equivalent demographic fields in order to explore the demographic characteristics of high seizure areas.

In seeking to determine which demographic fields are related to high seizure areas, one challenge is the very large number of possible demographic relationships. Therefore, a robust methodology was required to determine which were most closely related to increased seizures. A flexible model was also needed, as these relationships are not necessarily linear nor consistent across demographic fields. Tree-based models are valuable in exploratory cases such as these, because they allow for complex relationships between predictors and outcomes, including interactions between predictors in the way they affect the outcome (Hastie et al. 2009). We can create a tree that progressively splits each group of observations into two subgroups, based on the level of individual demographics fields. The choice of splits can be optimised automatically in software to best separate those with high and low risk of seizures. Note that the seizure outcomes used in these models was the rate per 100,000 of population, to allow for differences in the population in each SA2. The choices of splits yield useful information about the importance of the demographic fields in determining those areas with high levels of seizures. Complex interaction patterns are also naturally incorporated into tree models.

It is possible to construct one such tree model for these data but over fitting is a genuine risk. As there is a relatively large number of potential demographic predictors available when fitting the model, it is possible to estimate the risk for the existing data very well, but the model may not be generalisable to other pathways. For example, it may highlight spurious relationships due to only a few unusual observations in this particular data set. One way to mitigate this problem is to create many trees based on random samples of the data and average the results of these, so the results are not sensitive to the specific data available. These kinds of models are called random forests, as they involve many trees generated from random samples of the data. These kinds of models have been shown to have very high performance (Fernández-Delgado et al. 2014) and the number of SA2s makes this approach feasible in this case. The longitude and latitude of the centroid of each SA2 polygon was also included to allow for and assess spatial correlations (Mascaro et al. 2014). The particular random forests used for this analysis involved 5000 trees, with the minimum size at each terminal node set to 5.

All analysis used the open-source statistical environment R (R Core Team 2013), along with contributed packages randomForest (Liaw and Wiener 2002), sp (Pebesma and Bivand 2005), ggmap (Kahle and Wickham 2013), rgdal (Bivand et al. 2014), gstat (Pebesma 2004) and ggplot2 (Wickham 2009).

Total seizure numbers for the period are reported for each of 12 broad commodity categories and counts of these seizures are given in Tables 2 and 3, separated by state and region. Region in this case is defined using the ABS remoteness structures, which are based on a measure of relative access to services (Australian Bureau of Statistics 2016b). In order to identify patterns of dependence between the address locations and category (that is, whether certain states/regions are associated with more or less of certain seizure types) the cell deviations that contribute to a χ2-test were considered (Greenwood and Nikulin 1996). A cell deviation less than -2 or greater than 2 was considered evidence for a lack of independence; here we observe some deviations of 10 times this amount which would be highly unexpected if the categories and locations were truly independent. Based on the magnitude of these cell deviations, given in brackets, there were proportionally fewer Animal Products destined for Victoria (-28.9) compared with New South Wales (NSW) (10.1) and Queensland (29.3) and greater amounts of Plant/Plant Products (15.0), also compared with these states (-10.0, -12.8). Overall, there were greater amounts of Human Therapeutics destined for Victoria (21.9), Vegetable/Vegetable Products destined for NSW (8.5) and Contaminated Goods/Footwear/Packaging destined for Western Australia (16.1). There are increased numbers of Plant/Plant Product seizures destined for regional areas (19.4) compared with metropolitan areas (-8.6) and increased numbers of Fruit/Fruit Product, Vegetable/Vegetable Products and Mushroom/Fungi seizures destined for metropolitan areas (4.4, 4.6, 5.2) compared to regional areas (-9.4, -10.4, -10.8).

As the random forest approach involves the averaging of many tree models, no single tree or diagram can be used to display the relationship between the demographic data and seizure rates. However, it is possible to consider the effect of each key demographic field, averaged over the other fields, using overall measures of importance. These are a summary of impact of the splits based on each field. These values are only relative and have been standardised such that the largest value is set to 100. This gives an indication of the relative importance of each demographic field in the prediction of seizures, as compared with the relative importance values for the other demographic fields reported for the same model. Table 4 gives the relative importance values for the top 10 fields in a model based on all the available data. In each case the relationship was positive, with increasing seizure rate as the demographic field increased. Note that some language and ancestry fields have been deliberately obscured: Country A refers to an East Asian country for which Language A is the official language. These have been obscured in order to avoid potential labelling of social groups in society on the basis of race or ethnicity, gender, age etc. The point of the study was not to profile groups pejoratively, but to explore a method that can be used to assist in designing campaigns to raise awareness.

Overall, no household income, some relationship to Country A, recent arrival and university study dominate this list of important fields. These fields are also highly related, as international students are often recent arrivals from Country A without any household income. However, it is worth emphasising that the model did not consider counts of individuals with all of these characteristics at the SA2 level, only the counts of each within each SA2.

Given this relationship with students, it is worth considering the destinations of seizures in relation to universities. Figure 1 represents the seizure rates for each SA2 for the two capital cities with the most seizures, Sydney and Melbourne. In many cases the red hotspots are close to the locations of key university campuses.

Figure 1.

Seizure rate patterns by SA2 polygon for Sydney and Melbourne, per 100,000 people from 2008-2011. Green squares correspond to locations of university campuses.

As the units of study (SA2s) were spatially related, it is appropriate to comment on any spatial patterns observed. While considered in the model to assess spatial correlation, the longitude and latitude of the SA2 polygons did not feature highly in variable importance, with relative importance values of 0.4 and 0.2, respectively. This indicates that those SA2s with similar location did not, in fact, have similar seizures frequencies, after adjustment for demographic fields. A variogram of the residuals of the random forest model has been provided in the suppl. material 1 and also shows no evidence of spatial correlation.

This could be due, in part, to the spatial patterns in the seizure counts being strongly related to that of the predictors in the model, leaving no detectable residual spatial correlation. A simple linear regression analysis between seizure rate and the latitude and longitude of the SA2 indicated an association between seizure rates and longitude (but not latitude), suggesting this might be happening in part at least.

While the overall seizure rate is of primary concern, different kinds of seizures pose different kinds of risks and require different mitigation policies in response. As indicated in Tables 2 and 3, there are 12 broad commodity categories into which seizures can be classed. Excluding live animals due to small counts, separate random forests were constructed for the seizures in each category, to compare the important demographic fields. The methods used were otherwise identical to those used for the analysis of the overall seizure rate, with the same relative variable importance measures available for interpretation.

Fields relating to no household income featured highly for all commodity categories, with university study and recent arrival also dominating most categories. East Asian Language A and/or some Country A ancestry were related to all the food categories (Fruit products, Legumes, Herbs, Mushrooms and Vegetable products), whereas another East Asian country was related to the destinations of animal products, an Eastern European ancestry and language were related to biologicals (this general class includes animal or microbial derived products such as foods, therapeutics, laboratory materials and vaccines) and language and ancestry from the Horn of Africa were related to human therapeutics. Seizures due to contaminated materials or soil/mineral samples were conspicuous for having no language or ancestry fields identified as important, which may be to do with the broad nature of these categories, with the biosecurity risk typically unrelated to the imported item itself and, instead, related to the presence of additional contaminants. Contamination, soil/mineral presence and plant products were all related to attendance of a university or other tertiary institution, with plant products also associated with those with high English language proficiency in particular.