sfc软件包中文说明
27 Aug 2016sfc软件包于2016年8月25日正式出现在CRAN上,其主要功能是进行大规模物质流核算和不确定性分析。与传统物质流核算方法不同的是,sfc软件包是基于数据和模型分离的思想构建的。这种方法的优势在于能够把模型从复杂的数据结构中抽离出来,从而使研究者能够把注意力集中在模型的构建上,以避免复杂数据结构对建模过程的干扰。
1 安装sfc
sfc软件包依赖于dplyr、tidyr、triangle、zoo、sna等R软件包,因此在R平台上安装sfc之前,最好先安装好这些依赖的软件包。sfc软件包可以用如下方式安装:
install.packages("sfc")
如果要安装sfc软件包的最新开发版,则可以通过以下方式:
install.packages("devtools")
devtools::install_github("ctfysh/sfc")
2 数据文件
物质流核算涉及的数据的类型多且量大,同时数据存在着一定的不确定性,所以需要有一种既简单直观又易于被模型所有调用的数据格式,来管理这些数据。表1就是sfc软件包为了物质流核算所设计的一种数据文件格式,以.csv的形式存储,或者直接以数据框(data frame)的形式存在于内存中,供sfc软件包中的sfc
函数调用。
表1 数据文件
NAME | NOTE | UNIT | TIME | DATA | CV | DIST | P1 | P2 | P3 |
---|---|---|---|---|---|---|---|---|---|
DRP | P ore production | kt | Y2000 | 2.0e+01 | 0.1 | unif | 16.536 | 23.464 | NA |
PRP | P content in ore | 1 | all | 5.0e-01 | 0.1 | triangle | 0.378 | 0.622 | 0.500 |
PRR | P ore production ratio | 1 | all | 8.0e-01 | 0.2 | beta | 4.200 | 1.050 | 0.000 |
DFP | Fertilizer production | kt | Y2000 | 1.5e+01 | 0.1 | unif | 12.402 | 17.598 | NA |
PFP | P content in fertilizer | 1 | all | 4.0e-01 | 0.2 | triangle | 0.204 | 0.596 | 0.400 |
DCP | Daily use chemical production | kt | Y2000 | 5.0e+00 | 0.1 | unif | 4.134 | 5.866 | NA |
PCP | P content in daily use chemical | 1 | all | 3.0e-01 | 0.3 | triangle | 0.080 | 0.520 | 0.300 |
PCR | Waste chemical recycle rate | 1 | all | 5.0e-01 | 0.2 | beta | 12.000 | 12.000 | 0.000 |
DMP | Crop production | kt | Y2000 | 1.0e+03 | 0.1 | unif | 826.795 | 1173.205 | NA |
PMP | P content in crop | 1 | all | 5.0e-03 | 0.2 | triangle | 0.003 | 0.007 | 0.005 |
PHR | Proportion of feed | 1 | all | 4.0e-01 | 0.0 | none | 0.400 | NA | NA |
PNR | Proportion of food | 1 | all | 4.5e-01 | 0.0 | none | 0.450 | NA | NA |
DAP | Animal production | kt | Y2000 | 8.0e+02 | 0.1 | unif | 661.436 | 938.564 | NA |
PAP | P content in animal | 1 | all | 2.0e-03 | 0.2 | triangle | 0.001 | 0.003 | 0.002 |
PAR | Meat production ratio | 1 | all | 7.5e-01 | 0.2 | beta | 5.500 | 1.833 | 0.000 |
DRP | P ore production | kt | Y2005 | 4.0e+01 | 0.1 | unif | 33.072 | 46.928 | NA |
DFP | Fertilizer production | kt | Y2005 | 3.0e+01 | 0.1 | unif | 24.804 | 35.196 | NA |
DCP | Daily use chemical production | kt | Y2005 | 1.0e+01 | 0.1 | unif | 8.268 | 11.732 | NA |
DMP | Crop production | kt | Y2005 | 2.0e+03 | 0.1 | unif | 1653.590 | 2346.410 | NA |
DAP | Animal production | kt | Y2005 | 1.6e+03 | 0.1 | unif | 1322.872 | 1877.128 | NA |
表1中用于sfc
函数的字段有:NAME、TIME、DATA、DIST、P1、P2、P3。这里NAME、TIME、DATA是用于物质流核算的字段,用来表征变量名、变量对应的时间和变量的取值。这3个字段可以看成是一种函数关系:NAME(TIME)=DATA。其中NAME和DATA是必需的信息而TIME则根据实际情况决定是否需要。如果涉及到空间信息,还可以添加SITE字段。如果模型数据结构很复杂,既包括时间信息,又包括空间信息,还包括其他一些数据信息,比如省、市、县这些信息,可以用自定义字段C1、C2、C3等来对应这些信息。我们定义TIME、SITE、C1、C2、C3这些字段共同决定的一个核算单元为最小核算单元。对于最小核算单元,我们有一个要求,就是它必须包含一整套变量名(NAME)和变量的取值(DATA)。因为我们在核算的时候,假设每个最小核算单元之间是独立的,也就是说我们是对每个最小核算单元进行一一计算的。如果哪个最小核算单元所给的数据不齐全,那自然就无法进行核算。当然,我们也会遇到这样一种情况,就是某个数据在所有年份的取值都是一样的。如果我们每年都复制一下这个数值,会比较麻烦。所以,我们在TIME以及最小核算单元的其他字段中设计了一个特殊的字符”all”,当填入这个字符的时候,就表示对NAME所表示的变量,在该字段下各种情况下的取值都是一样的。而对于TIME这个字段,一般情况下都是以年为单位的,建议在年的前面加上”Y”(如表1所示),这样可以防止某些出错的可能1。另外,TIME这个字段是有顺序的,但有些时候我们的数据可能是分段的,即在某些时段的变量的取值是一样的。这时我们只需要给出分段点处的那些变量取值即可,其他年份的数据则可通过算法自动插补进去。在sfc包中所采用的算法是LOCF(时间序列中缺失的数据用前一个值代替),那么我们只需要给出那些取值一样的时段中的第一个变量及其取值即可。采用以上两种方式能够减少我们的工作量和数据文件的大小,能够规避我们在重复输入一个值的时候出现的笔误,但是却可能增大我们检查每个最小核算单元数据完整性的难度。在整理数据的时候,建议一个变量一个变量的整理,而不是按一个最小核算单元接着一个最小核算单元的方式整理。
DIST、P1、P2、P3是用于不确定性分析的字段,分别表示分布函数和3个分布参数。在R中,常用的分布函数及参数设定如表2所示。由于不同的分布函数的参数个数不同,sfc软件包为了统一,分别用P1、P2、P32的顺序对应分布函数参数的顺序。如果不足3个,缺的地方不填或者用NA代替。
表2 R中常见的分布函数及参数设定
概率分布 | DIST参数 | 分布参数P1、P2、P3的含义 |
---|---|---|
β分布 | beta | shape1, shape2, ncp |
二项式分布 | binom | size, prob |
Cauchy分布 | cauchy | location, scale |
卡方分布 | chisq | df, ncp |
指数分布 | exp | rate |
F分布 | f | df1, df1, ncp |
Γ分布 | gamma | shape, scale |
几何分布 | geom | prob |
超几何分布 | hyper | m, n, k |
对数正态分布 | lnorm | meanlog, sdlog |
logistic分布 | logis | location, scale |
负二项式分布 | nbinom | size, prob |
正态分布 | norm | mean, sd |
Poisson分布 | pois | lambda |
t分布 | t | df, ncp |
均匀分布 | unif | min, max |
Weibull分布 | weibull | shape, scale |
Wilcoxon分布 | wilcox | m, n |
无分布 | none | DATA中的数值 |
需要说明的是,除了以上提到的那些变量名外,其他的任何变量名原则上(除非出了bug)是不会参与计算的,所以可以随便增加这些字段,比如表1中的NOTE、UNIT、CV这些字段。为了区别R对象,字段名推荐全部用大写字母表示。
3 模型文件
模型是物质流核算的最核心部分,它其实就是每一条流计算方法的集合。由于我们是在R平台上开发的sfc软件包,所以模型文件其实就是符合R编码规则的代码块,如下所示:
## P rock production
PF["R","Ch",]=DRP*PRP*PRR
PF["R","W",]=DRP*PRP-PF["R","Ch",]
## P chemical production
PF["Ch","Cr",]=DFP*PFP
PF["Ch","H",]=DCP*PCP
PF["Ch","W",]=(PF["R","Ch",]-PF["Ch","Cr",]-PF["Ch","H",])/(1-PCR)
PF["W","Ch",]=PF["Ch","W",]*PCR
## Crop production
PF["Cr","A",]=DMP*PMP*PHR
PF["Cr","H",]=DMP*PMP*PNR
PF["Cr","W",]=PF["Ch","Cr",]-PF["Cr","A",]-PF["Cr","H",]
## Animal production
PF["A","H",]=DAP*PAP*PAR
PF["A","W",]=PF["Cr","A",]-PF["A","H",]
## Human consumption
PF["H","W",]=PF["Ch","H",]+PF["Cr","H",]+PF["A","H",]
以”=”为界(当然也可以用”<-“),”=”左边表示物质流,它用一个三维数组来进行存储,其中第三个维度为空,表示的是最小核算单元。前两个维度分别表示物质流的始端(START)和末端(END)。这里PF是数组变量名,也可以使用别的名称。”=”右边则是左边物质流的计算方法,每一条流都是由数据文件中的变量计算出来的。当然,我们也可以在前面已经计算的物质流的基础上进行其他物质流的计算。直接用数据文件中的变量计算物质流的方法我们称为”独立型”核算,利用前面已经计算出来的物质流来计算的方法我们称为”依附型”核算。而在依附型核算中有一个特殊的算法,它是根据物质守恒的原理进行计算的,我们称为”平衡型”核算。
考虑到我们比较习惯于用MS Excel来整理数据,在sfc软件包里推出了一款简化版的模型文件(表3)。与数据文件类似,这种模型文件也是.csv格式的,它包括四个关键字段:NAME、START、END、FUN。这四个关键字段分别表示物质流的名称、始端节点、末端节点和核算公式。其中,对于依附型核算,物质流的名称可以作为一个变量使用,这个与数据文件中的变量的作用是一致的(见表3)。需要说明的一点是:sfc软件包只识别四个关键字段都存在的行。所以我们可以在NAME那个字段里随意地写说明性文字,包括空行等。
表3 简化版sfc模型文件
NAME | START | END | FUN |
---|---|---|---|
P rock production | |||
PF_01 | R | Ch | DRP*PRP*PRR |
PF_02 | R | W | DRP*PRP-PF_01 |
P chemical production | |||
PF_03 | Ch | Cr | DFP*PFP |
PF_04 | Ch | H | DCP*PCP |
PF_05 | W | Ch | PF_06*PCR |
PF_06 | Ch | W | (PF_01-PF_03-PF_04)/(1-PCR) |
Crop production | |||
PF_07 | Cr | A | DMP*PMP*PHR |
PF_08 | Cr | H | DMP*PMP*PNR |
PF_09 | Cr | W | PF_03-PF_07-PF_08 |
Animal production | |||
PF_10 | A | H | DAP*PAP*PAR |
PF_11 | A | W | PF_07-PF_10 |
Human consumption | |||
PF_12 | H | W | PF_04+PF_08+PF_10 |
4 计算方法
有了数据文件和模型文件之后,我们就可以用sfc
函数来进行计算了。sfc
函数包括六个主要参数,分别是data
、model
、sample.size
、rand.seed
、check
和inner
。其中前两个参数是必填参数,而后四个参数都有默认值的设置。如果不需要做不确定性分析,则只要填data
和model
这两个参数即可,因为后四个变量都是针对不确定性分析的。这里,sample.size
指的是样本量,如果样本量为1,那么就执行普通单次计算;如果样本量大于1,那么就执行Monte Carlo不确定性分析。rand.seed
设置的是随机数种子,为任意整数,设置后能够保证每次不确定性分析的结果都是一样的。如果随机数种子取值为NULL,那么每次不确定性分析的结果就会不一样。check
判断是否执行模型检查的开关,如果设置为TRUE,那么就执行;反之设置为FALSE,那么就不执行。设置这个开关的目的在于进行有约束的Monte Carlo不确定性分析。因为原则上我们在验证了普通单次计算结果的有效性后(包括物质流的大小和正负号),至少我们要保证Monte Carlo不确定性分析过程中所抽取的每个样本计算的物质流的方向(正负号)与普通单次计算结果是一致的,否则抽样无效。这样做的结果是有效的样本量少于我们设定的样本量。inner
判断是否返回Monte Carlo不确定性分析的中间结果,如果设置为TRUE,那么就返回;反之设置为FALSE,那么就不返回。这种设置的最大目的是降低计算过程中的内存消耗,特别是在大规模计算的时候。sfc
函数调用方式如下:
sfc(data, model, sample.size, rand.seed, check, inner)
5 计算结果
采用sfc
函数计算的结果可能返回四个值:result
、sample.size
、sample
和inner.result
。其中,result
表示物质流核算结果,如果sample.size = 1
,则可以得到普通单次计算结果(表4),同时sample.size
的返回值为1;如果sample.size = n
(n>1),则可以得到Monte Carlo不确定性分析计算结果(表5),同时sample.size
的返回值为m,这里m为有效的样本量。另外两个结果是在inner = TRUE
时才返回的,其中sample
表示每个变量的抽样结果,inner.result
表示抽样后计算出来的物质流的结果。这两个结果一般没有必要返回,只有在我们要对计算进行过程进一步分析的时候才有可能用到这两个结果。
表4 普通单次计算结果
TIME | START | END | FLOW |
---|---|---|---|
Y2000 | R | Ch | 8.00 |
Y2000 | W | Ch | 0.50 |
Y2000 | Ch | Cr | 6.00 |
Y2000 | Ch | H | 1.50 |
Y2000 | Cr | H | 2.25 |
Y2000 | A | H | 1.20 |
Y2000 | Cr | A | 2.00 |
Y2000 | R | W | 2.00 |
Y2000 | Ch | W | 1.00 |
Y2000 | Cr | W | 1.75 |
Y2000 | H | W | 4.95 |
Y2000 | A | W | 0.80 |
Y2005 | R | Ch | 16.00 |
Y2005 | W | Ch | 1.00 |
Y2005 | Ch | Cr | 12.00 |
Y2005 | Ch | H | 3.00 |
Y2005 | Cr | H | 4.50 |
Y2005 | A | H | 2.40 |
Y2005 | Cr | A | 4.00 |
Y2005 | R | W | 4.00 |
Y2005 | Ch | W | 2.00 |
Y2005 | Cr | W | 3.50 |
Y2005 | H | W | 9.90 |
Y2005 | A | W | 1.60 |
表5 Monte Carlo不确定性分析计算结果
TIME | START | END | MEAN | MEDIAN | SD | CV | Q05 | Q25 | Q75 | Q95 |
---|---|---|---|---|---|---|---|---|---|---|
Y2000 | R | Ch | 8.966 | 8.841 | 1.361 | 0.152 | 6.821 | 8.056 | 9.775 | 11.211 |
Y2000 | R | W | 1.259 | 1.035 | 1.065 | 0.846 | 0.099 | 0.432 | 1.859 | 2.977 |
Y2000 | Ch | Cr | 5.718 | 5.744 | 1.058 | 0.185 | 4.079 | 4.960 | 6.406 | 7.517 |
Y2000 | Ch | H | 1.457 | 1.468 | 0.476 | 0.327 | 0.726 | 1.085 | 1.832 | 2.212 |
Y2000 | Ch | W | 3.701 | 2.823 | 3.067 | 0.829 | 0.329 | 1.450 | 5.235 | 9.581 |
Y2000 | Cr | H | 2.213 | 2.180 | 0.404 | 0.183 | 1.621 | 1.899 | 2.496 | 2.874 |
Y2000 | Cr | A | 1.967 | 1.938 | 0.359 | 0.183 | 1.441 | 1.688 | 2.219 | 2.555 |
Y2000 | Cr | W | 1.538 | 1.298 | 1.113 | 0.723 | 0.176 | 0.626 | 2.265 | 3.798 |
Y2000 | H | W | 4.874 | 4.912 | 0.781 | 0.160 | 3.599 | 4.325 | 5.406 | 6.032 |
Y2000 | A | H | 1.204 | 1.193 | 0.362 | 0.301 | 0.660 | 0.937 | 1.425 | 1.837 |
Y2000 | A | W | 0.763 | 0.744 | 0.459 | 0.601 | 0.081 | 0.384 | 1.093 | 1.557 |
Y2000 | W | Ch | 1.910 | 1.358 | 1.776 | 0.930 | 0.140 | 0.641 | 2.707 | 5.112 |
Y2005 | R | Ch | 18.555 | 18.456 | 2.694 | 0.145 | 14.340 | 16.723 | 20.011 | 23.212 |
Y2005 | R | W | 2.271 | 1.755 | 1.892 | 0.833 | 0.199 | 0.808 | 3.144 | 5.985 |
Y2005 | Ch | Cr | 11.775 | 11.620 | 2.157 | 0.183 | 8.468 | 10.257 | 13.127 | 15.675 |
Y2005 | Ch | H | 2.860 | 2.809 | 0.985 | 0.344 | 1.345 | 2.097 | 3.583 | 4.545 |
Y2005 | Ch | W | 8.264 | 6.171 | 6.766 | 0.819 | 0.470 | 3.116 | 11.692 | 20.485 |
Y2005 | Cr | H | 4.502 | 4.402 | 0.887 | 0.197 | 3.242 | 3.874 | 5.116 | 6.074 |
Y2005 | Cr | A | 4.002 | 3.912 | 0.789 | 0.197 | 2.882 | 3.444 | 4.547 | 5.399 |
Y2005 | Cr | W | 3.271 | 2.969 | 2.104 | 0.643 | 0.314 | 1.567 | 4.582 | 7.390 |
Y2005 | H | W | 9.630 | 9.474 | 1.512 | 0.157 | 7.208 | 8.730 | 10.667 | 12.144 |
Y2005 | A | H | 2.268 | 2.241 | 0.652 | 0.287 | 1.274 | 1.849 | 2.696 | 3.378 |
Y2005 | A | W | 1.734 | 1.714 | 0.919 | 0.530 | 0.221 | 1.078 | 2.333 | 3.114 |
Y2005 | W | Ch | 4.344 | 2.971 | 4.416 | 1.016 | 0.210 | 1.464 | 6.125 | 12.563 |
6 图形化界面
由于sfc软件包是基于R开发的,所以对于非R用户则不容易使用。为此,我们基于shiny平台,开发了一个sfc的图形化界面。界面的链接为https://ctfysh.shinyapps.io/sample_app/。 非R用户只需要知道数据文件和简化版的模型文件如何整理即可使用这个平台进行计算。这也说明了,如何有效管理我们核算数据,以及我们核算方法,是物质流核算关键之所在。至于如何计算,以及在什么平台上,用什么软件来计算,这些都是次要的。而sfc软件包的最大突破口,恰好是建立在对数据与模型的分开管理上的。至于核算本身,未来还有很多问题等着我们去解决,需要大家共同的努力。