以下内容主要译自《The R Book》，部分来自Stata软件的 help 文档。

仅用于学习、交流。疏漏和错误在所难免，欢迎指正。

 

问题：有sand, clay, loam三类土壤，从每个类型的庄稼地中随机抽取十块样地，计算样地的平均单位面积产量。如果每块样地栽培的作物类型、管理方式、病虫害防治投入等其他因素都相同，想知道土壤类型对单产的影响是否显著？如果显著，达到什么程度？再者，三类土壤两两比较，结果又如何？（书上没问这个）

## 土壤是分类变量，分sand，clay，loam三类。或者说，土壤这个因子有三个水平

 

Crop yields per unit area were measured from 10 randomly selected fields on each of three soil types.

sand	clay	loam
6	17	13
10	15	16
8	3	9
6	11	12
14	14	15
17	12	16
9	12	17
11	8	13
7	10	18
11	13	14

All fields were sown with the same variety of seed and provided with the same fertilizer and pest control inputs.

The question is whether soil type significantly affects crop yield, and if so, to what extent.

 

（一） 用 R 处理之

步骤：（红色=命令，蓝色=结果，绿色=注）

1 读入数据

> results <- read.csv("g:/ex/r/oneway.csv")

> attach(results)

> names(results)

[1] "sand" "clay" "loam"

> results

    sand  clay  loam

1     6   17    13

2    10   15    16

3     8    3     9

4     6   11    12

5    14   14    15

6    17   12    16

7     9   12    17

8    11    8    13

9     7   10    18

10   11   13    14

 

2 计算均值

> sapply(list(sand,clay,loam),mean)

[1]  9.9  11.5  14.3

 

3 计算方差

> sapply(list(sand,clay,loam),var)

[1]  12.544444  15.388889  7.122222 

 

 

 

4 等方差检验

> y<-c(sand,clay,loam)

> soil<-factor(rep(1:3,c(10,10,10)))

> soil

 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3

Levels: 1 2 3

> fligner.test(y~soil)

 

        Fligner-Killeen test of homogeneity of variances

 

data:  y by soil

Fligner-Killeen:med chi-squared = 0.3651, df = 2, p-value = 0.8332

 

 ##  there is no evidence of any significant difference in variance across the three samples，所以接下来可以进行one-way anova

 

> bartlett.test(y~soil)

 

## 实际上 Bartlett 检验是否属于正态分布的作用强于等方差检验

 

        Bartlett test of homogeneity of variances

 

data:  y by soil

Bartlett's K-squared = 1.2764, df = 2, p-value = 0.5283

 

5 方差分析：函数aov, summary

> summary(aov(y~soil))

           Df   Sum Sq    Mean Sq    F value    Pr(>F) 

soil         2    99.2       49.60      4.245     0.025 *

Residuals   27    315.5      11.69                

---

Signif. codes:  0   ‘***’ 0.001   ‘**’ 0.01   ‘*’ 0.05   ‘.’ 0.1   ‘ ’ 1

 

# 手工制表，一般是下表的模样。显然 summary.aov 没有给出 Total

Source	Sum of squares	Degrees of freedom	Mean square	F ratio	p value
Soil type	99.2	27	49.6	4.24	0.025
Error	315.5	2	11.685
Total	414.7	29

 

 

 

 

 

6 影响大小：函数aov, summary.lm

> summary.lm(aov(y~soil))

 

Call:

aov(formula = y ~ soil)

 

 

Residuals:

   Min   1Q   Median   3Q   Max

   -8.5  -1.8     0.3     1.7   7.1

 

 

 

# intercept指产量均值for sand，sand是第一个因子水平；soil2 = mean(clay)－mean(sand)； # soil3 = mean(loam)－mean(sand)。可见第二、三行均分别以第一个因子水平做基准进行计算。

 Coefficients:

            Estimate   Std. Error   t value      Pr(>|t|)   

(Intercept)     9.900     1.081      9.158      9.04e-10 ***

soil2         1.600     1.529      1.047       0.30456   

soil3         4.400     1.529      2.878       0.00773 **

---

Signif. codes:  0   ‘***’ 0.001   ‘**’ 0.01   ‘*’ 0.05   ‘.’ 0.1   ‘ ’ 1

 

 

Residual standard error: 3.418 on 27 degrees of freedom

Multiple R-squared: 0.2392,     Adjusted R-squared: 0.1829

F-statistic: 4.245 on 2 and 27 DF, p-value: 0.02495

 

# 首先来关注Coefficients: 跟拟合线性模型时给出的表相似，但含义大不一样。

1） Intercept 表示均值。本例中代表sand上栽培作物的平均产量，因为sand是第一个因子水平。应注意，通常情况下intercept是按字母顺序排在最低位置的那个因子水平。

2） soil2 表示均值之差。这里指mean(clay)－mean(sand)。clay是第二个因子水平。

3） 同理soil3指mean(loam)－mean(sand)。loam是第三个因子水平。

 

# 显然，没有给出soil2和soil3的比较结果，但是可以口算出来，mean(loam)－mean(clay)=4.4－1.6=2.8，则 t 值=2.8/1.529=1.83 < qt(0.975, 18) = 2.100922，所以二者差异不显著。根据 t 值，可以计算出

 p-value = 2*(1－pt(1.83,df=18)) = 0. 0838617，这是双尾检验(two-tailed test)。

 

********************************************************************************

ANOVA with aov or lm

函数aov和lm的主要差别在于输出结果形式上的不同。

但是，如果存在多重误差项（multiple error terms），则必须使用aov，因为 lm 不支持 Error 项。

********************************************************************************

 

以上是书中实例，实际应用时，通过调整数据存储样式，可简化操作：调整后的“.csv”格式文件见下面的数据表（只列出了一小部分）。

其中，soil = 1，代表土壤类型sand；

     soil = 2，代表clay；

     soil = 3，代表loam。

 

yield	soil
6	1
10	1
8	1
6	1
14	1
···	···

 

 

 

 

 

 

 

 

用上表进行分析时，操作如下：

 

> mydata <- read.csv("g:/ex/r/oneway2.csv")

> attach(mydata)

> myaov <- aov(yield~as.factor(soil))

## 与myaov2 <- aov(yield~soil)完全相同，why？

> summary(myaov)

## 得到anova报表，跟上面的anova分析结果完全相同

> TukeyHSD(myaov)

## anova 结果表明土壤类型不同，单产有显著差异，进一步做两两比较，找出哪些类型之间差异显著，结果见下图。此外，TukeyHSD is conservative without being ridiculously.

## 土壤类型 1 和 3 差异显著（ p = 0.020  ）

 

 

http://s15/mw690/b5c8908cgd3ec0671280e&690ANOVA）：R and Stata" TITLE="单向方差分析（One-way ANOVA）：R and Stata" />

 

（二）用 Stata 处理之

 

    用 R 很简单，用 stata 软件也是如此：

 

    读入文件：  insheet using "G:\ex\r\oneway2.csv"

 

    单向方差分析：  anova yield soil

 

    专用 one-way ANOVA 命令：  oneway yield soil

 

    ## 两种命令得到的 anova 结果相同，不同的是，oneway 同时给出了 Bartlett’ test 结果

 

 

Stata基于Tukey方法的多重比较（见下面的命令和运行后的结果图）：

 

pwcompare soil, effects sort mcompare(tukey)

 

http://s4/mw690/b5c8908cg7b979b4ae593&690ANOVA）：R and Stata" TITLE="单向方差分析（One-way ANOVA）：R and Stata" />

 

（三）Technique Notes  

 

1 R technique note

1.1 Syntax for aov

Description：

通过对每一层调用函数lm，拟合方差分析模型。

Usage：

aov(formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, ...)

Arguments：

•   formula       公式，用于指明模型

•   data          数据框，用于指定formula用到的数据。如果丢失，按标准方法搜索。

•   projections    逻辑标志：是否返回预测（projections）？

•   qr           逻辑标志：是否返回QR分解值？

•   contrasts     公式中做比较的数据清单。不能用于误差项，只在误差项中提供因子的比较会报错。

•   ...           lm 函数中的参数均可用，例如subset 或者na.action，详情见weights。

Details：

•   封装了根据 balanced 或者 unbalanced 实验设计拟合线性模型的lm 函数。

•   与lm 函数的主要区别是print, summary, 等的样式：主要用ANOVA的传统表达而不是线性模型的。

•   如果只包括单一误差项，用于指定误差层，并在每一误差层（error stratum）内拟合恰当的模型。

•   formula 可以指定多重响应（multiple response）。

Note：

•   aov 函数为 balanced designs 而设计，如果不平衡，结果就难以解释：当心响应数据丢失引起不平衡。如果有两个或者更多的误差层（error strata），不平衡该方法在统计学上是无效的，可能用nlme软件包里的lme函数更好。

•   可以用replications函数检查balance。

1.2 Syntax for TukeyHSD

Description：

Create a set of confidence intervals on the differences between the means of the levels of a factor with the specified family-wise probability of coverage. The intervals are based on the Studentized range statistic, Tukey's ‘Honest Significant Difference’ method.

Usage：

TukeyHSD(x, which, ordered = FALSE, conf.level = 0.95, ...)

Arguments：

•   x          拟合的模型对象，通常是用函数aov 拟合而成。

•   which      字符向量，列出拟合的模型中的项目，将要计算的是它们的区间。默认值为全部项目。

•   ordered     逻辑值，表示在取差值之前，因子的水平是否按令样本均值增加的顺序排列。如果为真，

计算得到的均值之差全为正值。显著性差异将是那些lwr end point 为正的。

•   conf.level   数值，介于零和1之间，给出供使用的the family-wise confidence level

•   ...         指其它供选参数，目前不可用。

Details：

在方差分析过程中，当按某个因子的水平计算均值时，简单的用 t 检验比较会令判定差异显著的概率放大，这样会把不显著的判定为显著。原因在于用给定的概率范围计算区间作为每一个区间，但是在解释范围的时候经常用的是关于区间的整个“族”。

John Tukey 基于样本均值的范围而不是个体差异引入区间。本函数给出区间返回值是建立在学生氏化的幅度统计量（Studentized range statistics）的基础上。

用该方法建立的区间只能应用于balanced designs，balanced designs 的因子，每个水平的观测数是相等的。本函数包含对样本大小的调整，对于略微 unbalanced designs，也能生成合理的区间。

如果指定了非因子项，就会报错：if no terms are left this is a an error.

 

2 Stata technique note

2.1 Syntax for anova

anova  varname  [termlist]  [if]  [in]  [weight]  [ , options]          分析方差和协方差

其中：termlist 是因子变量列表，具有以下性质

    假定变量是分类的；使用c. 因子变量操作可以修改变量的类型

 符号“ | ”（nesting的标志）可以换成符号“ # ”（interaction的标志）

 符号“ / ”用在项目之后，表示接下来的项目是关于预测项的误差项

 

http://s15/mw690/b5c8908cg7b9993ab6d3e&690ANOVA）：R and Stata" TITLE="单向方差分析（One-way ANOVA）：R and Stata" />

2.2 Syntax for pwcompare

Menu: Statistics > Postestimation > Pairwise comparisons

Syntax：pwcompare marginlist [, options]

其中 marginlist 是因子变量或者交互作用的清单，它们出现在当前估计结果或者指向参考方程的_equns 中。可以录入变量，或者说不带有前缀i.，可以用任意因子-变量语法：

pwcomare i.sex i.group i.sex#i.group

pwcompare sex group sex#group

pwcompare sex##group

 

http://s4/mw690/b5c8908cgd3ffcc692293&690ANOVA）：R and Stata" TITLE="单向方差分析（One-way ANOVA）：R and Stata" />
* tukey, snk, duncan, and dunnett are only allowed with results from anova, manova, regress, and mvreg.  tukey, snk, duncan, and dunnett are not allowed with results from svy.

* Time-series operators are allowed if they were used in the estimation.

 

 

2.3 Syntax for margins

 

margins 是一个 postestimation 命令，在使用 regress 或者 logistic 等估计命令拟合模型之后使用margins命令。

 

margins 命令估计和报告响应变量的margins，响应变量的derivatives的margins，也称marginal effects. margin 是一个基于拟合的模型的统计量，在模型中一些或者所有协变量（covariates）是固定的。marginal effects 是响应随协变量发生变化时的变化，可以报告为导数（derivative），弹性（elasticity），或者半弹性（semi-elasticity）。

 

要理解 margin 是什么，margins命令怎么用，最好的办法是看些实例：

 

例1：输入下面三行命令，得到结果如下图。

use http://www.stata-press.com/data/r12/margex

regress y i.sex i.group

margins sex

 

http://s13/bmiddle/b5c8908cg047fe7289b9c&690ANOVA）：R and Stata" TITLE="单向方差分析（One-way ANOVA）：R and Stata" />

        结果解读：Margin 那一列是 y 的均值。基于线性回归 y~sex+group 可知，如果数据中的每个人都当作 male，y 的均值约为60.6，如果都当作 female，则均值为78.9.

 

例2：还是上面的数据集，经 logistic 回归后的 margins 应用。

logistic outcome i.sex i.group

margins sex

http://s6/mw690/b5c8908cgd3ffd3fca6b5&690ANOVA）：R and Stata" TITLE="单向方差分析（One-way ANOVA）：R and Stata" />

 

    结果解读：Margins 列是预测概率的均值。对outcome~sex+group 做 logistic 回归后，如果把数据中人都当作 male，outcome 的平均概率为0.13；若为 female，则 outcome 的平均概率为0.19.