免费试用
banner

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

作者: afenxi来源: afenxi时间:2017-03-28 11:53:080

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

y_i = sum_ + 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)。 如果两项都是 因子,那么将产生“子类”因子1。 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 里面的模型公式 可以产生统计专家所期望的各种模型1。 例如,利用关联项而非主要效应的模型拟合 常常会产生令人惊讶的结果, 不过这些仅仅为统计专家们设计的。

原创文章,作者:古思特,如若转载,请注明出处:《用R语言进行数据分析:定义统计模型的公式》https://www.afenxi.com/post/3411

banner
看过还想看
可能还想看
热点推荐
Yonghong的价值观:以卓越的数据技术为客户创造价值,实现客户成功。
免费试用