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 数据文件

NAMENOTEUNITTIMEDATACVDISTP1P2P3
DRPP ore productionktY20002.0e+010.1unif16.53623.464NA
PRPP content in ore1all5.0e-010.1triangle0.3780.6220.500
PRRP ore production ratio1all8.0e-010.2beta4.2001.0500.000
DFPFertilizer productionktY20001.5e+010.1unif12.40217.598NA
PFPP content in fertilizer1all4.0e-010.2triangle0.2040.5960.400
DCPDaily use chemical productionktY20005.0e+000.1unif4.1345.866NA
PCPP content in daily use chemical1all3.0e-010.3triangle0.0800.5200.300
PCRWaste chemical recycle rate1all5.0e-010.2beta12.00012.0000.000
DMPCrop productionktY20001.0e+030.1unif826.7951173.205NA
PMPP content in crop1all5.0e-030.2triangle0.0030.0070.005
PHRProportion of feed1all4.0e-010.0none0.400NANA
PNRProportion of food1all4.5e-010.0none0.450NANA
DAPAnimal productionktY20008.0e+020.1unif661.436938.564NA
PAPP content in animal1all2.0e-030.2triangle0.0010.0030.002
PARMeat production ratio1all7.5e-010.2beta5.5001.8330.000
DRPP ore productionktY20054.0e+010.1unif33.07246.928NA
DFPFertilizer productionktY20053.0e+010.1unif24.80435.196NA
DCPDaily use chemical productionktY20051.0e+010.1unif8.26811.732NA
DMPCrop productionktY20052.0e+030.1unif1653.5902346.410NA
DAPAnimal productionktY20051.6e+030.1unif1322.8721877.128NA

表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的含义
β分布betashape1, shape2, ncp
二项式分布binomsize, prob
Cauchy分布cauchylocation, scale
卡方分布chisqdf, ncp
指数分布exprate
F分布fdf1, df1, ncp
Γ分布gammashape, scale
几何分布geomprob
超几何分布hyperm, n, k
对数正态分布lnormmeanlog, sdlog
logistic分布logislocation, scale
负二项式分布nbinomsize, prob
正态分布normmean, sd
Poisson分布poislambda
t分布tdf, ncp
均匀分布unifmin, max
Weibull分布weibullshape, scale
Wilcoxon分布wilcoxm, n
无分布noneDATA中的数值

需要说明的是,除了以上提到的那些变量名外,其他的任何变量名原则上(除非出了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模型文件
NAMESTARTENDFUN
P rock production
PF_01RChDRP*PRP*PRR
PF_02RWDRP*PRP-PF_01
P chemical production
PF_03ChCrDFP*PFP
PF_04ChHDCP*PCP
PF_05WChPF_06*PCR
PF_06ChW(PF_01-PF_03-PF_04)/(1-PCR)
Crop production
PF_07CrADMP*PMP*PHR
PF_08CrHDMP*PMP*PNR
PF_09CrWPF_03-PF_07-PF_08
Animal production
PF_10AHDAP*PAP*PAR
PF_11AWPF_07-PF_10
Human consumption
PF_12HWPF_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 普通单次计算结果

TIMESTARTENDFLOW
Y2000RCh8.00
Y2000WCh0.50
Y2000ChCr6.00
Y2000ChH1.50
Y2000CrH2.25
Y2000AH1.20
Y2000CrA2.00
Y2000RW2.00
Y2000ChW1.00
Y2000CrW1.75
Y2000HW4.95
Y2000AW0.80
Y2005RCh16.00
Y2005WCh1.00
Y2005ChCr12.00
Y2005ChH3.00
Y2005CrH4.50
Y2005AH2.40
Y2005CrA4.00
Y2005RW4.00
Y2005ChW2.00
Y2005CrW3.50
Y2005HW9.90
Y2005AW1.60

表5 Monte Carlo不确定性分析计算结果

TIMESTARTENDMEANMEDIANSDCVQ05Q25Q75Q95
Y2000RCh9.1138.9251.4450.1596.8238.17210.00211.663
Y2000RW1.4441.1371.2720.8800.1120.5072.0094.243
Y2000ChCr5.8615.7571.0450.1784.2405.1646.4207.734
Y2000ChH1.4591.4670.4020.2750.7941.1171.7802.062
Y2000ChW3.7032.9932.9500.7970.3381.3805.5299.179
Y2000CrH2.2012.1900.4280.1941.5721.8492.5052.948
Y2000CrA1.9561.9460.3800.1941.3971.6432.2262.620
Y2000CrW1.7041.4501.2160.7140.1820.7992.5733.801
Y2000HW4.8114.8170.7310.1523.5634.3325.3835.881
Y2000AH1.1511.1740.3390.2950.6160.8761.3671.722
Y2000AW0.8060.7740.4630.5740.0810.4731.1361.546
Y2000WCh1.9111.3861.7060.8930.1210.6342.7185.207
Y2005RCh18.15317.7952.8050.15513.82716.19519.97422.930
Y2005RW2.7342.3851.9210.7030.3081.1454.0376.299
Y2005ChCr11.39311.2571.9910.1758.46210.04212.41114.984
Y2005ChH2.8782.7850.8920.3101.5572.1973.4584.447
Y2005ChW8.0176.8296.0490.7540.8063.07010.98920.803
Y2005CrH4.2674.2500.8530.2002.9573.6564.8805.717
Y2005CrA3.7933.7780.7580.2002.6283.2504.3375.082
Y2005CrW3.3332.9512.2540.6760.3081.4534.7607.880
Y2005HW9.3489.3221.4120.1516.9838.56910.24311.586
Y2005AH2.2022.1600.6280.2851.2881.7412.6273.273
Y2005AW1.5911.5090.8260.5200.3830.9532.2022.810
Y2005WCh4.1363.2043.6540.8840.4111.4895.70011.671

6 图形化界面

由于sfc软件包是基于R开发的,所以对于非R用户则不容易使用。为此,我们基于shiny平台,开发了一个sfc的图形化界面。界面的链接为https://ctfysh.shinyapps.io/sample_app/。 非R用户只需要知道数据文件和简化版的模型文件如何整理即可使用这个平台进行计算。这也说明了,如何有效管理我们核算数据,以及我们核算方法,是物质流核算关键之所在。至于如何计算,以及在什么平台上,用什么软件来计算,这些都是次要的。而sfc软件包的最大突破口,恰好是建立在对数据与模型的分开管理上的。至于核算本身,未来还有很多问题等着我们去解决,需要大家共同的努力。


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

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