R中的统计模型.docx

1、R中的统计模型R中的统计模型这一部分假定读者已经对统计方法，特别是回归分析和方差分析有一定的了解。后面我们还会假定读者对广义线性模型和非线性模型也有所了解。R已经很好地定义了统计模型拟合中的一些前提条件，因此我们能构建出一些通用的方法以用于各种问题。R提供了一系列紧密联系的统计模型拟合的工具，使得拟合工作变得简单。正如我们在绪论中提到的一样，基本的屏幕输出是简洁的，因此用户需要调用一些辅助函数来提取细节的结果信息。1定义统计模型的公式下面统计模型的模板是一个基于独立的方差齐性数据的线性模型 用矩阵术语表示，它可以写成其中y是响应向量，X是模型矩阵（model matrix）或者设计矩阵（des

2、ign ma-trix）。X的列是决定变量（determining variable）。通常，列都是1，用来定义截距（intercept）项。例子在给予正式的定义前，举一些的例子可能更容易了解全貌。假定y,x,x0,x1,x2,.是数值变量，X是一个矩阵，而A,B,C,.是因子。下面的例子中，左边给出公式，右边给出该公式的统计模型的描述。yxy1+x二者都反映了y对x的简单线性模型。第一个公式包含了一个隐式的截距项，而第二个则是一个显式的截距项。y0+xy-1+xyx-1 y对x过原点的简单线性模型(也就是说，没有截距项)。log(y)x1+x2 y的变换形式log(y)对x1和x2进行的多重

3、回归(有一个隐式的截距项)。ypoly(x,2)y1+x+I(x2) y对x的二次多项式回归。第一种是正交多项式（orthogonal polynomial），第二种则显式地注明各项的幂次。yX+poly(x,2) y利用模型矩阵X和二次多项式项x进行多重回归。yA y的单因素方差分析模型，类别由A决定。yA+x y的单因素协方差分析模型，类别由A决定，协方差项为x。yA*ByA+B+A:ByB%in%AyA/B y对A和B的非可加两因子方差分析模型（two factor non-additive model）。前两个公式表示相同的交叉分类设计（crossed classification），

4、后两个公式表示相同的嵌套分类设计(nested classification)。抽象一点说，这四个公式指明同一个模型子空间。y(A+B+C)2yA*B*C-A:B:C 三因子实验。该模型包括一个主效应（main effects）和两个因子的交互效应（interactions）。这两个公式等价。yA*xyA/xyA/(1+x)-1 在A的各个水平独立拟合y对x的简单线性回归。三个公式的编码不一样。最后一个公式会对A各个水平分别估计截距项和斜率项的。yA*B+Error(C) 一个实验设计有两个处理因素A和B以及因子C决定的误差分层（error strata）。如在裂区实验设计（split plo

5、texperiment）中，所有区组（还包括子区组）都由因子C决定的。操作符用来定义R的模型公式（model formula）。一个普通的线性模型式可以表示为responseop 1 term 1 op 2 term 2 op 3 term 3.其中response是一个作为响应变量的向量或者矩阵，或者是一个值为向量/矩阵的表达式。op i是一个操作符。它要么是+要么是-，分别表示在一个模型中加入或者去掉某一项(公式第一项的操作符可选）。term i可以(1)是一个向量，矩阵表达式或者1，(2)因子，(3)是一个由因子，向量或矩阵通过公式操作符连接产生的公式表达式（formula expres

6、sion）。基本上，公式中的项决定了模型矩阵中的列要么被加入要么被去除。1表示截距项，并且默认就已加入模型矩阵，除非显式地去除这一选项。公式操作符（formula operators）在效果上和用于程序Glim和Genstat中的Wilkinson&Rogers标记符(notation)相似。一个不可避免的改变是操作符.在R里面变成了:,因为点号在R里面是合法的命名字符。这些符号总结如下(参考Chambers&Hastie,1992,p.29):YM 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的张量积（te

7、nsor product）。如果两项都是因子，那么将产生“子类”因子(subclasses factor)。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%M1.Mn M的所有各项以及所有到n阶为止的“交互作用”项I(M) 隔离M。M内所有操作符当一般的运算符处理。并且该项出现在模型矩阵中。注意，在常常用来封装函数参数的括弧中的操作符按普通的四则运算法则解释。I()是一个恒等函数（identity function），它使得常规的算术运算符可以用在模型公式中。还要特别注意模型公式仅仅指定了模型矩

8、阵的列项，暗含了对参数项的指定。在某些情况下可能不是这样，如非线性模型的参数指定。1对照我们至少要知道模型公式是如何指定模型矩阵的列项的。对于连续变量这是比较简单的，因为每一个变量对应于模型矩阵的一个列(如果模型中包含截距，会在矩阵中列出值都是1的一列)。对于一个k-水平的因子A该如何处理呢？无序和有序因子给出的结论是不一样的。对于无序因子，因子第2，.，第k不同水平的指标产生k?1列。(因此隐含的参数设置就是把其他水平和第一个水平的响应程度进行比较)。对于有序因子，k-1列是在1,.,k上的正交项（orthogonal polynomial），并且忽略常数项。尽管这里的回答有点复杂，但这不是

9、事情的全部。首先在含有一个因子项的模型中忽略截距项，这一项将会被编入指示所有因子水平的k列中。其次整个行为可以通过options设置参数contrasts而改变。R的默认设置为options(contrasts=c(contr.treatment,contr.poly)提这些内容的主要原因是R和S对无序因子采用不同的默认值。S采用Helmert对照。因此，当你需要比较你的结果和某本书上或论文上用SPLUS代码的结果时，你必须设置options(contrasts=c(contr.helmert,contr.poly)这是一个经过认真考虑的改变。因为处理对照（treatment contrast

10、）(R默认)对于新手是比较容易理解的。这还没有结束，因为在各个模型的各个项中对照方式可以用函数contrasts和C重新设置。我们还没有考虑交互作用项：这些交互作用项将会产生各分量项的乘积。尽管细节是复杂的，R里面的模型公式在要求不是太离谱的情况下可以产生统计专家所期望的各种模型。提供模型公式的各种扩展特性是让R更灵活。例如，利用关联项而非主要效应的模型拟合常常会产生令人惊讶的结果，不过这些仅仅为统计专家们设计的。2线性模型对于常规的多重模型（multiple model）拟合，最基本的函数是lm()。下面是调用它的方式的一种改进版：fitted.modelfm2fmanova(fitted.

11、model.1,fitted.model.2,.)结果将是一个方差分析表以显示依次加入的拟合模型的差异。需要比较的拟合模型常常是等级序列（hierarchical sequence）。这个和默认的设置实际上没有差别，只是使它更容易理解和控制。5更新拟合模型函数update()是一个非常便利的函数。它允许拟合一个比原先模型增加或减少一个项的模型。它的形式是new.modelfm05fm6smf6fmfullfitted.modelfmfmkalythoskalythos$Ymatfmpfmlsummary(fmp)summary(fml)两种模型都拟合的很好。为了计算LD50，我们可以利用一个简单ld50ldp-ld50(coef(fmp);ldlfmodnlfitxyfn-function(p)sum(y-(p1*x)/(p2+x)2)为了进行拟