Prose_如何用Vensim做人口预测？

Posted on In , Waline: Views: Word count in article: 949 Reading time ≈ 3 mins.

最近几天的主要精力是整理广东某市的统计年鉴与人口普查的数据，年鉴侧重于户籍人口，普查侧重于常住人口。这里用基于JavaScript的Echarts做了男女性别比，用VensimPLE做了一次系统动力学的常住人口预测示例；尔后用R语言的animation包与ggplot2包做了人口金字塔的动图，最终效果如下图。

Fig.1 Population Pyramid，2010-2035

这一常住人口模型是十分粗糙的，它假定了常住人口流动幅度很小（模型中设定为0）、死亡率变化幅度很小（模型中不变）、总和生育率恒定在1.6左右，乃及不同年龄别的男女比为计算方便，假定为1：1。如果需要精细预测，就要进一步考虑到年龄别的男女死亡率、外出务工人员的老年返乡等因素。
2010年，该市的常住人口结构很有意思，我们以性别比为例。男性数量在同期人群占比上存在两个峰值，一个是1948年前后出生的人口，一个是2003年前后出生的人口。当然，1913年前后出生的人口考虑到女性寿命长于男性，可能也会低估同期男性占比。但总得来讲，性别比应该是受到同期的社会环境影响，比如中间年龄段的洼地可能是由于该市是劳动力的输出大市。

同时，诚如Fig1的人口金字塔，2010年前后，该市人口结构是增长性人口结构；预期在2035年前后转为稳定型人口结构；如果进一步预测，该市可能2050年前后开始迅速老龄化。但相比于东北城市，该市的生育率仍然是相当高的，城市人口也十分年轻。
最后，本文使用的系统动力学的人口模型如下图；人口金字塔的制图参考了EasyCharts2017年5月的文章，代码如下。

Fig.2 System Dynamics Graph
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
library(ggplot2)
library(animation)
library(dplyr)
library(tidyr)
library(xlsx)
library(ggthemes)

male<-read.xlsx("population simulation.xlsx",sheetName="male2",header=T,encoding='UTF-8',check.names = FALSE)
female<-read.xlsx("population simulation.xlsx",sheetName="female2",header=T,encoding='UTF-8',check.names = FALSE)
female<-female%>%gather(Year,Population,-1)
male<-male%>%gather(Year,Population,-1)
female$Population<-female$Population*-1
male$sex<-"male";female$sex<-"female"
maoming<-rbind(male,female)
maoming$Population <- round(maoming$Population)
maoming$abs_pop <- abs(maoming$Population)
maoming$agegroup<-factor(maoming$agegroup,
levels=c("0-4","5-9","10-14","15-19","20-24","25-29","30-34","35-39","40-44","45-49",
"50-54","55-59","60-64","65-69","70-74","75-79","80-84",'85-90','91-94','95-99'),order=T)
m<-seq(2010,2035,by=1)

saveGIF({
  for (i in m) {
    title <- as.character(i)
    year_data <- filter(maoming,Year==i)
    g1<-ggplot(year_data,aes(x =agegroup,y=Population,fill=sex,width=1)) +
      coord_fixed()+
      coord_flip() +
      geom_bar(data=subset(year_data,sex=="female"),stat = "identity") +
      geom_bar(data=subset(year_data,sex=="male"), stat = "identity") +
      scale_y_continuous(breaks = seq(-400000,400000,length=9),
                         labels = paste0(as.character(c(abs(seq(-40,40,length=9)))), "w"),
                         limits = c(-450000, 450000)) +
      theme_economist(base_size = 14) +
      scale_fill_manual(values = c('#D40225', '#374F8F')) +
      labs(title=paste0("Population structure of Maoming:", title),
           caption="Data Source: Census data of Maoming City, the 2010 Revision"
           ,y="Population",x="Age") +
      guides(fill=guide_legend(reverse = TRUE))+
      theme(
        legend.position =c(0.8,0.9),
        legend.title = element_blank(),
        plot.title = element_text(size=20),
        plot.caption = element_text(size=12,hjust=0),
      )
    print(g1)
  }
},movie.name='maoming_population.gif',interval=0.5,ani.width=700,ani.height=600)