denver
期待楼主带领我们继续学习
areg
由于昨天去研究SPSS保存格式*.sav的导入问题，没有来得及帧出后帧，对不起啊！今天帧出
#############################################################
#
# CHAPTER 1
# An Introduction to R
# 1.6 Basic Data Manipulation
#
############################################################
############################################################
首先把该包相关材料载入，把修改部分也重新执行，学习中，我们不可能一次把所有的练习就学完，假如一口气学完，太多了，也记不住。当再次进入学时，当你重启R来学习时，以前执行过的命令，已经不而内存中。如果不这样做，往往会在后面学习过程中，命令执行不出结果，这是许多初学者在学习时常见的问题，以为学习资料中的材料错误啦，这是许多人常犯的错误。
############################################################
library("lattice"); library("MASS"); library("scatterplot3d"); library("HSAUR")
data("Forbes2000", package = "HSAUR")
Forbes2000[32,2]<-"FMISN1"
Forbes2000[54,2]<-"FMISN2"
Forbes2000[883,2]<-"FMISN3"
# 下面接着学习新内容的示例操作
companies <- Forbes2000[, "name"] # 把数据集Forbes2000的变量name赋予向量companies
companies[1] # 提取新向量companies的第一个元素
[1] "Citigroup"
>1:3 # 显示向量数据1至3
[1] 1 2 3
>companies[1:3] # 提取向量companies中，即name中，的第1至3行下标的元素，与上条命令比较，明白？
[1] "Citigroup" "General Electric" "American Intl Group"
> companies[-(4:2000)] # 提取向量companies中不包括第4至2000个的元素。
[1] "Citigroup" "General Electric" "American Intl Group"
> companies[-(4:3000)] # 补充，所用下标，不能超出向量的长度。
错误: 下标出界
> length(companies) # 补充，复查向量长度。
[1] 2000
# 用下标提取所需要目标的相关信息
> Forbes2000[1:3, c("name", "sales", "profits", "assets")]
name sales profits assets
1 Citigroup 94.71 17.85 1264.03
2 General Electric 134.19 15.59 626.93
3 American Intl Group 76.66 6.46 647.66
# 如果你在提取变量名单中，所列变量在数据中没有，如现在添加一个没有的变量mystudy，显示如下或显示
> Forbes2000[1:3, c("name", "sales", "profits", "assets"，"mytudy")]
错误: "Forbes2000[1:3, c("name", "sales", "profits", "assets"
> name sales profits assets
错误: " name sales"里有语法错误
> 1 Citigroup 94.71 17.85 1264.03
错误: "1 Citigroup"里有语法错误
> 2 General Electric 134.19 15.59 626.93
错误: "2 General"里有语法错误
> 3 American Intl Group 76.66 6.46 647.66
# 下面两个命令提取的功能相同。
> companies <- Forbes2000$name
> companies <- Forbes2000[, "name"]
>order_sales <- order(Forbes2000$sales) # 把变量sales升序排列，注意order_sales中的下划线，它也是一个赋值符号
> companies[order_sales[1:3]]
[1] "Custodia Holding" "Central European Media" "Minara Resources"
# 提取销售额前三强的相关信息,注意c(2000, 1999, 1998)，按降序排列提取。
> Forbes2000[order_sales[c(2000, 1999, 1998)], c("name","sales", "profits", "assets")]
name sales profits assets
10 Wal-Mart Stores 256.33 9.05 104.91
5 BP 232.57 10.27 177.57
4 ExxonMobil 222.88 20.96 166.99
# 可用逻辑表达式来解析下标。
>Forbes2000[Forbes2000$assets > 1000, c("name", "sales","profits", "assets")]
name sales profits assets
1 Citigroup 94.71 17.85 1264.03
9 Fannie Mae 53.13 6.48 1019.17
403 Mizuho Financial 24.40 -20.11 1115.90
# 逻辑运算，2000个大公司中，仅有3个的assets超过1000
> table(Forbes2000$assets > 1000)
FALSE TRUE
1997 3
# 判断没有上报盈利的企业数
> na_profits <- is.na(Forbes2000$profits)
> table(na_profits)
na_profits
FALSE TRUE
1995 5
# 提取没有上报盈利的企业相关信息
> Forbes2000[na_profits, c("name", "sales", "profits","assets")]
name sales profits assets
772 AMP 5.40 NA 42.94
1085 HHG 5.68 NA 51.65
1091 NTL 3.50 NA 10.59
1425 US Airways Group 5.50 NA 8.58
1909 Laidlaw International 4.48 NA 3.98
# 列出表中记录不完整的数目，即记录中缺失值
> table(complete.cases(Forbes2000))
FALSE TRUE
5 1995
# 提取Forbes2000中英国公司。
> UKcomp <- subset(Forbes2000, country == "United Kingdom")
> dim(UKcomp)
[1] 137 8 # 显示结果为英国公司有137个，每个有8个变量指标
# 提取中国公司，此部分由人areg示例，中国人当然关心中国的，如下，
> CHcomp<-subset(Forbes2000, country == "China")
> dim(CHcomp)
[1] 25 8
>CHcomp
rank name country category sales profits assets marketvalue
55 55 PetroChina China Oil & gas operations 29.53 5.67 58.36 90.49
81 81 China Petroleum & Chemical China Oil & gas operations 39.16 1.94 45.32 50.09
202 202 China Telecom China Telecommunications services 9.12 2.04 24.85 29.92
…………………… # 节省显示页面，此部分由本人省略
1973 1973 China Southern Airlines China Transportation 2.18 0.07 4.49 3.16
areg
#############################################################
#
# CHAPTER 1
# An Introduction to R
# 1.7 Simple Summary Statistics
#
############################################################
library("lattice"); library("MASS"); library("scatterplot3d"); library("HSAUR")
data("Forbes2000", package = "HSAUR")
Forbes2000[32,2]<-"FMISN1"
Forbes2000[54,2]<-"FMISN2"
Forbes2000[883,2]<-"FMISN3"
# 上面几个命令，由于重新启动R系统，只好重复再次操作，以保障后续操作得以完成。
>summary(Forbes2000) # 对数字变量给出“五点描述”，对字符变量指出其是字符，对计算变量给出频数，显示略
>lapply(Forbes2000, summary) # 更具体而敏感地给出每一因子水平的描述性统计。显示略。
# 分类因子各水平的中位数，注意na.rm=TRUE，因为在利润(profits)列中有缺失值，否则给出结果为NA。在SPSS中，算术表达式计算中有缺失，其结果为NA，而函数表达式则正常运算，是可能其在函数表达中默认设置有na.rm=TRUE，许多时候，我们编写的函数或表达式不能正常进行，往往就是缺失值设置不当或没有设置造成的，这在S-PLUS中同理。
>mprofits <- tapply(Forbes2000$profits, Forbes2000$category,median, na.rm = TRUE)
# 结果为各行业利润的中位数，显示略。在上面tapply()函数中，我们也可以尝试求平均数，而平均数是对正态分布描述集中趋势的有效指标,而本例是非正态分布，选用中位数更好恰当；平均数对极端正值敏感，而中位数不敏感，在这种极不对称的情况下，中位数指标更为恰当。
> median(Forbes2000$profits)
[1] NA # 缺失值造成不运算
> median(Forbes2000$profits,na.rm=TRUE)
[1] 0.2 # 不按行业分，只求这个数据集利润变量的中位数，排除了缺失值。
# 上面的命令由areg添加，在学习过程中，大家多试试各种与当前命令相关的命令，以达到对某个参量或参数的熟悉。
> rev(sort(mprofits))[1:3]
Oil & gas operations Drugs & biotechnology Household & personal products
0.35 0.35 0.31
# sort()是升序排列，rev()是隆序排列，往往原文资料说成是升序排列的反向排列，上面命令是从上面对象mprofits所求得的各行业中位数提取利润最高的三个行业。即把降序排列的前三位列出。
#############################################################
#
# 1.7 Simple Summary Statistics
# 1.7.1 Simple Graphics
#
############################################################
> fm <- marketvalue ~ sales # ? 加权关系？
> class(fm)
[1] "formula"
# 上面把一段长时间内之平均市价值(marketvalue)与销售值的关系赋予对象fm，即组成两个变量间的公式，由于class(fm)给出它的属性归类。上面给出是一个公式。
> layout(matrix(1:2, nrow = 2)) # 将调出一个图层面板，此面反板上可以同时列出两个图形，如果不设，第二个图形的出现，将覆盖第一个图形，即每一个图形消失，只留下第二个由当前命令得到的图形。因为nrow=2，如果你接着作第三个图形，前两个同时消失，只出现第三个，还可以再同时作出第四个。这样，关于layout()中的参数设置应该明白了吧？修改参数试试。
> hist(Forbes2000$marketvalue) # 给出marketvalue为横坐标的图形，显示与手册的一致，是一个极偏态分布，表现为marketvalue主要由排名前几位的公司占有，垄断！见手册。
> hist(log(Forbes2000$marketvalue)) # 前一个命令的横坐标用原始值来划分，刻度很大，并且不易得出准确的信息，或给出的信息不直观。现在把marketvalue进行对数转换，得出的图形更有参考价值。是否需要转换要根据数据分布特点来决定。
> plot(log(marketvalue) ~ log(sales), data = Forbes2000,pch = ".") # 得到以对数转换后，marketvalue为y轴、sales为x轴的图形。如图不转换，如下条由我添加的命令，得到的图形，表现出的信息很低。
> plot(marketvalue~sales, data = Forbes2000,pch = ".")
denver
今天学完了，非常感谢楼主，这种学习方式对我而言效果非常好。
areg
#######################################################
#
# CHAPTER 1
# An Introduction to R
# 1.8 Organising an Analysis
#
#######################################################
#####如果你的R系统还没有重启，下面几条命令可以不重设#####
library("lattice"); library("MASS"); library("scatterplot3d"); library("HSAUR")
data("Forbes2000", package = "HSAUR")
Forbes2000[32,2]<-"FMISN1"
Forbes2000[54,2]<-"FMISN2"
Forbes2000[883,2]<-"FMISN3"
# 上面几个命令，由于重新启动R系统，只好重复再次操作，以保障后续操作得以完成。此部分内容，以后各节不再重复，望学习者注意。
# 在R的命令窗口直接执行进行运算能完成任务，但是如果我们把要分析某数据库的执行步骤单独建个文件来保存，以便下次使用，这种方式更为方便可取。例如建个以.R为结尾的源文件analysis.R，要用时可以用源函数source来打开。可以把你进行的数据分析的步骤如数据进程、转换、描述性统计、构建模型、报表等步骤保存起来。
# 例如，本手册学习中，我们每次启动支要载入相关的包，要把几个错误之处用命令修正，我们就可以这些步骤以Forbes2000.R来保存，即从R窗口"文件_建立新的脚本程序"来打开脚本编辑器，例如把下面的步骤复制到脚本窗口中
library("lattice"); library("MASS"); library("scatterplot3d"); library("HSAUR")
data("Forbes2000", package = "HSAUR")
Forbes2000[32,2]<-"FMISN1"
Forbes2000[54,2]<-"FMISN2"
Forbes2000[883,2]<-"FMISN3"
然后从“文件_保存”弹出保存窗口，指定到想保存的文件夹中去，如本例保存为Forbes2000,文件类型选择为.R，其它类型以后再学习。当你要用到这些操作步骤的时候，可用source()函数去调出来，也可以在窗口中从“文件_输入R代码”弹出的窗口中去找你想输入的.R脚本文件，如我在D:/Rstudy文件夹中找到Forbes2000.R，一旦点击文件，.R的程序及步骤载入，并且在控制台窗口出现下面的
信息：
> source("D:\\Rstudy\\Forbes2000.R")
# 实际上，以后操作熟练啦，假如不想用窗口调入*.R，就用上面的这个命令source()去调吧。
areg
#########################################################
#
# CHAPTER 1
# An Introduction to R
# 1.9 Summary
#
#########################################################
# 文本小总结，读大家参阅原文。把数据读入R的有很多不同的方式，包括直接从数据库读入。关于多图形和数据的处理随后的章节将进一步介绍，此处只继续手册的示例。
>boxplot(log(marketvalue) ~ country, data = subset(Forbes2000,country %in% c("United Kingdom", "Germany","India", "Turkey")), ylab = "log(marketvalue)",varwidth = TRUE)
# 此命令得到的图是以每个国家x轴的盒形图，所列的四个国家盒形图很窄,并且因国家名太多，这四个国家名没有显示出来。
# %in% 表示从某个向量提取子向量或元素。前面命令是从数据集中提取四个国家的全部公司，如果把它反过来呢？
>boxplot(log(marketvalue) ~ country, data = subset(Forbes2000,c("United Kingdom", "Germany","India", "Turkey") %in% country), ylab = "log(marketvalue)",varwidth = TRUE)
# 得到很多盒形图！
##########################################################
与手册显示不同，原因尚未找到。先做练习，以后回头解决问题#经多次尝试，感觉是在此包中对country进行了永久赋值
##########################################################
三月小鱼
LZ真热心!
大大好人啊 !
areg
#########################################################
#
# CHAPTER 1
# An Introduction to R
# Exercises
#
#########################################################
载入相关命令后，进行练习,请先自己练习后，再参考本例
##############################
########## Ex. 1.1：##########
> tapply(Forbes2000$profits, Forbes2000$country,median, na.rm = TRUE)
Africa Australia
-0.005 0.280
…………………………………………# 省略 by areg
Finland France
0.200 0.190
France/ United Kingdom Germany
-2.830 0.230
…………………………………………# 省略 by areg
Turkey United Kingdom
0.180 0.205
…………………………………………# 省略 by areg
United Kingdom/ South Africa United States
-0.100 0.240
Venezuela
0.120
> # 上面操作是按书中例题演变而来，目的是验证下面自己的操作是否正确。提取子集用$符号在此题中比用下标方便。
> tapply(Forbes2000$profits, Forbes2000$country== "United States",median, na.rm = TRUE)
FALSE TRUE
0.17 0.24
> tapply(Forbes2000$profits, Forbes2000$country== "United Kingdom",median, na.rm = TRUE)
FALSE TRUE
0.200 0.205
> tapply(Forbes2000$profits, Forbes2000$country== "France",median, na.rm = TRUE)
FALSE TRUE
0.20 0.19
> tapply(Forbes2000$profits, Forbes2000$country== "Germany",median, na.rm = TRUE)
FALSE TRUE
0.20 0.23
############################
######### Ex. 1.2: #########
> Gercomp<-subset(Forbes2000, country == "Germany")
> table(complete.cases(Gercomp))
TRUE
65
> Gercomp[Gercomp$profits < 0, c("name","country", "sales","profits", "assets")]
name country sales profits assets
350 Allianz Worldwide Germany 96.88 -1.23 851.24
364 Deutsche Telekom Germany 56.40 -25.83 132.01
397 E.ON Germany 37.95 -0.73 115.57
431 HVB-HypoVereinsbank Germany 40.52 -0.87 705.36
500 Commerzbank Germany 22.43 -0.31 437.86
798 Infineon Technologies Germany 7.18 -0.51 11.79
869 BHW Holding Germany 7.46 -0.38 117.96
926 Bankgesellschaft Berlin Germany 9.43 -0.74 182.69
1034 W&W-W黶tenrot Germany 7.57 -0.08 56.44
1187 mg technologies Germany 8.54 -0.13 6.45
1477 N黵nberger Beteiligungs Germany 3.00 -0.03 15.97
1887 SPAR Handels Germany 6.84 -0.40 1.64
1994 Mobilcom Germany 2.16 -3.62 8.67
> # 德国公司名中的乱码，不知R系统怎么搞的，俺不知道:(
##########或者用下面方法来提取#######################
> Forbes2000[Forbes2000$profits < 0 & Forbes2000$country == "Germany", c("name","country", "sales","profits", "assets")]
name country sales profits assets
350 Allianz Worldwide Germany 96.88 -1.23 851.24
364 Deutsche Telekom Germany 56.40 -25.83 132.01
397 E.ON Germany 37.95 -0.73 115.57
431 HVB-HypoVereinsbank Germany 40.52 -0.87 705.36
500 Commerzbank Germany 22.43 -0.31 437.86
798 Infineon Technologies Germany 7.18 -0.51 11.79
869 BHW Holding Germany 7.46 -0.38 117.96
926 Bankgesellschaft Berlin Germany 9.43 -0.74 182.69
1034 W&W-W黶tenrot Germany 7.57 -0.08 56.44
1187 mg technologies Germany 8.54 -0.13 6.45
1477 N黵nberger Beteiligungs Germany 3.00 -0.03 15.97
1887 SPAR Handels Germany 6.84 -0.40 1.64
1994 Mobilcom Germany 2.16 -3.62 8.67
#########方法多种，自己试吧#########################
############################
######### Ex. 1.3: #########
这第三练习题是道很好题，首先你需要了解百慕大的相关情况，如下：
Industries: tourism, finance, insurance, structural concrete products, paints, perfumes, pharmaceuticals, ship repairing
Exports:
total value: $56 million (2000 est.)
commodities: semitropical produce, light manufactures, reexports of pharmaceuticals
partners : UK 29.5%, US 9.8% (1997)
Imports:
total value: $739 million (2000 est.)
commodities: miscellaneous manufactured articles, machinery and transport equipment, food and live animals, chemicals
partners: US 34%, UK 9%, Mexico 8% (1997)
通过以上信息，应该初步判断出该地区是以旅游业、运输业（航海、航空）、保险业务、通信等主要产业，而那些什么重工业不是这个地区的发展对象，再从这两千个公司去选择，应该建议哪些公司去投资。这是一个较为主观的题，没有唯一答案。如本人推荐旅游业、运输业、保险业务、银行、通信业务的公司可以考虑到此地区去设分公司。
> levels(Forbes2000[, "category"]) # 首先例出这2000个公司分别属于哪些类！
[1] "Aerospace & defense" "Banking"
[3] "Business services & supplies" "Capital goods"
[5] "Chemicals" "Conglomerates"
[7] "Construction" "Consumer durables"
[9] "Diversified financials" "Drugs & biotechnology"
[11] "Food drink & tobacco" "Food markets"
[13] "Health care equipment & services" "Hotels restaurants & leisure"
[15] "Household & personal products" "Insurance"
[17] "Materials" "Media"
[19] "Oil & gas operations" "Retailing"
[21] "Semiconductors" "Software & services"
[23] "Technology hardware & equipment" "Telecommunications services"
[25] "Trading companies" "Transportation"
[27] "Utilities"
>
# 找出选择的业务类型："Banking" ,"Business services & supplies" ,"Drugs & biotechnology" ,"Food drink & tobacco" ,"Hotels restaurants & leisure" , "Insurance" ,"Telecommunications services" ,"Transportation"等等。
> myidea<-Forbes2000[Forbes2000$category == c("Banking" ,"Business services & supplies" ,"Drugs & biotechnology" ,"Food drink & tobacco" ,"Hotels restaurants & leisure" , "Insurance" ,"Telecommunications services" ,"Transportation"), c("name","country", "category", "sales","profits", "assets")]
> myidea[myidea$profits>1,c("name","country", "category", "sales","profits", "assets")]
name country category sales profits assets
1 Citigroup United States Banking 94.71 17.85 1264.03
14 Berkshire Hathaway United States Insurance 56.22 6.95 172.24
25 Wells Fargo United States Banking 31.80 6.20 387.80
41 Wachovia United States Banking 24.47 4.25 400.87
57 BBVA-Banco Bilbao Vizcaya Spain Banking 24.10 2.81 288.80
62 Allstate United States Insurance 32.15 2.73 134.14
138 Canon Japan Business services & supplies 24.76 1.61 23.34
175 BCE Canada Telecommunications services 14.70 1.40 30.35
185 SunTrust Banks United States Banking 7.07 1.33 125.39
192 Union Pacific United States Transportation 11.55 1.31 33.46
207 Sprint FON United States Telecommunications services 14.19 1.62 21.86
>
# 上面的关于选择哪些公司适合去百慕大群岛开展业务，在此仅是演示练习的思路:)
############################
######### Ex. 1.4: #########
# 由于初学，尽可能用短命令，长命令，一旦出错，不好找错误之处。
> ForbesTopp <- subset(Forbes2000, profits != "NA") # 先提取利润排50名的公司子集
> length(ForbesTopp$profits) # 查看子集profits长度,为下面确定降序提取利润前50名的最大值。
> order_profits <- order(ForbesTopp$profits) # 用利润变量对数据框升序排列，而sort排序是只对向量。
> Forbes50<-ForbesTopp[order_profits[c(1995:1945)], c("name", "country", "sales", "profits", "assets")]
> Forbes50 # 查看提取的公司是否是符合要求的
> layout(matrix(1:2, nrow = 2)) # 给图层指定每页显示的图形行列数
> hist(Forbes50$assets) # 此表达式结果图形不是正态分布
> hist(Forbes50$sales) # 此表达式结果图形不是正态分布
> hist(log(Forbes50$assets)) # 此表达式结果图形仍然不是正态分布
> hist(log(Forbes50$sales)) # 此表达式结果图形近似正态分布
# 关闭图层面板，执行下表达式，以使图形美观一些。
> plot(log(sales) ~ log(assets), data = Forbes50,pch = ".")
> plot(log(assets) ~ log(sales), data = Forbes50,pch = ".")
# 下面的命令作出的图形很差，留得后而专门章节再深究。猜测是在这个数据包中存在有country是永久赋值！
boxplot(log(变量名) ~ country, data = subset(Forbes50,country %in% c("某个国家名", "某个国家名", "某个国家名" )),
ylab = "log(变量名)",varwidth = TRUE)
############################
######### Ex. 1.5: #########
> companies_mean<-tapply(Forbes2000$sales, Forbes2000$country,mean, na.rm = TRUE)
> companies_mean
Africa Australia
6.820000 5.244595
…………………………………………# 省略 by areg
Finland France
10.291818 20.102063
France/ United Kingdom Germany
1.010000 20.781385
Greece Hong Kong/China
2.528333 2.044000
…………………………………………# 省略 by areg
Italy Japan
10.213902 10.190633
Jordan Kong/China
1.330000 5.717500
…………………………………………# 省略 by areg
Turkey United Kingdom
4.713333 10.445109
United Kingdom/ Australia United Kingdom/ Netherlands
10.010000 7.540000
United Kingdom/ South Africa United States
2.060000 10.058256
Venezuela
0.980000
>
> summary(Forbes2000[Forbes2000$profits > 5, c("name","country", "profits")])
name country profits
Length:37 United States :20 Min. : 5.120
Class :character Switzerland : 3 1st Qu.: 5.965
Mode :character United Kingdom: 3 Median : 6.845
China : 1 Mean : 8.203
France : 1 3rd Qu.: 8.922
(Other) : 4 Max. :20.960
NA's : 5 NA's : 5.000
>
#########################################################
请对R熟悉的朋友贡献出更简洁的算法，下次介绍第二章
谢谢初学者之间互相交流
##########################################################
areg
######################################################
#
# CHAPTER 2
# Simple Inference: Guessing Lengths,Wave Energy, Water Hardness, Piston Rings, and Rearrests of Juveniles
#
######################################################
2.1 Introduction # 手册中没有带，统计分析的基础理论要自己去补。
2.2 Statistical Tests # 手册中没有带
2.3 Analysis Using R # 从手册中R应该开始，以后各章都是同样，本人areg也很无柰。
######################################################
#
# CHAPTER 2
# 2.3.1 Estimating the Width of a Room
#
######################################################
# 把下面表达式复制到一个新脚本中，并命名为Tiro2.R，以方便本章学习。
library("lattice"); library("MASS"); library("scatterplot3d"); library("HSAUR")
data("roomwidth", package = "HSAUR")
#######################################################
> roomwidth # 先看看数据集是什么样的。
unit width
1 metres 8
2 metres 9
3 metres 10
………# 省略by areg
44 metres 40
45 feet 24
………# 省略by areg
113 feet 94
>
convert <- ifelse(roomwidth$unit == "feet", 1, 3.28) # 如果单位nuit是英尺（即不是米），返回值1，如果不是英尺，返回3.28。在此表达式为yes=1,no=3.28
> convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
> convert
[1] 3.28 3.28 3.28 3.28 3.28 3.28 3.28 3.28 3.28 3.28
……………………# 省略by areg
[41] 3.28 3.28 3.28 3.28 1.00 1.00 1.00 1.00 1.00 1.00
[51] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
……………………# 省略by areg
[101] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[111] 1.00 1.00 1.00
> # 回头看看数据集，转换正确，同时本个也理解为这是一种替换或变相赋值，在提取分类中用处很大，逻辑转换。
tapply(roomwidth$width * convert, roomwidth$unit, summary) # 统一单位为feet。下面是执行结果。
> tapply(roomwidth$width * convert, roomwidth$unit, summary)
$feet
Min. 1st Qu. Median Mean 3rd Qu. Max.
24.0 36.0 42.0 43.7 48.0 94.0
$metres
Min. 1st Qu. Median Mean 3rd Qu. Max.
26.24 36.08 49.20 52.55 55.76 131.20
> # 换种方法试
> summary(roomwidth)
unit width
feet :69 Min. : 8.00
metres:44 1st Qu.:16.00
Median :35.00
Mean :32.92
3rd Qu.:44.00
Max. :94.00
> # 因单位不统一，正确地给出了实际上错误的信息，下面新一个数据集myroom
########################################################
########################################################
convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
unit<-roomwidth$unit
width<-roomwidth$width * convert
myroom<-data.frame(unit, width)
> summary(myroom)
unit width
feet :69 Min. : 24.00
metres:44 1st Qu.: 36.08
Median : 44.00
Mean : 47.15
3rd Qu.: 51.00
Max. :131.20
>
Max. :131.20
> tapply(myroom$width, myroom$unit, summary)
$feet
Min. 1st Qu. Median Mean 3rd Qu. Max.
24.0 36.0 42.0 43.7 48.0 94.0
$metres
Min. 1st Qu. Median Mean 3rd Qu. Max.
26.24 36.08 49.20 52.55 55.76 131.20
> # 请与tapply(roomwidth$width * convert, roomwidth$unit, summary)比较，同时也比较summary(roomwidth)与summary(myroom),为继续操作练习，请删除内存中我们自己多加进去的多余的对象，以免与后面手册中的操作相混。
> rm(convert,myroom,unit,width) # 如果不想用rm()，关闭R，重启载入Tiro2.R,回到本章开始时节。把相关步骤也进行重新进行.
> ls()
[1] "convert" "roomwidth"
######################################################
######################################################
> ls() # 查看当前操作窗口中的对象是否存在。
[1] "convert" "roomwidth"
> layout(matrix(c(1, 2, 1, 3), nrow = 2, ncol = 2,byrow = FALSE))
# 决定图层面板的摆放，此由一个矩阵来决定
> matrix(c(1, 2, 1, 3), nrow = 2, ncol = 2,byrow = FALSE)
[,1] [,2]
[1,] 1 1
[2,] 2 3
# 由矩阵决定第1个图在第1行占两列，第2图在第2行、第1列，第3图在第2行、第2列。
# byrow = FALSE是按列排列，=TRUE是按行。
> boxplot(I(width * convert) ~ unit, data = roomwidth, ylab = "Estimated width (feet)", varwidth = TRUE, names = c("Estimates in feet", "Estimates in metres (converted to feet)"))
# 绘制两种尺度推断精确程序的盒形图，注意横坐标认定。
> feet <- roomwidth$unit == "feet" # 对尺度因子分水平（英尺和米）
> qqnorm(roomwidth$width[feet], ylab = "Estimated width (feet)") # 画用英尺为单位的QQ图
> qqline(roomwidth$width[feet]) # 画用英尺为单位的QQ拟合线
> qqnorm(roomwidth$width[!feet], ylab = "Estimated width (metres)") # 画用米（非英尺）为单位的QQ图
> qqline(roomwidth$width[!feet]) # 画用米为单位的QQ拟合线
> # 上面各步骤均能顺利操作，关于图形的解释，请见相关统计教材，这本手册相关部分没有提供，毕竟免费的包啊！
######################################################
## 补充：手册上命令能进行正常操作，请把补充的与手册相比较，手册上的更为简洁，命令操作跟随在图形后面。
######################################################
对两个尺度进行独立样本t检验
> convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
> unit<-roomwidth$unit
> width<-roomwidth$width * convert
> t.test(width~unit, alternative='two.sided', conf.level=.95, var.equal=FALSE, data=roomwidth)
Welch Two Sample t-test
data: width by unit
t = 14.9557, df = 109.909, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
24.00599 31.33986
sample estimates:
mean in group feet mean in group metres
43.69565 16.02273
># 结果错在哪里?在t检验的公式中的变量虽然表面上看，我们已经进行赋值了，但是检验公式引用的数据集是roomwidth里的原始数据，它们的单位并不统一。看下面：
> convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
> unit<-roomwidth$unit
> width<-roomwidth$width * convert
> myroom<-data.frame(unit, width)
> t.test(width~unit, alternative='two.sided', conf.level=.95, var.equal=FALSE, data=myroom)
Welch Two Sample t-test
####################假设方差不等#########################
data: width by unit
t = -2.3071, df = 58.788, p-value = 0.02459
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-16.54308 -1.17471
sample estimates:
mean in group feet mean in group metres
43.69565 52.55455
####################假设方差相等#########################
> t.test(width~unit, alternative='two.sided', conf.level=.95, var.equal=TRUE, data=myroom)
Two Sample t-test
data: width by unit
t = -2.6147, df = 111, p-value = 0.01017
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-15.572734 -2.145052
sample estimates:
mean in group feet mean in group metres
43.69565 52.55455
> # 从t检验结果看出，这两种尺度用在测量标度上会带显著性差异(P=0.05)
##############回到手册中的操作#######################
########### 假设方差相等 #############
> convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
> t.test(I(width * convert) ~ unit, data = roomwidth, var.equal = TRUE)
Two Sample t-test
data: I(width * convert) by unit
t = -2.6147, df = 111, p-value = 0.01017
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-15.572734 -2.145052
sample estimates:
mean in group feet mean in group metres
43.69565 52.55455
> ########### 假设方差不等 #############
> t.test(I(width * convert) ~ unit, data = roomwidth, var.equal = FALSE)
Welch Two Sample t-test
data: I(width * convert) by unit
t = -2.3071, df = 58.788, p-value = 0.02459
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-16.54308 -1.17471
sample estimates:
mean in group feet mean in group metres
43.69565 52.55455
> ########### Wilcoxon秩和t检验 #############
> wilcox.test(I(width * convert) ~ unit, data = roomwidth, conf.int = TRUE)
Wilcoxon rank sum test with continuity correction
data: I(width * convert) by unit
W = 1145, p-value = 0.02815
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-9.3599953 -0.8000423
sample estimates:
difference in location
-5.279955
>
下节介绍 2.3.2 Wave Energy Device Mooring
denver
请问楼主，boxplot中varwidth=TRUE能够实现什么功能？我查了帮助，没有查到。将其去掉，输出的图形又没有看出有什么区别。请楼主明示，谢谢
areg
## Default S3 method:
boxplot(x, ..., range = 1.5, width = NULL, varwidth = FALSE,
notch = FALSE, outline = TRUE, names, plot = TRUE,
border = par("fg"), col = NULL, log = "",
pars = list(boxwex = 0.8, staplewex = 0.5, outwex = 0.5),
horizontal = FALSE, add = FALSE, at = NULL)
width a vector giving the relative widths of the boxes making up the plot.
varwidth if varwidth is TRUE, the boxes are drawn with widths proportional to the square-roots of the number of observations in the groups.
############################################################
上面是它的帮助文件的说明，我把varwidth=TRUE，=FALSE都试过，没有看出区别。但在有的数据集的图中，又历区别，好象是如果=FALSE，各个变量在横坐标上的图等宽，这是许多统计软件中的默认值，如果=TRUE，观测值数量明显多的所占比例就宽，比如盒形图。同一图中，各变量的盒形图本来是等宽的（X轴），长度由数据变异程度决定，如果=TRUE，而几个变量指标间的观测值数量差异太多，很可能有的是较窄的盒形图，有的胖扁的盒形图。这样的图我见过，只是当时没有去深入研究。我的猜测对还是不对，你自己试试吧？OK？
初学者，找到原因真难啊！
areg
######################################################
#
# CHAPTER 2
# 2.3.2 Wave Energy Device Mooring
#
######################################################
# 上面介绍的是独立样本t检验，就是参与检验的样本数据不同，如metres为44个，feet为69个，
> data("waves", package = "HSAUR")
> edit(waves) # 补充by areg
method1 method2
1 2.23 1.82
2 2.55 2.42
3 7.99 8.26
4 4.09 3.46
5 9.62 9.77
6 1.59 1.40
7 8.98 8.88
8 0.82 0.87
9 10.83 11.20
10 1.54 1.33
11 10.75 10.32
12 5.79 5.87
13 5.91 6.44
14 5.79 5.87
15 5.50 5.30
16 9.96 9.82
17 1.92 1.69
18 7.38 7.41
>
################配对t检验 补充by areg##################
> t.test(waves$method1, waves$method2, alternative='two.sided', conf.level=.95, paired=TRUE)
Paired t-test
data: waves$method1 and waves$method2
t = 0.9019, df = 17, p-value = 0.3797 # 对waves的两种方法间没有显著性差异(P=0.05)
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.08258476 0.20591810
sample estimates:
mean of the differences
0.06166667
###############回到手册#############################
> boxplot(waves$method1, ylab="method1")
> mooringdiff <- waves$method1 - waves$method2
> layout(matrix(1:2, ncol = 2))
[1] 2
> boxplot(mooringdiff, ylab = "Differences (Newton metres)",main = "Boxplot")
> abline(h = 0, lty = 2)
> qqnorm(mooringdiff, ylab = "Differences (Newton metres)")
$x # 下面的x和y按本手册在R控制台并不出现，这是在载入Rcmdr包得到的。
[1] 1.08532491 0.06968492 -1.08532491 1.91450583 -0.86163412 0.35549042
[7] -0.06968492 -0.35549042 -1.38299413 0.67448975 1.38299413 -0.67448975
[13] -1.91450583 -0.50848806 0.50848806 0.21042839 0.86163412 -0.21042839
$y
[1] 0.41 0.13 -0.27 0.63 -0.15 0.19 0.10 -0.05 -0.37 0.21 0.43 -0.08
[13] -0.53 -0.08 0.20 0.14 0.23 -0.03
> qqline(mooringdiff) # 图形比原手册精美
areg
######################################################
#
# CHAPTER 2
# 2.3.3 Mortality and Water Hardness
#
######################################################
> source("D:\\Rstudy\\Tiro2.R") # 自制脚本，提供大家以这种方式进行，自己把本章学习要用到的包和数据集调入。
> ls()
[1] "roomwidth" "water" "waves"
>
########## 下面三条命令是Rcmdr包中得到的#############
Hist(water$hardness, scale="frequency", breaks="Sturges", col="darkgray")
scatterplot(mortality~hardness, reg.line=lm, smooth=TRUE, labels=FALSE, boxplots='xy', span =0.5, data=water)
boxplot(water$mortality, ylab="mortality")
###############试试玩，下面回到手册##################
> nf <- layout(matrix(c(2, 0, 1, 3), 2, 2, byrow = TRUE),c(2, 1), c(1, 2), TRUE)
########下面是对上一条命令各组分的分析###############
> matrix(c(2, 0, 1, 3), 2, 2, byrow = TRUE)
[,1] [,2]
[1,] 2 0
[2,] 1 3
##矩阵由正整数组成，byrow = TRUE把C中的元素按行排列，= FALSE，把元素按列秩序排列；矩阵位置中的元素值就是指定所制第某个图形在面板的上摆放位置。如本例中，第[1,2]位置为0，该位置就不放图。
##矩阵后面的第一个向量c(2,1)，表示每列的相对宽度，本例中，第1列是2，第2列是1，意味着第1列宽度是第2列的两倍；矩阵后的第二个向量(1,2)，表示每行的相对或绝对高度，本例中，第1行高度为1，第2行高度为2。这样就把各图的大小和基本位置确定了。
##最后的TRUE表示按矩阵的位置来摆放各个图，如果为FALSE，从右上第1个位置摆放图，向左摆；摆满第一行后，换行，还是从右向左摆。注意啊，这是个细节，别小鞋套大脚！
> psymb <- as.numeric(water$location)
##把数据集中提取字符转化成数字并赋予向量psymb
> plot(mortality ~ hardness, data = water, pch = psymb)
##pch等于单个字符或整数，不同的整数代表不同的符号，如三角、星星、方块、圆圈等等，不过有的值不能用或有专门用途。
> abline(lm(mortality ~ hardness, data = water))
##对散点图加拟合线。
> legend("topright", legend = levels(water$location),pch = c(1, 2), bty = "n")
##在图形的右上区域加上继承水平层次的示意图，在此示意符号用1和2分别代表的圆卷和三解表示，bty有两个取值字母o和n，n是不把示意图区域用框框起来，而o是把该区域框起来。
> hist(water$hardness)
##画水硬度的直方图，就是水硬度各层次出现频率的图，注意与条形图的区别，条形图是各层不相连的。
> boxplot(water$mortality)
##水矿化程度的盒形图
areg
######################################################
#
# CHAPTER 2
# 2.3.4 Piston-ring Failures
#
######################################################
### 注意，下个命令原手册中没提示，你先把数据载入
>data("pistonrings", package = "HSAUR")
> ls()
[1] "pistonrings" "roomwidth" "water" "waves"
> edit(pistonrings) # 查看原数据，了解其列表结构是如何布置的，同时也想想假如这是自己的数据，那该如何分析，再继续看原作者的分析。这是一个很基础性的问题，许多人试验研究数据排列无章，无法分析。任何数据分析，必需对原始数据进行正确的组织结构放置。统计软件的学习，分析操作技能是基础，分析思维才是精华。本例打开数据一看，是一个R*C列联表(4*3)，用拟合优度检验中的卡方检验，数据不是成对的，与成对出现的卡方检验不同。
North Centre South
C1 17 17 12
C2 11 9 13
C3 11 8 19
C4 14 7 28
>
> chisq.test(pistonrings)$residuals
leg
compressor North Centre South
C1 0.6036154 1.6728267 -1.7802243
C2 0.1429031 0.2975200 -0.3471197
C3 -0.3251427 -0.4522620 0.6202463
C4 -0.4157886 -1.4666936 1.4635235
>
> library(grid); library(colorspace); library(vcd) # 手册中载入一个，实际中需要前两个的支持
> assoc(pistonrings) # 画残差列联图
> ?assoc()
############ assoc()帮助文件 ######################
# assoc()
Description
Produce an association plot indicating deviations from a specified independence model in a possibly high-dimensional
contingency table.
Association plots have been suggested by Cohen (1980) and extended by Friendly (1992) and provide a means for
visualizing the residuals of an independence model for a contingency table.
## Default S3 method:
assoc(x, row_vars = NULL, col_vars = NULL, compress = TRUE,
xlim = NULL, ylim = NULL,
spacing = spacing_conditional(sp = 0), spacing_args = list(),
split_vertical = NULL, keep_aspect_ratio = FALSE,
xscale = 0.9, yspace = unit(0.5, "lines"), main = NULL, sub = NULL,
..., residuals_type = "Pearson", gp_axis = gpar(lty = 3))
## S3 method for class 'formula':
assoc(formula, data = NULL, ..., subset = NULL, na.action = NULL, main = NULL, sub = NULL)
################### 残差 #####################
残差是指实际观察值与回归估计值的差。有多少对数据，就有多少个残差。残差分析就是通过残差所提供的信息，分析出数据的可靠性、周期性或其它干扰 。异常数据是指与其它数据产生的条件有明显不同的数据，因此异常数据的残差会特别的大。一旦发现异常数据应及时剔除，用剩余数据重新建立回归方程，以提高回归方程的质量 。
areg
######################################################
#
# CHAPTER 2
# 2.3.5 Rearrests of Juveniles
#
######################################################
> data("rearrests", package = "HSAUR")
> rearrests
Juvenile court
Adult court Rearrest No rearrest
Rearrest 158 515
No rearrest 290 1134
> # 这是一个四格表，并且数据是成对出现，常用卡方检验McNemar's test，函数是mcnemar.test。
> mcnemar.test(rearrests, correct = FALSE)
McNemar's Chi-squared test
data: rearrests
McNemar's chi-squared = 62.8882, df = 1, p-value = 2.188e-15
> ## P值极小，强烈证实在青少年被重捕过的，成年后仍然很可能再犯罪遭重捕。本性难移?少管所造成的?
> binom.test(rearrests[2], n = sum(rearrests[c(2,3)])) # 精确二项式检验
Exact binomial test
data: rearrests[2] and sum(rearrests[c(2, 3)])
number of successes = 290, number of trials = 805, p-value = 1.918e-15
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.3270278 0.3944969
sample estimates:
probability of success
0.3602484
>
## 本章后面，我没有找到手册中的练习题，可能是没有吧！下节介绍第三章的R分析
denver
> library(colorspace)
错误在library(colorspace) : 不存在叫'colorspace'这个名称的程辑包
> library(vcd)
错误在library(vcd) : 不存在叫'vcd'这个名称的程辑包
请问areg老师，我怎么解决此类问题？
areg
[quote]引用第35楼denver于2006-11-20 00:18发表的“”:
> library(colorspace)
错误在library(colorspace) : 不存在叫'colorspace'这个名称的程辑包
> library(vcd)
错误在library(vcd) : 不存在叫'vcd'这个名称的程辑包
请问areg老师，我怎么解决此类问题？[/quote]
你在R控制台窗口“程序包_选择镜像”中选择一个网速快一点，更新也即时的点，我常选的是Korea韩国，选择安装程度包，把需要的这几个包点上，很快就完成的。
areg
#####################################################
#
# CHAPTER 3
# Conditional Inference: Guessing Lengths, Suicides, Gastrointestinal Damage, and Newborn Infants
#
######################################################
3.1 Introduction
3.2 Conditional Test Procedures
3.3 Analysis Using R
#####################################################
#
# CHAPTER 3
# 3.3 Analysis Using R
# 3.3.1 Estimating the Width of a Room Revised
#
######################################################
## 每章学习，首先把相关的包及数据集调用制成个脚本。打开一个脚本，输入下列表达式，并保存为Tiro3.R (只要以.R结尾)
library("lattice"); library("MASS"); library("scatterplot3d"); library("HSAUR")
data("roomwidth", package = "HSAUR")
## 前面的第2章讲到在测量房间宽度时，用米为单位来推断比用英尺大了一此误差，在此我们在一定条件下来分析原因。首先把米为转换成英尺，并把它的观测值赋予为变量y。
> source("D:\\Rstudy\\Tiro3.R")
> ls()
[1] "roomwidth"
>
## 为方便复制到R控制台进行学习，连续条件的命令，我尽量不写提示符。
convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
feet <- roomwidth$unit == "feet"
metre <- !feet
y <- roomwidth$width * convert
T <- mean(y[feet]) - mean(y[metre])
T
## 得到下面的结果
> convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
> feet <- roomwidth$unit == "feet"
> metre <- !feet
> y <- roomwidth$width * convert
> T <- mean(y[feet]) - mean(y[metre])
> T
[1] -8.858893
>
## 只是简单地从平均数来看，两者相差为8.858893英尺。
## 为得到检验T(两种尺度平均数之差)的分布条件，我们用9999个值来赋予y.可以用抽样函数来置换y向量。
meandiffs <- double(9999)
for (i in 1:length(meandiffs))
{
sy <- sample(y)
meandiffs[ i ] <- mean(sy[feet]) - mean(sy[metre])
}
## sample(x, size, replace = FALSE, prob = NULL) 随机置换，## range(y)为此处的x，即(24:131.2), length(y)结果是113，此处为size=113。那么此处的sy<-sample(y),就是每次抽取一个在(24:131.2)范围内，由113个元素组成的样本。
## meandiffs <- double(9999)表示有9999个双精度数值，即它的长度为9999；i in 1:length(meandiffs),即i的取值范围，此为c(1:9999)，循环次数；每循环一次，对sy <- sample(y)进行一次随机置换抽样，进行第i次抽样得到的mean(sy[feet]) - mean(sy[metre])结果赋予meandiffs，共进行了9999次抽样，得到9999个meandiffs. > hist(meandiffs) ## 绘制这9999个meandiffs的直方图。
> abline(v = T, lty = 2) ## 添加线，abline(v=, untf = FALSE, ...)，此处添加的值v=-8.858893;lty线的类型，它的不同数据定义为各种各样的线，如此处为2，代表的是虚线，可以换了试试。下个命令同理，两值组成一个绝对值区域。
> abline(v = -T, lty = 2)
> greater <- abs(meandiffs) > abs(T)
## 把9999个meandiffs中绝对值大于8.858893的元素赋予对象greater。我们可以看看greater的情况，如下
################################
# > summary(greater)
# Mode FALSE TRUE
#logical 9915 84
################################
> mean(greater)
[1] 0.00840084
> binom.test(sum(greater), length(greater))$conf.int ## 提取精确二项检验的临界值
[1] 0.006706214 0.010390397
attr(,"conf.level") ## 提取精确二项检验的置信水平区间
[1] 0.95
############# 附结果 ###########################
> binom.test(sum(greater), length(greater))
Exact binomial test
data: sum(greater) and length(greater)
number of successes = 84, number of trials = 9999, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.006706214 0.010390397
sample estimates:
probability of success
0.00840084
################################################
################################################
## 载入下列包
library("survival")
library("splines")
library("mvtnorm")
library("modeltools")
library("coin")
################################################
> independence_test(y ~ unit, data = roomwidth, distribution = "exact")
Exact General Independence Test
data: y by groups feet, metres
Z = -2.5491, p-value = 0.008492
alternative hypothesis: two.sided
###############################################
> wilcox_test(y ~ unit, data = roomwidth, distribution = "exact")
Exact Wilcoxon Mann-Whitney Rank Sum Test
data: y by groups feet, metres
Z = -2.1981, p-value = 0.02763
alternative hypothesis: true mu is not equal to 0
##############################################
证明这两个尺度用在测量房间时得到的准确程度不相同！与第二章中对房间的所下结论近似。
areg
#####################################################
#
# CHAPTER 3
# 3.3.2 Crowds and Threatened Suicide ## 手册中没有
# 3.3.3 Gastrointestinal Damages
#
#####################################################
> data("Lanza", package = "HSAUR")
> edit(Lanza)
study treatment classification
1 I Misoprostol 1
2 I Misoprostol 1
……………………………## 略
30 I Placebo 1
……………………………## 略
60 II Misoprostol 1
……………………………## 略
90 II Placebo 1
……………………………## 略
120 III Misoprostol 1
……………………………## 略
150 III Placebo 2
……………………………## 略
179 IV Misoprostol 1
……………………………## 略
189 IV Placebo 4
……………………………## 略
198 IV Placebo 5
>
> xtabs(~treatment + classification + study, data = Lanza)
, , study = I
classification
treatment 1 2 3 4 5
Misoprostol 21 2 4 2 0
Placebo 2 2 4 9 13
, , study = II
classification
treatment 1 2 3 4 5
Misoprostol 20 4 6 0 0
Placebo 8 4 9 4 5
, , study = III
classification
treatment 1 2 3 4 5
Misoprostol 20 4 3 1 2
Placebo 0 2 5 5 17
, , study = IV
classification
treatment 1 2 3 4 5
Misoprostol 1 4 5 0 0
Placebo 0 0 0 4 6
>
# xtabs() # 从因子组合来建列联表
# xtabs(formula = ~., data = parent.frame(), subset, na.action, exclude = c(NA, NaN), drop.unused.levels = FALSE)
# 从第一个小组研究结果来看，零假设处理组不是独立的，检验米索前列醇的处理效果。
> library("coin")
> cmh_test(classification ~ treatment, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)), subset = Lanza$study == "I")
Asymptotic Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
chi-squared = 28.8478, df = 1, p-value = 7.83e-08
>
# classification = c(0, 1, 6, 17,30) 没有看明白是怎么来的?
# 极低的P值，显示出强烈的处理效果
# 第二个研究小组
> cmh_test(classification ~ treatment, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)), subset = Lanza$study == "II")
Asymptotic Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
chi-squared = 12.0641, df = 1, p-value = 0.000514
> # 从第二结果的精确P值，得到的置信区间是正确的。
##########
> p <- cmh_test(classification ~ treatment, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)), subset = Lanza$study == "II",distribution = approximate(B = 19999))
> p
Approximative Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
chi-squared = 12.0641, p-value = 5e-05
> pvalue(p)
[1] 5.00025e-05
99 percent confidence interval:
2.506396e-07 3.714653e-04
##########
> cmh_test(classification ~ treatment, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)), subset = Lanza$study == "III")
Asymptotic Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
chi-squared = 28.1587, df = 1, p-value = 1.118e-07
##########
> cmh_test(classification ~ treatment, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)), subset = Lanza$study
== "IV")
Asymptotic Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
chi-squared = 15.7414, df = 1, p-value = 7.262e-05
> # 最后，这样对每一个研究小组的结果进行分析是不满意的。因为试验设计时，把四个小组设为区组变量，应该从全局来研究它们线性关系。
> cmh_test(classification ~ treatment | study, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)))
Asymptotic Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
stratified by study
chi-squared = 83.6188, df = 1, p-value < 2.2e-16
> # 假如不考虑各个小组，直接用呢？补充by areg
#####################################################
>cmh_test(classification ~ treatment, data = Lanza,scores = list(classification = c(0, 1, 6, 17,30)))
Asymptotic Linear-by-Linear Association Test
data: classification (ordered) by groups Misoprostol, Placebo
chi-squared = 78.6296, df = 1, p-value < 2.2e-16
> #虽然在此结果差不多，实际算法还是有区别的，在分析中要注意。
#####################################################
areg
#####################################################
#
# CHAPTER 3
# 3.3.4 Teratogenesis
#
#####################################################
# 矩阵构建时，默认是按列把数据读入。
> anomalies <- as.table(matrix(c(235, 23, 3, 0, 41,35, 8, 0, 20, 11, 11, 1, 2, 1, 3, 1), ncol = 4,dimnames = list(MD = 0:3,
RA = 0:3)))
> anomalies
RA
MD 0 1 2 3
0 235 41 20 2
1 23 35 11 1
2 3 8 11 3
3 0 0 1 1
>
###############################################
> mh_test(anomalies)
Asymptotic Marginal-Homogeneity Test
data: response by groups MD, RA
stratified by block
chi-squared = 21.2266, df = 3, p-value = 9.446e-05
>
###############################################
> mh_test(anomalies, scores = list(c(0, 1, 2, 3)))
Asymptotic Marginal-Homogeneity Test for Ordered Data
data: response (ordered) by groups MD, RA
stratified by block
chi-squared = 21.0199, df = 1, p-value = 4.545e-06
>
###############################################
## 第四章 会晚点一些再帧出，我要先分析一下自己的试验数据。大家先学，如不明白的问题，提出来，共同学习探讨
#####################################################
#
# CHAPTER 4
# Analysis of Variance: Weight Gain,Foster Feeding in Rats, WaterHardness and Male Egyptian Skulls
#
#####################################################
