设计矩阵的理解

发表于 | 分类于
| 字数统计: 1,925 | 阅读时长 ≈ 8

设计矩阵

在统计学中，一个设计矩阵(design matrix)也被称为模型矩阵(model matrix)或回归矩阵(regressor matrix)，它是一系列解释变量(explanatory variable)值的矩阵，通常命名为X。每行代表一个独立的观测对象，每列代表这个对象相应的变量。一些统计模型中常常会用到设计矩阵，例如广义线性模型(generalized linear models)。这种模型中通常会包含一些哑变量（例如0或1），这些哑变量用于区分不同的分组。设计矩阵的数据是一些解释变量的数据，解释变量用于解释应变量。

设计矩阵举例

在一个设计矩阵$X$中，$X_{ij}$表示矩阵$X$矩阵中第$j$列，第$i$行的一个变量值，在回归模型中，它可以写为$y=X\beta$这样的形式，其中$X$就是一个设计矩阵，$\beta$就是一个模型系数的向量，而$y$就是预测值的一个向量。我的理解就是，这个向量含有n值的组合，如果是简单线性回归，那么这个向量就是一个值，如果是多元线性回归，那么这个$\beta$就是多个值组成一个向量，在此，多元线性回归与简单线性回归在本质上是一样的。

一个$n \times p$的矩阵指的就是n个观测样本，p个变量（也可以理解为特征），如下所示：

1
2
3
4
matrix(1:12,nrow=2)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    1    3    5    7    9   11
# [2,]    2    4    6    8   10   12

在上面的这个矩阵里，n就是2，6就是p，这个矩阵就是$2\times6$矩阵。而在统计学中，当用设计矩阵（我觉得这样命名，主要是因为这个矩阵是与实验设计有关的）来表示一个实验时，这个n值就是在一个实验中的不同的重复。而列则表示一个样本中要检测的不同变量。例如当我们进行一个问卷调查时，问了10个人4个问题，那么最终在处理数据时，用设计矩阵来表示就是$10 \times 4$矩阵，这个10表示10个人，4表示4个问题，如果用i表示行（一共有10行），j表示列（一共有4列），那么这个矩阵中的某个值就可以看作是我问了第i个人的第j个问题。

再来看一个案例，这个案例是简单线性回归模型。在这个简单线性回归模型中，它的解释变量只有1个，也就是说，用设计矩阵来表示时，它只有1列，我们先假设它有7行，也就是i=7，用公式来表示这个简单线性模型就是下面的这个样子：

mark

其中$\beta_{0}$表示上面这个方程在y轴上的截矩，而$\beta_{1}$表示斜率，上面的这个方程用设计矩阵来表示就是下面的这个样子：

mark

再来看一个多元线性回归的案例，它有两个协变量，分别为w和x，假设我们有7个观测值，用方程来表示这个多元线性回归就是如下所示：

mark

用设计矩阵来表示就是下面的这个样子：

mark

这个矩阵就是一个$7\times 3$矩阵。

来看one-way ANOVA的案例。在这个案例中，我们有3个实验组，7个观测值。其中前3个观测值属于1组，第4和第5这两个观测值属于1组，第6和第7这两个观测值属于一组，那么这个模型可以写为下面的公式：

mark

用设计矩阵来表示就是下面的这个样子：

mark

其中$\mu_{i}$表示第$i$组的均值。

ANOVA模型通常可以写为每组的参数$\tau_{i}$与某些参考值的差值。通常这个参考值指提某个组的组。这样在多组进行比较时，例如多个实验组与一个对照组比较时就比较方便。在下面的案例中，我们选择group 1作为参考组，那么上面的这个公式就可以写为下面的这个样子：

mark

如果写为设计矩阵，就是下面的这个样子：

mark

在这个模型中，$\mu$是参考组的均值，而$\tau_{i}$则是group i与参考组的差值，上面的矩阵中不含有$\tau_{1}$是因为它与参考组的差值为0。

当我们使用R中的lm或aov或glm来进行线性描述时，就需要从formula和data参数中来创建一个模型矩阵。

看一个R中的简单线性回归的矩阵，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Formaldehyde
# carb optden
# 1  0.1  0.086
# 2  0.3  0.269
# 3  0.5  0.446
# 4  0.6  0.538
# 5  0.7  0.626
# 6  0.9  0.782
summary(fm1 <- lm(optden ~ carb, Formaldehyde))
# 
# Call:
#   lm(formula = optden ~ carb, data = Formaldehyde)
# 
# Residuals:
#   1        2        3        4        5        6 
# -0.00671  0.00103  0.00277  0.00714  0.00751 -0.01174 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.00509    0.00783    0.65     0.55    
# carb         0.87629    0.01353   64.74  3.4e-07 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.00865 on 4 degrees of freedom
# Multiple R-squared:  0.999,	Adjusted R-squared:  0.999 
# F-statistic: 4.19e+03 on 1 and 4 DF,  p-value: 3.41e-07
model.matrix(fm1)
# (Intercept) carb
# 1           1  0.1
# 2           1  0.3
# 3           1  0.5
# 4           1  0.6
# 5           1  0.7
# 6           1  0.9
# attr(,"assign")
# [1] 0 1

我们还可以从一个formula和data中直接创建一个模型矩阵，如下所示：

1
2
3
4
5
6
7
8
9
10
model.matrix(~ carb, Formaldehyde)
# (Intercept) carb
# 1           1  0.1
# 2           1  0.3
# 3           1  0.5
# 4           1  0.6
# 5           1  0.7
# 6           1  0.9
# attr(,"assign")
# [1] 0 1

在上面的这个公式中，~只有右边的变量，没有左边的变量，因为我们并不需要应变量来研究模型矩阵。双边公式(a two-sided formula)也可以写为单边公式，

再看一个单因素方差分析模型，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
test_aov <- aov(iris$Sepal.Length~iris$Species)
model.matrix(test_aov)
# (Intercept) iris$Speciesversicolor iris$Speciesvirginica
# 1             1                      0                     0
# 2             1                      0                     0
# ....
# 69            1                      1                     0
# 70            1                      1                     0
# 71            1                      1                     0
# 72 ...
# 110           1                      0                     1
# 111           1                      0                     1
# 112           1                      0                     1
# 113           1                      0                     1
# ...
# attr(,"assign")
# [1] 0 1 1
# attr(,"contrasts")
# attr(,"contrasts")$`iris$Species`
# [1] "contr.treatment"

线性回归的截矩

在R中，formula是一个符号表达式对象，它表示的一个或多个自变量（右边）与应变量（左边）的关系，在经典回归或ANOVA中，lm()函数能够闺怨一个formula对象转化为一个设计矩阵X和一个系数向量B，然后通过一定的统计学方法得到向量B的最佳估计。通常公式中的每项都会被转换为X上的一个或多个列。而model.matrix()函数则能将参数转换为设计矩阵X。看一个案例，如下所示：

1
Fuel ~ Weight + Displacement

其中Fuel是一个应变量，Weight和Displacement是自变量。这个公式的拟合过程就是找到下面统计模型的参数，如下所示：

$Fuel = \alpha + \beta_{1}Weight + \beta_{2}Displacement + \epsilon$

需要注意的是，这个公式并没有明确地包含α或ε，但是这个截矩在上面的公式中已经隐含了，如果我们添加了1，那么就明确表示截矩存在，如下所示：

1
2
Fuel ~ 1 + Weight + Displacement    
#  明确含有截矩这一项

如果我们要拟合的线性模型中不含截矩，那么此时要添加上0，或者是-1，如下所示：

1
2
3
4
Fuel ~ 0 + Weight + Displacement      
# 0表示公式中明确不含截矩
Fuel ~ -1 + Weight + Displacement    
# -1也有同样的效果

参考资料

  1. Design matrix
  2. Statistical Formulas