下载 - 景观生态与可持续性科学研究中心

Download Report

Transcript 下载 - 景观生态与可持续性科学研究中心

张海成
李希希
张艳聪
主要内容
1
基本概念原理
2
生态学问题
3
统计学解答
2
一 基本概念原理
1.1 引入
1.2 一元线性回归
1.3 多元线性回归
1.4 非线性回归
张海成
3
1.1
引入
考察变量之间关系的统计方法:
自变量
连续
分类
混合
连续
回归分析
方差分析
(ANOVA)
协方差分析
(ANCOVA)
分类
逻辑斯蒂回归
分对数模型
逻辑斯谛回归
或分对数模型
响应变量
连续变量: 相关分析、回归分析
4
回归
相关
 不必区分自变量与因变量
 区分自变量和因变量
 描述两变量数量上联系的密  自变量对因变量数量上的影
切程度
响程度
 所研究变量属于对等的关系  所研究变量属于非对等关系
5
回归分析法:
通过实验和观测建立变量之间关系的一种数理统计方法
线性回归分析
一元回归分析
回归分析
多元回归分析
非线性回归分析
6
主要目的:
在于了解自变量对因变量影响
内容:
 探索和确定变量之间的相关关系和相关程度
 建立回归模型,检验变量之间的相关程度
 应用回归模型进行估计和预测等
7
1.2 一元线性回归
Y   0  1x  
Y因变量; x自变量;β0, β1回归系数,ε随机误差
适用条件:
无偏性: E(ε) = 0
正态性, ε 为服从正态分布的随机变量
方差齐性:各组ε的方差都相同,即 ε 与自变量、因变量之间相互独立
独立性:各组ε之间相互独立,且COV(ε i, ε j) = 0, (i ≠ j)
8
参数β0, β1 的普通最小二乘估计
基本思想:回归直线与所有样本数据点距离最近
n
Q (  0 , 1 ) 

( yi   0   1 xi )
2
i 1
 y  ˆ
i

2
n
0
 ˆ1 x i
2
n
 min Q (  0 ,  1 )  min
i 1
  y  
i
0
  1 x i 
i 1
体重—身高
9
n
 Q
  2   y i    0   1 x i   0

  0
i 1

n
Q

  2   y i    0   1 x i x i  0
   1
i 1

 n  0   x i  1   y i

2

  x i  0  (  x i )  1 
ˆ1 
ˆ 0 
n  xi yi 
n x
y
n
i
2
i
x
x
  x 
 ˆ1
i
2
i
x
n
yi
对β0,β1分别求偏导
且令两偏导等于0
整理后得到
i
yi

 ( x  x )( y 
 (x  x)
i
i
y)
2
i
求解
i
 y  ˆ1 x
10
回归方程拟合优度的评价
拟合优度: 样本观测值聚集在样本回归线周围的紧密程度
决定系数:R²

n
SST 
 (y
i
 y)
2
i 1
 y
2
n
i 1
i
 y

2
n

 y
i
 yˆ i  
i 1
SST
SSE
总离差平方和 残差平方和
  yˆ
2
n
i
 y

i 1
SSR
回归平方和
11
R 
2
SSR
SST
 yˆ


 y
SSE 
i

 y
 y
2
2
i
yˆ
2
ˆ
(
y

y
)

SSR 
 ( yˆ  y )
2
SST 
 ( y  y)
2
y
R²表示全部偏差中有百分之多少可由x和y
的回归关系来解释
12
回归方程的显著性检验
1
• 对整个回归方程的显著性检验(F检验)
2
• 对各回归系数的显著性检验(t检验)
13
对整个回归方程的显著性检验
① 提出假设:H0:回归方程不显著(βi = 0); H1: 回归方程显著(βi不全为0)
② F-检验(ANOVA)
方差来源 平方和 自由度
回归
SSR
1
均方
MSR 
F值
SSR
1
误差
SSE
n-2
总计
SST
n-1
MSE 
SSE
n2
F 
MSR
MSE
③ 给出显著性水平α,确定临界值F(1, n-2)
④ 若F ≥ Fα (1, n-2), 拒绝H0,说明回归方程总体显著;反之
14
对回归系数的显著性检验
① 提出建设:H0: 回归系数不显著(β1 = 0 ); H1:回归系数显著(β1 ≠0)
② t-检验
^
t
1
t 统计量
S1
S1 
S yx
Var ( ˆ1 ) 

S yx  ˆ 
ˆ 
2
y
2

(x  x)
2
 ˆ0  y  ˆ1  xy
n2
 (y
^
i
 yi )
n2
2

e
2
i
n2
S1:回归系数 ˆ1 的标准差
Syx, σ: 估计标准误差
σ²: 误差项方差
③ 给定显著性水平α,确定临界值tα/2(n-2)
④ 若 |t|≥ tα/2(n-2), 则拒绝H0,认为回归系数显著;反之
15
两条或多条回归方程的差异性比较
例:1)A,B两样地观测数据得到的回归方程是否有显著性差异
2)男、女生身高与体重的回归方程是否有显著性差异
Y = a 1X + b 1 + ε 1
Y = a 2X + b 2 + ε2
…
Y = a n X + b n + εn
原理: 检验斜率、残差、截距三项
1)残差ε : F-检验
2)斜率a,截距b使用其他复杂统计量,SPSS以
协方差分析进行检验(见ANCOVA-63)
SPSS具体处理步骤:
多用协方差分析检验截距和斜率的差异,以SPSS为例:1.先重新整理数据,将y2数据列加到y1下面,变成一个变量y;将x2数据列加到
x1下面,变成一个变量x;然后再设定一个新的分组变量group,原来第1组值为1,第2组值为2.
2.进行协方差分析(第一步分析斜率是否无差异).
Analyze->General Linear Model->Univariate
Dependent List:填入y---------将y做为因变量
Fixed Factor:填入group
Covaraites:填入x--------将x做为协变量
Model:选Custom
Model:填入 x group x*group---------注意如果变量填入顺序不一样,结果也会不一样. (注意:x*group是通过双选得到)
Sum of squares下拉列表框:选TypeI
然后点击ok,看结果里x*group这一行的Sig.P值,若大于0.05,则接受原假设,即两条回归直线的斜率无差异,否则拒绝.
3.再来进行截距的无差异分析
其实过程跟上面一样,只是Model里去掉了x*group交叉项.
Analyze->General Linear Model->Univariate
Dependent List:填入y---------将y做为因变量
Fixed Factor:填入group
Covaraites:填入x--------将x做为协变量
Model:选Custom
Model:填入 x group ---------注意如果变量填入顺序不一样,结果也会不一样.
Sum of squares下拉列表框:选TypeI
点击ok后,看group一行的Sig.P值,若P值大于0.05说明两条回归直线截距也无差异,若小于0.05说明截距是有差异的.
PS:paired test是student-t检验的一种,但检验的结论是两组数据所来源总体的差异。所以判断两回归曲线的差异,需用协方差分
析或单独的检验统计两来判断截距和斜率的差异。
回归方程的估计与预测
点估计
回归方程估计
置信区间估计
19
回归方程的估计与预测
点估计
对于给定的 X 值,求出 Y 平均值的一个估计值或
Y 的一个个别值的预测值。
y

x0
yˆ  a  bx
x
20
回归方程的估计与预测
利用估计的回归方程,对于自变量 x 的一个给
置信区间估计 定值 x0 ,求出因变量 y 的平均值E(y0)的估计区
间 ,这一估计区间称为置信区间
对于给定的 x = x0 ,Y 的1-置信区间为:
yˆ 0  t 2
yˆ
自由度为n-2的 t 分布
的  水平双侧分位数
21
回归方程的估计与预测
置信区间估计
Sy: 回归估计标准差
22
1.3 多元线性回归
Y   0   1 x 1   2 x 2  ......   p x p  
(假定:x1,x2…xp无交互作用)



Y 




y1 

y2 
 

y n 
1

1
X 


1

x11
x12

x 21
x 22




x n1
xn 2

x1 p 

x2 p 
 

x np 
 0 



 1
 
 


 
 p
ε=
 1 
 
2 
  
 
 
 n
Y:因变量;
ε :随机误差
x1, x2, … xp: 自变量; Β0, β1, … βp: 回归系数
适用条件:
无偏性E(ε) = 0
服从正态分布
方差齐性
独立性
23
多元线性回归模型的参数估计
基本思想:使估计值
yˆ
方法:最小二乘法
与观测值y之间的残差在所有样本点上达到最小
(详见路径分析)
Q最小
步骤:
 分别对β0, β1…… βp求偏导
 令各偏导等于0
 求得β0, β1…… βp
24
拟合优度检验
总偏差 (SST) = 回归偏差 (SSR) + 剩余偏差 (SSE)
R 
2
SSR
SST

  yˆ
 y
i

 y
 y
2
2
i
25
回归方程显著性检验
对整个回归模型的显著性检验(似一元回归模型)
① 提出假设,H0:回归模型不显著(β1 = β2……= βp = 0);H1:回归
模型显著( β1 , β2…… βp 不全为0)
② 建立F统计量
方差来源
平方和
自由度
回归
SSR
p
误差
SSE
n-p-1
总计
SST
n-1
均方和
MSR 
F值
SSR
p
MSE 
SSE
n  p 1
F 
MSR
MSE
③ 给出显著性水平α,查F分布表,得临界值Fα (p, n-p-1)
④ F ≥ Fα (p, n-p-1), 则拒绝H0,回归模型显著
26
回归系数显著性检验
① 提出假设: H0: βi = 0; H1:βi = 0;
② 建立t检验计算公式:③ ④
^
t i 
Si 
i
Si
Var ( ˆ i ) 
c ii ˆ
Si:回归系数标准差;Cii是(XTX)1中第i+1个主对角线元素
③ 给出显著性水平α,查t分布表,得临界值tα / 2 (n-p-1)
④ |t| ≥ tα / 2 (n-p-1), 则拒绝H0,回归系数显著
27
1.4 非线性回归模型
线性模型:
y   0   1 x 1   2 x 2  ......   p x p  
yi 是解释变量xi 的线性函数
yi 是相应参数βi 的线性函数
28
非线性模型:
 建立因变量与自变量之间的非线性关系
 无假设条件,可建立任意形式的模型
本质线性模型
非线性模型
本质非线性模型
29
本质线性模型
30
模型的选取
 了解数据特点或可辨别变量间的相关关系
——直接选择函数模型
B = pDq
 数据认识不充分,不能辨别变量间的相关关系
—— 绘制散点图,根据图形特点选择模型形式
 对于无法直接从散点图看出应选取何种回归模型的
——尝试多种可能的模型,并通过R²的大小、图形
本身的比较选择出最优模型
31
函数转换法不足:
 可能导致更复杂的计算
 可能导致数据失真
 不是所有非线性模型均可转换为线性模型
32
本质非线性模型
无法通过变量变换转化为线性模型
33
非线性回归模型解法

最小二乘法(模型主体部分占优势)
Gauss-Newton算法
Newton-Raphson算法
 最大似然法(模型随机误差部分占优势)
1) 估算模型中参数的起始值及其取值范围
2) 利用迭代算法寻找使得残差平方和最小的参数估计值
Newton-Raphson算法
计数算法
BHHH算法
34
二 生态学问题
2.1 生态学问题
2.2 数学模型
李希希
35
动物生态学的基本知识
动物生态学(animal ecology):是从生物种群和
群落的角度研究动物与其周围环境相互关系的科学,
是生态学的分支,是由动物学与生态学等交叉形成
的学科。
36
研究内容
(1)阐明动物与生存条件的关系,生存条件的变化对动
物的生理结构、形态特征和行为方式的影响;
(2)研究在一定的生存条件下各种动物种群的数量关系,
出生率和死亡率的变化,种群密度和年龄分布;
(3)研究一定的环境条件下种内和种间关系以及它们对
动物进化的意义,种内与种间的合作与竞争,捕食-被捕
食,种间各种共生关系,以及动物种群的结构和演化;
(4)研究不同生态条件下动物种群和群落的形成、适应
性和演化;
(5)人类对动物资源开发利用和动物遗传资源的保护等。
37
2.1 生态学问题
自然界中有很多例子能证明捕食者具有
控制其捕食猎物密度的功能。
例如:1、1950年~
1960年间北美五大洲
中的七鳃鳗几乎导致
湖红点鲑的灭绝;
2、在人类开
始捕杀狼之前,黄石
地区的狼使得当地麋
鹿的密度受到控制。
38
功能响应(Functional Response):一个捕食者个
体所能猎杀的猎物数量(或者一个寄生物能攻击的寄
主数量)是关于猎物密度的函数。
一定时间内被猎杀的猎物数随着猎物密度的增加
逐渐接近一个渐近线。
39
第Ⅰ类功能响应
※捕食者随机搜寻猎物
※搜寻时间为一常数
※食欲有限
※不考虑猎物被捕食者处理所
需的时间。
被捕食猎物数在一定范围内随
猎物种群的增加而直线增加,
超过则不再增加。
被猎杀猎物比例保持最大(非
密度依赖型),继而降低。
其中,Ne是被捕食猎物数;
N是猎物密度。
40
第Ⅱ类功能响应
※捕食者随机搜寻猎物
※搜寻时间与处理时间成反比
※捕食者食欲有限
捕食的猎物数随猎物密度的增加,
以双曲线型接近渐近线,超过这
一范围将不再增加
被猎杀猎物的比例渐进型减少
(反密度依赖型)。
T 表示总时间。
Th表示每个捕食者处理每个猎物所需的时
间(其中包括用于吞食猎物和无法攻击其
他猎物的所有时间)。
41
第Ⅲ类功能响应
※取食猎物的效率
~食量限制
~猎物种群低密度下的敏感
程度
※猎物数与猎物种群密度为
“S”型曲线
猎物被猎杀比例逐渐增加至S型
曲线的拐点(密度依赖型),然
后开始减少。
42
研究功能响应涉及的相关问题
1. 功能响应的类型是什么?(比较C和E、D和F)
2. 功能响应的机制模型中系数的最佳估计是什么?
3. 描述两种功能响应的模型中系数是否有显著的不同?
43
• 模型选择——包括使用对被捕杀猎物与猎物总数之
比的逻辑斯蒂回归方法
• 假设检验——包括对被杀猎物数量和所提供猎物数
量比值的非线性最小二乘回归来估计功能响应的参
数以及比较不同功能响应的参数。
44
2.2 功能响应的数学模型
• 与功能响应相关的因素
1. 被捕食量,用Ne表示。
2. 攻击系数(瞬时的搜寻速率)使得猎物的遇
到概率和猎物密度相关联。用a表示。
3. 猎物密度,用N表示。
4. 从开始捕食到捕到猎物所用的时间。用Ts表示。
5. 总时间,用T表示。
6. 每个捕食者处理每个猎物所需的时间(其中包
括用于吞食猎物和无法攻击其他猎物的所有时
间),用Th表示。
Ne = a*TsN
①
45
Ne = a*TsN
• 将Ts考虑为一个变量,又根据:
Ts =T - ThNe
• 将 ② 代入
①
①
②
中,变形可得:
公式2: Ne = aNT/(1+aNTh)
46
第Ⅰ类功能响应
公式2:
Ne = aNT/(1+aNTh)
当Th=0, 即不考虑猎物的处理时间时:
公式1:
Ne = aNT
47
Disc方程模拟第Ⅱ类功能响应
模型假定:
1.渐近线由时限决定;
2.捕食者遇到猎物概率和猎物密度线性相关;
3.当捕食者在处理猎物时不再进行其他捕获;
4.猎物密度恒定。
48
第Ⅱ类功能响应
公式2:
Ne = aNT/(1+aNTh)
1. 被捕食量,用Ne表示。
2. 攻击系数(瞬时的搜寻速率)使得
猎物的遇到概率和猎物密度相关联。
用a表示。
3. 猎物密度,用N表示。
4. 总时间,用T表示。
5. 每个捕食者处理每个猎物所需的时
间(其中包括用于吞食猎物和无法
攻击其他猎物的所有时间),用Th
表示。
49
存在猎物消耗时描述功能响应的适合模型
ɑ关于N的双曲线型函数
a 
d  bN
1  cN
③
b,c和d均为常数
50
式③产生的各种形式的ɑ和N之间的关系
a 
等于零的系数
d  bN
1  cN
③
ɑ的计算公式
ɑ与N之间的关系
无
(d+bN)∕(1+cN)
增加至渐近线;截距为d,不等于
零
c
d+bN
线性增加;截距为d,不等于零
d
bN∕(1+cN)
增加至渐近线;截距等于零
d和c
bN
线性增加;截距等于零
ɑ 假定功能响应已知是Ⅲ型;因而至少有b>0
51
猎物消耗的功能响应
③
a 
d  bN
1  cN
代入
aN T
Ne 
1  aN Th
dN T  bN T
2
Ne 
1  cN  d N T b  b N T h
2
52
第Ⅲ类功能响应
当a是一个N的增函数时,就会出现第Ⅲ类功能响应:
dN T  bN T
2
Ne 
1  cN  d N T b  b N T h
2
53
上述模型假定猎物密度N恒定,而实
际上N往往随实验的进行而降低。当N
减少时,描述功能响应的合适模型应是
对公式2和3的时间积分式,从而可应对
变化的猎物密度。
54
Ne 
aN T
1  aN Th
N e  N 0 1 exp  a (T h N e  T ) 
55
对于第Ⅲ类功能响应,当ɑ如式 ③ 中所示是初始密
度的函数时为最简单的形式:
③
a 
d  bN 0
1  cN 0
N e  N 0 1  exp  a (T h N e  T ) 
N e  N 0 1  exp  ( d  bN 0 )(T h N e  T ) / (1  cN 0 ) 
56
总结
1.功能响应(Functional Response):一个捕食
者个体所能猎杀的猎物数量(或者一个寄生物能攻
击的寄主数量)是关于猎物密度的函数。
实验时猎物密度N恒定
类型II
类型III
Ne 
aNT
1  aNT h
dNT  bN T
2
Ne 
1  cN  dNT h  bN T h
2
a: 攻击系数
b,c,d: 决定攻击系数a
Th:处理每个猎物所需的时 Th:处理每个猎物所需的时间
间
实验时猎物密度N不恒定(初始N0)
类型II
类型III
Ne  N 0 {1  exp[ a (Th N e  T )]} N e  N 0 {1  exp[( d  bN 0 )( Th N e  T ) /(1  cN 0 )]}
a: 攻击系数
Th:处理每个猎物所需的
时间
b,c,d: 决定攻击系数a
Th:处理每个猎物所需的时间
59
3.1 统计学问题解答
3.2 实例分析
张艳聪
60
1. 功能响应的类型是什么?(比较C和E 或 D和F)
2. 功能响应的机制模型中系数的最佳估计是什么?
3. 描述两种功能响应的模型中系数是否有显著的不同?
61
62
3.1.1 模型选择:确定功能响应形状
4
1) 求解过程(最大似然法)
2) 求解过程(区分相应曲线)
3.1.2 假设检验:参数估计和比较
1)
方法(最小二乘法)
2) 求解过程(牛顿迭代法)
63
变量特征:
因变量Y取值:Y=1 死亡,Y=0 活着, 0-1 变量
自变量N:
N0 猎物初始密度,
为连续变量
方法:被杀猎物和现有猎物比例的逻辑斯蒂回归
描述离散型因变量和连续型自变量的关系的回归模型
自变量(N)
因变量(Y)
连续型
N∈[0,+∞)
离散型
Y=0,1
Y=1的概率
(μY·N)
连续型
μY·N=P(Y=1)
∈[0,1]
μY·N的Logistic
变换
连续型
Q∈(-∞,+∞)
64
Logistic变换关系
被捕食的概率
猎物密度
P=f(N)
P  [ 0 ,1]
N  [ 0 ,  )
?
Logistic 变换
Q  Logit ( p )  ln
P
1 P
 (  ,  )
N  [ 0 ,  )
Q  f ( N )  a 0  a1 N  a 2 N  a 3 N    a z N
2
Q  ln
P
e
Q

z
P
1 P
exp( a 0  a1 N  a 2 N  a 3 N    a z N )
2
Q
1 e
3
3
z
1  exp( a 0  a1 N  a 2 N  a 3 N    a z N )
2
3
z
65
最大似然法:
(1)最大似然估计是一种统计方法,它用来求一个样本
集的相关概率密度函数的参数。
(2)当从模型总体随机抽取n组样本观测值后,最合理
的参数估计量应该使得从模型中抽取该n组样本观测值的
概率最大。
66
Logistic回归方程:
被捕食: Pr ob {Yi  1} 
存活:
Pr ob {Yi  0}  1 -
e
1 e
e
Q
Q
1 e
Q
[exp( a 0  a1 N  a 2 N    a z N )]
2
Q


z
Yi
1  exp( a 0  a1 N  a 2 N    a z N )
2
z
[exp( a 0  a1 N  a 2 N    a z N )]
2
1
1 e
Q

z
Yi
1  exp( a 0  a1 N  a 2 N    a z N )
2
z
(其中a0, a1, a2 … az是需要估计的系数)
因变量Y为(0, 1)变量
n
L
[exp( a 0  a1 N  a 2 N  a 3 N    a z N )]
 1  exp( a
i 1
利用迭代法寻找参数数值
使得似然函数L最大化
2
3
Yi
 a1 N  a 2 N  a 3 N    a z N )
2
0
z
3
z
67
例题1:假设B组每只老鼠被捕食的概率为P,
则同时出现以下10项观察值的概率L是多少?
e
P ( Y i  1) 
Q
1 e
Q

)
Yi
1 e
Q
(e
Q
P (Y i  0 )  1  P (Y i  1) 
1
1 e
Q

)
Yi
1 e
Q
(e
Q
n
 L  P (Y1 )  P (Y 2 )    P (Y n ) 
 P (Y )
i
i 1
n

Q Yi
(e )
 1 e
Q
i 1
n

[exp( a 0  a1 N  a 2 N  a 3 N    a z N )]
 1  exp( a
i 1
2
3
Yi
 a1 N  a 2 N  a 3 N    a z N )
2
0
z
3
z
68
n
L 
[exp( a 0  a1 N  a 2 N
 1  exp(
i 1
2
 a3 N
a 0  a1 N  a 2 N
2
3
   a z N )]
 a3 N
z
3
Yi
   az N )
z
对于给定的自变量N值,当参数a0, a1, a2, …, az 取不同
的值时,L大小会发生变化,表明不同参数估计对实际结果
的拟合情况。L越大,说明产生实际结果的概率越大,即与
实际结果拟合得越好。
当参数a0, a1, a2, …, az 取某些特定值,使得L最大时,表
明这些参数使得实际结果发生的概率最大,说明这些参数是
拟合模型的最佳值。使用得到实际观测值的概率来作为拟合
参数优劣标准的方法,就是最大似然估计。
69
使用Logistic回归区分II、III型功能响应
Ne ~ N
难区分
70
使用Logistic回归区分II、III型功能响应
P=Ne/N ~ N
P ( 0 )  0
P ( 0 )  0
容易区分
71
Ne
 P
N
e
exp( a 0  a1 N  a 2 N  a 3 N    a z N )
2
Q
1 e
Q

2
( e )   (1  e )  e  (1  e ) 
Q
Q
(1  e )
Q
2

2
Q ( 0 ) e
(1  e
Q (0)
Q (0)
)
2

a1e
3
Q e
Q
z
Q
(1  e )
Q   a 1  2 a 2 N  3 a 3 N    za z N
 P ( 0 ) 
z
1  exp( a 0  a1 N  a 2 N  a 3 N    a z N )
Q
P 
3
Q
2
z 1
a0
(1  e )
a0
2
 a1  0  P ( 0 )  0  III 型曲线
a1  0  P ( 0 )  0  II 型曲线
72
3.1.1 模型选择:确定功能响应形状
4
1) 求解过程(最大似然法)
2) 求解过程(区分相应曲线)
3.1.2 假设检验:参数估计和比较
1)
方法(最小二乘法)
2) 求解过程(牛顿迭代法)
73
方法:最小二乘法
寻找特定的参数使得模型预测值与观测值的离差平
方和最小,以此作为拟合参数优劣的标准。
L
 (y
i
 yˆ i )
实际观测值
2
理论预测值
74
实验时猎物密度N恒定
类型II
类型III
Ne 
aNT
1  aNT h
dNT  bN T
2
Ne 
1  cN  dNT h  bN T h
2
a:攻击系数
b,c,d: 决定攻击系数a
Th:处理每个猎物所需的时 Th:处理每个猎物所需的时间
间
实验时猎物密度N不恒定(初始N0)
类型II
类型III
Ne  N 0 {1  exp[ a (Th N e  T )]} N e  N 0 {1  exp[( d  bN 0 )( Th N e  T ) /(1  cN 0 )]}
a: 攻击系数
Th:处理每个猎物所需的
时间
b,c,d: 决定攻击系数a
Th:处理每个猎物所需的时间
75
例子:II型曲线
Nˆ e i 
(1)猎物密度N恒定
L

i
2
( Ne i  Nˆ e i ) 

i
aN i T
1  aN i T h
( Ne i 
aN i T
1  aN i T h
)
2
(2)猎物密度N不恒定 Nˆ ei  N 0 {1  exp[ a (T h Nˆ ei  T )]}
i
Nˆ e i  ?
76
使用牛顿迭代法求
Nˆ e  ?
解关于 Nˆ e 的方程:
Nˆ e  N 0 {1  exp[ a (T h Nˆ e  T )]}
令 f ( N e )  N 0  N 0 exp[ a (T h N e  T )]  N e,当函数
f ( N e )  0 时的 N e 值即为需要求的
Nˆ e
77
k  f ( Ne 1 ) 
( Ne 2 , f ( Ne 2 ))
Nˆ e
( Ne 2 , 0 )
0  f ( Ne 1 )
Ne 2  Ne 1
k  f ( Ne 2 ) 
0  f ( Ne 2 )
Ne 3  Ne 2
( Ne 3 , 0 )
( Ne 1 , 0 )
( Ne 1 , f ( Ne 1 ))
78
最小二乘法
最大似然法
拟合回归模型的参数
目的
模型自变量(X)
连续变量
连续变量
离散变量
因变量(Y)
(0-1型分布)
连续变量,一般 Y=1的概率,∈[0,1]
条件平均数 
Y X
∈[0,+∞]
函数
L

( y i  yˆ i )
2
L

Pˆ ( y  y i )
最优标准
参数取值使得L最小
参数取值使得L最大
统计意义
参数取值使得观测值与 参数取值使得利用模型预测
模型预测值的离差平方 产生实际观测值的概率最大
和最小,即为最优参数 ,即为最优参数,模型与实
,模型与实际最符合。
际最符合。
79
Notonecta捕食Asellus
捕食者:Notonecta glauca 中国大仰蝽属
猎物:Asellus aquaticus 栉水虱
1)猎物密度从5,7,8…..100等11个取值;
2)每个密度至少八个重复;
3)总共有89个观测值
4)持续72h。
Q1:功能响应的曲线是那种,Ⅱ或者Ⅲ?
Q2:最佳的参数估计是什么?
80
原始数据如右表:
N0 = initial number of prey
REP = replicate number,
FATE: 0 = prey eaten
1 = prey alive
NE = count of prey in each FATE
FATE=0时,
NE:为被捕食数目,即Ne
FATE=1时,
NE:为活着的数目
对一个重复,N0=ΣNE
N0
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
5
7
7
7
REP FATE NE
1
0
0
1
1
5
2
0
1
2
1
4
3
0
1
3
1
4
4
0
2
4
1
3
5
0
2
5
1
3
6
0
2
6
1
3
7
0
2
7
1
3
8
0
3
8
1
2
1
0
0
1
1
7
2
0
0
81
第一步:模型选择
(1)进行Logistic回归分析,建立回归方程
Ne
2
= P 
N0
3
exp (P0 + P1 N 0 + P2 N 0 + P3 N 0 )
2
3
1 + exp (P0 + P1 N 0 + P2 N 0 + P3 N 0 )
(2)根据回归方程估计相应曲线的形状
P ( 0 ) 
Q ( 0 ) e
(1  e
Q (0)
Q (0)
)
2

p1 e
p0
(1  e
p0
)
2
p1  0  P ( 0 )  0  III 型曲线
p1  0  P ( 0 )  0  II 型曲线
82
SAS界面
run
log
program
output
log
editor
83
84
SAS procedure CATMOD
/*对参数等于0进行最大似然检验,同时对似然比率进行检验来
检验模型的整体拟合度*/
DATA NOTONECT; /* name the database*/
INPUT N0 REP FATE NE;
N02=N0**2; /* initial number of prey squared */
N03=N0**3; /*initial number of prey cubed */
CARDS; /*data followed, ends with “;”*/
85
SAS procedure CATMOD
PROC CATMOD DATA=NOTONECT;
DIRECT N0 N02 N03;
不使用广义最
小二乘法估计
MODEL FATE = N0 N02 N03 / ML NOGLS NOPROFILE ;
POPULATION N0 REP;
最大似然法
WEIGHT NE; /*指定权重变量*/
不显示样本值和响应值
对一个重复,N0=ΣNE
86
The result of CATMOD procedure
Maximum Likelihood Analysis of Variance
Source
DF Chi-Square Pr > ChiSq
-------------------------------------------------Intercept
1
28.68
<.0001
N0
1
6.66
0.0098
N02
1
7.52
0.0061
N03
1
5.67
0.0172
Likelihood Ratio 85
206.82
<.0001
Analysis of Maximum Likelihood Estimates
Standard
ChiParameter Estimate
Error Square Pr > ChiSq
---------------------------------------------------------Intercept
-1.3457
0.2513
28.68
<.0001
N0
0.0478
0.0185
6.66
0.0098
N02
-0.00103 0.000376
7.52
0.0061
N03
5.267E-6 2.211E-6
5.67
0.0172
P ( 0 )  0
P ( 0 )  0
所有参数显著不等于0,但是模型整体拟合不佳
一次项系数大于0,可以判断应属于III类功能相应
87
DATA NOTONEC2; /* 修正数据,去除猎物存活的数据行,添加捕食比
例数据列*/
SET NOTONECT;
IF FATE= 0; PROPEAT= NE/N0;
观
测
值
估
计
值
PROC MEANS DATA=NOTONEC2; /*计算基本统计量:均值,SE等*/
BY N0 NOTSORTED; /*按照N0值分组且不排序*/
VAR PROPEAT; /*指定被食比例变量,计算该变量的均值、方差等*/
OUTPUT OUT=NOTOMEAN MEAN=MEANPROP; /*输出结果到
NOTOMEAN数据集中,输出每个N0对应的均值*/
DATA NOTONEC3; /* 产生估计值数据集,包括每组N0对应的预测值*/
SET NOTOMEAN;
K=EXP(-1.3457+(0.0478*N0)-(0.00103*N0**2)+(0.000005267*N0**3));
PRED=K/(1+K);
PROC PLOT DATA=NOTONEC3; /*绘制平均观测值与估计值散点图 */
PLOT PRED*N0='P'
MEANPROP*N0='*'/OVERLAY;
88
每个点代表在给定猎物初始密度下的被捕食概率
* : 实际观测值
P :拟合估计值
89
第二步:拟合数据并估计参数
非线性最小二乘法( SAS procedure NLIN)
应用条件:
(1)正态分布
(2) 方差齐性
选用存在猎物消耗的III类功能响应模型:
N e  N 0 {1  exp[( d  bN 0 )( Th N e  T ) /( 1  cN 0 )]}
H0:b=0 c=0 d=0 Th=0
N0: 猎物初始密度 Ne: 被食数量
T:可利用时间 Th:处理每个猎物所需时间
90
第二步:假设检验
N e  N 0 {1  exp[( d  bN 0 )( Th N e  T ) /( 1  cN 0 )]}
初始值的设定:
d: 许多情况下比较小,一般初值设置为0
b和c: 通常是选择覆盖两个数量级的不同初始值,
因为当攻击系数a (瞬时搜索速率)是一个N的增函数
时,出现第三类功能相响应,所以参数 b>0 ,c>=0,d>=0
Th: 猎物处理时间 Th>0
N0: 猎物初始密度
T:可利用时间 (72h)
Ne: 被食数量
Th:处理每个猎物所需时间
91
完整模型:
N e  N 0  N 0 exp[( d  bN 0 )( Th N e  T ) /(1  cN 0 )]
PROC NLIN DATA=NOTONEC2;
PARMS BHAT= 0.001 0.01 0.1
CHAT= 0.001 0.01 0.1
DHAT= 0 THHAT=3.0;
BOUNDS BHAT>0,CHAT>=0,
THHAT>0;
T=72;
X=NE;
A=(DHAT+BHAT*N0)/(1+CHAT*N0);
C1=EXP(-A*T);
C2=A*THHAT;
H=N0*C1*EXP(C2*X)+X-N0;
ITER=0;
待续……
92
完整模型:
续上一页
牛顿迭代法
DO WHILE(ABS(H)>0.0001 AND ITER<50);
X=X-H/(N0*C1*C2*EXP(C2*X)+1);
H=N0*C1*EXP(C2*X)+X-N0;
/*新的预测值Nei+1*/
/* 将Nei+1 带入函数求取f(Nei+1) */
ITER=ITER+1;
END;
MODEL NE=X;
93
Procedure NLIN: Full model
Nonlinear Least-Squares Summary Statistics
Dependent Variable NE
Sum of
Mean
Approx
Source
DF Squares
Square
F Value Pr > F
Model
3
9526.8
3175.6 157.57 <.0001
Error
86
1733.2
20.1535
Uncorrected Total
89 11260.0
Approx
Approximate 95%
Parameter Estimate
Std Error Confidence Limits
Label
BHAT
0.000416
0.000228
-0.00004
0.000869
CHAT
0
0
0
0
DHAT
0.000832
0.00382
-0.00676
0.00842
THHAT
4.0640
0.3575
3.3534
4.7747
Bound0
107.9
315.6
-509.1
725.0
0 <= CHAT
N e  N 0 {1  exp[( d  bN 0 )( Th N e  T ) /( 1  cN 0 )]}
N e  N 0 {1  exp[( d  bN 0 )( Th N e  T )]}
94
简化模型1:
PROC NLIN DATA=NOTONEC2;
PARMS BHAT=0.001 0.01 0.1
DHAT=0
THHAT=3.0;
BOUNDS BHAT>0,THHAT>0;
T=72;
X=NE;
A=(DHAT+BHAT*N0);
C1=EXP(-A*T);
C2=A*THHAT;
H=N0*C1*EXP(C2*X)+X-N0;
ITER=0;
待续……
95
简化模型1:
续上一页
DO WHILE(ABS(H)>0.0001 AND ITER<50);
X=X-H/(N0*C1*C2*EXP(C2*X)+1);
H=N0*C1*EXP(C2*X)+X-N0;
牛顿迭代法
ITER=ITER+1;
END;
MODEL NE=X;
96
Procedure NLIN: Reduced model 1
Nonlinear Least-Squares Summary Statistics
Dependent Variable NE
Sum of
Mean
Approx
Source
DF Squares
Square F Value
Pr > F
Model
3
9526.8
3175.6
157.57
<.0001
Error
86
1733.2
20.1535
Uncorrected Total
89 11260.0
Parameter
BHAT
DHAT
THHAT
Approx
Estimate Std Error Approximate 95% Confidence Limits
0.000416 0.000228 -0.00004
0.000869
0.000831 0.00382 -0.00676
0.00842
4.0641
0.3575
3.3534
4.7747
N e  N 0 {1  exp[( d  bN 0 )( Th N e  T )]}
N e  N 0 {1  exp[ bN 0 ( Th N e  T )]}
97
简化模型2:
PROC NLIN DATA=NOTONEC2;
PARMS BHAT= 0.001 0.01 0.1
THHAT= 3.0;
BOUNDS BHAT>0,THH>0;
T=72;
X=NE;
A=BHAT*N0;
C1=EXP(-A*T);
C2=A*THHAT;
H=N0*C1*EXP(C2*X)+X-N0;
ITER=0;
待续……
98
简化模型2:
续上一页
DO WHILE(ABS(H)>0.0001 AND ITER<50);
X=X-H/(N0*C1*C2*EXP(C2*X)+1);
H=N0*C1*EXP(C2*X)+X-N0;
牛顿迭代法
ITER=ITER+1;
END;
MODEL NE=X;
99
Procedure NLIN: Reduced model 2
Nonlinear Least-Squares Summary Statistics
Dependent Variable NE
Sum of
Mean
Approx
Source
DF Squares
Square
F Value
Pr > F
Model
2
9525.8
4762.9
238.93 <.0001
Error
87
1734.2
19.9338
Uncorrected Total
89
11260.0
Parameter
BHAT
THHAT
Estimate
0.000460
4.1043
Approx
Std Error Approximate 95% Confidence Limits
0.000104 0.000254
0.000666
0.2892
3.5294
4.6792
N e  N 0 {1  exp[( d  bN 0 )( Th N e  T )]}
N e  N 0 {1  exp[ 0.0005 N 0 ( 4.1043 N e  T )]}
100
小结
统计方法: 回归分析
线性回归
非线性回归
生态学问题:
应
用
捕食和功能响应
确定功能
响应类型
Logistic回归
功能响应参数
的估计和比较
非线性最小二乘法
101
作业
102
A林场有100公顷落叶松林地,林分密度为1000棵/公顷,整个林地落叶
松的平均胸径D=10cm,现随机伐取30棵落叶松样本,测得各样本胸径(D)及
材积(V)的数据如下表,
(1) 求胸径(D)和材积(V)的最佳回归关系(取线性、幂、指数、对数方程中
最优的)
(2) 求置信度α=0.05,D=10cm时的点估计,置信区间估计 (提示:可先转
化为线性,tα/2(28)=2.048)
(3)落叶松价格为800元/m3, 求当前A林场落叶松总经济价值。
样本N
胸径D(cm)
材积V(m3)
样本N
胸径D(cm)
材积V(m3)
1
8.3
0.0300
16
8.5
0.0388
2
6.7
0.0133
17
3.8
0.0035
3
4.2
0.0035
18
10.2
0.0703
4
9
0.0333
19
3.4
0.0008
5
8.6
0.0390
20
13.4
0.1388
6
8.1
0.0233
21
7.9
0.0368
7
4.5
0.0035
22
10.8
0.0708
8
6.8
0.0115
23
4.5
0.0105
9
4.8
0.0043
24
13.8
0.1550
10
8.3
0.0278
25
6.1
0.0200
11
4.5
0.0035
26
11.9
0.1055
12
9.4
0.0463
27
8.9
0.0400
13
6.6
0.0183
28
5.2
0.0113
14
3.9
0.0048
29
10.6
0.0703
15
9.9
0.0515
30
11.2
0.0925
103
猎物消耗的功能响应
等于零的系数
无
a
c
ɑ的计算公式
ɑ与N之间的关系
增加至渐近线;截距为d,不等于
(d+bN)∕(1+cN)
d  bN
a

代入
零N T
1  cN
d+bN
Ne 
线性增加;截距为d,不等于零
1  aN Th
d
bN∕(1+cN)
增加至渐近线;截距等于零
d和c
bN
线性增加;截距等于零
dN T  bN T
2
Ne 
1  cN  dN T b  bN T h
2
104