用R语言进行数据分析：定义统计模型的公式

下面的统计模型的模板是一个基于 独立的方差齐性数据的线性模型

     y_i = sum_{j=0}^p beta_j x_{ij} + e_i,     i = 1, ..., n,

其中 e_i 属于 NID(0, sigma^2)。 用矩阵格式表示，它可以写为

     y = X  beta + e

其中 y 是响应向量，X 是模型 矩阵（model matrix）或者设计矩阵（design matrix）。X 的列 x_0, x_1, …, x_p是 决定变量（determining variable）。通常，x_0 列就是截距（intercept）项。

例子

在给予正式的定义前，下面的例子可能会给出 一些不错的描述。

假定 y, x, x0, x1, x2, … 是数值 变量，X 是一个矩阵，而A, B, C, … 是因子。下面的例子中，左边给出公式， 右边给出该公式的统计模型的描述。

y ~ x

y ~ 1 + x

二者都反映了y 对 x 的简单线性模型。第一个包含了一个隐式的截距项， 而第二个则是一个显式的截距项。

y ~ 0 + x

y ~ -1 + x

y ~ x - 1

y 对 x 的通过原点的简单线性模型 (也就是没有截距项)。

log(y) ~ x1 + x2

y 的变换形式 log(y) 对 x1 和 x2 进行的多重回归 (有一个隐式的截距项)。

y ~ poly(x,2)

y ~ 1 + x + I(x^2)

y 对 x 的二阶多项式回归。第一种是 正交多项式（orthogonal polynomial），第二种是显 式地注明各项的幂次。

y ~ X + poly(x,2)

y 利用模型矩阵 X 和二阶多项式项x 进行多重回归。

y ~ A

y 的单一分类因子的方差分析模型。分类因子 由 A 决定。

y ~ A + x

y 的单一分类因子的方差分析模型。因子项为 A，共变项为 x。

y ~ A*B

y ~ A + B + A:B

y ~ B %in% A

y ~ A/B

基于 A 和 B的非可加两因子方差分析模型y。 前两个公式表示相同的交叉分类设计， 后两个公式表示相同的嵌套分类设计。概略地说， 四个公式指明同一个模型子空间。

y ~ (A + B + C)^2

y ~ A*B*C - A:B:C

三因子实验。该模型包括一个主要影响因子，另外两个因子 有相互作用。这两个公式等价。

y ~ A * x

y ~ A/x

y ~ A/(1 + x) - 1

A的每一水平拟合 y 对 x 的简单线性回归。 三个公式的编码不一样。最后一个公式会对 A 每一水平分布估计截距项 和斜率项的。

y ~ A*B + Error(C)

一个实验设计有两个处理因素 A 和 B以及 因子 C 决定的误差分层（error strata）。例如在一个 区组由因子 C 决定的裂区实验设计（split plot experiment）。

操作符 ~ 用于定义R 的模型公式。 一个普通线性模型的公式可以表示为

     response ~ op_1 term_1 op_2 term_2 op_3 term_3 ...

其中

response

是一个作为响应变量的向量或者矩阵， 或者是一个值为向量/矩阵的表达式。

op_i

是一个操作符。它要么是 + 要么是 -，分别表示在一个模型中加入 或者去除某一项(默认值为前者）。

term_i

可以

• 是一个向量，矩阵表达式或者1，
• 因子
• 是一个由因子，向量或矩阵通过公式操作符（formula operators） 连接产生的公式表达式（formula expression）。

基本上，公式中的项决定了模型矩阵中列要么被加入要么被去除。 1 表示截距项，并且 默认加入模型矩阵除非 显式地去除这一选项。

公式操作符在效果上和用于程序 Glim 和 Genstat 中的 Wilkinson & Rogers 标记符相似。一个不可避免的改变是 操作符 . 变成了 : 因为点号在 R 里面是合法的名字字符。

这些符号总结如下(参考 Chambers & Hastie, 1992, p.29):

Y ~ M

Y 由模型 M解释。

M_1 + M_2

同时包括 M_1 和 M_2项。

M_1 - M_2

包括 M_1 但排除 M_2项。

M_1 : M_2

M_1 和 M_2 的张量积（tensor product）。 如果两项都是 因子，那么将产生“子类”因子。

M_1 %in% M_2

和 M_1:M_2 类似，但编码方式不一样。

M_1 * M_2

M_1 + M_2 + M_1:M_2.

M_1 / M_2

M_1 + M_2 %in% M_1.

M^n

M 的所有各项以及所有到n阶为止的“交互作用”项

I(M)

隔离 M。M内所有操作符当一般的运算符处理。 该项出现在模型矩阵中。

注意，在常常用来封装函数参数的括弧中 的操作符按普通的四则运算法则解释。 I() 是一个恒等函数（identity function），它使得 常规的算术运算符可以用在模型公式中。

还要注意模型公式仅仅指定了模型矩阵的 列项，参数则没有显式地指定。 在某些情况下可能不是这样，如 非线性模型。

对照

我们至少要知道模型公式是如何指定模型矩阵的列项的。 对于连续变量这是比较简单的，因为每一个变量对应于 模型矩阵的一个列(如果模型中包含截距， 会在矩阵中列出这一列的)。

对于一个 k-水平的因子 A 该如何处理呢？无序和有序因子的 结论是不一样的。对于无序因子 k – 1 列用于作为因子第二，…， 第 k 水平的指标(因此隐含的参数就是 其他水平和第一个水平的响应程度的比较)。对于 有序因子，k – 1 列是 在 1, …, k 上的正交项，并且忽略常数项。

尽管这里的回答有点复杂，但这不是事情的全部。 首先在含有一个因子项的模型中忽略截距项， 这一项将会被编入指示所有因子水平的 k 列中。 其次整个行为可以通过 options 设置参数 contrasts 而改变。 R 的默认设置为

     options(contrasts = c("contr.treatment", "contr.poly"))

提这些内容的主要原因是 R 和 S 对无序因子 采用不同的默认值。S采用 Helmert 对照。因此，当你需要比较某本书上或者论文上用 S-Plus 写的代码，你必须设置

     options(contrasts = c("contr.helmert", "contr.poly"))

这是一个经过认真考虑的改变。因为处理对照（treatment contrast）(R默认) 对于新手是比较容易理解的。

这还没有结束，因为在各个模型的各个项中对照方式可以用 函数 contrasts 和 C 重新设置。 我们还没有考虑交互作用项：这将会产生 各分量项的复合项。

尽管细节是复杂的，R 里面的模型公式 可以产生统计专家所期望的各种模型。 例如，利用关联项而非主要效应的模型拟合 常常会产生令人惊讶的结果， 不过这些仅仅为统计专家们设计的。