CN107391877A

CN107391877A - 一种基于r软件的入侵植物频度数量变化调查分析方法

Info

Publication number: CN107391877A
Application number: CN201710695009.2A
Authority: CN
Inventors: 黄小荣; 侯远瑞; 彭玉华
Original assignee: Guangxi Zhuang Autonomous Region Forestry Research Institute
Current assignee: Guangxi Zhuang Autonomous Region Forestry Research Institute
Priority date: 2017-08-15
Filing date: 2017-08-15
Publication date: 2017-11-24
Anticipated expiration: 2037-08-15
Also published as: CN107391877B

Abstract

一种基于R软件的入侵植物频度数量变化调查分析方法，使用沿样线规则等距间隔6～7 m设置小样方块，记录某一入侵植物的有无，调整小样方块尺寸使频度数据在40%～60%的范围；利用pwr、TrialSize、MESS和R自带功效分析功能，计算达到给定显著性水平和功效高低所需要的样本数，并对临时样地和固定样地入侵植物的频度数据变化进行显著性检验和统计功效分析；计算时只需代入实测数值、复制粘贴代码。本发明的入侵植物频度数据调查不受季节限制，调查覆盖面积大，数据重复性高，方法简单易行，为入侵植物频度数据的小样方块尺寸调整、独立两样本差异卡方检验、成对两样本McNemar检验提供了范例。

Description

一种基于R软件的入侵植物频度数量变化调查分析方法

技术领域

本发明属于植物生态学领域，涉及入侵植物频度数量变化的监测技术和应用R软件进行统计分析的技术方法，具体涉及入侵植物频度变化的调查分析方法。

背景技术

生物入侵是威胁林地、宜林地、荒山、灌草丛等生态系统的重大安全隐患。入侵植物争夺土壤养分、光照和水分，排挤乡土植物，改变群落结构，降低单位面积载畜量，导致林木生产、畜牧生产的巨大经济损失，给人类与家畜带来过敏、中毒、哮喘等威胁。监测入侵植物分布和动态变化，是有效预防和控制入侵植物扩散蔓延的前提条件。

频度数据调查只与小样方块中有没有某一植物有关，不记录植物丰度和物种个数。频度数据调查不受季节限制，不同调查人员之间数据重复性高。频度数据主要受小样方块的尺寸制约。监测入侵植物的频度(frequency)变化，是发现入侵威胁最快、成本最低的方法。入侵植物频度数据调查，一般是用分散均匀、数量较多的小样方块调查某一入侵植物的有或无(present/absent)，可以记录为目标植物出现的小样方块占总样方块的百分数(0～100)，也可以记录为目标植物出现的比例数(0～1)，还可记录为有或无目标植物的样方个数的计数(count)；文献中对频度数据有多种称谓，如比例数据(proportional)、二进制数据(binary，1和0)或类型数据(categorical)等等。

频度数据为伯努利分布(Bernoulli)，不符合正态分布，检验频度数据变化是否显著不能用传统的方差分析、t检验或wilcox检验，而要用比例检验、卡方检验、fisher检验或McNemar检验。目前统计教科书和说明书中提供的频度数据检验方法，绝大多数为医学、农学和经济学方面的例子，比如McNemar检验孪生兄弟肺癌与抽烟的关联，卡方检验金发和碧眼是否连锁、灌溉与不灌溉对作物黄叶枯叶数量的影响，比例检验掷多少次骰子才有中奖机会等等。频度数据统计分析五花八门的称谓，计数、百分数、0～1比例频度的交叉使用，加上入侵植物频度统计案例的缺乏，使从业人员在分析入侵植物频度数量变化是否显著时无从下手；甚至有人避简求繁，用排序分类来分析频度数据，把本来简单的问题复杂化了。梳理总结一套适用于入侵植物频度数据调查和变化是否显著的统计分析方法，具有重要的现实意义。

当监测目的为探索某一入侵植物频度数量变化是否“显著”时，样本大小(samplesize)和统计功效(power analysis)是必需权衡考虑的。样本大，则统计功效高、精确度高、调查成本高；样本小，则省时省力，但统计结果说服力不强、可检测最小变化(minimumdetectable change)较大。恰当地选择样本大小、显著水平和功效高低，才能以有限的成本产生显著水平和功效两方面均较满意的调查结果。本发明研究对比了多个样本功效分析程序，如MESS、TrialSize、pwr和R自带功效分析程序，根据我们监测入侵植物的经验，参考Elzinga(2001，Monitoring Plant and Animal Population)和Bonham(2013，Measurements forTerrestrial Vegetation)提供的相关技术参数，梳理出一套适用于入侵植物频度数量变化监测的样本和功效计算方法。本方法需在免费软件R中运行(https://www.R-project.org)。

发明内容

本发明的目的是针对目前对入侵植物频度数量变化监测中存在的技术问题，旨在提供一套包括样地预调查、小样方块尺寸调整、样本大小计算、频度数据变化显著性检验、统计功效检验的入侵植物数量变化测量分析方法；利用R软件作为计算器，直接输入代码实现样本计算、显著性检验和功效分析。

为了实现上述目的，本发明的技术方案如下：

一种基于R软件的入侵植物频度数量变化调查分析方法，其特征在于：主要操作步骤如下：

(1)在计算机上安装R软件及所需软件包：在https://www.R-project.org下载R3.4.1或更高版本，文件大小约173M；根据提示直接运行安装；在64位系统可安装成中文视窗界面，在32位系统只呈现全英文界面；安装完成后，双击桌面R 3.4.1图标打开RGui视窗，从“程序包”菜单的下拉菜单中选择“安装程序包”，在弹出的“Secure CRAN mirrors”中选择“China”后“确定”，再在弹出的“Packages”中选择“MESS”后“确定”，安装好MESS程序包；使用类似过程安装“pwr”、“TrialSize”程序包；

(2)样地预调查：选择典型样地，以系统抽样法(systematic sampling)进行入侵植物频度数据的初步调查；调查以样线沿山脚——山顶方向设置，详见图1，样线与样线平行且间隔7m；并以1m×1m小样方块沿样线分布，小样方块随机起点、等距间隔6m；记录64个1m×1m小样方块有无某一入侵植物之后，检查入侵植物出现的频度是否在20％～80％的范围内，如果频度低于或高于此范围，则需调整小样方块尺寸后重新开始预调查；

(3)调整小样方块的尺寸：频度数据与小样方块尺寸关系密切；小样方块的面积越小，植物的频度越低；可以通过增减小样方块的面积，使某一植物的频度控制在40％～60％，为观察其后续的频度变化预留较大的空间；以p1代表用1m×1m小样方块预调查所得的某一植物百分频度，p1＝((有/(有+无))×100)，若要将这一植物的百分频度调整为p2，则调整后小样方块面积为：

(log(100)-log(100-p2))/(-log(1-p1/100))

将p1和p2数值代入上式后，将整条代码复制粘贴到RGui，回车直接得到调整后的小样方块面积(m²)；计算公式参照Bonham(2013，p107)，以初始面积为1m²简化后编成R代码；R代码编写过程必需在英文输入法状态下进行，“()”和“～”的输入需特别注意，区分大小写、上下划线；

(4)显著水平和功效大小选择：一般统计分析的α＝0.05,β＝0.20，功效＝1-β；在入侵植物频度变化的监测中，误报与漏报后果同样严重，因此，采用α＝0.10,β＝0.10是比较合理的；

(5)单样本监测与达到给定精确度所需样本数计算：利用pwr程序包(Clay Ford，2017)进行计算：

library(pwr)

pwr.p.test(h＝ES.h(p1,p2),sig.level＝0.1,power＝0.9)

将预调查中的比例频度p1＝有/(有+无)和给定精确度|p1-p2|中的比例频度p2代入上式，将两条代码复制粘贴到RGui，回车得到样本数n；

(6)两样本临时样地频度比较的技术步骤：

a、样本数计算使用R软件自带功效分析函数，不需调用程序包，代码为：

power.prop.test(p1＝？,p2＝？,sig.level＝0.1,power＝0.9)

将预调查中两样本的入侵植物出现比例频度p1和p2代入上式，在RGui中运算得到每个组所需样本数n；

b、两样本临时样地入侵植物频度差异显著性检测；临时样地两样本的调查数据分布如表1，其中4个数值为样方个数计数；

表1：两样本临时样地入侵植物有无数据表

	第一年	第二年
			有某一入侵植物	C11	C21
无某一入侵植物	C10	C20

差异显著性检验代码为：

chisq.test(matrix(c(C11,C10,C21,C20),nrow＝2))

c、两样本临时样地入侵植物频度变化检验的功效分析使用如下代码：

power.prop.test(p1＝？,p2＝？,n＝？,sig.level＝0.1)

其中：p1、p2为最终实测两样本入侵植物的比例频度，n为最终实测的样本大小，这3个数值代入上式后，复制粘贴到RGui，回车得到统计功效power；

(7)两样本固定样地重复观测频度变化的技术步骤：

a、固定样地样本数使用TrialSize(Ed Zhuang,2013)计算较合理，代码为：

library(TrialSize)

(McNemar.Test(α,β,psai,paid))

其中psai＝p01/p10或psai＝p10/p01，psai≥1；paid＝p01+p10；数据收集方法见表2；将α,β,psai,paid代入上式后，两条代码复制粘贴至RGui中，回车得到所需样本大小n值；

表2：固定样地两次重复观测数据表

在表2中，每个观测小样方块可能有4种表现：两次记录都有某一入侵植物(p11)、第一次有但第二次无(p10)、第一次无但第二次有(p01)、两次记录都无该入侵植物(p00)；4个数值用0～1比例频度表示，4个数值之和等于1；

b、两样本固定样地重复观测入侵植物频度变化显著性检验，使用McNemar检验，代码如下：

mcnemar.test(matrix(c(p11,p10,p01,p00),nrow＝2))

另外，也可以在以上代码中采用按表2分列的计数数据；

C、两样本固定样地重复观测入侵植物频度变化检验的功效分析采用ClausEkstrom 2017年7月提供的MESS程序包，因为目前的TrialSize和R自带程序中没有成对样本频度检验的功效分析功能，代码为：

library(MESS)

power_mcnemar_test(n＝？,paid＝？,psi＝？,sig.level＝0.1)

其中：n为每组小样方块的个数；psi＝p01/p10或psi＝p10/p01，psi≥1；paid＝p01+p10；代入数值后，上面两条代码复制粘贴到RGui，回车得到统计功效power。

本发明具有的优点及有益效果如下：

1、本发明一种基于R软件的入侵植物频度数量变化调查分析方法，只记录小样方块中目标植物的有无，外业调查不受季节限制，不同调查人员之间数据重复性高，是监测入侵植物威胁的好方法。

2、使用沿样线规则等距间隔6～7m设置小样方块的频度调查方法，分散性好，调查覆盖面积大，对入侵植物空间分布变化很敏感。

3、采用预先进行初步调查，根据预调查结果调整小样方块的尺寸和选择样本大小，既不多调查也不少调查，以有限的成本产生显著水平和统计功效两方面均符合要求的调查结果。

4、根据实地调查结果和文献资料对比，从多种程序包选择适合入侵植物频度变化调查分析的方法；将只有医学试验案例的程序包用于解决植物生态学中频度数据成对样本卡方检验问题，有创新力度。

5、为入侵植物频度数据的小样方块尺寸调整、单样本精确度计算、独立两样本差异卡方检验、成对两样本McNemar检验提供了范例。

6、所有计算均编辑成R代码，使用时只需代入实测数值、复制粘贴代码，不需掌握高深的R语言；本发明技术含量高，使用简便。

附图说明

图1：为入侵植物频度数据调查样方块设置图。

具体实施方式

下面结合具体实施例对本发明做进一步说明。

实施例1：

小样方块尺寸的调整。

在平果柚木山调查入侵植物假臭草频度数量变化，预调查64个1m²小样方块(图1)，发现假臭草频度为8％，为了预留空间观测假臭草后续数量变化，需要将小样方块尺寸调整到多大才能使假臭草频度控制为50％，通过以下代码计算：

(log(100)-log(100-50))/(-log(1-8/100))

计算结果为8.31m²,需将小样方块面积调整为8.31m²。

实施例2：

监测入侵植物单样本变化，计算达到给定精确度所需的样本数。

入侵植物假臭草在平果陇兀屯垭口山分布较多，以系统抽样法设置1m×1m小样方块(图1)，预调查64个小样方块，发现38个小样方块有假臭草，26个小样方块无假臭草，假臭草比例频度为38/64＝0.59。如果期望能够以80％的把握、在0.05的显著水平上监测到10％的假臭草频度变化，需要调查多少个小样方块？

此例中，因p1＝0.59，|p1-p2|＝0.1，所以p2＝0.49；另外，α＝0.05,power＝0.80。

打开R 3.4.1软件，在其RGui视窗的“>”提示符后以英文状态输入以下两条代码(或在word中编辑后直接复制粘贴到RGui)：

library(pwr)

pwr.p.test(h＝ES.h(0.59,0.49),sig.level＝0.05,power＝0.8)

计算结果为n＝194.2975，需要调查195个1m×1m小样方块才能以0.8的功效、在α＝0.05显著水平上监测假臭草小至10％的频度变化。

实施例3：

检验入侵植物临时样地两样本频度差异显著性。

在平果陇兀任豆山，以系统抽样法预调查64个1m×1m小样方块，发现有假臭草样方30个，无假臭草样方34个。我们期望第二年能够以90％的把握检测假臭草10％的频度变化，可以接受10％的误报率，需要调查多少个小样方块？

预调查的p1＝30/64＝0.47；p2＝0.47+0.10＝0.57，α＝0.10,power＝0.90。将以下代码复制粘贴到RGui得到样本数n＝426。

power.prop.test(p1＝0.47,p2＝0.57,sig.level＝0.1,power＝0.9)

因此，我们在该地点继续以系统抽样法扩大范围调查，第一年调查426个1m×1m小样方块，其中205个有假臭草，221个无假臭草，样地没有打桩作永久标记。第二年在原地点重新拉样线，以系统抽样法调查426个1m×1m小样方块，其中247个有假臭草，179个无假臭草。在这两年内，该地点的假臭草数量是否发生显著变化？将以下代码输入RGui：

chisq.test(matrix(c(205,221,247,179),nrow＝2))

卡方检验得到χ²＝7.92，p＝0.005。再用以下代码检验统计功效：

power.prop.test(p1＝205/426,p2＝247/426,n＝426,sig.level＝0.1)

得到power＝0.89。结论：我们有89％的把握在α＝0.1的显著水平上认为陇兀任豆山的假臭草数量在这两年里极显著增加。

实施例4：

检验入侵植物固定样地成对样本频度差异显著性。

一年生入侵植物适合用临时样地调查，多年生入侵植物可以用固定样地调查。本例用固定样地调查平果陇兀顶果山入侵植物飞机草2015～2016频度数量变化；样方分布同图1。2015年5月调查64个1m×1m小样方块，有飞机草样方23个，无飞机草样方41个；我们可以接受10％的误报率和10％的漏报率；估计原来有飞机草样方的1/6在第二年无飞机草，而原来无飞机草样方的1/4在第二年有飞机草；需要设置多少个固定小样方块？

p01＝41/4/64＝0.16；p10＝23/6/64＝0.06；paid＝p01+p10＝0.16+0.06＝0.22；

输入以下代码计算：

library(TrialSize)

(McNemar.Test(0.10,0.10,2.674,0.22))

计算得到所需小样方块为184。用PVC管打桩标示样线位置和样方起始点，2015年5月调查184个小样方块，2016年5月再次调查同样位置的184个小样方块。两次调查都有飞机草的小样方块64个，第一年有第二年无的小样方块3个，第一年无第二年有的小样方块8个，两次调查都无飞机草的小样方块109个。用McNemar检验飞机草数量变化是否显著：

mcnemar.test(matrix(c(64,3,8,109),nrow＝2))

运算得：χ²＝1.4545，p＝0.2278。再进行功效检验：

library(MESS)

power_mcnemar_test(n＝184,paid＝(3/184+8/184),psi＝8/3,sig.level＝0.1)

运算得power＝0.898。结论：我们有90％的把握在α＝0.1的显著水平上认为2015年至2016年陇兀顶果山飞机草数量增加不显著。

Claims

1.一种基于R软件的入侵植物频度数量变化调查分析方法，其特征在于：主要操作步骤如下：

(1)在计算机上安装R软件及所需软件包：在https://www.R-project.org下载R 3.4.1或更高版本，根据提示直接运行安装，安装完成后，双击桌面R 3.4.1图标打开RGui视窗，从“程序包”菜单的下拉菜单中选择“安装程序包”，在弹出的“Secure CRAN mirrors”中选择“China”后“确定”，再在弹出的“Packages”中选择“MESS”后“确定”，安装好MESS程序包；使用类似过程安装“pwr”、“TrialSize”程序包；

(2)样地预调查：选择典型样地，以系统抽样法进行入侵植物频度数据的初步调查；调查以样线沿山脚——山顶方向设置，样线与样线平行且间隔7m；并以1m×1m小样方块沿样线分布，随机起点、等距间隔6m；记录64个1m×1m小样方块有无某一入侵植物之后，检查入侵植物出现的频度是否在20％～80％的范围内，如果频度低于或高于此范围，则需调整小样方块尺寸后重新开始预调查；

(3)单样本监测与达到给定精确度所需样本数计算：利用pwr程序包(Clay Ford，2017)进行计算：

library(pwr)

pwr.p.test(h＝ES.h(p1,p2),sig.level＝0.1,power＝0.9)

将预调查中比例频度p1＝有/(有+无)和给定精确度|p1-p2|中的比例频度p2代入上式，之后，将两条代码复制粘贴到RGui，回车得到样本数n；

(4)两样本临时样地频度比较所需样本数计算：使用R软件自带功效分析函数计算，代码为：

power.prop.test(p1＝？,p2＝？,sig.level＝0.1,power＝0.9)

将预调查中两样本入侵植物的比例频度p1和p2代入上式，在RGui中运算得到每个组所需样本数n；

(5)两样本临时样地入侵植物频度差异显著性检测。临时样地两样本的调查数据分布如表1，采用卡方检验：

chisq.test(matrix(c(C11,C10,C21,C20),nrow＝2))

两样本临时样地入侵植物有无计数列表中的4个数值为样方个数计数；

(6)两样本临时样地入侵植物频度变化检验的功效分析使用如下代码：

power.prop.test(p1＝？,p2＝？,n＝？,sig.level＝0.1)

其中，p1、p2为两样本入侵植物比例频度，n为最终实测的样本大小，这3个数值代入上式后，复制粘贴到RGui，回车得到统计功效power。

(7)两样本固定样地重复观测频度变化所需样本数的计算中，选择临床医学程序包TrialSize，计算得到的样本数最小也最合理。

(8)两样本固定样地重复观测入侵植物频度变化显著性检验，使用McNemar检验，数据收集方法见表2：

表2：固定样地两次重复观测数据表

两样本固定样地频度差异检验的代码如下：

mcnemar.test(matrix(c(p11,p10,p01,p00),nrow＝2))

另外，也可以采用按表2分列的计数数据。

2.根据权利要求1所述的一种基于R软件的入侵植物频度数量变化调查分析方法，其特征在于：步骤(2)所述调整小样方块尺寸是使入侵植物频度控制在40％～60％范围内，为观察其后续变化预留较大的空间；以p1代表用1m×1m小样方块预调查所得的某一植物百分频度，p1＝((有/(有+无))×100)，若要将这一植物的百分频度调整为p2，则调整后样方面积为：

(log(100)-log(100-p2))/(-log(1-p1/100))。

3.根据权利要求1所述的一种基于R软件的入侵植物频度数量变化调查分析方法，其特征在于：步骤(3)、(4)、(5)、(7)所述的样本计算或功效分析中，是选择α＝0.10、β＝0.10作为一般常用的显著水平和功效大小，兼顾入侵植物监测中的误报与漏报两个方面；但上述步骤也可用于另外指定的α和β值。

4.根据权利要求1所述的一种基于R软件的入侵植物频度数量变化调查分析方法，其特征在于：步骤(8)所述两样本固定样地重复观测频度变化显著性检验中，是选择MESS程序包来进行两样本固定样地重复观测入侵植物频度变化检验的功效分析，代码为：

library(MESS)

power_mcnemar_test(n＝？,paid＝？,psi＝？,sig.level＝0.1)

其中：n为每组小样方块的个数；psi＝p01/p10或psi＝p10/p01，psi≥1；paid＝p01+p10。