areg
#####################################################
#
# 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[ i ]，共进行了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
