Prose_如何用Vensim做人口预测?

最近几天的主要精力是整理广东某市的统计年鉴与人口普查的数据,年鉴侧重于户籍人口,普查侧重于常住人口。这里用基于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)