问题

如题,在进行回归分析时(本例还是个平衡分组的数据),R语言的输出结果与JMP不一致。

JMP的分析结果如下
Regression Analysis in JMP

R的分析结果如下:
Regression Analysis in R

很奇怪,Batch的主效应和JMP不一样,JMP中的Effect Tests也在R中不同,那么,如何得到相同的结果呢?

因为大多统计软件Minitab, JMP,SAS默认使用TypeIII进行离差平方和计算。而R默认是TypeI,所以本例建模的时候设置contrastscontr.sum

td<-read.csv("regression_interactive.csv")

td_fit<-lm(PotencyBatch*Time,data=td,contrasts=list(Batch="contr.sum"))
summary(td_fit)
library(car)
Anova(td_fit,type="III")

也可以在建模的时候,不使用contrasts参数,而是通过运行options(contrasts=c("contr.sum","contr.sum"))设置R的环境选项。

文件下载:

JMP 文件
CSV 数据文件

环境说明:

JMP:没什么好说的,就是商业软件
R:

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Simplified)People's Republic of China.936
[2] LC_CTYPE=Chinese (Simplified)
People's Republic of China.936

[3] LC_MONETARY=Chinese (Simplified)People's Republic of China.936
[4] LC_NUMERIC=C

[5] LC_TIME=Chinese (Simplified)
People's Republic of China.936

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] car_3.0-2 carData_3.0-2 RJSONIO_1.3-1.1 XML_3.98-1.16 RCurl_1.95-4.11 bitops_1.0-6 bindrcpp_0.2.2
[8] dplyr_0.7.7 stringr_1.3.1

loaded via a namespace (and not attached):
[1] zip_1.0.0 Rcpp_0.12.19 cellranger_1.1.0 pillar_1.3.0 compiler_3.5.1 bindr_0.1.1

[7] forcats_0.3.0 tools_3.5.1 tibble_1.4.2 pkgconfig_2.0.2 rlang_0.3.0.1 openxlsx_4.1.0

[13] rstudioapi_0.8 curl_3.2 haven_2.0.0 rio_0.5.16 hms_0.4.2 tidyselect_0.2.5
[19] glue_1.3.0 data.table_1.12.0 R6_2.3.0 readxl_1.2.0 foreign_0.8-70 purrr_0.2.5

[25] magrittr_1.5 abind_1.4-5 assertthat_0.2.0 stringi_1.2.4 crayon_1.3.4 rstudio_0.98.1103

因为 JMP 悄悄地给你把 time 减掉了 7.5 做了个中心化而 R 没这么做,你运行下下面代码就会得到 JMP 的结果了。

td$time <- td$Time-7.5
td_fit<-lm(Potency~Batch*time,data=td,contrasts=list(Batch="contr.sum"))
summary(td_fit)

你的这组数据对中心化是很敏感的,JMP很大胆地替你做了个默认设定,谨记商业软件也是有免责条款, 类似这种差异如果你自己不做说明出问题锅永远是你自己的,不论开源软件还是商业软件。

    yufree 谢谢大佬的指点,果然 ,中心化了。。。可是,中心化以后,截距变了……

    另外,大神,我想再问一下contrasts对方差分析表和系数表会有哪些影响啊,能大概说一下吗?付费也行。我整理了一个我的笔记,不过感觉还是不够专业,不够全。

    R的输出粘过来太乱了,所以截图了。

    -----------------------------------------我是分割线:(2小时以后)自问自答--------------------------------------------

    我明白了, slope=tan(A)=a/b, 其中a是截距变化,b是横坐标变化。

    所以,中心化以后,横坐标减掉了7.5,所以截距减少(a)=斜率(slope)横坐标变化(b)=-0.227.5=-1.65,所以真实截距=98.88333-(-1.65)=100.53333。

    好像JMP虽然在计算系数的时候中心化了,但是计算截距的时候,还原成了真实截距。

    方差分析里不同模型效应对比方法可以看这里这里的解释,我没把握说的比他们更清楚。

    R的结果输出问题可以用 broom 包来整理输出。

      yufree 我翻译了您给的2个参考资料之一,如下 :

      方差分析(ANOVA,以下称方差分析)是用来分析不同因子对样本的方差贡献水平的统计学过程。它源起于R。A.Fisher in 1925,关于平衡数据(每个因子水平有相同的观测数)的一个例子。

      当数据是非平衡的时候,有几种不同的方法来计算方差分析中的平方和(sums of suqare)。至少有三种方法,通常称为Type I, II, III 平方和(这个名称最初好像是从SAS包中被传播到统计界世界的,并且现在已被广泛传播和应用)。空间使用哪个类型的平方和被统计学家广泛讨论 。然而,对数据使用不同的假设,比如Type I, II III平方和,来进行检验,是非常有必要的。

      想象一下,假设有一个模型,含有两个因子A和B,因此,会有两个主效应和一个交互效应AB。整个模型应该是SS(A, B, AB)

      其它的模型可能类似:SS(A, B)表示 模型没有交互效应。SS(B, AB)表示模型的因子A的效应不明显。诸如此类。

      特定的因子(包括交互效应)的影响,可以通过对不同的模型进行检查来检验出来。比如,为了确定交互作用是否显著,将会对模型SS(A, B, AB)和SS(A, B)两个模型进行F检验。

      非常容易通过平方和的增加来表示这些模型的不同。让:
      SS(AB | A, B) = SS(A, B, AB) – SS(A, B)
      SS(A | B, AB) = SS(A, B, AB) – SS(B, AB)
      SS(B | A, AB) = SS(A, B, AB) – SS(A, AB)
      SS(A | B) = SS(A, B) – SS(B)
      SS(B | A) = SS(A, B) – SS(A)
      这些表达式表明了平方和增加的不同。比如SS(AB | A, B)表示“交互作用在主效应之后”,以及SS(A | B)表示“A主效应的平方和在B的主效应平方和之后,并且忽略交互作用”。

      不同的平方和类型基于在执行的时候,模型减少的阶段(可以理解为逐步法的减少)。尤其是:
      Type I,通常被称为“序贯”平方和。
      SS(A) for factor A.
      SS(B | A) for factor B.
      SS(AB | B, A) for interaction AB.

      这个用来检验因子A的主效应,在A的主效应之后是B的主效应,然后在所有主效应后是AB的交互作用。

      因为序贯的天然属性,以及两个主效应的检验是按一定顺序进行的这个事实,对于非平衡数据,这个类型的平方和会根据主效应的顺序给出不同的结果。

      对于非平衡数据,这个方法用来检验边际均值效应。对于特定的项,这意味着结果会依赖于实际的样本量,也就是说,特定样本的比例。换句话说,它检验第一个因子时并不考虑控制其它的因子。(对于以后的讨论和实例,可以看Zahn )

      注意,在处理非平衡数据时,这经常不是一个有趣的假设。
      Type II:
      SS(A | B) for factor A.
      SS(B | A) for factor B.

      这个类型在其它主效应后检验每个主效应

      注意
      注意没有明确假设交互项(换句话说,你应该先检验交互作用(SS(AB | A, B)),并且仅当交互作用AB不显著的时候,继续分析主效应。)

      如果实际上没有交互项,Type II在统计上比Type III有更高的功效。(详细情况参考Langsrud )

      计算上,这等于使用不同的因子顺序,执行Type I分析,并且得到合适的结果(然后,在其它效应之后运行一个主效应,在上面的例子中)
      Type III:
      SS(A | B, AB) for factor A.
      SS(B | A, AB) for factor B.

      这个类型在其它主效应和交互作用之后检验一个主效应。这个方法因此可以在存在交互项的时候使用。

      然而,如果存在交互作用,解释一个主效应的意义通常也不是非常有意思。(总的来说,如果存在交互作用,主效应应该在后面进行分析)

      如果交互项不显著,Type II的功效更高。
      注意:当数据是平衡的时候,因子是正交的,Type I, II, III给出的结果是一样的。总结来讲,是控制其它因子水平的情况下,检验一个因子的显著性。这等于使用Type II或III 平方和。大体上,如果没有交互作用,Type II功效更高,遵守边际原则。如果存在交互作用,Type II就不太合适了,不过仍然可以使用Type III。但是结果解释应该小心(存在交互作用的时候,主效应很难解释)。

      方差分析和R中的aov函数
      方差分析和R中的aov函数执行序贯平方和(Type I)。就像上面提到的,这个很难检验假设的效应,因为,计算一个因子的效应时,其它因子的水平是变化的。这意味着结果仅可以解释在(非平衡)数据样本中出现的特定水平的观察值。幸运的是,基于以上讨论,这可以很清楚的直接与Type II的结果做对比。

      因为Type II平方和在其它所有主效应之后检验一个主效应,并且假设没有交互作用。可以使用anova()函数得到正确的平方和,并且可以变化因子顺序。

      比如,假设一个数据框。(search) 响应变量是用户使用信息检索系统找到相关答案的时间(time),用户被分成2组使用不同的系统进行搜索(sys)。他们也被要求使用一些不同的索引条件(topic).

      Type I SS in R:
      anova(lm(time ~ sys * topic, data=search))
      如果数据是非平衡的,你使用下面的命令会得到稍微不同的结果。
      anova(lm(time ~ topic * sys, data=search))
      使用上面每一条命令,Type II平方和通过每个结果的第二行得到。(因为Type I平方和,第二个元素将会在第一个因子后成为第二个因子)。也就是,你通过第一条命令,得到了topic的Type II平方和结果。通过第二条命令,得到了sys的Type II平方和结果。

      Type III SS in R:
      这与Type II的结果比稍微更复杂一点。首先,应该在R中设定contrasts选项,因为多因子方差分析模型是过度参数化的,有必要选择一个对比对象,设定平方和是0,否则方差分析将会在假设检验上给出错误结果。(默认的对比类型不满足这个要求)
      options(contrasts = c(“contr.sum”,”contr.poly”))
      然后建模:
      model <- lm(time ~ topic * sys, data=search)
      最后,对每个模型元素调用drop1函数.
      drop1(model, .~., test=”F”)
      这个结果给出的就是Type III的平方和,包括了F检验的P值。
      使用car包计算Type II和Type III平方和。
      使用其它更方便的方法得到Type II, III平方和结果是使用car包,它定义了一个新的函数,Anova(),可以直接计算type II和Type III平方和。
      Type II,使用上面相同的数据。
      Anova(lm(time ~ topic * sys, data=search, type=2))
      Type III:
      Anova(lm(time ~ topic * sys, data=search, contrasts=list(topic=contr.sum, sys=contr.sum)), type=3))
      注意:再次强调,因为SS是根据交互作用计算的,对于Type III你必须指定contrasts 选项来得到更合理的结果(这里是解释)
      背景: 平方和的由来
      源头
      如果有一个类型I平方和和类型II平方和, 就得到了一个问题。有更多的吗?很幸运,实际上,有4个平方和类型,每个都有他们相关的统计方法。这4个类型,称为Type I, II, III, IV(Goodnight 1978)。虽然我们本课程只使用类型I和类型II,这里是4个类型的介绍。
      Type I (sequential or incremental SS/序贯平方和)
      Type I计算平方和时,按顺序考虑每个因子,按他们在模型中列出的顺序。类型I平方和对非平衡数据,多因子结构的计算可能没有太大用处,但是可能对于平衡数据的嵌套模型会非常有用。类型I平方和对于简单多项式模型同样比较有用(比如回归分析)。允许在使用更复杂模型(二次,三次等)前,使用简单因素(比如线性)来解释同样变量数的变异。并且,通过将类型I与其它类型比较,提供了更多的非平衡数据差异的信息。类型II和类型III平方和被大家熟知为部分平方和,每个效应是根据其它效应进行调整的。

      Type III
      类型III是部分平方和方法,但是它比类型II更容易解释,所以我们从这里开始。在这种模型里,每个效应是根据所有其它效应来调整的。对于某个数据集,如果缺失值被最小二乘解的值代替,类型III平方和将会产生和类型I一样的平方和。类型III平方和与耶茨加权平方和分析方法一致。如果使用这个平方和方法,要求在有交互项的时候对比主效应。(类型II并不做这些,顺便,很多统计学家说不根本不应该做)。特别是,主效应A和B根据交互作用A*B调整,并且所有这些项都在模型中。如果模型只包括主效应,你会发现类型II和类型III分析结果一样。
      Type II
      类型II部分平方和有些难理解。大体上,类型II平方和对于可能是主效应或交互作用的效应U,根据效应V来调整,并且只有当V不包括U的时候。特别是,对于有交互作用的二因子结构,主效应A和B不是根据AB交互项调整的,因为这个交互项包括了A和B。因子A根据B调整,因为B不包括A,类似的,B根据A调整。最终,AB交互作用根据两个主效应的每一个来调整,因为没有一个主效应同时包括A和B。换句话说,类型II平方和根据所有其它不包括完整效应项目的因子进行调整。在一些情况下,你可以把它想象成一个序贯部分平方和,在里面,你允许低阶项解释尽可能多的变量,在让高阶项进入考虑之前,根据另一个进行调整。
      Type IV
      类型IV函数最初是在考虑到一些空格的情况下被设计出来的,同样被熟知为"radical"数据丢失。类型IV平方和的原则是非常复杂的,并且仅可以在一般结构方程的框架下进行讨论。应该注意的是类型IV函数不是当有空格的时候唯一的选择,但是等同于没有空格的类型III的结果。

        Jonie_Y 问题问得这么整齐细致,还能整理答案;这位客官若不嫌弃,不妨来加入论坛版主队伍吧。

          Jonie_Y 感谢整理!有一点需要澄清下:这里 powerful 是统计学意义上的,最好不要理解为字面上的更有效,应该是显著性检验的功效更高,也就是发现真实差异的能力更强。一般来说功效强的推断模型的模型假设更多(例如符合xx分布之类),也因此鲁棒性往往不太理想。但如果模型假设基本符合事实,功效越强的假设检验也可以理解为更有效。

            yihui 哇塞,我竟然得到了谢大神肯定和邀请,让我这个30多岁的老男人也有点激动啊。不过很少有时间,就清明节这两天好不容易自己独处,有时间看看书,整理整理材料,问问问题。虽然我很乐意,不过怕是版主的话会没时间弄。不过还是谢谢大神了。

            至于问问题的方式呢呢,其实也是一个过程吧。其中两个事情对我有深远影响。

            1,2011年开始学习统计(应该是重拾),那会儿已经忘光了,在QQ群里问了一句别人,“什么是假设检验”。那个人回答“去看书”。我觉得挺有道理,后来看了书果然理解了。所以,我明白其实好多问题,都是看书能解决的,而且更加系统。
            2,我是2013年左右开始学习R语言,因为第一件事,很少问问题,往往都是自己学。一般都是自己查也查不到,自己问研究了好几天的时候去问,其实也是问了也没人回答。然后就是自己继续研究。所以总得来说,问问题不多。
            3,另外一件对我有深远影响的就是大神在另一个网站的“论程序员的自我修养”,我深以为然。
            4,这次也是研究了好多天,搜索引擎各种关键字,往后翻n页,还是没找到答案(其实有一篇帖子提到了中心化,可是我用他的代码和数据没能重现他贴的JMP的结果。如果当时自己用JMP也重做一遍就好了)
            5,这次@yufree给了这么明确的答案,让我豁然开朗。我理应分享一下学习结果。而且就算是google,也只搜到了一篇stackoverflow里面的相关关于具体问题的帖子。我们身为统计之都的国内论坛,我想贡献一下力量(国内资料太少了)。
            6,之前发帖前也是认真看了公告才发纠正过,算是学习过程吧。就是我可以提供更多信息不一定有用,但是不能缺少需要的信息。(大神有一句话这个意思)
            7,看过大神的博客,感觉大神在某些方面还是有点完美主义的,我也怕被批评,哈哈哈。

            总之,谢谢大神。

            ------------------------我是分割线(又是2小时以后)----------------------
            刚才其实我只是随便点了一下不同的标签,然后点到版主版块,发现说明写着普通用户不可见,但是我能看见啊,第一反应就是屁颠屁颠去报告BUG,贴子都写好了,突然反应过来,去看了一下级别,发现变版主了~~

            如果后面老大不移除权限我就去学习一下版主需知,如果时间不能保证没关系的话。

            再次感谢谢大神~~ 这次算是近距离接触了~

            yufree 嗯嗯,谢谢大神,确实powerful翻译有点问题,后面我更新一下。

            ---------我是分割线(2小时以后)-----------
            已更新