本文主要介紹以下內容：1.基本統計方法 2頻數表和列聯表 3.相關 4.t檢驗 5.組間差異的非參數檢驗
1.基本統計方法
簡單實例 mycars<-c("mpg","hp","wt") summary(mtcars[mycars]) 結果： mpg hp wt Min. :10.40 Min. : 52.0 Min. :1.513 1st Qu.:15.43 1st Qu.: 96.5 1st Qu.:2.581 Median :19.20 Median :123.0 Median :3.325 Mean :20.09 Mean :146.7 Mean :3.217 3rd Qu.:22.80 3rd Qu.:180.0 3rd Qu.:3.610 Max. :33.90 Max. :335.0 Max. :5.424
summary()函數提供了最小值、最大值、四分位數、均值，另外還可以因子向量和邏輯型向量的頻數統計。另有apply()函數和sapply()函數計算所選擇的任意描述性統計量。對於sapply(x)函數，使用格式如下：
sapply(x,FUN,options)
其中x是數據框，FUN作為一個任意函數。如果指定了options,他們將被傳遞給FUN。可以插入的函數有mean()、sd()、var()、min()、max()、median()、length()、range()、quantile()。函數fivenum()可返回圖基五數總括(最小值，下四分位數、中位數、上四分位數、最大值)
mystats <- function(x, na.omit=FALSE){ if (na.omit) x <- x[!is.na(x)] m <- mean(x) n <- length(x) s <- sd(x) skew <- sum((x-m)^3/s^3)/n kurt <- sum((x-m)^4/s^4)/n - 3
return(c(n=n, mean=m, stdev=s, skew=skew, kurtosis=kurt)) } myvars<-c("mpg","hp","wt") sapply(mtcars[myvars],mystats) #輸出值
mpg hp wt n 32.000000 32.0000000 32.00000000 mean 20.090625 146.6875000 3.21725000 stdev 6.026948 68.5628685 0.97845744 skew 0.610655 0.7260237 0.42314646 kurtosis -0.372766 -0.1355511 -0.02271075
其他描述性統計函數，包含Hmisc、pastecs、psych。
library(Hmisc) myvars<-c("mpg","hp","wt") describe(mtcars[myvars]) #輸出結果
3 Variables 32 Observations -------------------------------------------------------------------- mpg n missing distinct Info Mean Gmd .05
32 0 25 0.999 20.09 6.796 12.00 .10 .25 .50 .75 .90 .95
14.34 15.43 19.20 22.80 30.09 31.30 lowest : 10.4 13.3 14.3 14.7 15.0, highest: 26.0 27.3 30.4 32.4 33.9
-------------------------------------------------------------------- hp n missing distinct Info Mean Gmd .05
32 0 22 0.997 146.7 77.04 63.65 .10 .25 .50 .75 .90 .95
66.00 96.50 123.00 180.00 243.50 253.55 lowest : 52 62 65 66 91, highest: 215 230 245 264 335
-------------------------------------------------------------------- wt n missing distinct Info Mean Gmd .05
32 0 29 0.999 3.217 1.089 1.736 .10 .25 .50 .75 .90 .95
1.956 2.581 3.325 3.610 4.048 5.293 lowest : 1.513 1.615 1.835 1.935 2.140, highest: 3.845 4.070 5.250 5.345 5.424
--------------------------------------------------------------------
library(pastecs) myvars<-c("mpg","hp","wt") stat.desc(mtcars[myvars]) #輸出結果
mpg hp wt nbr.val 32.0000000 32.0000000 32.0000000 nbr.null 0.0000000 0.0000000 0.0000000 nbr.na 0.0000000 0.0000000 0.0000000 min 10.4000000 52.0000000 1.5130000 max 33.9000000 335.0000000 5.4240000 range 23.5000000 283.0000000 3.9110000 sum 642.9000000 4694.0000000 102.9520000 median 19.2000000 123.0000000 3.3250000 mean 20.0906250 146.6875000 3.2172500 SE.mean 1.0654240 12.1203173 0.1729685 CI.mean.0.95 2.1729465 24.7195501 0.3527715 var 36.3241028 4700.8669355 0.9573790 std.dev 6.0269481 68.5628685 0.9784574 coef.var 0.2999881 0.4674077 0.3041285
library(psych) myvars<-c("mpg","hp","wt") describe(mtcars[myvars]) #輸出結果
vars n mean sd median trimmed mad min max range mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40 33.90 23.50 hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00 335.00 283.00 wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51 5.42 3.91 skew kurtosis se mpg 0.61 -0.37 1.07 hp 0.73 -0.14 12.12 wt 0.42 -0.02 0.17
分組描述性統計量：在比較多組個體或觀測時，關注的焦點經常是各組的描述性統計信息，而不是樣本整體的描述性統計信息，
myvars<-c("mpg","hp","wt") aggregate(mtcars[myvars],by=list(am=mtcars$am),mean) #輸出結果 am mpg hp wt 1 0 17.14737 160.2632 3.768895 2 1 24.39231 126.8462 2.411000
#aggregate()僅允許在每次調用中使用平均值、標准差這樣的單返回值函數
#by(data,INDICES,FUN)
dstats<-function(x){sapply(x,mystats)
myvars<-c("mpg","hp","wt")
by(mtcars[myvars],mtcars$am,dstats)
mtcars$am: 0
mpg hp wt
n 19.00000000 19.00000000 19.0000000
mean 17.14736842 160.26315789 3.7688947
stdev 3.83396639 53.90819573 0.7774001
skew 0.01395038 -0.01422519 0.9759294
kurtosis -0.80317826 -1.20969733 0.1415676
--------------------------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.00000000 13.0000000 13.0000000
mean 24.39230769 126.8461538 2.4110000
stdev 6.16650381 84.0623243 0.6169816
skew 0.05256118 1.3598859 0.2103128
kurtosis -1.45535200 0.5634635 -1.1737358
分組計算的擴展，doBy包和psych包提供了分組計算的描述性統計量的函數，doBy包中的summaryBy()函數使用的基本格式：
#doBy()包中summaryBy()函數的使用格式： #summaryBy(formula,data=dataframe,FUN=function) #formula接受以下格式： #var1+var2+var3+var4+……+varN~groupvar1+groupvar2+……#+groupvarN #在~左側的變量師需要分析的數值型變量，而在右側的變量是類別型的分組變#量。
library(doBy) summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats) #輸出結果
am mpg.n mpg.mean mpg.stdev mpg.skew mpg.kurtosis hp.n hp.mean hp.stdev hp.skew 1 0 19 17.14737 3.833966 0.01395038 -0.8031783 19 160.2632 53.90820 -0.01422519
2 1 13 24.39231 6.166504 0.05256118 -1.4553520 13 126.8462 84.06232 1.35988586 hp.kurtosis wt.n wt.mean wt.stdev wt.skew wt.kurtosis 1 -1.2096973 19 3.768895 0.7774001 0.9759294 0.1415676
2 0.5634635 13 2.411000 0.6169816 0.2103128 -1.1737358
library(psych) myvars<-c("mpg","hp","wt") describeBy(mtcars[myvars],list(am=mtcars$am)) #輸出結果
Descriptive statistics by group am: 0 vars n mean sd median trimmed mad min max range skew kurtosis se mpg 1 19 17.15 3.83 17.30 17.12 3.11 10.40 24.40 14.00 0.01 -0.80 0.88 hp 2 19 160.26 53.91 175.00 161.06 77.10 62.00 245.00 183.00 -0.01 -1.21 12.37 wt 3 19 3.77 0.78 3.52 3.75 0.45 2.46 5.42 2.96 0.98 0.14 0.18
-------------------------------------------------------------------- am: 1 vars n mean sd median trimmed mad min max range skew kurtosis se mpg 1 13 24.39 6.17 22.80 24.38 6.67 15.00 33.90 18.90 0.05 -1.46 1.71 hp 2 13 126.85 84.06 109.00 114.73 63.75 52.00 335.00 283.00 1.36 0.56 23.31 wt 3 13 2.41 0.62 2.32 2.39 0.68 1.51 3.57 2.06 0.21 -1.17 0.17
2.頻數表和列聯表
#一維列聯表
#table()函數生成簡單的頻數統計表
library(vcd) mytable<-with(Arthritis,table(Improved)) mytable #輸出結果
Improved None Some Marked 42 14 28
#prop.table()將頻數轉化成比例值：
prop.table(mytable) Improved None Some Marked 0.5000000 0.1666667 0.3333333
#prop.table()*100轉化成百分比
prop.table(mytable)*100 Improved None Some Marked 50.00000 16.66667 33.33333
#二維列聯表 #table()函數的使用格式 #mytable<-table(A,B) A為行變量 B為列變量 #xtabs()函數還可以使用公式風格的輸入創建列聯表 #mytable<-xtabs(~A+B,data=mydata) #其中mydata是一個矩陣或者數據框
library(vcd) mytable<-xtabs(~Treatment+Improved,data=Arthritis) mytable #輸出結果
Improved Treatment None Some Marked Placebo 29 7 7 Treated 13 7 21
#margin.table()和prop.table()函數生成頻數和比例
margin.table(mytable,1) #輸出結果
Treatment Placebo Treated 43 41 prop.table(mytable,1) #輸出結果
Improved Treatment None Some Marked Placebo 0.6744186 0.1627907 0.1627907 Treated 0.3170732 0.1707317 0.5121951
#1代表table()語句中的第一個變量 #列和列比例的計算
margin.table(mytable,2) Improved None Some Marked 42 14 28 prop.table(mytable,2) Improved Treatment None Some Marked Placebo 0.6904762 0.5000000 0.2500000 Treated 0.3095238 0.5000000 0.7500000
#addmargins()函數為表增加邊際和
addmargins(mytable) Improved Treatment None Some Marked Sum Placebo 29 7 7 43 Treated 13 7 21 41 Sum 42 14 28 84 addmargins(prop.table(mytable)) Improved Treatment None Some Marked Sum Placebo 0.34523810 0.08333333 0.08333333 0.51190476 Treated 0.15476190 0.08333333 0.25000000 0.48809524 Sum 0.50000000 0.16666667 0.33333333 1.00000000 addmargins(prop.table(mytable,1),2) #只增加了行的和
Improved Treatment None Some Marked Sum Placebo 0.6744186 0.1627907 0.1627907 1.0000000 Treated 0.3170732 0.1707317 0.5121951 1.0000000 addmargins(prop.table(mytable,2),1) #只增加列的和
Improved Treatment None Some Marked Placebo 0.6904762 0.5000000 0.2500000 Treated 0.3095238 0.5000000 0.7500000 Sum 1.0000000 1.0000000 1.0000000
#使用CrossTable()創建二維列表
library(gmodels) CrossTable(Arthritis$Treatment,Arthritis$Improved) Cell Contents |-------------------------|
| N |
| Chi-square contribution |
| N / Row Total |
| N / Col Total |
| N / Table Total |
|-------------------------| Total Observations in Table: 84
| Arthritis$Improved Arthritis$Treatment | None | Some | Marked | Row Total |
--------------------|-----------|-----------|-----------|-----------| Placebo | 29 | 7 | 7 | 43 |
| 2.616 | 0.004 | 3.752 | |
| 0.674 | 0.163 | 0.163 | 0.512 |
| 0.690 | 0.500 | 0.250 | |
| 0.345 | 0.083 | 0.083 | |
--------------------|-----------|-----------|-----------|-----------| Treated | 13 | 7 | 21 | 41 |
| 2.744 | 0.004 | 3.935 | |
| 0.317 | 0.171 | 0.512 | 0.488 |
| 0.310 | 0.500 | 0.750 | |
| 0.155 | 0.083 | 0.250 | |
--------------------|-----------|-----------|-----------|-----------| Column Total | 42 | 14 | 28 | 84 |
| 0.500 | 0.167 | 0.333 | |
--------------------|-----------|-----------|-----------|-----------|
#多維列聯表
mytable<-xtabs(~Treatment+Sex+Improved,data=Arthritis) mytable #輸出結果
, , Improved = None Sex Treatment Female Male Placebo 19 10 Treated 6 7 , , Improved = Some Sex Treatment Female Male Placebo 7 0 Treated 5 2 , , Improved = Marked Sex Treatment Female Male Placebo 6 1 Treated 16 5
# ftable(mytable) Improved None Some Marked Treatment Sex Placebo Female 19 7 6 Male 10 0 1 Treated Female 6 5 16 Male 7 2 5
#邊際頻數
margin.table(mytable,1) Treatment Placebo Treated 43 41 margin.table(mytable,2) Sex Female Male 59 25 margin.table(mytable,3) Improved None Some Marked 42 14 28
#治療情況*改善情況的邊際頻數
margin.table(mytable,c(1,3)) Improved Treatment None Some Marked Placebo 29 7 7 Treated 13 7 21
#治療情況*性別的各類改善情況說明
ftable(prop.table(mytable,c(1,2))) Improved None Some Marked Treatment Sex Placebo Female 0.59375000 0.21875000 0.18750000 Male 0.90909091 0.00000000 0.09090909 Treated Female 0.22222222 0.18518519 0.59259259 Male 0.50000000 0.14285714 0.35714286 ftable(addmargins(prop.table(mytable,c(1,2)),3)) Improved None Some Marked Sum Treatment Sex Placebo Female 0.59375000 0.21875000 0.18750000 1.00000000 Male 0.90909091 0.00000000 0.09090909 1.00000000 Treated Female 0.22222222 0.18518519 0.59259259 1.00000000 Male 0.50000000 0.14285714 0.35714286 1.00000000
下面介紹三種檢驗：1.卡方獨立檢驗 2.Fisher精確檢驗 3.Cochran-Mantel-Haenszel檢驗
#卡方獨立性檢驗 library(vcd) #治療情況與改善情況 mytable<-xtabs(~Treatment+Improved,data=Arthritis) chisq.test(mytable) #結果 data: mytable X-squared = 13.055, df = 2, p-value = 0.001463 #性別和改善情況 mytable<-xtabs(~Improved+Sex,data=Arthritis) chisq.test(mytable) #結果 data: mytable X-squared = 4.8407, df = 2, p-value = 0.08889 #(p<0.01說明治療和改善有某種關系，p>0.05說明性別和改善效果沒有關系) #Fisher精確檢驗 mytable<-xtabs(~Treatment+Improved,data=Arthritis) fisher.test(mytable) #結果 data: mytable p-value = 0.001393 alternative hypothesis: two.sided #備注：fisher.test()函數可以在任意行列數大於2的二維列聯表中使用，但是不能用於2*2的列聯表 #Cochran-Mantel-Haenszel檢驗 #mantelhaen.test()函數用於Cochran-Mantel-Haenszel卡方檢驗 #檢驗治療情況和GIA是你情況在性別的每一水平下是否獨立 mytable<-xtabs(~Treatment+Improved+Sex,data=Arthritis) mantelhaen.test(mytable) data: mytable Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647 #結果表明，患者接受的治療和得到的改善在性別的每一個水平小並不獨立(分性別來看，用葯治療的患者較接受安慰劑的患者有了更好的改善) #相關性的度量 library(vcd) mytable<-xtabs(~Treatment+Improved,data=Arthritis) assocstats(mytable) #結果 X^2 df P(> X^2) Likelihood Ratio 13.530 2 0.0011536 Pearson 13.055 2 0.0014626 Phi-Coefficient : NA Contingency Coeff.: 0.367 Cramer's V : 0.394
3.相關
相關系數可以用來描述定量變量之間的關系。相關系數的符號(+-)表明關系的方向，其值的大小表示關系的強弱程度。
R可以計算多種相關系數，包括Pearson相關系數、Spearman相關系數、Kendall相關系數、偏相關系數、多分格相關系數和多系列相關系數
#1.Pearson、Spearman和Kendal相關 #Pearson積差相關系數衡量兩個變量之間的線型相關程度 #Spearman等級相關系數則衡量分級定序變量之間的相關程度 #Kendal's Tau相關系數也是一種非參數等級相關度量 #cor()函數計算這三種相關系數，而cov()函數可以用來計算協方差 #cor(x,use=,method=) #x 矩陣或者數據框 #use 指定趨勢值的處理方式 all.obs-遇到缺失值報錯、everything-遇到確實數據時,相關系數的計算結果將被設定為missing、complete.obs-行刪除、pairwise.complete.obs-成對刪除
states<-state.x77[,1:6] #計算方差和協方差
cov(states) #輸出結果
Population Income Illiteracy Life Exp Murder HS Grad Population 19931683.7588 571229.7796 292.8679592 -407.8424612 5663.523714 -3551.509551 Income 571229.7796 377573.3061 -163.7020408 280.6631837 -521.894286 3076.768980 Illiteracy 292.8680 -163.7020 0.3715306 -0.4815122 1.581776 -3.235469 Life Exp -407.8425 280.6632 -0.4815122 1.8020204 -3.869480 6.312685 Murder 5663.5237 -521.8943 1.5817755 -3.8694804 13.627465 -14.549616 HS Grad -3551.5096 3076.7690 -3.2354694 6.3126849 -14.549616 65.237894
#計算Pearson積差相關系數
cor(states) #輸出結果
Population Income Illiteracy Life Exp Murder HS Grad Population 1.00000000 0.2082276 0.1076224 -0.06805195 0.3436428 -0.09848975 Income 0.20822756 1.0000000 -0.4370752 0.34025534 -0.2300776 0.61993232 Illiteracy 0.10762237 -0.4370752 1.0000000 -0.58847793 0.7029752 -0.65718861 Life Exp -0.06805195 0.3402553 -0.5884779 1.00000000 -0.7808458 0.58221620 Murder 0.34364275 -0.2300776 0.7029752 -0.78084575 1.0000000 -0.48797102 HS Grad -0.09848975 0.6199323 -0.6571886 0.58221620 -0.4879710 1.00000000
#計算Spearman等級相關系數
cor(states,method="spearman") Population Income Illiteracy Life Exp Murder HS Grad Population 1.0000000 0.1246098 0.3130496 -0.1040171 0.3457401 -0.3833649 Income 0.1246098 1.0000000 -0.3145948 0.3241050 -0.2174623 0.5104809 Illiteracy 0.3130496 -0.3145948 1.0000000 -0.5553735 0.6723592 -0.6545396 Life Exp -0.1040171 0.3241050 -0.5553735 1.0000000 -0.7802406 0.5239410 Murder 0.3457401 -0.2174623 0.6723592 -0.7802406 1.0000000 -0.4367330 HS Grad -0.3833649 0.5104809 -0.6545396 0.5239410 -0.4367330 1.0000000
#結論 #收入和高中畢業率之間有很強的正相關 #文盲率和預期壽命之間存在很強的負相關
#偏相關
#控制一個或多個定量變量時，另外兩個定量變量之間的相互關系。
#ggm中的pcor()函數用於計算偏相關數
pcor(u,s)
#u是一個數值向量，其那兩個值表示要計算的相關系數的變量下標，其余的數值為條件變量的下標。s為變量的協方差陣。
library(ggm)
colnames(states)
pcor(c（1,5,2,3,6），cov(states))
#控制收入、文盲率和高中畢業率，人口和謀殺率之間的相關系數為0.346。
#其他相關
polycor()包的hetcor()函數可以計算一種混合的相關矩陣，其中包含數值型變量的Person積差相關系數、數值變量和有序變量之間的多系列相關系數、有序變量之間的多分格相關系數以及二分變量的四分相關系數。
相關性檢驗：使用cor.test()函數對單個的Person、Spearman、Kendall相關系數進行檢驗，簡化后的格式為：
cor.test(x,y,alternative,method=)
x,y為要檢驗的相關性的變量，alternative則用來指定進行雙側檢驗或單側檢驗，而method用以指定要計算的相關模型("person","kendall","spearman")。當研究的假設為總體的相關系數小於0時，要使用alternative="less";在研究的假設為總體的相關系數大於0時，要使用alternative="greater"。在默認情況下，alternative="two.side"
states<-state.x77[,1:6] cor.test(states[,3],states[,5]) #輸出結果 data: states[, 3] and states[, 5] t = 6.8479, df = 48, p-value = 1.258e-08 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.5279280 0.8207295 sample estimates: cor 0.7029752
#檢驗了預期壽命和謀殺率的Pearson相關系數為0的原假設。假設總體的相關度為0，則預計在一千萬次中只會右少於一次的機會見到0.703這樣打的樣本相關度。由於這種情況幾乎不可能發生，可以拒絕原來的假設
#對劇中進行顯著性檢驗
library(psych)
corr.test(states,use="complete")
Correlation matrix
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.00 0.21 0.11 -0.07 0.34 -0.10
Income 0.21 1.00 -0.44 0.34 -0.23 0.62
Illiteracy 0.11 -0.44 1.00 -0.59 0.70 -0.66
Life Exp -0.07 0.34 -0.59 1.00 -0.78 0.58
Murder 0.34 -0.23 0.70 -0.78 1.00 -0.49
HS Grad -0.10 0.62 -0.66 0.58 -0.49 1.00
Sample Size
[1] 50
Probability values (Entries above the diagonal are adjusted for multiple tests.)
Population Income Illiteracy Life Exp Murder HS Grad
Population 0.00 0.59 1.00 1.0 0.10 1
Income 0.15 0.00 0.01 0.1 0.54 0
Illiteracy 0.46 0.00 0.00 0.0 0.00 0
Life Exp 0.64 0.02 0.00 0.0 0.00 0
Murder 0.01 0.11 0.00 0.0 0.00 0
HS Grad 0.50 0.00 0.00 0.0 0.00 0
人口數量和高中畢業率的相關系數(-0.10)並不顯著地不為0(p=0.5)
#其他顯著性檢驗
psych包中的pcor.test()函數可以用來檢驗在控制一個或多個額外變量時兩個變量之間的條件獨立性。使用格式為：
pcor.test(r,q,n)
其中r是由pcor()函數計算得到的偏相關系數，q為要控制的變量數，n為樣本大小。
4.t檢驗
研究中最常見的行為就是對兩組進行比較。接受某種新葯治療的患者是否較使用某種現有葯物的患者表現處了更大的改善？某種工藝是否較另外一種工藝制品的不合格率更少？兩種教學手法，哪一種更有效？這里我們將關注結果為連續型的組間比較，並假設其呈正態分布。
#t檢驗調用格式：
t.test(y~x,data) 其中y是一個數值型變量，x是一個二分變量。調用格式為： t.test(y1,y2) 其中用y1和y2為數值型向量。可選參數data的取值為一個包含了這些變量的矩陣或數據框 #比較南方group1和非南方group0各州的監禁概率
library(psych) states<-state.x77[,1:6] corr.test(states,use="complete") library(MASS) t.test(Prob~So,data=UScrime) Welch Two Sample t-test data: Prob by So t = -3.8954, df = 24.925, p-value = 0.0006506 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.03852569 -0.01187439 sample estimates: mean in group 0 mean in group 1
0.03851265 0.06371269 p<0.001 南方各州和非南方各州監禁的概率不同
非獨立樣本的t檢驗
例子：較年輕(14~24)男性的失業率是否和年長的(35~39）男性的失業率更高?在這種情況下，這兩組數據並不獨立。你不能說亞拉巴州的年輕男性和年長男性的失業率沒有關系。在兩組的觀測之間相關時，你獲得的是一個非獨立設計。前后測設計或者重復測量設計同樣產生非獨立的組。
非獨立樣本的t檢驗假定組間的差異呈正態分布。調用的格式為：
t.test(y1,y2,paired=TRUE)
其中y1和y2為兩個非獨立組的數值向量。結果如下：
library(MASS) sapply(UScrime[c("U1","U2")],function(x)(c(mean(x),sd=sd(x))))
with(UScrime,t.test(U1,U2,paried=TRUE)) #輸出結果 U1 U2 95.46809 33.97872 sd 18.02878 8.44545
#輸出結果
data: U1 and U2
t = 21.174, df = 65.261, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
55.69010 67.28862
sample estimates:
mean of x mean of y
95.46809 33.97872
#差異的均值(61.5）足夠大，可以保證拒絕年終和年輕的平均失業率相同的假設。年輕男性的失業率更高一些。
5.組間差異的非參數檢驗
若兩組數據獨立，可以使用Wilcoxon秩和檢驗來評估觀測是否是從相同的概率分布中抽得的，調用格式為：
wilcox.test(y~x,data)
wilcox.test(y1,y2)
其中y1和y2為各組的結果變量。可選參數data的取值為一個包含了這些變量的矩陣或者數據框。默認進行一個雙側檢驗。可以添加參數excat來進行精確檢驗，指定alternative="less"或alternative="greater"進行方向的檢驗。
#使用Mann-Whitney U檢驗監禁率的問題
with(UScrime,by(Prob,So,median)) So: 0 [1] 0.038201
-------------------------------------------------------------------- So: 1 [1] 0.055552 wilcox.test(Prob~So,data=UScrime) data: Prob by So W = 81, p-value = 8.488e-05 alternative hypothesis: true location shift is not equal to 0 #p<0.001 說明南方各州和非南方各州監禁率相關的假設不成立
sapply(UScrime[c("U1","U2")],median) U1 U2 92 34 with(UScrime,wilcox.test(U1,U2,paired=TRUE)) data: U1 and U2 V = 1128, p-value = 2.464e-09 alternative hypothesis: true location shift is not equal to 0
多於兩組的比較，如果各組不獨立，Friedman檢驗會更合適。
Kruskal-Wallis檢驗的調用格式為：
kruskal.test(y~A,data)
其中y是一個數值型結果變量，A是一個擁有兩個或更多水平的分組變量。而Friedman檢驗的調用格式為：
friedman.test(y~A|B,data)
其中A是一個分組變量，而B是一個用以認定匹配觀測的區間變量。
以上為基本的統計分析。