R語言實戰學習筆記-基本統計分析

時間 2019-11-26

標籤語言實戰學習筆記基本統計分析简体版

原文原文鏈接

本文主要介紹如下內容：1.基本統計方法 2頻數表和列聯表 3.相關 4.t檢驗 5.組間差別的非參數檢驗 app

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)函數，使用格式以下：less

sapply(x,FUN,options)ide

其中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。spa

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.1737358code

分組計算的擴展，doBy包和psych包提供了分組計算的描述性統計量的函數，doBy包中的summaryBy()函數使用的基本格式：orm

#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檢驗blog

#卡方獨立性檢驗
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.相關

相關係數能夠用來描述定量變量之間的關係。相關係數的符號(+-)代表關係的方向，其值的大小表示關係的強弱程度。ip

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