1. 程式人生 > >R語言:常用統計檢驗方法

R語言:常用統計檢驗方法

轉自http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=240107

正態總體均值的假設檢驗

t檢驗
單個總體
例一 
某種元件的壽命X(小時),服從正態分佈,N(mu,sigma^2),其中mu,sigma^2均未知,16只元件的壽命如下:問是否有理由認為元件的平均壽命大於255小時。
命令:
X<-c(159, 280, 101, 212, 224, 379, 179, 264,
222, 362, 168, 250, 149, 260, 485, 170)
t.test(X, alternative = "greater", mu = 225)

兩個總體

例二
X為舊鍊鋼爐出爐率,Y為新鍊鋼爐出爐率,問新的操作能否提高出爐率
命令:
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)

t.test(X, Y, var.equal=TRUE, alternative = "less")

成對資料t檢驗

例三
對每個高爐進行配對t檢驗
命令:
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
t.test(X-Y, alternative = "less")

正態總體方差的假設檢驗

例四
從小學5年級男生中抽取20名,測量其身高(釐米)如下:
問,在0.05顯著性水平下,
平均值是否等於149
sigma^2 是否等於 75
命令:
X<-scan()
136 144 143 157 137 159 135 158 147 165
158 142 159 150 156 152 140 149 148 155
var.test(X,Y)

例五

對鍊鋼爐的資料進行分析
命令:
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
var.test(X,Y)

二項分佈的總體檢驗

例六

有一批蔬菜種子的平均發芽率為P=0.85,現在隨機抽取500粒,用種衣劑進行浸種處理,結果有445粒發芽,問種衣劑有無效果。
命令:
binom.test(445,500,p=0.85)

例七

按照以往經驗,新生兒染色體異常率一般為1%,某醫院觀察了當地400名新生兒,有一例染色體異常,問該地區新生兒染色體是否低於一般水平?
命令:
binom.test(1,400,p=0.01,alternative="less")
非引數檢驗
#資料是否正態分佈的Neyman-Pearson 擬合優度檢驗-chisq

例八 

5種品牌啤酒愛好者的人數如下
A 210
B 312
C 170
D 85 
E 223
問不同品牌啤酒愛好者人數之間有沒有差異?
命令:
X<-c(210, 312, 170, 85, 223)
chisq.test(X)

例九 

檢驗學生成績是否符合正態分佈
命令:
X<-scan()
25 45 50 54 55 61 64 68 72 75 75
78 79 81 83 84 84 84 85 86 86 86
87 89 89 89 90 91 91 92 100
A<-table(cut(X, br=c(0,69,79,89,100)))
p<-pnorm(c(70,80,90,100), mean(X), sd(X))
p<-c(p[1], p[2]-p[1], p[3]-p[2], 1-p[3])
chisq.test(A,p=p)
# cut 將變數區域劃分為若干區間
# table 計算因子合併後的個數
# 均值之間有無顯著區別
大麥的雜交後代芒性狀的比例 無芒:長芒: 短芒=9:3:4,而實際觀測值為335:125:160 ,檢驗觀測值是否符合理論假設?
命令:
chisq.test(c(335, 125, 160), p=c(9,3,4)/16)

例十

# 現有42個數據,分別表示某一時間段內電話總機借到呼叫的次數,
# 接到呼叫的次數 0   1   2   3   4   5   6
# 出現的頻率     7   10  12  8   3   2   0
# 問:某個時間段內接到的呼叫次數是否符合Possion分佈?
命令:
x<-0:6
y<-c(7,10,12,8,3,2,0)
mean<-mean(rep(x,y))
q<-ppois(x,mean)
n<-length(y)
p[1]<-q[1]
p[n]<-1-q[n-1]
for(i in 2:(n-1))
p<-q-q
chisq.test(y, p=p)
Z<-c(7, 10, 12, 8)
n<-length(Z); p<-p[1:n-1]; p[n]<-1-q[n-1]
chisq.test(Z, p=p)


內容來自 

薛毅 陳立萍 《統計建模與R軟體》 清華大學出版社 2006

-----------------------------------

轉自http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=240107

# 理論分佈依賴於若干未知引數時
# Kolmogorov-Smirnov 檢驗
# ks.test() 
例一
對一臺裝置進行壽命檢驗,記錄十次無故障操作時間,並按從小到大的次序排列如下,
# 用ks檢驗方法檢驗此裝置無故障工作時間是否符合rambda=1/1500的指數分佈
# 命令:
X<-c(420, 500, 920, 1380, 1510, 1650, 1760, 2100, 2300, 2350)
ks.test(X, "pexp", 1/1500)
 
例二
假設從分佈函式F(x)和G(x)的總體中分別隨機抽取25個和20個觀察值樣本,檢驗F(x)和G(x)是否相同。
# 命令
X<-scan()
0.61 0.29 0.06 0.59 -1.73 -0.74 0.51 -0.56 0.39
1.64 0.05 -0.06 0.64 -0.82 0.37 1.77 1.09 -1.28
2.36 1.31 1.05 -0.32 -0.40 1.06 -2.47
Y<-scan()
2.20 1.66 1.38 0.20 0.36 0.00 0.96 1.56 0.44
1.50 -0.30 0.66 2.31 3.29 -0.27 -0.37 0.38 0.70
0.52 -0.71
ks.test(X, Y)
# ks多樣本檢驗的侷限性,只用在理論分佈為一維連續分佈,且分佈完全已知的情形。ks檢驗可用的情況下,功效一般優於Pearson chisq檢驗
 
# 列聯表(contingerncy table)的獨立性檢驗
# Pearson chisquare 進行獨立性檢驗
例三
為了研究吸菸是否與肺癌有關,對63位患者及43名非肺癌患者調查了其中的吸菸人數,得到2*2列聯表
# 資料     肺癌     健康      合計
# 吸菸     60        32          92
# 不吸菸    3        11          14
# 合計     63        43         106
# 命令
x<-c(60, 3, 32, 11)
dim(x)<-c(2,2)
chisq.test(x,correct = FALSE) # 不帶連續校正的情況
chisq.test(x) # 帶連續校正的情況
 
例四
在一次社會調查中,以問卷方式調查了901人的年收入,及其對工作的滿意程度,其中年收入A分為四檔:小於6000元,6000-15000元,15000元至25000元,超過25000元。對工作的滿意程度B 分為 很不滿意,較不滿意,基本滿意和很滿意四檔,結果如下
#                     很不滿意  較不滿意   基本滿意   很滿意        合計
# < 6000               20        24         80           82            206
# 6000 ~15000          22        38         104        125           289
# 15000~25000          13        28         81          113           235
# > 25000               7          18         54           92            171
# 合計                 62        108       319        412           901
# 命令如下
x<-scan()
20 24 80 82 22 38 104 125
13 28 81 113 7 18 54 92
dim(x)<-c(4,4)
chisq.test(x)

# Fisher 精確的獨立檢驗
# 試用條件 樣本數小於4
例五
某醫師研究乙肝免疫球蛋白防止子宮內胎兒感染HBV的效果,將33例HBsAg陽性孕婦隨機分為預防注射組和對照組,結果由下表所示,兩組新生兒HBV總體感染率有無差別
# 組別           陽性    陰性    合計   感染率
# 預防注射組     4       18      22     18.8
# 對照組         5       6       11     45.5
# 命令如下
x<-c(4,5,18,6); dim(x)<-c(2,2)
fisher.test(x)
# 對前面提到的肺癌進行檢驗
x<-c(60, 3, 32, 11);
dim(x)<-c(2,2)
fisher.test(x)
 
# McNemar檢驗
# McNemar檢驗不是獨立性檢驗,但是是關於列連表的檢驗
例六
甲乙兩種方法檢測細菌的結果
#          乙方法          
#                            合計
# 甲方法   +        -            
# +        49      25         74
# -        21      107        128
# 合計     70      132        202
# 命令
X <- c(49, 21, 25, 107);
dim(X) <- c(2,2)
mcnemar.test(X,correct=FALSE)
 
# 符號檢驗
# 1 假設一個樣本是否來自某個總體
例七
聯合國人員在世界上66個大城市的生活花費指數(以紐約1996年12月為100),按照從小到大的次序排列如下,其中北京的指數為99。假設這個樣本是從世界大城市中隨機抽樣得到的。用符號檢驗分析,北京是在中位數之上,還是中位數之下。
X<-scan()
66 75 78 80 81 81 82 83 83 83 83
84 85 85 86 86 86 86 87 87 88 88
88 88 88 89 89 89 89 90 90 91 91
91 91 92 93 93 96 96 96 97 99 100
101 102 103 103 104 104 104 105 106 109 109
110 110 110 111 113 115 116 117 118 155 192
binom.test(sum(X>99), length(X), al="l")

# 2 用成對樣本檢驗兩總體間是否有差異
例八
兩種不同飼料,對豬增重情況如下,分析兩種飼料養豬有無差異
# 命令
x<-scan()
25 30 28 23 27 35 30 28 32 29 30 30 31 16
y<-scan()
19 32 21 19 25 31 31 26 30 25 28 31 25 25
binom.test(sum(x<y), length(x))
 
例九
某飲料店為調查了顧客對飲料的愛好情況,某日隨機調查了13為顧客,喜歡奶茶超過咖啡用-表示,喜歡咖啡超過奶茶用+表示,兩者都喜歡用0表示,結果如下,分析顧客是更喜歡咖啡開始奶茶。
# 顧客編號 1 2 3 4 5 6 7 8 9 10 11 12 13
# 喜歡咖啡 1   1 1 1 0 1   1  1  1    1
# 喜歡奶茶   1           1          1
binom.test(3,12,p=1/2, al="l", conf.level = 0.90)

# 秩統計量 
# Spearman秩相關檢驗
例十
一項有六人蔘加表演的競賽,有兩人進行評定,評定結果用如表所示,試用Spearman秩相關檢驗方法檢驗這兩個評定員對於等級評定有無相關性
# 選手編號 1 2 3 4 5 6
# 甲的打分 4 2 2 4 5 6
# 乙的打分 5 3 4 3 2 5
x<-c(4,2,2,4,5,6); y<-c(5,3,4,3,2,5)
cor.test(x, y, method = "spearman")
 
# Kendall相關檢驗
例十一
某幼兒園對9對雙胞胎的智力進行測驗,並按照百分制打分,試用Kendall相關檢驗方法檢驗雙胞胎的智力是否相關。
# 1   2  3  4  5  6  7  8   9
# 86  77 68 91 70 71 85 87 63
# 88  76 64 96 65 80 81 72 60
X<-c(86, 77, 68, 91, 70, 71, 85, 87, 63)
Y<-c(88, 76, 64, 96, 65, 80, 81, 72, 60)
cor.test(X, Y, method = "kendall")

# Wilcoxon秩檢驗—— 考慮了樣本觀察值月總體中位數的差。
# 1 對於來自同一個總體樣本的檢驗
例十二
某電池廠稱其生產的某種電池,中位數為140安培小時,現隨機從其新生產的電池中抽取20個,檢驗其壽命,137.0 140.0 138.3 139.0 144.3 139.1 141.7 137.3 133.5 138.2 141.1 139.2 136.5 136.5 135.6 138.0 140.9 140.6 136.3 134.1
# 用Wilcoxon符號檢驗分析該廠生產的電池是否符合標準
X<-scan()
137.0 140.0 138.3 139.0 144.3 139.1 141.7 137.3 133.5 138.2
141.1 139.2 136.5 136.5 135.6 138.0 140.9 140.6 136.3 134.1
wilcox.test(X, mu=140, alternative="less",
exact=FALSE, correct=FALSE, conf.int=TRUE)

# 該方法也可用於成對樣本的檢驗 
例十三
為檢驗某種新肥料,將現有麥地分為十塊,再將每一塊分為兩部分,一半施普通肥料,一半兒施新肥料,用Wilcoxon符號檢驗法檢驗新複合肥能否顯著提高小麥產量。
# 1   2   3   4   5   6   7   8   9   10
# 459 367 303 392 310 342 421 446 430 412
# 414 306 321 443 281 301 353 391 405 390
#   
x<-c(459, 367, 303, 392, 310, 342, 421, 446, 430, 412)
y<-c(414, 306, 321, 443, 281, 301, 353, 391, 405, 390)
wilcox.test(x, y, alternative = "greater", paired = TRUE)
wilcox.test(x-y, alternative = "greater")
binom.test(sum(x>y), length(x), alternative = "greater")
 
# 非成對樣本的秩次和檢驗
# Wilcoxon-Mann-Whitney 統計量 U
例十四
測量了10名不同作業組的工人血鉛含量,分析兩組之間是否有差別。
# 非鉛作業組 24 26 29 34 43 58 63 72 87 101
# 含鉛作業組 82 87 97 121 164 208 213
x<-c(24, 26, 29, 34, 43, 58, 63, 72, 87, 101)
y<-c(82, 87, 97, 121, 164, 208, 213)
wilcox.test(x,y,alternative="less",exact=FALSE,correct=FALSE)
wilcox.test(x, y, alternative="less", exact=FALSE)

例十五
學生數學能力排序
新方法 3 5 7 9 10
原方法 1 2 4 6 8
問新舊方法之間是否有差別。
x<-c(3, 5, 7, 9, 10)
y<-c(1, 2, 4, 6, 8)
wilcox.test(x, y, alternative="greater")

例十六
檢驗一種藥物對於慢性支氣管炎有沒有效果,抽取了216個病例,治療效果。分析該藥物對兩種慢性支氣管炎的治療效果是否相同。
#       控制 顯效 進步 無效
# 單純型 62   41   14   11
# 喘息型 20   37   16   15
x<-rep(1:4, c(62, 41, 14,11))
y<-rep(1:4, c(20, 37, 16, 15))
wilcox.test(x, y, exact=FALSE)