sfc软件包中文说明

sfc软件包于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函数包括六个主要参数,分别是datamodelsample.sizerand.seedcheckinner。其中前两个参数是必填参数,而后四个参数都有默认值的设置。如果不需要做不确定性分析,则只要填datamodel这两个参数即可,因为后四个变量都是针对不确定性分析的。这里,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函数计算的结果可能返回四个值:resultsample.sizesampleinner.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软件包的最大突破口,恰好是建立在对数据与模型的分开管理上的。至于核算本身,未来还有很多问题等着我们去解决,需要大家共同的努力。

  1. 实际上最小核算单元对应的字段所起的作用是标识,为此最保险的做法就是全部用字符型变量。

  2. 常用的分布函数参数一般都不超过3个(见表2),当然如果要超过3个也可以继续编号下去。