Economia Regional: Lei de Zipf em R

Abstract

This is an undergrad student level instruction for class use.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Economia Regional: Lei de Zipf em R. Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em http://rpubs.com/amrofi/Regional_Economics_zipfms.

1 Introdução

As hierarquias de centros comerciais são refletidas na complexidade de suas funções e em seu tamanho em relação a outros centros, conforme Christaller (1933) e Losch (1940). Para construir uma “hierarquia de tamanho urbano”, os centros urbanos devem ser classificados por tamanho de população.

O esforço mais conhecido para criar tal hierarquia é a regra de tamanho de rank desenvolvida por GK Zipf em 1949, citado por Ruiz (2005).

Segundo Ruiz (2005, p.716), na economia urbana, quando as cidades são ordenadas de forma decrescente a partir de sua população, há uma relação peculiar entre população e o ranking das cidades.

2 Método

Partindo da relação entre o rank (\(R_i\)) e o tamanho da população (\(N_i\)) da cidade \(i\), é possível escrever uma relação não linear tal que

\[ R_i = \frac {\alpha}{{N_i}^\beta} \]

Observando as estimativas de \(\beta\) em valores absolutos, pode-se inferir sobre a simetria e a polarização da estrutura urbana. Para o caso de \(\beta = 1\) existirá uma distribuição do tipo Pareto e a validade da Lei de Zipf.

Uma estrutura urbana com grande assimetria entre cidades (tamanhos muito diferentes) e poucas cidades grandes (indicando maior polarização), devem apresentar estimativas de \(\beta\) menores que 1 em valor absoluto (\(\lvert {\beta} \rvert< 1\)). Observar que na saída da estimação deve aparecer \(\beta\) negativo.

De outro lado, valores de \(\beta\) maiores que 1 em valor absoluto (\(\lvert {\beta} \rvert> 1\)) indicarão vários grandes centros urbanos (e portanto menor polarização ou uma estrutura descentralizada) e simetria nos tamanhos das cidades (portes similares).

2.1 Dados básicos

Neste exemplo, utilizarei os dados populacionais dos \(i\) municípios de Mato Grosso do Sul (MS) em 2015. Os passos envolvem a ordenação (ranqueamento) dos municípios por população e a regressão é feita no formato log-log, portanto:

\[ R_i = \frac {\alpha}{{N_i}^\beta} \]

e,

\[ \ln{R_i} = \ln{\alpha} - \beta \ln{N_i} \]

em que: \(N_i\) é a população de cada município \(i\) em 2015, e \(R_i\) é a posição no ranqueamento de cada município. Os parâmetros \(\alpha\) e \(\beta\) são estimados por mínimos quadrados ordinários no modelo de regressão.

Conforme Rosen e Resnick (1980, p.165), a regra do tamanho-rank implica em um expoente de Pareto unitário, ou seja, \(\beta = 1\), na também chamada rank-size rule.

A ordenação municipal em Mato Grosso do Sul gerará algo como:

1. Campo Grande (rank = 1);
2. Dourados (rank = 2);
3. Três Lagoas (rank = 3);
4. Corumbá (rank = 4);
5. … .
   .
   .
   
6. menor população!

Vamos olhar os dados:

# script de Adriano M R Figueiredo, ESAN-UFMS
# adaptado para estudar a primazia e beta especializacao em economia regional
# variaveis do dataset, a partir do shapefile do IBGE e adaptado
# 79 municipios de Mato Grosso do Sul
# ordem - identificador do municipio
# codmun - codigo do municipio no IBGE 7 digitos
# nomemun - nome do municipio
# pop - populacao municipal estimada pelo IBGE para o TCU de 2002 a 2015 e rank para 2015
#
# Primazia das S maiores
# Ps = N1/sum(N1+N2+N3+N4+N5)
# s = limite do tamanho da populacao - 5, 10, 30 maiores (conforme Ruiz, 2005)
# ou por pop> 20000, 50000, 100000 (conforme Aveni et al, 2015)
# são 79 municipios ao todo

dados.ord2 <- structure(list(ordem = c(20, 32, 78, 26, 62, 53, 73, 56, 9, 50, 59, 
    3, 65, 28, 18, 51, 45, 70, 4, 8, 14, 64, 42, 40, 23, 48, 12, 22, 16, 41, 77, 
    55, 27, 67, 34, 52, 74, 63, 38, 25, 1, 54, 19, 60, 29, 33, 17, 10, 13, 75, 30, 
    72, 6, 37, 36, 5, 7, 43, 15, 61, 39, 69, 49, 44, 11, 71, 47, 79, 21, 31, 24, 
    68, 58, 2, 66, 57, 46, 76, 35), codmun = c("5002704", "5003702", "5008305", "5003207", 
    "5006606", "5005707", "5007901", "5006200", "5001102", "5005400", "5006309", 
    "5000609", "5007208", "5003306", "5002407", "5005608", "5005004", "5007695", 
    "5000708", "5001003", "5002100", "5007109", "5004700", "5004502", "5002951", 
    "5005202", "5001904", "5002902", "5002209", "5004601", "5008008", "5006002", 
    "5003256", "5007406", "5003801", "5005681", "5007935", "5006903", "5004304", 
    "5003157", "5000203", "5005806", "5002605", "5006358", "5003454", "5003751", 
    "5002308", "5001243", "5002001", "5007950", "5003488", "5007703", "5000856", 
    "5004106", "5004007", "5000807", "5000906", "5004809", "5002159", "5006408", 
    "5004403", "5007554", "5005251", "5004908", "5001508", "5007802", "5005152", 
    "5008404", "5002803", "5003504", "5003108", "5007505", "5006275", "5000252", 
    "5007307", "5006259", "5005103", "5007976", "5003900"), nomemun = c("Campo Grande - MS", 
    "Dourados - MS", "Três Lagoas - MS", "Corumbá - MS", "Ponta Porã - MS", "Naviraí - MS", 
    "Sidrolândia - MS", "Nova Andradina - MS", "Aquidauana - MS", "Maracaju - MS", 
    "Paranaíba - MS", "Amambai - MS", "Rio Brilhante - MS", "Coxim - MS", "Caarapó - MS", 
    "Miranda - MS", "Jardim - MS", "São Gabriel do Oeste - MS", "Anastácio - MS", 
    "Aparecida do Taboado - MS", "Bela Vista - MS", "Ribas do Rio Pardo - MS", "Ivinhema - MS", 
    "Itaporã - MS", "Chapadão do Sul - MS", "Ladário - MS", "Bataguassu - MS", "Cassilândia - MS", 
    "Bonito - MS", "Itaquiraí - MS", "Terenos - MS", "Nova Alvorada do Sul - MS", 
    "Costa Rica - MS", "Rio Verde de Mato Grosso - MS", "Fátima do Sul - MS", "Mundo Novo - MS", 
    "Sonora - MS", "Porto Murtinho - MS", "Iguatemi - MS", "Coronel Sapucaia - MS", 
    "Água Clara - MS", "Nioaque - MS", "Camapuã - MS", "Paranhos - MS", "Deodápolis - MS", 
    "Eldorado - MS", "Brasilândia - MS", "Aral Moreira - MS", "Batayporã - MS", "Tacuru - MS", 
    "Dois Irmãos do Buriti - MS", "Sete Quedas - MS", "Angélica - MS", "Guia Lopes da Laguna - MS", 
    "Glória de Dourados - MS", "Anaurilândia - MS", "Antônio João - MS", "Japorã - MS", 
    "Bodoquena - MS", "Pedro Gomes - MS", "Inocência - MS", "Santa Rita do Pardo - MS", 
    "Laguna Carapã - MS", "Jaraguari - MS", "Bandeirantes - MS", "Selvíria - MS", 
    "Juti - MS", "Vicentina - MS", "Caracol - MS", "Douradina - MS", "Corguinho - MS", 
    "Rochedo - MS", "Paraíso das Águas - MS", "Alcinópolis - MS", "Rio Negro - MS", 
    "Novo Horizonte do Sul - MS", "Jateí - MS", "Taquarussu - MS", "Figueirão - MS"), 
    pop2015 = c(853622, 212870, 113619, 108656, 86717, 51535, 51355, 50893, 47162, 
        43078, 41495, 37590, 34776, 33139, 28437, 27104, 25473, 24982, 24748, 24414, 
        24113, 23167, 22928, 22896, 22620, 21860, 21775, 21622, 21047, 20162, 19914, 
        19656, 19508, 19462, 19220, 17884, 17483, 16514, 15637, 14815, 14474, 14233, 
        13731, 13494, 12650, 12128, 11903, 11399, 11208, 11035, 10965, 10832, 10149, 
        10136, 9992, 8844, 8679, 8567, 7898, 7794, 7664, 7633, 7017, 6860, 6771, 
        6455, 6399, 6027, 5838, 5723, 5513, 5252, 5150, 5038, 4910, 4306, 4038, 3570, 
        3012), rank = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
        18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 
        37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 
        56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 
        75, 76, 77, 78, 79)), row.names = c(NA, -79L), class = c("tbl_df", "tbl", 
    "data.frame"))
knitr::kable(dados.ord2)

ordem	codmun	nomemun	pop2015	rank
20	5002704	Campo Grande - MS	853622	1
32	5003702	Dourados - MS	212870	2
78	5008305	Três Lagoas - MS	113619	3
26	5003207	Corumbá - MS	108656	4
62	5006606	Ponta Porã - MS	86717	5
53	5005707	Naviraí - MS	51535	6
73	5007901	Sidrolândia - MS	51355	7
56	5006200	Nova Andradina - MS	50893	8
9	5001102	Aquidauana - MS	47162	9
50	5005400	Maracaju - MS	43078	10
59	5006309	Paranaíba - MS	41495	11
3	5000609	Amambai - MS	37590	12
65	5007208	Rio Brilhante - MS	34776	13
28	5003306	Coxim - MS	33139	14
18	5002407	Caarapó - MS	28437	15
51	5005608	Miranda - MS	27104	16
45	5005004	Jardim - MS	25473	17
70	5007695	São Gabriel do Oeste - MS	24982	18
4	5000708	Anastácio - MS	24748	19
8	5001003	Aparecida do Taboado - MS	24414	20
14	5002100	Bela Vista - MS	24113	21
64	5007109	Ribas do Rio Pardo - MS	23167	22
42	5004700	Ivinhema - MS	22928	23
40	5004502	Itaporã - MS	22896	24
23	5002951	Chapadão do Sul - MS	22620	25
48	5005202	Ladário - MS	21860	26
12	5001904	Bataguassu - MS	21775	27
22	5002902	Cassilândia - MS	21622	28
16	5002209	Bonito - MS	21047	29
41	5004601	Itaquiraí - MS	20162	30
77	5008008	Terenos - MS	19914	31
55	5006002	Nova Alvorada do Sul - MS	19656	32
27	5003256	Costa Rica - MS	19508	33
67	5007406	Rio Verde de Mato Grosso - MS	19462	34
34	5003801	Fátima do Sul - MS	19220	35
52	5005681	Mundo Novo - MS	17884	36
74	5007935	Sonora - MS	17483	37
63	5006903	Porto Murtinho - MS	16514	38
38	5004304	Iguatemi - MS	15637	39
25	5003157	Coronel Sapucaia - MS	14815	40
1	5000203	Água Clara - MS	14474	41
54	5005806	Nioaque - MS	14233	42
19	5002605	Camapuã - MS	13731	43
60	5006358	Paranhos - MS	13494	44
29	5003454	Deodápolis - MS	12650	45
33	5003751	Eldorado - MS	12128	46
17	5002308	Brasilândia - MS	11903	47
10	5001243	Aral Moreira - MS	11399	48
13	5002001	Batayporã - MS	11208	49
75	5007950	Tacuru - MS	11035	50
30	5003488	Dois Irmãos do Buriti - MS	10965	51
72	5007703	Sete Quedas - MS	10832	52
6	5000856	Angélica - MS	10149	53
37	5004106	Guia Lopes da Laguna - MS	10136	54
36	5004007	Glória de Dourados - MS	9992	55
5	5000807	Anaurilândia - MS	8844	56
7	5000906	Antônio João - MS	8679	57
43	5004809	Japorã - MS	8567	58
15	5002159	Bodoquena - MS	7898	59
61	5006408	Pedro Gomes - MS	7794	60
39	5004403	Inocência - MS	7664	61
69	5007554	Santa Rita do Pardo - MS	7633	62
49	5005251	Laguna Carapã - MS	7017	63
44	5004908	Jaraguari - MS	6860	64
11	5001508	Bandeirantes - MS	6771	65
71	5007802	Selvíria - MS	6455	66
47	5005152	Juti - MS	6399	67
79	5008404	Vicentina - MS	6027	68
21	5002803	Caracol - MS	5838	69
31	5003504	Douradina - MS	5723	70
24	5003108	Corguinho - MS	5513	71
68	5007505	Rochedo - MS	5252	72
58	5006275	Paraíso das Águas - MS	5150	73
2	5000252	Alcinópolis - MS	5038	74
66	5007307	Rio Negro - MS	4910	75
57	5006259	Novo Horizonte do Sul - MS	4306	76
46	5005103	Jateí - MS	4038	77
76	5007976	Taquarussu - MS	3570	78
35	5003900	Figueirão - MS	3012	79

# é possivel detectar 30 municipios com pop>20000 em dados.ord2 fazendo para as
# 5, 10 e 30 maiores para dataset das 5 maiores
options(scipen = 999)  # desativando notacao cientifica
data5 <- dados.ord2[1:5, c("ordem", "nomemun", "pop2015")]
data10 <- dados.ord2[1:10, c("ordem", "nomemun", "pop2015")]
data30 <- dados.ord2[1:30, c("ordem", "nomemun", "pop2015")]

# Para todos os municipios estabelece o rank do maior para menor em populacao
# rank.pop=1 para o maior e 79 para o menor
rank.pop <- rank(-dados.ord2$pop2015)

3 Regressões

Observam-se assimetrias e poucas cidades grandes, \(|β|<1\) para as regressões (1 : -0.664), (2 : -0.768) e (4 : -0.932). Na regressão (3 : -1.001), a estrutura urbana é descentralizada com vários grandes centros, pois \(|β|>1\).

regressao <- lm(log(rank.pop) ~ log(dados.ord2$pop2015))
summary(regressao)

Call:
lm(formula = log(rank.pop) ~ log(dados.ord2$pop2015))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.59702 -0.14191  0.01674  0.17532  0.31619 

Coefficients:
                        Estimate Std. Error t value            Pr(>|t|)    
(Intercept)             12.43208    0.24846   50.04 <0.0000000000000002 ***
log(dados.ord2$pop2015) -0.93199    0.02554  -36.49 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2156 on 77 degrees of freedom
Multiple R-squared:  0.9453,    Adjusted R-squared:  0.9446 
F-statistic:  1332 on 1 and 77 DF,  p-value: < 0.00000000000000022

# para os 5 maiores
rank.pop5 <- rank(-data5$pop2015)  # estabelece o rank do maior para menor em populacao
knitr::kable(rank.pop5)

x
1
2
3
4
5

regressao5 <- lm(log(rank.pop5) ~ log(data5$pop2015))
summary(regressao5)

Call:
lm(formula = log(rank.pop5) ~ log(data5$pop2015))

Residuals:
       1        2        3        4        5 
 0.07263 -0.15680 -0.16839  0.08962  0.16294 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)   
(Intercept)         8.99977    1.14895   7.833  0.00433 **
log(data5$pop2015) -0.66429    0.09468  -7.016  0.00595 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1759 on 3 degrees of freedom
Multiple R-squared:  0.9426,    Adjusted R-squared:  0.9234 
F-statistic: 49.23 on 1 and 3 DF,  p-value: 0.005947

# para os 10 maiores
rank.pop10 <- rank(-data10$pop2015)  # estabelece o rank do maior para menor em populacao
regressao10 <- lm(log(rank.pop10) ~ log(data10$pop2015))
summary(regressao10)

Call:
lm(formula = log(rank.pop10) ~ log(data10$pop2015))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26541 -0.14132  0.01299  0.14337  0.19403 

Coefficients:
                    Estimate Std. Error t value    Pr(>|t|)    
(Intercept)         10.30039    0.72305   14.25 0.000000574 ***
log(data10$pop2015) -0.76769    0.06296  -12.19 0.000001898 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1757 on 8 degrees of freedom
Multiple R-squared:  0.9489,    Adjusted R-squared:  0.9426 
F-statistic: 148.7 on 1 and 8 DF,  p-value: 0.000001898

# para os 30 maiores
rank.pop30 <- rank(-data30$pop2015)  # estabelece o rank do maior para menor em populacao
regressao30 <- lm(log(rank.pop30) ~ log(data30$pop2015))
summary(regressao30)

Call:
lm(formula = log(rank.pop30) ~ log(data30$pop2015))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.43197 -0.09932 -0.04755  0.13212  0.58559 

Coefficients:
                    Estimate Std. Error t value            Pr(>|t|)    
(Intercept)         13.08186    0.48549   26.95 <0.0000000000000002 ***
log(data30$pop2015) -1.00075    0.04573  -21.88 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2034 on 28 degrees of freedom
Multiple R-squared:  0.9448,    Adjusted R-squared:  0.9428 
F-statistic: 478.9 on 1 and 28 DF,  p-value: < 0.00000000000000022

library(stargazer)
resultados <- stargazer(list(regressao5, regressao10, regressao30, regressao), type = "text", 
    style = ("all"), omit.stat = ("f"), rownames = NULL)

==================================================================================
                                         Dependent variable:                      
                    --------------------------------------------------------------
                    log(rank.pop5) log(rank.pop10) log(rank.pop30)  log(rank.pop) 
                         (1)             (2)             (3)             (4)      
----------------------------------------------------------------------------------
pop2015)              -0.664***                                                   
                       (0.095)                                                    
                      t = -7.016                                                  
                      p = 0.006                                                   
pop2015)                              -0.768***                                   
                                       (0.063)                                    
                                     t = -12.193                                  
                                     p = 0.00001                                  
pop2015)                                              -1.001***                   
                                                       (0.046)                    
                                                     t = -21.884                  
                                                      p = 0.000                   
pop2015)                                                              -0.932***   
                                                                       (0.026)    
                                                                     t = -36.491  
                                                                      p = 0.000   
Constant               9.000***       10.300***       13.082***       12.432***   
                       (1.149)         (0.723)         (0.485)         (0.248)    
                      t = 7.833      t = 14.246      t = 26.946      t = 50.036   
                      p = 0.005      p = 0.00000      p = 0.000       p = 0.000   
----------------------------------------------------------------------------------
Observations              5              10              30              79       
R2                      0.943           0.949           0.945           0.945     
Adjusted R2             0.923           0.943           0.943           0.945     
Residual Std. Error 0.176 (df = 3) 0.176 (df = 8)  0.203 (df = 28) 0.216 (df = 77)
==================================================================================
Note:                                                  *p<0.1; **p<0.05; ***p<0.01

# exportando resultados para um arquivo texto
write.table(resultados, "resultados.txt", sep = "\t")

4 Gráfico

O gráfico de dispersão permite observar a relação de Zipf.

library(ggplot2)
library(Hmisc)
ggplot() + geom_point(data = data5, aes(x = log(rank.pop5), y = log(sort(data5$pop2015, 
    decreasing = TRUE))))

ggplot() + geom_point(data = dados.ord2, aes(x = log(rank.pop), y = log(sort(pop2015, 
    decreasing = TRUE))))

# together with variable you want to plot against
# regressao<-lm(log(rank.pop)~log(dados.ord2$pop2015))
predicted_df <- data.frame(rank_pred = predict(regressao, dados.ord2), log(dados.ord2$pop2015))

# this is the predicted line of multiple linear regression
ggplot(data = dados.ord2, aes(x = log(rank.pop), y = log(dados.ord2$pop2015))) + 
    geom_point(color = "blue") + geom_line(color = "red", data = predicted_df, aes(x = rank_pred, 
    y = log(dados.ord2$pop2015)))

5 Primazia

A medida da primazia dá uma ideia do tamanho relativo da principal cidade em relação à soma das populações das 5, 10, 30 maiores ou até à soma de todas as demais cidades.
Portanto, a primazia para as 5 maiores será:

\[ {P_5} = \frac{{{N_1}}}{{\sum\limits_{i = 1}^5 {{N_i}} }} \]

Uma opção na literatura é adotar uma população de corte ou em alguns casos coloca-se o denominador exclusive a \(N_1\).

# Calculo da primazia conforme RUIZ (2005) e Rosen e Resnick (1979)
# http://www.scielo.br/pdf/ee/v35n4/v35n4a05.pdf para P5 - cinco maiores
P5 = (data5$pop2015[1]/sum(data5$pop2015))
P5

[1] 0.6205975

# para P10 - dez maiores
n = length(data10$pop2015)
P10 = (data10$pop2015[1]/sum(data10$pop2015))
P10

[1] 0.5270876

# para P30 - trinta maiores
n = length(data30$pop2015)
P30 = (data30$pop2015[1]/sum(data30$pop2015))
P30

[1] 0.3981715

# para todos
P79 = (dados.ord2$pop2015[1]/sum(dados.ord2$pop2015))
P79

[1] 0.3219715

Referências

CHRISTALLER, Walter. Die zentralen Orte in Süddeutschland (Central places in southern Germany). Jena: Gustav Fischer, 1933.

LOSCH, August. The Economics of Location. Jena: Fischer, 1940. English translation. New Haven, Conn.: Yale Univ. Press, 1954.

ROSEN, Kenneth T.; RESNICK, Mitchel. The size distribution of cities: an examination of the Pareto law and primacy. Journal of Urban Economics, v. 8, n. 2, p. 165-186, 1980.

RUIZ, Ricardo Machado. Estruturas Urbanas Comparadas: Estados Unidos e Brasil. Estudos econômicos. São Paulo: FEA-USP, 35(3):715-737, out-dez 2005. Disponível em: http://www.scielo.br/pdf/ee/v35n4/v35n4a05.pdf.

ZIPF, George Kingsley. Human behavior and the principle of least effort. Oxford, England: Addison-Wesley Press. 1949.