---
title: Economic Complexity of Australian Regions
author:
  - name: Hamish Gamble
    orcid: 0000-0003-1277-1227
    corresponding: true
    email: hamish.gamble@flinders.edu.au
    affiliations:
      - Flinders University
  - name: Sharif Rasel
    orcid: 
    corresponding: false
    roles: []
    affiliations:
      - Flinders University
keywords:
  - Economic Complexity
  - Advanced Manufacturing
abstract: |
plain-language-summary: |
key-points:
date: last-modified
bibliography: references.bib
citation:
  container-title: |
number-sections: true
---

## Introduction

## Literature Review

@hidalgo2009 introduced the concept of economic complexity as a means of quantifying and explaining differences in the economic development trajectory of different countries. Their method used bilateral trade data to identify the network structure of countries and the products they export and built on the concept of relatedness introduced in @hidalgo2007. Economic complexity has been shown to be a positive predictor of Gross Domestic Product (GDP), and GDP growth. Increasing economic complexity has also been shown to decrease unemployment and increase employment @adam2023, reduce green house gas emissions @romero2021 and reduce income inequality @hartmann2017.

Relatedness has since been applied across industry [@neffke2012], research areas [@guevara2016], occupation [@muneepeerakul2013] and technology (patents) [@kogler2013].

The relatedness approach has also been used to quantify economic complexity across cities, states, and regions, using employment data[@ecus, @ecnz, @ecmexico], business counts[@ecchina], patent classifications [@balland2021], and interstate and international trade data [@ecaus].

Despite differences in data sources, the method for calculating economic complexity in the literature is relatively standard. The presence of an activity in a region is often identified using a location quotient method, such that an activity is said to be present in a region if:

$$\frac{X_a^r/\sum_{a}X_a^r}{\sum_{r}X_a^r/\sum_{r,a}X_a^r} \geq 1$$

Where $X$ is the measure of an activity $a$ in region $r$ - such as the level of employment in an occupation in a city, or the number of businesses classified in an industry in a province, or the value of exports of a product from a country. The location quotient method creates a binary matrix $M$ with $a$ rows and $r$ columns.

### Regional economic complexity of small areas

The location quotient method can be unreliable due to the discontinuity at 1. This is especially relevant when economic complexity is calculated in regional areas where either $X_a^r$ or $\frac{\sum_{r}X_a^r}{\sum_{r,a}X_a^r}$ are small. In these cases, small changes, or measurement error in $X_a^r$ can significantly change the location quotient.

The choice of region size and activity classification is important. In a study of the economic complexity of US regions, @ecus use metropolitan areas as the basis for calculations. Metropolitan areas in the United States are defined such that jobs within a given area are held by residents who live in that area.Metropolitan areas have a population of at least 50,000 people. The smallest MSA was estimated to have a 2023 population of 57,700 (about 0.015% of US population). They find a poor correlation between ECI calculated at higher level aggregated industry classifications indicating the importance of a high degree of disaggregation to provide as much information to the model as possible @ecus.

In New Zealand, @ecnz use weighted correlations of local employment shares. Regions range from a population of 1,434 to 573,150 with a mean population of 29,947 and median population of 6,952. Employment is measured as an industry-occupation pair.

-   Differences in relationship between complexity and relatedness on indicators may be entirely context dependent.

## Data & Methods {#sec-data-methods}

### Data

-   Calculate economic complexity indicators for Australian regions using employment data from the 2021 Census.

-   Regions classified by Statistical Areas Level 3 (SA3)

-   Economic activity classified by ANZSIC industry division and ANZSCO major group

```{r}
#| label: setup
library(tidyverse)
library(strayr)
library(ecomplexity)
library(sf)
library(tmap)
library(spdep)
library(spatialreg)
library(huxtable)
library(flextable)

industry_occupation_sa3 <- read_complexitydata("sa3_indp1_occp1", export_dir = "data")

regions <- industry_occupation_sa3 |> 
  group_by(x = str_detect(sa3, "POW|Migratory")) |>  
  summarise(count = sum(count)) |> 
  mutate(p = count/sum(count),
         x = ifelse(!x, "included", "excluded")) |> 
  split(f = ~x)

sa3 <- read_absmap("sa32021", export_dir = "data")


industry_occupation_sa3 <- industry_occupation_sa3 |> 
  filter(!str_detect(sa3, "POW|Migratory")) |> 
  mutate(anzsic_division = fct_inorder(anzsic_division),
         anzsco_major = fct_inorder(anzsco_major), 
         industry_occupation = paste0(anzsic_division, " (", anzsco_major,")"),
         industry_occupation = fct_inorder(industry_occupation)) |> 
  group_by(sa3, industry_occupation) |> 
  summarise(count = sum(count),
            .groups = "drop") 
```

-   We exclude individuals who identify their place of work as a Migratory - Offshore - Shipping region or as No Fixed Address. Employment in these regions totals `r scales::comma(regions$excluded$count)` or about `r scales::percent(regions$excluded$p)` of the total sample.

-   Following @ecnz, employment is aggregated into industry-occupation pairs, allowing for differentiation between, for example, managers working in agriculture, forestry, and fishing, and managers working in retail trade.

-   Dataset covers `r length(unique(industry_occupation_sa3$sa3))` regions and `r length(unique(industry_occupation_sa3$industry_occupation))` industry-occupations. @fig-employment-density shows the presence of any level of employment within a region and industry-occupation. As can be seen, there is a high level of employment density across our data, with 88.5% of all combinations of region, industry, and occupation.

```{r}
#| label: fig-employment-density
#| fig-cap: "Presence of employment across regions and industry-occupations."

employment_density <- function(data, region, activity, flip = FALSE) {
  region_size = paste0(region, "_size")
  activity_size = paste0(activity, "_size")
  
  data |> 
    group_by(.data[[region]])  |> 
    mutate({{region_size}} := sum(count)) |> 
    group_by(.data[[activity]]) |> 
    mutate({{activity_size}} := sum(count)) |> 
    ungroup() |> 
    mutate(presence = count > 0,
           {{region}} := fct_infreq(.data[[region]], w = .data[[region_size]]),
           {{activity}} := fct_infreq(.data[[activity]], w = .data[[activity_size]])) |> 
    ggplot(aes(x = if(!flip) .data[[region]] else .data[[activity]],
               y = if(!flip) .data[[activity]] else .data[[region]],
               fill = presence)) + 
    geom_raster() + 
    scale_fill_manual(breaks = c(TRUE, FALSE),
                      values = c("#3182db", "white")) + 
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.border = element_rect(fill = NA),
        legend.position = "none") + 
  coord_equal()
}

employment_density(industry_occupation_sa3, "sa3", "industry_occupation") + 
  labs(x = "Region",
       y = "Industry-Occupation")

```

### Method

This section follows the method of @ecnz using correlations of employment shares rather than a location quotient method.

#### Relatedness

Activities are related based on the weighted correlation between the local activity share of employment, weighted by each regions share of total employment.

-   First calculate the weighted covariance

$$
cov_{aa} = \sum_{c \in C} (\frac{E_c^{a_i}}{E_c}-\frac{E^{a_i}}{E})(\frac{E_c^{a_j}}{E_c}-\frac{E^{a_j}}{E})
$$

-   Divide the weighted covariance by the city share-weighted standard deviations of the local activity shares to get the weighted correlation.

-   Map the correlation to the interval $[0,1]$ such that:

    $$
    r_{aa} = \frac{1}{2}(cor(a_i, a_j) + 1)
    $$

City relatedness is calculated symmetrically such that:

$$
r_{cc} = \frac{1}{2}(cor(c_i, c_j) + 1)
$$

#### Complexity

Activity complexity is defined by the second eigenvector of the matrix $r_{aa}$ and city complexity is defined by the second eigenvector of the matrix $r_{cc}$. The sign of activity complexity is set such that it is positively correlated with the weighted mean size of cities that contain activity $a$, and the sign of city complexity is set such that it is positively correlated with the local share-weighted mean complexity of activities in city $c$

```{r}
#| label: complexity-function
# function to calculate complexity from a region * activity matrix m

complexity_nz <- function(m, base_data) {
  
  activity_share <- m / rowSums(m)
  national_share_employment <- rowSums(m) / sum(m)
  
  region_name <- colnames(base_data)[1]
  activity_name <- colnames(base_data)[2]
  
  # Relatedness of activities is the weighted covariance between the local share vectors for activities i and j, 
  # weighted by each regions share of total employment.
  
  
  r_aa <- 0.5*(cov.wt(x = activity_share, wt = national_share_employment, cor = TRUE)$cor + 1)
  
  
  #Complexity of activity a is the element a of the standardized second eigenvector of the row-standardized relatedness matrix r.
  complexity <- list()
  complexity$activity <- Re(eigen(r_aa/rowSums(r_aa))$vector[,2])
  complexity$activity  <- (complexity$activity  - mean(complexity$activity ))/sd(complexity$activity)
  names(complexity$activity) <- colnames(m)
  
  #Complexity is positively correlated with the weighted mean size of cities that contain activity a
  wmsc <- colSums(m/colSums(m)*rowSums(m))
  
  if (cor(complexity$activity , wmsc) < 0) {
    complexity$activity  = -1*complexity$activity
  } else complexity$activity  = complexity$activity 
  
  message(glue::glue("most complex activity: {names(complexity$activity[complexity$activity == max(complexity$activity)])}")) 
  
  
  # Relatedness of cities is symmetric to activities.
  city_share <- t(m / rowSums(m))
  national_share_activity <- colSums(m)/sum(m)
  r_cc <- 0.5*(cov.wt(x = city_share, wt = national_share_activity, cor = TRUE)$cor + 1)
  
  complexity$city <- Re(eigen(r_cc/rowSums(r_cc))$vector[,2])
  complexity$city <- (complexity$city- mean(complexity$city))/sd(complexity$city)
  names(complexity$city) <- rownames(m)
  
  #City complexity is positively correlated with the local share weighted mean complexity of activities in city c.
  wmc <- rowSums(m/colSums(m)*complexity$activity)
  
  if (cor(complexity$city, wmc) < 0 || complexity$city["Brisbane Inner"] < 0) {
    complexity$city = -1 * complexity$city
    } else complexity$city = complexity$city
  
  message(glue::glue("most complex city: {names(complexity$city[complexity$city == max(complexity$city)])}")) 
  
  df.complexity <- inner_join(base_data, enframe(complexity$city, name = region_name, value = "city_complexity")) |> 
    inner_join(enframe(complexity$activity, name = activity_name, value = "activity_complexity"))
  return(df.complexity)
}

```

## Results

```{r}
#| label: complexity-place-of-work
#|
m <- industry_occupation_sa3 |> 
  pivot_wider(id_cols = sa3,
              names_from = industry_occupation,
              values_from = count) |> 
  column_to_rownames(var = "sa3") |> 
  as.matrix()

complexity <- list()

complexity$indp_occp <- complexity_nz(m, base_data = industry_occupation_sa3)
complexity$indp_occp_rca <- calculate_complexity(industry_occupation_sa3, 
                                                 region = "sa3",
                                                 product = "industry_occupation",
                                                 value = "count", 
                                                 year = 2021) |> 
  select(sa3, industry_occupation, count, city_complexity = country_complexity_index, activity_complexity = product_complexity_index)
  
```

@fig-gcc-complexity shows the regional complexity of SA3 regions in Australian Greater Capital City Areas based on 2021 Census data. Complexity is highest in capital cities and surrounding regions.

```{r}
#| label: fig-gcc-complexity
#| fig-cap: Complexity of Australian Greater Capital City Areas

complexity_map <- function(data, layer) {
  data <- data |> 
    pluck(layer) |> 
    distinct(sa3, city_complexity) |> 
    left_join(sa3, by = c("sa3" = "sa3_name_2021")) |> 
  st_as_sf() 
  
  data |> 
  filter(str_detect(gcc_name_2021, "Australian Capital Territory|Greater")) |> 
  tm_shape() +
  tm_polygons(fill = "city_complexity",
              fill.scale = tm_scale_continuous(values = "greek", ),
              fill.legend = tm_legend(title = "City Complexity Index",
                                      orientation = "landscape")) +
  tm_facets("gcc_name_2021") +
  tm_place_legends_bottom()
  
  }


complexity_map(complexity, "indp_occp")

```

### Regression

$ECI_r = CC + log(employment) + log(population) + log(income_{hh}) + log(business_{entries}) + business_{share}$

```{r}
#| label: city-complexity-ols
load("data/data_by_region.rda")


complexity.reg <- complexity$indp_occp |> 
  distinct(sa3, city_complexity) |> 
  inner_join(data_by_region, by = c("sa3" = "region")) |> 
  inner_join(sa3, by = c("sa3" = "sa3_name_2021")) |> 
  group_by(state_name_2021) |> 
  mutate(business_share = total_businesses_no/sum(total_businesses_no, na.rm = TRUE)) |>
  ungroup() |> 
  mutate(turnover_greater_2m = total_businesses_no - (businesses_with_turnover_of_10m_or_more_no + businesses_with_turnover_of_5m_to_less_than_10m_no + businesses_with_turnover_of_2m_to_less_than_5m_no),
         turnover_greater_2m_share = turnover_greater_2m/total_businesses_no,
         across(c(total_businesses_no, 
                  total_business_entries_no,
                  total_business_exits_no,
                  total_persons_employed_aged_15_years_and_over_no,
                  estimated_resident_population_persons_no,
                  median_equivalised_total_household_income_weekly), 
                .fn = ~log(.x),
                .names = "log_{.col}")) |> 
  drop_na(city_complexity,
          log_total_businesses_no,
          log_total_business_entries_no,
          log_total_business_exits_no,
          turnover_greater_2m_share,
          log_estimated_resident_population_persons_no,
          business_share,
          log_median_equivalised_total_household_income_weekly) |> 
  st_as_sf() 

regression_spec <- parse(text = "city_complexity ~ log_total_businesses_no + log_total_business_entries_no + log_total_business_exits_no + turnover_greater_2m_share + log_total_persons_employed_aged_15_years_and_over_no + log_estimated_resident_population_persons_no + business_share + log_median_equivalised_total_household_income_weekly")[[1]]

fit.ols <- lm(regression_spec,
              data = complexity.reg)


```

### Spatial Correlation

```{r}
#| label: spatial-setup
nb <- poly2nb(complexity.reg, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
```

```{r}
#| label: global-moran


moranI <- moran(complexity.reg$city_complexity, listw = lw, n = length(nb), S0 = Szero(lw))

moranIp <- moran.mc(complexity.reg$city_complexity, lw, nsim = 1000)
```

Based on @fig-gcc-complexity, it looks like there are clusters of complexity, centred around capital cities.

```{r}
#| label: fig-residuals-map
#| fig-cap: Residuals from linear regression

complexity.reg <- complexity.reg |>
  mutate(olsresid = resid(fit.ols)) 

tm_shape(complexity.reg) + 
  tm_polygons(fill = "olsresid",
              fill.scale = tm_scale_continuous(values = "reds"),
              fill.legend = tm_legend(outside = TRUE,
                                      title = "OLS Residuals",
                                      text.size = 0.8),
              col = "grey80",
              lwd = 0.5) +
  tm_layout(frame = FALSE)
```

The residuals from the linear regression are shown in @fig-residuals-map which also shows that the distribution of the residuals appears non-random.

```{r}
#| label: fig-moran-plot
#| fig-cap: Moran Scatterplot for City Complexity in Australia

moran.plot(as.numeric(scale(complexity.reg$city_complexity)),
           listw = lw,
           xlab = "City Complexity",
           ylab = "Neighbours City Complexity")
```

The correlation between complexity and lagged complexity is shown in @fig-moran-plot which also shows a dependency. Finally, we observe a global Moran's I of `r moranI$I` with a p.value of `r round(moranIp$p.value)`. As such, the data appear to be spatially autocorrelated, so a lagged AR or lagged error model should be estimated instead.

```{r}
#| label: city-complexity-spatialreg


fit.lag <- lagsarlm(regression_spec,
                    data = complexity.reg,
                    listw = lw, 
                    zero.policy = TRUE)

fit.err <- errorsarlm(regression_spec,
                      data = complexity.reg,
                      listw = lw, 
                      zero.policy = TRUE)



```

```{r}
#| label: regression-results
huxreg("OLS" = fit.ols, "Spatial AR" = fit.lag, "Spatial Error" = fit.err)

```

## Hot spots

```{r}
#| label: fig-complexity-hot-spots
#| fig-cap: City Complexity hot spots (based on local Moran's I p.values)

MC.i <- localmoran_perm(complexity.reg$city_complexity, lw, nsim = 9999) 

MC.i.df <- MC.i |> 
  as.data.frame() 

complexity.reg$p <- MC.i.df$`Pr(folded) Sim`
complexity.reg$Ii <- hotspot(MC.i, 
                              Prname = "Pr(folded) Sim", 
                              cutoff = 0.05, 
                              p.adjust = "fdr")
complexity.reg$Ii <- factor(complexity.reg$Ii, 
                             levels = c("High-High",
                                        "Low-Low",
                                        "Low-High", 
                                        "High-Low",
                                        ">0.05"))
complexity.reg$Ii[is.na(complexity.reg$Ii)] <- ">0.05"

pal2 <- c( "#FF0000", "#0000FF", "#a7adf9", "#f4ada8","#ededed")

tm_shape(complexity.reg) + 
  tm_polygons(fill = "Ii",
              fill.scale = tm_scale_categorical(values = pal2),
              fill.legend = tm_legend(outside = TRUE, 
                                      text.size = 0.8),
              col = "grey80",
              lwd = 0.5) +
  tm_layout(frame = FALSE) + 
  tm_facets(by = "gcc_name_2021")


```

## Conclusion

## Appendix

Regional economic complexity can be calculated using other data, including employment by industry and employment by occupation.

### Other Employment Indicators

```{r}
#| label: complexity-indp-class
industry_sa3 <- read_complexitydata("sa3_indp4", export_dir = "data") |> 
  mutate(year = 2021) |> 
  filter(!str_detect(indp, "[A-S]|&&&&|@@@@"))
occupation_sa3 <- read_complexitydata("sa3_occp4", export_dir = "data") |> mutate(year = 2021)


complexity$indp4 <- calculate_complexity(industry_sa3, region = "sa3", product = "indp", value = "count", year = 2021)  |> mutate(product = "indp4")
complexity$occp4 <- calculate_complexity(occupation_sa3, region = "sa3", product = "occp", value = "count", year = 2021) |> mutate(product = "occp4")

```

Employment density is much sparser when using 4-digit industry of employment data. Only `r industry_sa3 |> summarise(mean(count >0)) |> pull() |> scales::comma(scale = 100)` of the industry class-region combinations contain any level of employment.

```{r}
#| label: fig-employment-density-classifications
#| fig-cap: "Employment density, SA3 regions: (A) ANZSIC industry class, (B) ANZSCO occupation unit"

library(patchwork)

p1 <- employment_density(industry_sa3, "sa3", "indp") +
  labs(x = "Region",
       y = "Industry Class") 
p2 <- employment_density(occupation_sa3, "sa3", "occp") +
  labs(x = "Region",
       y = "Occupation Unit")

wrap_plots(p1, p2) + 
  plot_annotation(tag_levels = "A")

```

```{r}
#|label: city-complexity-industry

industry_sa3 <- industry_sa3 |> 
  group_by(indp) |> 
  mutate(indp_size = sum(count)) |> 
  filter(indp_size > 0) |> 
  ungroup()

m <- industry_sa3 |> 
  pivot_wider(id_cols = sa3,
              names_from = indp,
              values_from = count) |> 
  column_to_rownames(var = "sa3") |> 
  as.matrix()

complexity$indp4 <- complexity_nz(m, industry_sa3)

m <- occupation_sa3 |> 
  pivot_wider(id_cols = sa3,
              names_from = occp,
              values_from = count) |> 
  column_to_rownames(var = "sa3") |> 
  as.matrix()

complexity$occp4 <- complexity_nz(m, occupation_sa3)

full_join(distinct(complexity$indp4, sa3, city_complexity),
          distinct(complexity$indp_occp, sa3, city_complexity), by = "sa3") |>
  ggplot(aes(x = city_complexity.x, y = city_complexity.y)) + 
  geom_point() + labs(x = "City Complexity (Industry Division)", y = "Economic Complexity (Industry-Occupation)")



```

```{r}
#| label: regression-other-eci

complexity.reg.all <- complexity.reg |> 
  inner_join(distinct(complexity$indp4, sa3, indp_complexity = city_complexity), by = "sa3") |> 
  inner_join(distinct(complexity$occp4, sa3, occp_complexity = city_complexity))

reg_spec <- parse(text = "indp_complexity ~ log_total_businesses_no + log_total_business_entries_no + log_total_business_exits_no + turnover_greater_2m_share + log_total_persons_employed_aged_15_years_and_over_no + log_estimated_resident_population_persons_no + business_share + log_median_equivalised_total_household_income_weekly")[[1]]

fit.ols <- lm(reg_spec, data = complexity.reg.all)

```

### Specific Industries

Following @mealy2022 analysis of green economic complexity, we can apply an advanced manufacturing lens to determine advanced manufacturing complexities.

```{r}
#| label: advanced-manufacturing
#|
am <-  tibble::tribble(
   ~indp,                                                           ~description,
                    1811L,                                         "Industrial gas manufacturing",
                    1812L,                                 "Basic organic chemical manufacturing",
                    1813L,                               "Basic inorganic chemical manufacturing",
                    1821L,                   "Synthetic resin and synthetic rubber manufacturing",
                    1829L,                                    "Other basic polymer manufacturing",
                    1831L,                                             "Fertiliser manufacturing",
                    1832L,                                              "Pesticide manufacturing",
                    1841L,             "Human pharmaceutical and medicinal product manufacturing",
                    1842L,        "Veterinary pharmaceutical and medicinal product manufacturing",
                    1851L,                                      "Cleaning compound manufacturing",
                    1852L,                      "Cosmetic and toiletry preparation manufacturing",
                    1891L,                          "Photographic chemical product manufacturing",
                    1892L,                                              "Explosive manufacturing",
                    1899L,                    "Other basic chemical product manufacturing n.e.c.",
                    2311L,                                          "Motor vehicle manufacturing",
                    2312L,                         "Motor vehicle body and trailer manufacturing",
                    2313L,                        "Automotive electrical component manufacturing",
                    2319L,                              "Other motor vehicle parts manufacturing",
                    2391L,                                     "Shipbuilding and repair services",
                    2392L,                                     "Boatbuilding and repair services",
                    2393L,              "Railway rolling stock manufacturing and repair services",
                    2394L,                           "Aircraft manufacturing and repair services",
                    2399L,                       "Other transport equipment manufacturing n.e.c.",
                    2411L,         "Photographic, optical and ophthalmic equipment manufacturing",
                    2412L,                         "Medical and surgical equipment manufacturing",
                    2419L,            "Other professional and scientific equipment manufacturing",
                    2421L,               "Computer and electronic office equipment manufacturing",
                    2422L,                                "Communication equipment manufacturing",
                    2429L,                             "Other electronic equipment manufacturing",
                    2431L,                                "Electric cable and wire manufacturing",
                    2432L,                            "Electric lighting equipment manufacturing",
                    2439L,                             "Other electrical equipment manufacturing",
                    2441L,                                    "Whiteware appliance manufacturing",
                    2449L,                               "Other domestic appliance manufacturing",
                    2451L,                                    "Pump and compressor manufacturing",
                    2452L, "Fixed space heating, cooling and ventilation equipment manufacturing",
                    2461L,                   "Agricultural machinery and equipment manufacturing",
                    2462L,                      "Mining and construction machinery manufacturing",
                    2463L,                           "Machine tool parts and parts manufacturing",
                    2469L,              "Other specialised machinery and equipment manufacturing",
                    2491L,                         "Lifting and handling equipment manufacturing",
                    2499L,                          "Other machinery and equipment manufacturing"
   ) |>
  mutate(indp = as.character(indp))






```

```{r}
#| label: city-complexity-advanced-manufacturing-complexity

amci <- calculate_complexity(industry_sa3, region = "sa3", product = "indp", value = "count", year = 2021) |>
  mutate(m = as.numeric(rca >= 1),
         is_am = as.numeric(indp %in% am$indp),
         pci = scales::rescale(product_complexity_index, to = c(0, 1)),
         p = m * is_am) |>
  group_by(sa3) |>
  summarise(amci = sum(pci * p),
            amp = sum((1-p)*density*pci)/sum(1-p),
            amce = sum(count*is_am)) |>
  mutate(across(c(amci, amp), ~ (.x - mean(.x))/sd(.x)))

complexity.reg <- complexity.reg |>
  inner_join(amci)


lm.amci <- lm(city_complexity ~ median_age_persons_years + estimated_resident_population_persons_no + total_business_entries_no + total_business_exits_no + business_share + amci + amp, data = complexity.reg)


```

### Smaller Areas

```{r}
#| label: fig-employment-density-sa2
#| fig-cap: Employment density, Industry-occupation, SA2 regions


industry_occupation_sa2 <- read_complexitydata("sa2_indp1_occp1", export_dir = "data")|>
  filter(!str_detect(sa2, "POW|Migratory")) |>
  mutate(anzsic_division = fct_inorder(anzsic_division),
         anzsco_major = fct_inorder(anzsco_major),
         industry_occupation = paste0(anzsic_division, " (", anzsco_major, ")"),
         industry_occupation = fct_inorder(industry_occupation)) |>
  group_by(sa2, industry_occupation) |>
  summarise(count = sum(count), .groups = "drop") |>
  group_by(sa2) |>
  mutate(sa2_size = sum(count)) |>
  group_by(industry_occupation) |>
  mutate(industry_occupation_size = sum(count)) |>
  ungroup() |>
  filter(sa2_size > 0,
         industry_occupation_size > 0)

employment_density(industry_occupation_sa2, "sa2", "industry_occupation", flip = TRUE) +
  labs(x = "Industry-Occupation",
       y = "Region") +
  coord_equal(ratio = 0.25)

```

```{r}
#| label: sa2-indp-occp-complexity
m <- industry_occupation_sa2 |>
  pivot_wider(id_cols = sa2,
              names_from = industry_occupation,
              values_from = count) |>
  column_to_rownames(var = "sa2") |>
  as.matrix()

complexity_small <- complexity_nz(m, industry_occupation_sa2)

complexity_small |>
  distinct(sa2, city_complexity) |>
  inner_join(read_absmap("sa22021", export_dir = "data"),
             by = c("sa2" = "sa2_name_2021")) |>
  st_as_sf() |>
  filter(str_detect(gcc_name_2021, "Australian Capital Territory|Greater")) |>
  tm_shape() +
  tm_polygons(fill = "city_complexity",
              lwd = 0.05,
              fill.scale = tm_scale_continuous(values = 'greek'),
              fill.legend = tm_legend(title = "City Complexity Index",
                                      orientation = "landscape"))+
  tm_facets("gcc_name_2021") +
  tm_place_legends_bottom()
```

## References {.unnumbered}

::: {#refs}
:::
