SAS系统Contrast 语句详解

十二月 18, 2009 · 科技翻译, 统计计算 · (No comments)

我在"【总结】方差分析系列<一>------用好ANOVA中的Contrast选项"中曾简单介绍了contrast选项，下面我翻译并整理了SAS中contrast语句常数向量的获取方法，也适用于SAS的estimate语句和SPSS软件。

本文主要内容整理自"A quick guide to writing contrast statements in SAS and SPSS"。

Contrast语句使你可以用自定义的方式进行假设检验，它必须出现在model语句之后。如果用到manova语句、repeated语句、random语句或test语句，contrast语句必须出现在这些语句之前。标记用来标识所进行的检验，用以标识的文字或符号需用单引号括起来。效应表达式用以指定假设检验的因素（组合），这些因素（组合）必须是model语句中出现过的。效应表达式后的常数向量用以指定相应因素（组合）各水平的值，在指定各水平的情况下进行相关因素的分析。

下面将举例说明效应表达式常数向量获取过程。

1. 影响因素的顺序及编码方式

在SAS中，class语句指定了因素的顺序，这决定着contrast语句中常数向量的写法。

在默认情况下，因素按字母顺序或数值大小进行编码。例如，如果变量"sex"有"男性"和"女性"两个值，并分别用M和F表示。那么，程序将先为F（女性）计算系数。

2. 数据的描述统计

在SAS中，使用proc univariate可以使我们对数据集有个粗略的认识。输出结果里还包括盒子图。

3. 确定感兴趣的差异比较

假设一个完全随机区组试验设计，有A、B和C三个因素。A因素有2个水平(1和2)，B因素有3个水平(1、2和3)，C因素有4个水平(1、2、3和4)。 我们设定一个简单比较：当b为2时，比较A 2个水平之间的差异。

4. 新建一个系数表，并关注影响因素顺序及其水平顺序。

根据A和B因素在class语句里的排列顺序，A因素出现在行，而B因素出现在列，并把A和B因素水平按数值大小排列（如下图）。在表格左侧和顶部的单元格标明需比较水平的系数（红色斜体表示）。

5. 表格左侧的系数与顶部的系数两两相乘，并在相应单元格里标明结果。那么，A*B交互项的系数为0 1 0 0 -1 0（红色表示）。

6. 把交互项系数横向和纵向相加，并写到表格右侧和底部。

这些值是A(表格右侧)和B(表格底部)主效应系数。由于B主效应的系数均为零，可以不出现在contrast语句。对于A，我们必需按顺序写为1 -1。

7. 表格右侧和底部系数相加，即得截距的系数，并写到右下角单元格。因为其值为零，截距不出现在contrast语句中。

8.最后，即得contrast语句。

对于SAS系统的proc glm 或proc mixed过程步，contrast语句可写为：contrast "A1 v A2 in B2" A 1 -1 A*B 0 1 0 0 -1 0;

在SAS系统里，小于1的至少4位十进制小数也可作为常数向量 （如：0.3333）。

现在，设想同一试验数据，但我们对一个更复杂的差异比较感兴趣，这涉及到A、B和C三个因素。设定C的4水平是对照，那么你可以让1、2和3水平与4水平进行比较，而且是在B为2和A为1条件下。

下面，我们将以在B为2和A为1条件下C的1水平与4水平比较为例。

1. 首先取A和B两个因素，建立参数表。

从上图可知，contrast语句需要保留A、B以及截距项，但这一步还不能下这个结论。

2. 现在，检查A*B交互项与C，以找到A*B*C交互项的系数。

从上图可知，contrast语句不需保留截距和A*B交互项，而保留C。这里才能确定C项需要保留。

3. 由于C保留，所以我们检查A与C，B与C，以确定contrast语句还需要保留哪些项。

A*C

这样，contrast语句不保留A和截距，而将保留C。

B*C

因此，contrast语句将不保留B和截距，但C将保留。

4. 最后，还需要B*C与A和A*C与B，以确定最后两项是否进入contrast语句。

B*C与A

我们已经知道contrast语句不保留A和截距项，但B*C将保留。

A*C与B

contrast语句不保留B和截距，但保留A*C。

5. 到此为止，我们已经确定了contrast语句里需指定的效应：C, A*C, B*C, and A*B*C。

contrast语句将写为：

contrast "C1 v C4 in A1 and B2"
C 1 0 0 -1
A*C 1 0 0 -1 0 0 0 0
B*C 0 0 0 0 1 0 0 -1 0 0 0 0
A*B*C 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;

注：于非零系数之后的零可以省略。所以，A*B*C 0 0 0 0 1 0 0 -1与上面的意义相同。

6. 我们对其它比较也可能感兴趣，如水平2，3与4比较。

对于C 2 v 4 in A1 and B2

contrast "C2 v C4 in A1 and B2" C 0 1 0 -1
A*C 0 1 0 -1 0 0 0 0
B*C 0 0 0 0 0 1 0 -1 0 0 0 0
A*B*C 0 0 0 0 0 1 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;

对于 C 3 V C4 in A1 and B2

contrast "C3 v C4 in A1 and B2" C 0 0 1 -1
A*C 0 0 1 -1 0 0 0 0
B*C 0 0 0 0 0 0 1 -1 0 0 0 0
A*B*C 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;

延伸资料：

How can I do ANOVA contrasts?

如何检验两组回归系数之间的差异显著性？

十一月 20, 2009 · R语言专栏, 统计计算 · (No comments)

在 Li Jingwen 的这篇文章中，图3中显示了两个生育时期的线性回归模型。随后作者比较了两个生育时期线性回归模型的回归系数(斜率)和截距，作者发现两个生育时期回归系数(斜率)差异不显著，而截距差异显著。

这种两组或多组回归系数之间的差异性如何检验？如何在R软件中实现？为此，我总结了回归系数的比较方法，如下。

回归系数的比较通常可以分为两类，线性回归模型回归系数比较和非线性回归模型回归系数比较。

我们先谈谈线性回归模型回归系数比较，而本帖只针对上面的文献讲解两组回归系数之间的比较。多组线性回归模型的回归系数比较与两组之间比较类似，只是多了几个虚变量，而非线性回归系统比较则使用的是残差平方和简化测验(sum of square reduction test, SSRT)，你可以参考"不同株型小麦干物质积累与分配对氮肥响应的动态分析"。

我们虚构有一个数据集，有gender、height和weight三个变量，文件名为new.csv。

# 设置工作目录
setwd("E:\\My Documents\\R\\data")
#读取外部csv格式数据
mydata <- read.table(file="new.csv", header=TRUE, sep=",")
# 查看数据集
mydata

这样我们首先得到了线性回归方程。现在假定零假设Ho：Bf-Bm=0，其中Bf为女性组的回归系数，Bm为男性组的回归系数。

我们需要定义两个虚变量，虚变量female的值为1表示女性，为0表示男性。虚变量femht为female与女性身高的乘积。


上面回归系数的统计学意义如下：

INTERCEP 5.601677 : 男性组线性回归截距
FEMALE -7.999147 : 女性组线性回归截距 - 男性组线性回归截距
HEIGHT 3.189727 : 男性组斜率，即Bm.
FEMHT -1.093855 : 女性组斜率 - 男性组斜率

FEMHT项对应的是零假设Ho：Bf-Bm=0，从P值=3.5456e-15可知，拒绝零假设，表明女性组斜率与男性组斜率之间存在显著差异。

R与通径分析

十一月 11, 2009 · R语言专栏, 统计计算 · 2 comments

笔者写的"SAS与通径分析"在所有帖子中搜索的次数比较多，可以看出大家对这种数据处理方法的认知度。由于统计知识掌握程度的差异，如何方便快捷计算出结果可能是大多数人的想法。

今天笔者给大家介绍另一种实现通径分析的方法，使用R语言及agricolae软件包进行通径分析。

首先介绍一下agricolae软件包，这个是秘鲁国立工程大学硕士学位论文"农业研究统计分析工具"中编制的R语言软件包。

现在agricolae软件包在国际马铃薯中心(CIP)、拉莫林国立农业大学 (UNALM-PERU)和Instituto Nacional de Innovacion Agraria(INIA-PERU)使用较广泛。

大家从软件包的名称就可以看出，agricolae软件包针对农业和植物育种田间试验设计的统计分析功能，这包括格子设计、析因设计、完全区组和不完全区组设计、拉丁方设计、希腊拉丁方设计、Alpha设计、Cyclic设计、条裂区设计、多点试验比较、处理比较、重采样、模拟、多样性指数和共识集群。而通径分析需要用到agricolae软件包里的path.analysis()和correlation()两个函数。

我们还是以" SAS与通径分析"中的数据为例子进行计算，大家可以把两个代码的计算结果进行对比。

实验数据，包括4个自变量和1个因变量

y,x1,x2,x3,x4
15.7,10,23,3.6,113
14.5,9,20,3.6,106
17.5,10,22,3.7,111
22.5,13,21,3.7,109
155,10,22,3.6,110
16.9,10,23,3.5,103
8.60,8,23,3.3,100
17.0,10,24,3.4,114
13.7,10,20,3.4,104
13.4,10,21,3.4,110
20.3,10,23,3.9,104
102,8,21,3.5,109
7.40,6,23,3.2,114
11.6,8,21,3.7,113
12.3,9,22,36,105

计算过程：

> mydata <- read.table(file="E:\\My
Documents\\R\\data\\pathanalysis.csv", header=TRUE, sep=",") #读取外部数据文件pathanalysis.csv
> x <- mydata[,-1] #把自变量x1、x2、x3和x4从数据框mydata中提出，赋值给x
 > y <- mydata[,1] #把因变量y从数据框mydata中提出，赋值给y
> cor.y <- correlation(y,x)$correlation #计算向量y与向量 x的相关系数

Correlation Analysis

Method : pearson
Alternative: two.sided

> cor.x <- correlation(x)$correlation #计算向量x与向量 x的相关系数

Correlation Analysis

Method : pearson
Alternative: two.sided

> path.analysis(cor.x,cor.y) #进行通径分析
 Correlations
============
x1 x2 x3 x4
x1 1.00 -0.14 -0.06 -0.09
x2 -0.14 1.00 0.01 0.12
x3 -0.06 0.01 1.00 -0.21
x4 -0.09 0.12 -0.21 1.00

============
y
x1 0.04
x2 -0.10
x3 -0.11
x4 0.11

Direct(Diagonal) and indirect effect path coefficients
======================================================
x1 x2 x3 x4
x1 0.029523150 0.015115246 0.0050697909 -0.009708187
x2 -0.004133241 -0.107966043 -0.0008449652 0.012944249
x3 -0.001771389 -0.001079660 -0.0844965151 -0.022652435
x4 -0.002657083 -0.012955925 0.0177442682 0.107868740

Residual Effect^2 = 0.9668623

附path.analysis()和correlation()的详细使用说明。

path.analysis 通径分析

描述

如果因果关系明确，就可能使用图表形式即通径分析表示整个变量系统。这个函数可以计算直接效应和间接效应，并使用变量相关性和协方差。

用法

path.analysis(corr.x, corr.y)

自变量

corr.x 自变量相关矩阵

corr.y 与每个自变量的向量相关系数

详细信息

首先必须计算出相关系数

自变量值

corr.x 数值型

corr.y 数值型

作者

Felipe de Mendiburu

参考文献

Biometrical Methods in Quantitative Genetic Analysis, Singh, Chaudhary. 1979

参见

相关

实例

# Path analysis. Multivarial Analysis. Anderson. Prentice Hall, pag 616

library(agricolae)

# Example 1

corr.x<- matrix(c(1,0.5,0.5,1),c(2,2))

corr.y<- rbind(0.6,0.7)

names<-c("X1","X2")

dimnames(corr.x)<-list(names,names)

dimnames(corr.y)<-list(names,"Y")

path.analysis(corr.x,corr.y)

# Example 2

# data of the progress of the disease related bacterial wilt to the ground

# for the component CE Ca K2 Cu

data(wilt)

data(soil)

x<-soil[,c(3,12,14,20)]

y<-wilt[,14]

cor.y<-correlation(y,x)$correlation

cor.x<-correlation(x)$correlation

path.analysis(cor.x,cor.y)

correlation 相关分析 皮尔逊法、斯皮尔曼法和Kendall法

描述

获得所有变量的的相关系数及P值。计算方法有皮尔逊法、斯皮尔曼和Kendall法。若未指定计算方法，皮尔逊法是默认值。计算结果与SAS计算结果相似。

用法

correlation(x,y=NULL, method = c("pearson", "kendall", "spearman") ,alternative="two.sided")

自变量

x 表、矩阵或向量

y 表、矩阵或向量

method "pearson", "kendall", "spearman"

alternative "two.sided", "less", "greater"

详细信息

参数与函数cor()相同

自变量值

table 数值型

作者

Felipe de Mendiburu

参见

correl

实例

library(agricolae)

# example 1

data(soil)

analysis<-correlation(soil[,2:8],method="pearson") analysis

# Example 2: correlation between pH, variable 2 and other elements from soil.

data(soil)

attach(soil)

analysis<-correlation(pH,soil[,3:8],method="pearson",alternative="less")

analysis

detach(soil)

# Example 3: correlation between pH and clay method kendall.

data(soil)

attach(soil)

correlation(pH,clay,method="kendall", alternative="two.sided")

detach(soil)

R语言AIS软件包详解－－look函数

十一月 8, 2009 · R语言专栏, 统计计算 · (No comments)

Micah Altman于2006年写了R语言软件包AIS，包括四个函数DWIMxtabs、look、HTML.DWIMxtabs和print.DWIMxtabs，利用该软件包可以对数据集进行初步的探索性分析。Agri521翻译了该包的手册，感兴趣的可以到 这里 下载。

本文将以look函数为例进行详细的讲解，并列出统计结果，和大家一同学习，共同进步。

look函数可以对多元结构数据集有个初步的认识，这包括数据集的最大值、最小值、均值、中位数、百分位数、主成分分析、相关矩阵、相关矩阵热力地图(a heatmap of the correlation matrix)、频数分布图、平行坐标图(parallel coordinates plot)。

代码及结果如下：

> data(swiss) #打开R软件随带的数据集swiss
> swiss #显示数据集swiss，若数据集较大，可以使用head(swiss) ，仅显著部分数据

> library(AIS)　#打开AIS包
> fl <- look(swiss) #look函数赋给fl
> print(fl)　#输出对象fl

Fertility Agriculture Examination Education
Min. :35.00 Min. : 1.20 Min. : 3.00 Min. : 1.00
1st Qu.:64.70 1st Qu.:35.90 1st Qu.:12.00 1st Qu.: 6.00
Median :70.40 Median :54.10 Median :16.00 Median : 8.00
Mean :70.14 Mean :50.66 Mean :16.49 Mean :10.98
3rd Qu.:78.45 3rd Qu.:67.65 3rd Qu.:22.00 3rd Qu.:12.00
Max. :92.50 Max. :89.70 Max. :37.00 Max. :53.00
Catholic Infant.Mortality
Min. : 2.150 Min. :10.80
1st Qu.: 5.195 1st Qu.:18.15
Median : 15.140 Median :20.00
Mean : 41.144 Mean :19.94
3rd Qu.: 93.125 3rd Qu.:21.70
Max. :100.000 Max. :26.60


Importance of components(成分重要性):
Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
Standard deviation 1.7887865 1.0900955 0.9206573 0.66251693 0.45225403
Proportion of Variance 0.5332928 0.1980514 0.1412683 0.07315478 0.03408895
Cumulative Proportion 0.5332928 0.7313442 0.8726125 0.94576729 0.97985624
Comp.6
Standard deviation 0.34765292
Proportion of Variance 0.02014376
Cumulative Proportion 1.00000000

(相关矩阵)
Fertility Agriculture Examination Education Catholic
Fertility 1.0000000 0.2426643 -0.66090300 -0.44325769 0.41364556
Agriculture 0.2426643 1.0000000 -0.59885994 -0.65046381 0.28868781
Examination -0.6609030 -0.5988599 1.00000000 0.67460383 -0.47505753
Education -0.4432577 -0.6504638 0.67460383 1.00000000 -0.14441631
Catholic 0.4136456 0.2886878 -0.47505753 -0.14441631 1.00000000
Infant.Mortality 0.4371367 -0.1521287 -0.05915436 -0.01898137 0.06611714
Infant.Mortality
Fertility 0.43713670
Agriculture -0.15212866
Examination -0.05915436
Education -0.01898137
Catholic 0.06611714
Infant.Mortality 1.00000000

> plot(fl) #绘制对象fl

# 若需保存图可以使用如下代码

png(filename="images/look_%03d.png" ,width=480, height=480)
plot(fl)
dev.off()

MS ACCESS与SAS Dataset数据传输的六种方法

十二月 5, 2008 · 统计计算 · 2 comments

Agri521最近对Scientific data management 非常感兴趣，这对将来的研究工作有很大的帮助。为此，本人设计了科研数据管理数据库，用于大量实验数据的存储；实现了与SAS Dataset的无缝传输，用于大量数据的统计分析。

下面列出了MS ACCESS与SAS DATA TRANSFER的几种方法。

方法 附加软件 说明
CSV None 1. Save each Access table as csv file, then read in by DATA step.
2. Input statement should match the column ordering and data type in database table.
DDE None 1. Type each Access table name in DDE statement.
2. An input statement be written to match the layout of columns.
3. Invoke MS Access and open the table file to be input.
4. MS Access cannot have a '.' as missing value for numerics.
PROC ACCESS SAS/Access to PC file formats 1. Need to save Access file as a Microsoft Excel 5.0/95 Workbook file, but doesn't need input statement.
2. Column limit is 256, long value may be truncated with an error message.
ODBC 1. SAS/Access interface to ODBC.
2. Microsoft Access ODBC driver.		 1. Get ODBC installed and configured.
2. Define DSN reference to database(.mdb file)
Third party solution 第三方软件，如 DBMS/COPY, DBMS/Engine. 1. Buy and get software installed and configured.
2. Define some parameters if necessary.
PROC IMPORT SAS/Access to PC file formats Type each table name in SAS progrma before data
transfer.
← Older entries
苏ICP备08103697号
Agri521'Blog is prowdly powered by Wordpress · © 2009 · All rights reserved.
OrangeJuice theme by Theme Museum. 登录