Small cities face greater impact from automation

The city has proved to be the most successful form of human agglomeration and provides wide employment opportunities for its dwellers. As advances in robotics and artificial intelligence revive concerns about the impact of automation on jobs, a question looms: how will automation affect employment in cities? Here, we provide a comparative picture of the impact of automation across US urban areas. Small cities will undertake greater adjustments, such as worker displacement and job content substitutions. We demonstrate that large cities exhibit increased occupational and skill specialization due to increased abundance of managerial and technical professions. These occupations are not easily automatable, and, thus, reduce the potential impact of automation in large cities. Our results pass several robustness checks including potential errors in the estimation of occupational automation and subsampling of occupations. Our study provides the first empirical law connecting two societal forces: urban agglomeration and automation's impact on employment.


Firm Size Increases with City Size
The U.S. Bureau of Labor Statistics (BLS) uses the annual tax filings of companies to produce a yearly census of those companies. Unfortunately, the data available to the public doesn't include the specific distribution of BLS jobs comprising each firm. Previous work has shown that firm sizes nation wide follow a Zipf distribution [2] indicating that a majority of firms are small, but surprisingly large ones also exist infrequently. Figure 1 shows that the average number of workers per firm increases logarithmically with city size. Larger firms have more capital with which to hire specialized workers along with organization/managerial staff to coordinate those workers. According to our theory, there exists a positive feedback loop where large firms provide demand for specialized workers and cities provide a richer market of skilled workers to meet that demand.
2 Measuring Labor Specialization

Normalized Shannon Entropy
We employ normalized Shannon entropy, as opposed to the standard Shannon entropy definition, to control for size effects on the distributions in cities. For example, it has been shown that the number of different occupations grows with city size (see [1], and SM Fig. 1B), but this result may be due to randomness as more people (e.g. sampled from a long-tailed distribution) are added to a city, and does not account for how workers are distributed amongst these occupations. Large cities may have only a few workers of otherwise absent occupations in small cities, but, perhaps, this distinction does not mean much qualitatively. This motivates us to consider the distribution of workers amongst different occupations, rather than only considering the number of occupations, and to apply relevant information theoretic methods to measure the diversity/specialization of these distributions. We normalize the standard Shannon entropy calculation to control for the number of different occupations in a city, or, equivalently, we normalize Shannon entropy by the maximum possible Shannon entropy given the number of different occupations in the city (i.e. given a number of occupations, Shannon entropy is maximized for the uniform distribution). This normalization is a standard practice for comparing the diversity or information of systems of different sizes. For a summary of normalized entropy, see [3]. In particular, normalized Shannon entropy has been used in a variety of fields, including virology [4], climatology [5], and city science [6]. To understand this normalization, consider that a sufficient number of roles of a fair 6-sided dice and, separately, of a fair 20-sided dice should each produce uniform distributions with maximized Shannon entropy. However, the Shannon entropy of the distribution for the 6-sided dice is − ∑ 1 6 · log( 1 6 ) = 1.79 and the entropy of the distribution for the 20-sided dice is − ∑ 1 20 · log( 1 20 ) = 3.00; specifically, they are not equivalent despite both being discrete uniform distributions because the distributions have a different number of bins. This is analogous to cities having a different number of unique occupations due, potentially, to randomness that occurs with increased city size. We control for this effect by normalizing Shannon entropy according to the maximum possible Shannon entropy given the number of bins in the discrete distribution. Specifically, given a discrete system with N bins (i.e. N-sided dice, or a city with N unique occupations), Shannon entropy is maximized when the distribution is uniform, and the maximum value is given by Therefore, to normalize Shannon entropy according to the maximum possible Shannon entropy, we divide the standard Shannon entropy calculation by log(N) to obtain where p i is the probability of bin i. This normalized Shannon entropy produces a value of 1 for discrete uniform distributions regardless of the number of bins (i.e. regardless of N). In particular, this normalization allows us to control for the number of unique occupations across cities of different sizes to determine the uniformity of job and skill distributions in cities. We present BLS jobs ordered by decreasing skill specialization in Table 6.2. We also provide the scaling exponent of each BLS job, along with the Pearson correlation of the relative abundance of each job to the expected job impact from automation (discussed below) across cities. Figure 2A shows that specialized jobs tend to have larger scaling exponents. Figure 2B shows the distribution of job specialization.

Characterizing Specialization through O * NET Skills
We want to understand how each O * NET skill contributes to the relationships we observe. We present our findings in Table 6.4. First, we compare the raw importance of a skill in each city by summing the raw importance of the skill across each job. We then measure the Pearson correlation of the sum of a given skill compared to the expected job impact of each city (denoted ρ E , second column of table), the skill entropy each city (denoted ρ H , third column of table), and the size of each city (denoted ρ log(N) , right-most column of table). The skills in the Table 6.4 are ordered according to their correlation with expected job impact in cities. For each column, the p-value for the correlation is presented in parentheses. Figure 3 allows us to understand how related each correlation is by taking the Pearson correlation of each ρ we described above. Figure 3A demonstrates that skills which indicate lower expected job impact in cities also indicate greater skills specializations in cities. Figure 3B demonstrates that skills which indicate lower expected job impact in cities also indicate larger cities. Interestingly, Figure 3C demonstrates that skills which indicate skill specialization in  Figure 3: Comparing the relationships of O * NET skills to city size, expected job impact, and labor specialization. (A) We plot Pearson correlation of raw skill importance to expected job impact (ρ E ) on the x-axis versus the Pearson correlation of raw skill importance to city skill entropy (ρ H ) on the y-axis. We see that which indicate job impact from automation also indicate decreased specialization in cities. (B) We plot Pearson correlation of raw skill importance to expected job impact (ρ E ) on the x-axis versus the Pearson correlation of raw skill importance to city size (ρ log(N) ) on the y-axis. We see that which indicate job impact from automation also indicate smaller city sizes. (C) We plot Pearson correlation of raw skill importance to city skill entropy (ρ H ) on the x-axis versus the Pearson correlation of raw skill importance to city size (ρ log(N) ) on the y-axis. The correlation between these two variables is not significant cities are not significantly related to the skills which indicate city size. This finding is surprising given the other panels of the figure, and motivates us to consider the relationship between occupational specialization and city size through the jobs in each city (see main text).

Estimating the Affects of Automation
Automation and its impact on labor are increasingly important topics to researchers [7][8][9]. Examples throughout history, such as the industrial revolution and the advent of computers, demonstrate how technological advancement can lead to both job loss and job creation [10,11]. However, it is extremely difficult to predict how quickly a seemingly imminent technology will reach maturity and what the impact of that technology will be. For example, it's currently topical to discuss self-driving cars, but, while autonomous-capable cars are available for purchase, no self-driving cars are currently operated on the mass market. On the other hand, early leaders in computer hardware famously offered pessimistic predictions on the impact of computing: • "There is no reason anyone would want a computer in their home" -Ken Olsen, founder of Digital Equipment Corporation (1977) • "I think there is a world market for maybe five computers" -Thomas Watson, former president of IBM (1943)

Estimating Automation Impact using Frey/Osborne Data
Frey and Osborne [12] produced probabilities of computerization for each BLS job. They convened a workshop of leaders in automation to identify which of the BLS jobs were certainly automatable and certainly not-automatable. They used the O * NET skills dataset to identify the raw importance of nine workplace skills to each BLS job: Manual Dexterity, Finger Dexterity, Cramped Workspace/Awkward Positions, Originality, Fine Arts, Social Perceptiveness, Negotiation, Persuasion, and Assisting & Caring for Others. These O * NET skills represent "known bottlenecks to computerization." Using the importance of these skills to the jobs whose automatability was clear, they used a Gaussian process classifier to produce a probability of computerization for each BLS job. Frey and Osborne used these probabilities to conclude that 47% of the current U.S. jobs are at "high risk" of computerization. n Several studies [13][14][15] utilize these same probabilities to investigate the impacts of automation, which highlights the utility of the probabilities despite the difficulty of the prediction undertaken in [12]. We use the same probabilities in combination with the distribution of BLS jobs across U.S. cities to add spatial resolution to their findings. For a city, m, the expected job impact from automation is calculated according to where p auto ( j) is the probability of computerization according to [12] and p m ( j) is the proportion of workers in city m with job j. Table 6.1 demonstrates the ordered list of cities according to expected job impact. As mentioned above, it's difficult to validate automation predictions. Nonetheless, our calculations for expected job impact represent an aggregate signal for the types of jobs in a city in relation to imminent automation technology. In the Table 6.1, we present the U.S. cities ordered by their expected job impact from automation. The list produces an ordering that appears to make sense; cities with technology companies and research institutes, such as Boston, M.A., and Boulder, C.O., have the lowest expected job impact, while cities relying on the tourist industry and agriculture, such as Myrtle Beach, S.C., and Napa, C.A., have the highest expected job impact. While the absolute proportions can only be validated with time, we believe the overall trend embodied in expected job impact in cities represents an underlying true signal.
To demonstrate the robustness of our results further, we perform two robustness checks to verify the negative trend between city size and expected job impact from automation (see Fig. 1B from the main text). The probability of computerization (i.e. p auto ( j)) from [12] are produced through a machine learning process applied to predictions of the automatability of jobs from experts. Therefore, we expect some errors in the predictions of these experts, and our task is to demonstrate that the error in the resulting p auto ( j) would need to be substantial to invalidate our finding. We perform this analysis by artificially adding random noise to each p auto ( j) according to where e j is chosen uniformly at random from the interval [−error, +error] for each occupation. For each choice of error, we perform 500 trials calculating a new p * auto ( j) for each occupation and recalculating the expected job impact from automation in each city according to similar to equation 2. Finally, we measure the Pearson correlation between log 10 the total employment in each city and E * m so that we can compare to the empirical relationship we observe in Figure 1B of the main text (Pearson ρ = −0.53, p val < 10 −28 ). Figure 5A demonstrates the results of this exercise. We find that substantial error (error ≈ 0.15) needs to be added to the empirical probabilities of computerization for each occupation before our result from the main text no longer represents the observed trend. Even if we make the extremely strong assumption of error = .5, we would still observe a strong negative trend, and we would still conclude that small cities face greater impact from automation.  [12], we perform 500 trials measuring the resulting Pearson correlation between log 10 city size and expected job impact from automation after the error has been added to each occupation's probability of computerization (y-axis). (B) After selecting a proportion of occupations (x-axis), we perform 500 trials of randomly selecting that many occupations to remove while measuring the resulting Pearson correlation between log 10 city size and the expected job impact from automation in cities (y-axis).
In the second robustness check, we test the robustness of our observed relationship between city size and expected job impact if a randomly selected subset of occupations are removed from the analysis. For each proportion of occupations to be removed, we perform 500 trials of randomly selecting occupations to be ignored and recalculate E m using the p auto ( j) presenting in [12]. We then measure the resulting Pearson correlation between these new E m and log 10 the total employment in each city. Figure 5B demonstrates that our empirical observation from Figure 1B in the main text holds even if very large proportions of occupations are ignored. In fact, only when we ignored half of all occupations did we observe any trials demonstrating a trend contrary to the one presented in the main text. Therefore, we conclude that small cities face greater impact from automation.

Estimating Automation Impact using OECD Data
The Organization for Economic Co-operation Development (OECD) released alternative estimates for the probability of job automation with a focus on job categories used by OECD countries [16]. Rather than the job-based approach used in [12], assessments on the automatability of workplace skills were derived. These skill assessments can be used in combination with government data relating the importance of skills to jobs to assess the likelihood of computerization for jobs. Contrary to the alarming 47% of jobs at "high risk of computerization" found by Frey and Osborne, these new probabilities produce a more mild conclusion of only 9%. These job probabilities were derived with OECD job definitions in mind, but collaborations between OECD and U.S. BLS have lead to an official mapping between the two job definitions. We utilize this mapping to assess the resilience of labor markets in cities as a function of city size in Figure 6. Despite the more conservative estimates in [16], our results remain; we again observe significantly decreased expected job impact in large cities (Fig. 6A).

Expected Job Impact & Labor Specialization in Cities
In Figure 7, we further characterize the relationship between a city's resilience to job impact from automation and labor specialization. We provide additional figures in the main text detailing how workplace skills explain the positive correlation we observe between labor specialization and resilience to job impact in cities. Here, we demonstrate that resilience to job impact is significantly correlated to the number of unique jobs in a city, and more weakly correlated to the Shannon entropy of job distributions. This weaker correlation motivates our investigation into workplace skills, in addition to the distribution of jobs, presented in the main text.

Explaining Differences in Expected Job Impact
From equation 2, we may observe that both the automatability and the employment share of individual occupations contribute to the expected job impact from automation in a city. Correspondingly, we measure the difference in expected job impact for cities m and n according to where we have utilized ∑ E n · (share m ( j) − share n ( j)) = 0. Here, we let Jobs denote the set of all occupation types across all cities, p auto ( j) denotes the probability of computerization of occupation j according to [12], and share m ( j) denotes the employment share of occupation j in city m.
denote the percent influence of occupation j on the difference in expected job impact for cities m and n. Figure 8 demonstrates a visualization of equation 5 that we call an "occupation shift." Correspondingly, if we add the employment distributions in the 50 largest cities and 50 smallest cities together (respectively), then we can quantify how each occupation contributes to the differential impact of automation on employment in large and small cities. We present this occupation shift in Figure 9 (also Figure 5 of the main text).
Referring to the job clusters from Figure 5 in the main text, we see that purple occupations and blue occupations contribute the most to the difference in expected job impact, while green and yellow occupation types effectively diminish the difference in both occupation shifts. However, certain occupations, such occupations of the green job cluster, can both increase and decrease the difference between resilient and susceptible cities. The occupation shift allows us to understand which occupations explain the overall trend and which occupations go against the overall trend. If we had only considered occupations that add to the difference (i.e. occupations corresponding to dark colored bars on the right side of the plot), then we may have incorrectly concluded that the differences in relatively susceptible occupations explain the difference we observe in these two examples. This transparency can help urban policy makers determine how labor shifts in different industries may effect their preparedness for the impact of new technology.

Employment Difference in
Occupations that are: The occupation title is provided next to the corresponding bar and colored according to its job cluster as identified in Figure 4 of the main text. Red bars represent occupations with higher risk of computerization compared to Boston's expected job impact. Blue bars represent occupations with lower risk of computerization compared to Boston's expected job impact. Dark colors represent occupations that effectively increase the difference, while pale colors represent occupations that effectively decrease the difference in expect job impact. Bars in each of the quadrants are vertically ordered according to |δ (Las Vegas,Boston) ( j)|. The inset in the bottom left of the plot summarizes the overall influence of resilient occupations compared to occupations that are at risk of computerization.

Employment Difference in
Occupations that are:

Labor Specialization as a Mediator for City Size and Automation Impact
Our analysis in the main text predominantly relies on linear regression to explore the relationship between city size, labor specialization, and the expected impact of automation in cities (see Figures 1,2,& 3). In particular, we find evidence that the relationship between city size and expected impact may be mediated by the labor specialization in cities. This conceptualization leads us to perform a formal mediation analysis [17] with city size as a treatment variable (i.e. log 10 total employment in cities, denoted size m ), our various measures for labor specialization (i.e. H job (m), H skill (m), and 1 − T m ) as independent mediators, and the expected impact from automation (i.e. E m ) as the outcome variable. We take the generic urban variables from Figure 3 in the main text as additional control variables; these variables include the median household income (income m ), the per capita GDP (GDP m ), the percent of population with a bachelor's degree (bachelor m ), and the number of unique occupations in each city ( jobs m ). The purpose of these control variables are to mitigate omitted variable bias, but the U.S. labor system is a sufficiently complicated system that omitted variable bias can never fully be controlled for. In our opinion, this observation limits the strength of any conclusion about causality [18].
This analysis is subject to further assumptions as well. Firstly, we are assuming that the effect of city size on job impact is constant across cities. Again, this assumption is extremely difficult to prove given the complexity of the U.S. labor system, regional geographies, regional politics, and economic trade. Secondly, the effects of unobserved causes for the mediator and outcome variables (denoted e 1 and e 2 , respectively) are uncorrelated. In the analysis below, we measure the Pearson correlation between e 1 and e 2 for each choice of labor specialization measure. −0.107 * 0.389 0.282 6.01 × 10 −17 * p-value < .1, * * p-value < .01, * * * p-value < .001 Table 1: The results of mediation analysis. All variables were standardized prior to analysis. Ignoring whether the necessary assumptions for mediation analysis are met, we find some evidence that job specialization (H job (m)) may act as a mediator for the effect of city size on the expected impact of automation.  Figure 11: To confirm the validity of the regression models, we perform 1,000 trials where half of the cities are randomly selected without replacement as training data and the remaining cites are used for validation. We under go this process for the regression model using only generic urban indicators (A) and the regression model using all variables (B). The resulting distributions of variance explained (R 2 ) when the trained models are applied to separate validation data confirms that the full regression model accounts for an additional 10% of variance on average. U p d a ti n g a n d U s in g R e le v a n t K n o w le d g e

Simplifying Jobs & Skills
In an effort to clearly identify how jobs contribute to labor specialization in larger cities, we identify aggregate job types based on common workplace skills. Previous studies, such as [1], examined the relationship between industry size and city size for various abstractions for industry according to NAICS. Here, we are seeking an organic representation of the forces in effect, and so we use K-means clustering based on the raw skill values for each job to identify five clusters of similar jobs (i.e. occupations are instances and the raw O * NET importance of each skill are features). These job groups represent collections of jobs which rely on similar skills for completion. The BLS jobs comprising each job type are shown in Section 6.3. Note that our results and interpretations are consistent for anywhere from three to seven clusters (see 6.3.1). This simplification of the space of jobs allows us to clearly understand which job types are disproportionately emphasized in large cities through the scaling behaviors of these job types.
We also seek to explain our results on the basis of workplace skills. To this end, we measure the correlation of raw skill values across all BLS jobs for each pair of O * NET skills and employ K-means clustering to identify ten groups of co-occurring skills (i.e. workplace skills are instances and the Pearson correlation of the raw O * NET importance of that skill to the importance of each other skill are the features). The complete lists of raw O * NET skills comprising each skill type are presented in Section 6.5. We summarize the skills comprising each skill type with the groups' titles. This simplification of the space of skills clarifies how different types of skills explain our results, and trends that we present using these aggregated skill groups are apparent when reproduced using the raw O * NET skills instead.

O * NET Task Groups
An alternative simplification of the raw O * NET skills is the O * NET Task Groups, which represent collections of similar work activities. We provide the definitions for these task groups in Table 3. These task groups have been used to investigate the task connectivity of urban labor markets in relation to employment growth [19]. In Figure 13, we use these groups as alternative skill aggregations and assess which tasks indicate resilience to job displacement from automation in cities (Fig. 13A), which tasks indicate occupational specialization in cities (Fig. 13B), and which tasks indicate superlinear scaling of job types (Fig. 13C). We find that Mental Process tasks are indicative of increased specialization in cities, increased resilience to job displacement in cities, and superlinear scaling of job types. On the other hand, Work Output tasks, which focuses on physical skills, indicate less specialization in cities, less resilience to job displacement in cities, and linear or sublinear scaling of job types. These findings are in agreement with the results in the main text.  Figure 13: The relationships between O * NET tasks expected job impact from automation, labor specialization, and the scaling of job types. (A) We bin cities according to their expected job impact from automation (x-axis). For each task (legend), we normalize the importance of that task across bins to a probability P(E m | task) representing how strongly that task indicates each level of job displacement (y-axis). (B) We bin cities according to their skill specialization (x-axis) and sum the importance of each task for each bin. For each task (leg-end), we normalize the importance of that task across bins to a probability P(H skill (m) | task) representing how strongly that task indicates each level of specialization (y-axis). (C) By summing the importance of each task to each job type, we assess how strongly a task indicates a scaling relationship according to its z score. For a given task, z scores are calculated according to the distribution of importance across job clusters.  Table 2: Summarizing the relationship between tasks, job impact from automation, and city size. In the middle (right) column, we present the Pearson correlation of the proportion of each task to the expected job impact (log 10 city size). We provide the associated p-values in parentheses

The Routineness of Tasks
Autor et al. [20,21] identify workplace tasks according to their type and how routine the task is. They find that nonroutine tasks are becoming increasingly important to workers relative to routine tasks. We provide the definitions for these task groups in Table 5. In Figure 14, we use these groups as alternative skill aggregations and assess which tasks indicate resilience to job impact from automation in cities (Fig. 14A), which tasks indicate occupational specialization  in cities (Fig. 14B), and which tasks indicate superlinear scaling of job types (Fig. 14C). We find that all non-routine tasks are indicative of increased specialization in cities and increased resilience to job impact in cities. Non-routine analytic tasks and non-routine interactive tasks are indicative of superlinear scaling of job types, while non-routine manual tasks indicate linear or sublinear scaling of job types. Routine tasks indicate less specialization in cities, less resilience to job impact in cities, and linear or sublinear scaling of job types. These findings are in agreement with the results in the main text.  Table 4: Summarizing the relationship between tasks, job impact, and city size. In the middle (right) column, we present the Pearson correlation of the proportion of each task to the expected job impact (log 10 city size). We provide the associated p-values in parentheses  Figure 14: The relationships between O * NET tasks expected job impact from automation, labor specialization, and the scaling of job types. (A) We bin cities according to their expected job impact from automation (x-axis). For each task (legend), we normalize the importance of that task across bins to a probability P(E m | task) representing how strongly that task indicates each level of job displacement (y-axis). (B) We bin cities according to their skill specialization (x-axis) and sum the importance of each task for each bin. For each task (leg-end), we normalize the importance of that task across bins to a probability P(H skill (m) | task) representing how strongly that task indicates each level of specialization (y-axis). (C) By summing the importance of each task to each job type, we assess how strongly a task indicates a scaling relationship according to its z score. For a given task, z scores are calculated according to the distribution of importance across job clusters.

Relating City Trends to BLS Jobs
We present BLS jobs ordered by decreasing skill specialization in Table 6.2. We also provide the scaling exponent of each BLS job, along with the Pearson correlation of the relative abundance of each job to the expected job impact from automation (discussed below) across cities. p-values for the correlations are presented in parentheses.

Job Groups
The O * NET skills database allows us to identify how important each of 230 workplace skills is to completing each of the BLS jobs. We use K-means clustering to group jobs into five groups according to the skills required to perform those jobs. The complete list of BLS jobs comprising each job group is presented in the table below. Our interpretation about the scaling behaviors of jobs, and how aggregate skills indicate those scaling behaviors, is the same if we use anywhere between three and seven job groups instead of five while computing the K-means clustering algorithm.

Alternative Job Groups using K-means
We demonstrate that our choice to focus on five groups of jobs according to skills produces results that are consistent for several alternative numbers of groups. Using K-means to identify between three and seven job groups continues demonstrate that computational/analytical and managerial skill are more indicative of super linear job growth, while physical skills are more indicative of linear or sub linear job growth. Likewise, our conclusions relating job scaling to expected job impact by comparing skills hold as well.
K-means clustering of similar jobs (k = 3)

Stability Testing for Job Groups
We want to test the stability of the scaling results we observe when using five job clusters obtained from k-means clustering. In particular, how robust to sub-sampling is our observation that one job clusters scales faster than the rest? For a single trial, we sub-sample from the complete list of BLS occupations (percent indicated in plot titles) to obtain a matrix where each row represents a single occupation which was sub-sampled and each column represents the raw O*NET importance of a skill to each occupation. We apply k-means to this occupation-skill matrix (i.e. occupations are instances and skills are features) to obtain five occupation clusters (note: examination of prescribing between three and seven clusters is discussed in the SM). We then measure the scaling exponent (β) of each occupation cluster and rank the occupation clusters according to scaling exponent (rank indicated by color in plots). We perform 100 independent trials for each sub-sampling proportion in {10%, 20%, . . . , 90%, 100%} and plot the resulting scaling exponent distributions. In good agreement with our original findings, we find that indeed one occupation cluster (indicated in purple) tends to grow much faster than the other occupation clusters despite varying sub-sampling of occupations.   Figure 15: Boot-strapping at various rates of sub-sampling demonstrates the stability of our result that one job cluster scales at a greater rate than the rest when using five clusters obtained from k-means clustering.

Checking the Statistical Robustness of Job Group Scaling
Readers who are familiar with the urban scaling literature may be aware of an ongoing debate about the statistical significance of exponent measurements and identification of underlying statistical models to explain that growth.
For example, what model should one assume to test if a trend is significantly superlinear? Rather than solving this ongoing and important problem, the goal of this study is only to understand the relationship between automation and urbanization. Our narrative requires only that highly specialized occupations (represented by purple dots in Figure 3A of the main text) exhibit superlinear growth and be notably different from the growth exhibited by other occupations. Recent work by Leitao et al. [22] proposes several statistical models that may explain urban scaling trends, and they apply them to a variety of datasets to test the models' ability to explain urban scaling. Here, we employ these same models to test if our requirements on the scaling of highly specialized occupations are met according to the five job groups discussed in the main text (i.e. K-means clustering with k = 5). As an example, Figure 16 provides estimates of the scaling exponent along with standard errors for the scaling of each job group according to the unconstrained logarithm model. Table 9 details the complete analysis in line with the methods in [22]. For each model tested, we find that the purple job group, which represents highly specialized occupations, exhibits significantly superlinear scaling, and, furthermore, consistently exhibits faster growth rates than other job groups.