Rolnik szuka żony
W 2023 roku, w artykule zatytułowanym „Migration and Skewed Subnational Sex Ratios among Young Adults”, Michał Gulczyński pokazał, że iloraz płci (inaczej wskaźnik maskulinizacji, czyli stosunek liczby mężczyzn do liczby kobiet) w Europie zależy od gęstości zaludnienia na danym obszarze. W regionach o większej gęstości jest, średnio rzecz biorąc, relatywnie więcej kobiet. Wynika to ze zwiększonej (w porównaniu do mężczyzn) migracji kobiet, w uproszczeniu: ze wsi do miasta.
Oto wykres z tego artykułu:
Wykres pokazuje, że o ile „sex ratio” (czyli relacja liczby mężczyzn do liczby kobiet, współczynnik maskulinizacji) w grupie dzieci 0–9 lat nie zależy od gęstości zaludnienia regionu, to dla młodych dorosłych (25–34 lata) widać wyraźną zależność pomiędzy gęstością zaludnienia (pokazaną na skali logarytmicznej) a tym wskaźnikiem.
Spróbujmy odtworzyć taki wykres dla Polski, w R.
W tym celu skorzystamy z Banku Danych Lokalnych (BDL) GUS i pakietu bdl zbudowanego dla R. Potrzebujemy informacji o gęstości zaludnienia oraz liczbie ludności w przedziałach 0–9 oraz 25–34 dla obu płci we wszystkich regionach Polski.
Dane możemy pobrać dla różnych poziomów agregacji: są one dostępne dla makroregionów (poziom 1), województw (poziom 2), regionów NUTS2 (poziom 3), podregionów NUTS3 (poziom 4), powiatów (poziom 5) czy gmin (poziom 6). Wykres przygotujemy dla powiatów dla roku 2024.
library(bdl)
lvl <- 5
yr <- 2024
Pakiet bdl wymaga podania numerów Id zmiennych. Gęstość zaludnienia ma Id 60559, zaś numery zmiennych zawierających liczbę ludności w przedziałach wiekowych osobno dla płci możemy znaleźć tutaj.
Kliknij, by zobaczyć kod
pop_density <- get_data_by_variable("60559", unitLevel = lvl, year = yr)
mal <- get_data_by_variable(c("72301","72302", "47736", "47724"),
unitLevel = lvl, year = yr)
fem <- get_data_by_variable(c("72296","72297", "47696", "47695"),
unitLevel = lvl, year = yr)
Ponieważ przedziały wiekowe w danych GUS są pięcioletnie, dokonujemy odpowiednich agregacji.
Kliknij, by zobaczyć kod
mal$m09 <- mal$val_72301 + mal$val_72302
mal$m2534 <- mal$val_47736 + mal$val_47724
fem$f09 <- fem$val_72296 + fem$val_72297
fem$f2534 <- fem$val_47696 + fem$val_47695
Przekształcamy dane, aby przygotować wykresy. Poniższy kod wykorzystuje funkcje z pakietów dplyr i tidyr: wybiera odpowiednie wartości z wcześniej utworzonych tabel, łączy je ze sobą, oblicza iloraz płci dla dwóch grup wiekowych, wreszcie przekształca tabelę z formatu szerokiego na długi – tak jak wymaga tego pakiet ggplot2.
Kliknij, by zobaczyć kod
library(dplyr)
library(tidyr)
pop_density %>%
select (c(id, name, pop_density=val)) %>%
filter(substring(id,12,12)<"4") %>%
merge(select(mal, c(id, m09, m2534))) %>%
merge(select(fem, c(id, f09, f2534))) %>%
mutate(sexratio09 = m09/f09, sexratio2534 = m2534/f2534) %>%
select(id, name, pop_density, sexratio09, sexratio2534) %>%
pivot_longer(4:5, names_to='group', values_to='sex_ratio') %>%
mutate(group = ifelse(group=='sexratio09', '0-9 years', '25-34 years')) -> sexratio
head(sexratio)
## # A tibble: 6 × 5
## id name pop_density group sex_ratio
## <chr> <chr> <dbl> <chr> <dbl>
## 1 011212001000 Powiat bocheński 165. 0-9 years 1.08
## 2 011212001000 Powiat bocheński 165. 25-34 years 1.07
## 3 011212006000 Powiat krakowski 247. 0-9 years 1.06
## 4 011212006000 Powiat krakowski 247. 25-34 years 1.01
## 5 011212008000 Powiat miechowski 68.5 0-9 years 1.03
## 6 011212008000 Powiat miechowski 68.5 25-34 years 1.07
Następnie budujemy dwa modele regresji, po jednym dla każdej grupy wiekowej. Zmienną objaśnianą jest współczynnik maskulinizacji (sex ratio) w danej grupie wiekowej (0–9 lub 25–34 lata), a zmienną objaśniającą – gęstość zaludnienia.
Kliknij, by zobaczyć kod
reg09 <- lm(sex_ratio~I(log10(pop_density)), data=sexratio[sexratio$group=="0-9 years",])
reg2534 <- lm(sex_ratio~I(log10(pop_density)), data=sexratio[sexratio$group=="25-34 years",])
summary(reg09)
##
## Call:
## lm(formula = sex_ratio ~ I(log10(pop_density)), data = sexratio[sexratio$group ==
## "0-9 years", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.092851 -0.015425 -0.000443 0.015577 0.096695
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.051195 0.005755 182.645 <2e-16 ***
## I(log10(pop_density)) 0.001556 0.002640 0.589 0.556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02741 on 378 degrees of freedom
## Multiple R-squared: 0.0009176, Adjusted R-squared: -0.001725
## F-statistic: 0.3472 on 1 and 378 DF, p-value: 0.5561
summary(reg2534)
##
## Call:
## lm(formula = sex_ratio ~ I(log10(pop_density)), data = sexratio[sexratio$group ==
## "25-34 years", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.138142 -0.034130 -0.003395 0.033375 0.141351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.275149 0.010534 121.05 <2e-16 ***
## I(log10(pop_density)) -0.089491 0.004832 -18.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05017 on 378 degrees of freedom
## Multiple R-squared: 0.4757, Adjusted R-squared: 0.4743
## F-statistic: 343 on 1 and 378 DF, p-value: < 2.2e-16
Jesteśmy teraz gotowi do przedstawienia danych na wykresie:
Kliknij, by zobaczyć kod
# Graphics packages
library(ggplot2)
library(plotly)
# Prepare summary labels for regression results
xmin <- min(sexratio$pop_density, na.rm = TRUE)
ymin <- min(sexratio$sex_ratio, na.rm = TRUE)
ymax <- max(sexratio$sex_ratio, na.rm = TRUE)
text_labels <- data.frame(
x = xmin,
y = c(ymin + 0.05*(ymax-ymin), ymin),
label = c(
paste0("Slope 0-9: ", sprintf("%.4f", reg09$coefficients[2]),
", p-value: ", sprintf("%.3f", summary(reg09)$coefficients[2,4])),
paste0("Slope 25-34: ", sprintf("%.4f", reg2534$coefficients[2]),
", p-value: ", sprintf("%.3f", summary(reg2534)$coefficients[2,4]))))
# ggplot chart
ggplot(sexratio, aes(x=pop_density, y=sex_ratio, col=group, label=name)) +
scale_x_continuous(trans='log10') +
geom_point(alpha=.4) +
stat_smooth(method='lm', se=FALSE) +
scale_color_manual(name='', values=c('blue', 'red'))+
xlab("population density (persons per km²)") +
ylab("sex ratio (# of males/# of females)\nin the age group") +
geom_text(
data = text_labels, aes(x = x, y = y, label = label),
color = c("blue", "red"), inherit.aes = FALSE,
hjust = 0,
vjust = 0,
size = 4
)-> G1
ggplotly(G1) %>%
style(textposition = "top right")
Analogiczne wykresy możemy, wykorzystując ten sam kod, wykonać dla poziomu województw (lvl <- 2):
dla poziomu regionów (lvl <- 3):
dla poziomu podregionów (NUTS3, lvl <- 4):
oraz gmin (lvl <- 6):
Linki do wykresów w większym rozmiarze: województwa, regiony, podregiony, powiaty, gminy.
