统计模型——施肥效果分析（一）

背景

某地区作物生长所需的营养主要是氮（N）、磷（P）、钾（K）。某作物研究所在该地区对土豆与生菜做了一定数量的实验，实验数据如表4-1记录了土豆的实验数据，表4-2记录了生菜的实验数据。其中hm2表示公顷，t表示吨，kg表示千克。当一个营养的施肥发生变化时，总将另外两个营养施肥量保持在第七个水平。如对土豆产量关于N的施肥做实验时，P和K的施肥量分别取196kg/hm2与372kg/hm2.

目标：试分析施肥量与产量之间的关系，并对所有结果从应用价值与如何改进方案方面做出估计。

表4-1土豆的N、P、K效果：

表4-2生菜的N、P、K效果：

>> tn=[0 34 67 101 135 202 259 336 404 471];
yn=[15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75];
tp=[0 24 49 73 98 147 196 245 294 342];
>> yp=[33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73];
tk=[0 47 93 140 186 279 372 465 558 651];
yk=[18.98 27.35 34.86 39.52 38.44 37.73 38.43 43.87 42.77 46.32];
>> sn=[0 28 56 84 112 168 224 280 336 392];
xn=[11.02 12.70 14.56 16.27 17.75 22.59 21.63 19.34 16.12 14.11];
sp=[0 49 98 147 196 294 391 489 587 685];
>> xp=[6.39 9.48 12.46 14.33 17.10 21.49 22.46 21.34 22.07 24.53];
sk=[0 47 93 140 186 279 372 465 558 651];
xk=[15.75 16.76 16.89 16.24 17.56 19.20 17.97 15.84 20.11 19.40];

1、土豆产量与氮肥的关系

1.1 绘制散点图

以土豆产量与氮肥施肥量为例，初步观察，用二次函数作为经验方程。

1.2 回归模型

【3.1】

1.3 模型参数求解与显著性检验

调用matlab的回归求解函数regress,下面逐步解释。

>> clear     %清除内存里的一切变量;
>> load d:\tudou  %调回存储的数据;
>>[b,bt,r,rt,st]=regress(yn‘,[ones(length(tn),1),tn’,(tn.^2)‘])  %调用格式

• b 回归系数，升幂排列
• bt 回归系数的置信区间，默认95%
• r 残差，即理论值与测量值的差
• rt 残差的95%的置信区间
• st 模型检验参数（R2,F,P,sig2).（R2表示方程因变量与自变量之间的相关性检验，R也称为可决系数；F，p是方程显著性检验；sig2是σ2的估计值。）

regress   matlaba回归命令

yn’ 因变量（列）

[ones(length(tn),1),tn’,(tn.^2)‘])  矩阵，形式是

回归系数

b =    

14.7416  a0    

0.1971  a1  

-0.0003  a2

回归系数的95%的置信区间

bt =  

12.6301   16.8532  

0.1736    0.2207

-0.0004   -0.0003

回归系数的置信区间不应包括0；若包括0，说明该系数不显著（100次实验，95次以上离0非常近，非常危险（或者贡献不大））

方程显著性指标F=251.7971>Fα(2,7)说明回归结果显著。

残差的置信区间以包含0为好，说明残差总在0附近，意味着残差较小。如果出现残差置信区间不包括零，这个样本观测值极可能是极端值。

残差平方和记为SSE,利用SSE/(n-p-1)可以作为σ2的无偏估计值。

在该问题里，R2=0.9863表示控制氮肥施肥量，在100次种植实验里，至少有98次土豆的产量可以用该经验公式解释

1.4 模型修正或改进

根据上面的计算，发现第10个样本的残差置信区间右端点小于0，即第10个样本可能不正常。去掉第10个样本后（可以再线性插值替换），回归得到结果为R2=0.9956;F=678.5246;p<0.000

回归方程极显著：

1.5 预测预报

根据上面的分析当土豆的施肥量tn给定时，土豆的产量Yn服从分布

其中σ2根据残差平方和作点估计

 那么，在土豆的种植中，给定每公顷施肥量为tn=380kg时，土豆产量Yn的分布为

即土豆产量的95%的置信区间为( 39.8725  , 42.9181)

2、土豆产量与磷肥的关系（1）

2.1 模型猜测

磷肥施肥量过了100以后，土豆的产量几乎没有怎么增加，即土豆与磷肥的关系里，有一个上界，下面几个经验函数可以试试：

1. （指数函数）
2. （双曲函数）
3. (Logistic曲线）
4. （对数函数）
5. y=pm (x)  (多项式函数)
6. （幂函数）

 不管用哪个，都是非线性回归，误差和可信度分析都不如线性回归那么成熟好用。

2.2 数据预处理

（1）阅读数据和散点图发现，磷肥施肥量为24kg/hm2时的土豆产量反而不如不施肥的产量，认为第二个数据被污染了，用第一和第三个数据线性插值替代。

上述几个非线性模型，不管用哪个，都会涉及到自变量取负指数函数，由于自变量太大，取负指数后在0附近，太敏感。所以先对磷肥施肥量标准化：

2.3 经验公式线性化

如上图所示，土豆的最大产量不超过44，而威布尔函数 的最大值就是A，不妨假定A=44，则

取Y=ln(44-y),a0=lnB,a1=-C,则韦布尔函数就线性划为 绘制(tp’,Y)散点图验证猜测。

如下图所示，猜测正确。