阅读提示
要更好的理解这篇文章，您可能需要对以下知识有一定的了解。
- 方差分析
- 矩阵的乘法
- 线性模型
- 回归分析
- 对比
- 正交对比
前言
如果阅读过通过线性模型详解方差分析中“对比”（Contrasts）的数理原理,同学们应该知道，通过设置处理效应的对比系数，可以对各处理效应进行比较并检验。例如令：
但是，在对比中，存在一个较为特殊的对比——正交多项式对比。其系数矩阵非常奇怪，虽然也是让各处理效应加加减减，但并不是比较处理效应间的差异，而是对处理效应的“趋势”进行检验。这里的“趋势”可以简单理解为：一次函数关系，二次函数关系，三次函数关系……
好奇的小伙伴可能会有如下问题：
- 正交多项式对比与其他的对比有什么区别吗？
- 为什么对处理效应加加减减就能对各种“趋势”进行分析？
- 什么叫做正交多项式？为什么一定要正交的多项式？
这篇文章将带领大家一起探讨上述问题
1. 对多项式进行换汤不换药的重组
多项式是非常基础的概念，一个一元N次多项式通常指如下形式(为了与后面的内容一致，此处自变量用
很显然，
令：
则
此外，还可以进行如下形式的重组：
同样
考虑到系数，还可以进行如下形式的重组：
则
可以看出，对于任意一个多项式
可以重组成
其中
如果将
将1.1式写为：
也就是将一元n次回归替换成了n元一次回归
从这个角度来看，将1.1式变为1.2式，也就是进行了换元
此时如果对
1.2式的回归系数是
可以得知：
- 可以对任一多项式进行重新组合，看作是**带系数的一族多项式（从
到
）**之和
- 重新组合后相当于对自变量进行了换元，可采用多元线性回归进行分析
但是，组合的方式非常多，到底怎么样组合比较好呢？其实，比较好的形式就是：正交多项式
2. 正交多项式
正交多项式指一族多项式，其中任意两个多项式在某个区间内的内积积分为0，所以叫做正交。例如一族多项式
其中
且
则上述就是一个正交多项式（族）
- 准确来说应该是
其中
称为权函数，即权重，此处取
正交多项式在方差分析中的作用
如果我们将1.2式的
3. 方差分析模型和多项式回归模型
在这一部分内容中，会很详细的介绍方差分析模型与多项式回归模型的相互推导。只需要满足：
- 自变量的各水平是等距数据（至少为顺序，即可以转换成体现大小的数字），而不仅仅是分类数据 例如，在一项心理旋转的研究中，自变量是旋转角度，包含四个水平（0°，90°，180°，270°），因变量是反应时。 这是一个单因素四水平的研究。 分别用
来表示四个水平，假设每个水平测量了5次，用
来表示测量的次数，结果如下表：
上式中
例如
对上式进行回归分析，相当于四元一次回归，可以估计四个回归系数，对应的就是各处理效应量的估计。 考虑各自变量
可以发现这是一个4×4的单位矩阵。 既然是单位矩阵，那各列向量肯定两两正交，自变量相互独立（向量点积为0） 实际上，上一句话体现了方差分析的前提条件之一：各处理观测值独立。正因为正交，所以没有相关，所以相互独立。
考虑到自变量的各水平实际上并不是分类数据，而是等间距的测量数据（0°，90°，180°，270°）。所以，也可以将自变量
其中
为了方便，我们将
其中
将其带入（2）式，则（2）式变为：
令
则
注意到
而且
代入(3)，得到
代回（1）
可以看到，此时相当于：
也就是说，可以通过调整（1）式中各
答案是：可以！
我们现在假设
其中
则
代回
上式中，令
此时
代回(1)式，得到
即
也就是说，按照上式进行回归分析的话，同样可以分析二次项，一次项，截距的回归系数，而且自变量中没有任何二次项
同理，按照上面的逻辑，我们可以将方差分析模型转换为任意次数的多项式模型，而且模型中只是对各虚拟变量
4 方差分析的正交多项式对比
有了上面的知识，我们可以进入正题了。 还是回到心理旋转的例子，我们现在假设反应时
其中
可以发现该式与方差分析的模型非常相似：
（这也是为什么我们假设自变量与因变量之间为三阶，而不是二阶，不是四阶。因为保证了待估参数的数量相同，而只有这样才能保证后面的正交多项式的构建）
区别在于方差分析模型各自变量是相互正交的，即相互独立的，所以可以通过平方和的分解来分析各处理效应（回归系数） 那么，如果能够构建出相互正交的
我们可以这样来构建一个正交多项式族，首先，令：
可以发现，因为
由于
将
则
用
其中
同理，我们可以把
实际上，这种有限取值的等间距变量，其对应的各阶正交多项式已经有现成的递推公式了。
其中n为自变量的水平数，即取值个数。z是标准化以后的自变量。
根据递推公式，可得到：
虽然还可以推出后面的多项式，但常用的也就这几项。
再次回到我们的例子，此时已经将原始的模型（3.1式）进行了正交多项式的重组(取
此时，关键的一步。用虚拟变量
从表得知，可进行如下替换：
将其代回，得到
将各虚拟变量
上述矩阵就是对比矩阵，下面我们通过一个数据例子进行验证
数据例子如下：
进行自定义的多项式对比
from statsmodels.formula.api import ols #拟合线性模型的包
from statsmodels.stats.anova import anova_lm#进行方差分析的包
import pandas as pd
data={
'A':[1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4],
'Y':[1,3,2,2,4,2,4,5,2,2,3,4,7,5,7,4,6,9,8,10]
}
df=pd.DataFrame(data)
poly_mat=[[1,-1.5,1,-0.3],[1,-0.5,-1,0.9],[1,0.5,-1,-0.9],[1,1.5,1,0.3]]
model_poly=ols('Y~C(A,poly_mat)+0',data=df).fit()
print(model_poly.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Y R-squared: 0.613
Model: OLS Adj. R-squared: 0.540
Method: Least Squares F-statistic: 8.434
Date: Thu, 02 Apr 2020 Prob (F-statistic): 0.00137
Time: 23:54:14 Log-Likelihood: -37.380
No. Observations: 20 AIC: 82.76
Df Residuals: 16 BIC: 86.74
Df Model: 3
Covariance Type: nonrobust
===========================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------------
C(A, poly_mat)[custom0] 4.5000 0.392 11.476 0.000 3.669 5.331
C(A, poly_mat)[custom1] 1.7200 0.351 4.904 0.000 0.977 2.463
C(A, poly_mat)[custom2] 0.4000 0.392 1.020 0.323 -0.431 1.231
C(A, poly_mat)[custom3] -0.2667 0.585 -0.456 0.654 -1.506 0.972
==============================================================================
Omnibus: 0.569 Durbin-Watson: 1.915
Prob(Omnibus): 0.752 Jarque-Bera (JB): 0.627
Skew: -0.174 Prob(JB): 0.731
Kurtosis: 2.205 Cond. No. 1.67
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.再进行内置的多项式对比，比较一下结果
model_poly=ols('Y~C(A,Poly)',data=df).fit()
#A是类别，所以要用操作符C()来定义，C表示Categories
print(model_poly.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Y R-squared: 0.613
Model: OLS Adj. R-squared: 0.540
Method: Least Squares F-statistic: 8.434
Date: Thu, 02 Apr 2020 Prob (F-statistic): 0.00137
Time: 23:54:23 Log-Likelihood: -37.380
No. Observations: 20 AIC: 82.76
Df Residuals: 16 BIC: 86.74
Df Model: 3
Covariance Type: nonrobust
========================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------------
Intercept 4.5000 0.392 11.476 0.000 3.669 5.331
C(A, Poly).Linear 3.8460 0.784 4.904 0.000 2.184 5.509
C(A, Poly).Quadratic 0.8000 0.784 1.020 0.323 -0.862 2.462
C(A, Poly).Cubic -0.3578 0.784 -0.456 0.654 -2.020 1.305
==============================================================================
Omnibus: 0.569 Durbin-Watson: 1.915
Prob(Omnibus): 0.752 Jarque-Bera (JB): 0.627
Skew: -0.174 Prob(JB): 0.731
Kurtosis: 2.205 Cond. No. 2.00
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.可以看到，除了截距之外，各回归系数的估计值均不相等，难道上面推导的全错了？？
实际上，这是因为各多项式的回归系数为
虽然系数不相同，但t检验结果完全一样，这也说明了
我们可以通过代码查看一下内置的多项式对比矩阵
#查看内置多项式对比矩阵
from patsy.contrasts import Poly
contrast = Poly().code_with_intercept(levels=[1,2,3,4])
print(contrast.matrix)
[[ 1. -0.67082039 0.5 -0.2236068 ]
[ 1. -0.2236068 -0.5 0.67082039]
[ 1. 0.2236068 -0.5 -0.67082039]
[ 1. 0.67082039 0.5 0.2236068 ]]
#计算其逆矩阵
import numpy as np
C_Poly=contrast.matrix
C_Poly_inverse=np.mat(C_Poly).I
print(C_Poly_inverse)
[[ 0.25 0.25 0.25 0.25 ]
[-0.67082039 -0.2236068 0.2236068 0.67082039]
[ 0.5 -0.5 -0.5 0.5 ]
[-0.2236068 0.67082039 -0.67082039 0.2236068 ]]内置的对比矩阵为：
其逆矩阵为：
有时候为了方便，我们也可以将各元素化为整数，例如：
同样也能得到相同的结果
小结
- 正交多项式对比也是对比的一种形式，前提是自变量水平是等距的
- k个处理的方差分析可以分析的最高次项为k-1次
- 其逻辑为：①假设自变量与因变量之间是k-1次多项式回归；②将该多项式重组成正交多项式之和；③利用虚拟变量替换各多项式，形成新变量；④对新变量的偏回归系数进行检验，体现了该新变量（多项式）的趋势是否显著。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
