forked from fanofmeasurement/Rlectures
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathR_meandiff_1.Rmd
301 lines (202 loc) · 8.41 KB
/
R_meandiff_1.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
---
title: "R的假设检验之一"
author: "李峰"
date: '`r strftime(Sys.time(), format = "%B %d, %Y")`'
output:
html_document: default
---
```{r,echo=F}
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
```
```{r}
library(sm)
```
---
#### 1. 独立样本的均值差异检验
##### 1.1 理论依据
* 若总体方差$\sigma_1^2$和$\sigma_2^2$未知，但由经验可知相等，则$\sigma^2_{x_1-x_2}$的理论估计为：
\begin{equation}\label{eq:1}
\sigma^2_{x_1-x_2} = \frac{S_p}{n_1}+ \frac{S_p}{n_2},
\end{equation}
其中，
\begin{equation}\label{eq:2}
S_p = \frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2},
\end{equation}
称为合并方差。
```{r}
par(mfrow=c(2,1),mar=c(4,4,4,4))
set.seed(12345)
Pop1<-rnorm(10000,mean=2,sd=2) ###两总体方差相等
Pop2<-rnorm(10000,mean=10,sd=2)
Diff<-vector()
Sdx1<-vector()
Sdx2<-vector()
for(i in 1:2000){
x1<-sample(Pop1,size=100,replace=TRUE)
x2<-sample(Pop2,size=120,replace=TRUE)
Diff<-c(Diff,(mean(x1)-mean(x2)))
Sdx1<-c(Sdx1,sd(x1))
Sdx2<-c(Sdx2,sd(x2))
}
plot(density(Diff),xlab="mean(x1)-mean(x2)",ylab="Density",main="均值差的抽样分布(等方差)",cex.main=0.7,cex.lab=0.7)
points(mean(Diff),sd(Diff),pch=1,col=1)
S1<-mean(Sdx1)
S2<-mean(Sdx2)
Sp<-((100-1)*S1^2+(120-1)*S2^2)/(100+120-2)
points((2-10),sqrt(Sp/100+Sp/120),pch=2,col=2)
```
红色圆圈为估计值，黑色三角为理论值，二者重合。
---
* 若总体方差$\sigma_1^2$和$\sigma_2^2$未知且不相等时，则$\sigma^2_{x_1-x_2}$的理论估计为：
\begin{equation}\label{eq:3}
\sigma^2_{x_1-x_2} = \frac{S_1}{n_1}+ \frac{S_2}{n_2},
\end{equation}
```{r}
set.seed(12345)
Pop1<-rnorm(10000,mean=2,sd=2) ###两总体方差不等
Pop2<-rnorm(10000,mean=10,sd=4)
Diff<-vector()
Sdx1<-vector()
Sdx2<-vector()
for(i in 1:2000){
x1<-sample(Pop1,size=100,replace=TRUE)
x2<-sample(Pop2,size=120,replace=TRUE)
Diff<-c(Diff,(mean(x1)-mean(x2)))
Sdx1<-c(Sdx1,sd(x1))
Sdx2<-c(Sdx2,sd(x2))
}
plot(density(Diff),xlab="mean(x1)-mean(x2)",ylab="Density",main="均值差的抽样分布(不等方差)",cex.main=0.7,cex.lab=0.7)
points(mean(Diff),sd(Diff),pch=1,col=1)
S1<-mean(Sdx1)
S2<-mean(Sdx2)
points((2-10),sqrt(S1^2/100+S2^2/120),pch=2,col=2)
```
红色圆圈为估计值，黑色三角为理论值，二者重合。
方差不等时，Welch提出仍采用$\frac{\bar{x}_1-\bar{x}_2}{\sigma_{x_1-x_2}}$作为统计量，但并不服从自由度为$n_1+n_2-2$的t分布，而是服从自由度为：
\begin{equation}\label{eq:4}
df = \frac{(\frac{S^2_1}{n_1}+\frac{S^2_2}{n_2})^2}{\frac{(\frac{S^2_1}{n_1})^2}{n_1-1}+\frac{(\frac{S^2_2}{n_2})^2}{n_2-1}},
\end{equation}
的t分布，再依据t分布进行统计判断，即Welch调整。
---
##### 1.2 R函数
> t.test(数值型变量名~因子，data=数据框名，paired=FALSE，var.equal=TRUE/FALSE，mu=检验值，alternative=检验方向)
```{r}
Forest<-read.table(file="ForestData.txt",header=TRUE,sep=" ")
Forest$month<-factor(Forest$month,levels=c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec"))
Tmp<-subset(Forest,Forest$month=="jan" | Forest$month=="aug")
t.test(temp~month,data=Tmp,paired=FALSE,var.equal=TRUE)
t.test(temp~month,data=Tmp,paired=FALSE,var.equal=FALSE)
```
对气温数据的分析显示，一月份和八月份森林的气温有显著差异，t(df=`r AA$parameter`)=`r AA$statistic`，一月份的均值是`r AA$estimate[1]`(SD=)，八月份的气温均值是
应采用以上哪个结果应禁言是否方差齐性，采用levene's方差同质检验。levene's是对离差的绝对值的均值差异进行判断，如果有显著差异即方差不齐，否则方差齐性。
```{r}
library("car")
leveneTest(Tmp$temp,Tmp$month, center=mean)
```
---
以上可以形成一个自动化的部分。
```{r}
varequ<-leveneTest(Tmp$temp,Tmp$month, center=mean)
varequp<-varequ$`Pr(>F)`[1]
if (varequp>=.05) t.test(temp~month,data=Tmp,paired=FALSE,var.equal=TRUE) else{
t.test(temp~month,data=Tmp,paired=FALSE,var.equal=FALSE)
}
```
#### 2. 配对样本的均值差异检验
配对样本的均值差异检验本质上是单样本的t检验问题。即，可以将一一对应的观测值求差，得到差值样本，若差值样本的均值和0有显著差异，则认为配对样本的两总体的均值差有显著差异，反之则不显著。
即，对于单样本的t检验：
\begin{equation}\label{eq:5}
t=\frac{\bar{x}}{S/\sqrt{n}},
\end{equation}
\begin{equation}\label{eq:6}
t=\frac{\bar{D}}{S_D/\sqrt{n}},
\end{equation}
可以先得到差值样本，然后用将其和0比较：
```{r}
Diff<-ReportCard$chi-ReportCard$math
t.test(Diff,mu=0)
```
也可以直接使用函数：
```{r}
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
ReportCard<-na.omit(ReportCard)
t.test(ReportCard$chi,ReportCard$math,paired=TRUE)
```
#### 3. 关于效应值
一类错误为弃真错误，即$\alpha$水平，二类错误为取伪错误，即$\beta$。$1-\beta$为正确否定原假设的概率，即虚无假设是错的，而且你正确否定了。
$\beta$。$1-\beta$的大小取决于几个因素：
1. $\alpha$越大，则$1-\beta$越大；
---
2. 均值差异的大小，若均值差异很大，则两个分布间重叠程度低，$1-\beta$会增大；
3. 样本大小也是很重要的因素，样本越大，方差减小，则$1-\beta$会增大。
4. 单侧还是双侧检验也有影响，因双侧检验时，是用$p/2$和$\alpha$进行比较。
---
给定$\alpha$和$n$，并且规定了检验方向时，可以通过判断均值差异的大小判断两个分部间的重叠程度，判定其统计功效，即$1-\beta$的大小。
在样本均值差异的功效分析中，R采用Cohen提出的效应量定义。
对独立样本t检验，其效应量为：
\begin{equation}\label{eq:7}
ES=\frac{|\bar{x}_1-\bar{x}_2|}{\sqrt{\frac{S^2_1+S^2_2}{2}}},
\end{equation}
配对样本均值t检验的效应量为：
\begin{equation}\label{eq:8}
ES=\frac{|\bar{D}|}{S_{\bar{D}}},
\end{equation}
单样本t检验的效应量为：
\begin{equation}\label{eq:9}
ES=\frac{|\bar{x}-\mu_0|}{S},
\end{equation}
根据Cohen的研究，ES的值0.2时效应量比较小，取0.5时中等，取0.8时比较大。
在R里，用*pwr*包来实现统计功效的计算。样本量相同时，用**pwr.t.test**，样本量不相同时，用**pwr.t2n.test**。
```{r}
# install.packages("pwr")
library("pwr")
pwr.t2n.test(n1=2,n2=184,d=4.8,sig.level=0.05,alternative="two.sided")
```
也可以让上述程序自动化
```{r}
Des.Fun<-function(x,...){
Av<-mean(x,na.rm=TRUE)
Var<-var(x,na.rm=TRUE)
N<-length(x[!is.na(x)])
result<-list(avg=Av,var=Var,n=N)
return(result)
}
temp.month<-tapply(Tmp$temp,INDEX=Tmp$month,FUN=Des.Fun,na.rm=TRUE)
ES<-abs(temp.month$jan$avg-temp.month$aug$avg)/(sqrt((temp.month$jan$var+temp.month$aug$var)/2))
pwr.t2n.test(n1=temp.month$jan$n,n2=temp.month$aug$n,d=ES,sig.level=0.05,alternative="two.sided")
```
对于学生成绩，前面配对样本t检验的结果是语文和数学成绩存在显著差异，如果这是一个正确的决策，那么在$\alpha=0.05$的情况下，统计功效达到0.8且n=58时，ES至少应多大？
```{r}
pwr.t.test(n=58,sig.level=0.05,power=0.8,type="paired",alternative="two.sided")
```
试着让上述程序自动化？判断在例子中，正确拒绝虚无假设的概率是否大于0.8？
* 相关系数的统计效应分析
```{r}
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
Tmp<-ReportCard[complete.cases(ReportCard),]
corr<-cor.test(Tmp[,5],Tmp[,7],alternative="two.side",method="pearson")
library("pwr")
pwr.r.test(r=corr$estimate,sig.level=0.05,n=(corr$parameter+2),alternative="two.sided")
```
* chi-square的功效分析
其效应量是：
\begin{equation}\label{eq:10}
ES=\sqrt{\sum_{i=1}^{rc}\frac{(p^o_i-p^e_i)^2}{p^e_i}}=\sqrt{\frac{\chi^2}{n}},
\end{equation}
即$\phi$系数。
```{r}
(CrossTable<-table(Tmp[,c(2,12)]))
(ResChisq<-chisq.test(CrossTable,correct=FALSE))
(phi<-sqrt(ResChisq$statistic/nrow(Tmp)))
library("pwr")
(chip<-pwr.chisq.test(sig.level=0.05,N=nrow(Tmp),power=0.9,df=ResChisq$parameter))
if (phi<chip$w) print("power小于0.9") else{
print("power大于0.9")
}
```