1. 程式人生 > >R語言實戰之基本統計分析

R語言實戰之基本統計分析

第7章 基本統計分析

在資料被組織成合適的形式後,可以使用圖形探索資料,接下來是使用數值描述每個變數的分佈,然後則是兩兩探索所選擇變數之間的關係。
本章將評述用於生成基本的描述性統計量和推斷統計量的R函式。

7.1 描述性統計分析

本節介紹分析連續型變數中心趨勢、變化性和分佈性的方法。
使用第1章中Motor Trend雜誌的車輛路試(mtcars)資料集。重點關注每加侖汽油行駛英里數(mpg)、馬力(hp)和車重(wt)。

> vars <- c("mpg", "hp", "wt")
> head(mtcars[vars])
> > #檢視三個變數的資料
mpg hp wt Mazda RX4 21.0 110 2.620 Mazda RX4 Wag 21.0 110 2.875 Datsun 710 22.8 93 2.320 Hornet 4 Drive 21.4 110 3.215 Hornet Sportabout 18.7 175 3.440 Valiant 18.1 105 3.460

將首先檢視所有32中車型的描述性統計量,然後按照變速箱型別(am)和汽缸數(cyl)考察描述性統計量。變速箱型別是以0表示自動檔、1表示手動檔來編碼的二分變數,汽缸數為4、5或6。

7.1.1 方法雲集

在描述性統計量的計算方面,R中的選擇非常多。從基礎安裝中包含的函式開始,然後檢視那麼些使用者貢獻包中的擴充套件函式。
集中安裝中,可用summary( )函式獲取描述性統計量。
程式碼清單7-1 通過summary( )計算描述性統計量

> #7-1 
> summary(mtcars[vars])
      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( )函式提供了最小值、最大值、四分位數和數值型變數的均值,以及因子向量和邏輯型向量的頻數統計。

第5章中的apply( )函式或sapply( )函式計算所選擇的任意描述性統計量。對於sapply( )函式,使用格式為:
sapply(x, FUN, options)
其中x是資料框(或矩陣),FUN為一個任意函式。如果指定了options,它們將被傳遞給FUN。可以在這裡插入的典型函式有mean、sd、var、min、median、length、range和quantile。函式fivenum( )可返回圖基五數總括(即最小值、下四分位數、中位數、上四分位數和最大值)。

程式碼清單7-2 通過sapply( )計算描述性統計量

> mystats <- function(x, na.omit=FALSE) { #na.omit為0不移除缺失觀測的整行
+              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))
+ }
> sapply(mtcars[vars], 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 #峰度

樣本中的車型,每加侖汽油行駛英里數的平均值為20.1,標準差為6.0。分佈呈現右偏(偏度+0.61),並且較正態分佈稍平(峰度-0.37)。
如果只希望單純地忽略缺失值,應當使用:
sapply(mtcars[vars], mystats, na.omit=TRUE)

Himsc包中的describe( )函式可返回變數和觀測的數量、缺失值和唯一值的數目、平均值、分位數,以及五個最大的值和五個最小的值。

程式碼清單7-3 通過Hmisc包中的describe( )函式計算描述性統計量

> library(Hmisc) 
> describe(mtcars[vars])
mtcars[vars]
3 Variables 32 Observations —————————————————————————– mpg n missing distinct Info Mean Gmd .05 .10 32 0 25 0.999 20.09 6.796 12.00 14.34 .25 .50 .75 .90 .95 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 .10 32 0 22 0.997 146.7 77.04 63.65 66.00 .25 .50 .75 .90 .95 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 .10 32 0 29 0.999 3.217 1.089 1.736 1.956 .25 .50 .75 .90 .95 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 —————————————————————————–

pastecs包中有一個名為stat..desc( )的函式,可以計算種類繁多的描述性統計量。使用格式為:
stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
其中的x是一個數據框或時間序列。若basic=TRUE(預設值),則計算其中所有值、空值、缺失值的數量,以及最小值、最大值、值域和綜合。若desc=TRUE(預設值),則計算中位數、平均數、平均數的標準誤差,平局數置信度為95%的置信區間、方差、標準差以及變異係數。若norm=TRUE(非預設),則返回正態分佈的統計量,包括偏度和峰度(以及它們的統計顯著程度)和Shapiro - Wilk正態檢驗結果。這裡使用了p值來計算平均數的置信區間(預設為0.95)。

程式碼清單7-4 通過pastecs包中的stat.desc( )函式計算描述性統計量

> library(pastecs)
> stat.desc(mtcars[vars])
                     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

psych包中的describe( )函式,可以計算非缺失值的數量、平均數、標準差、中位數、截尾均值、絕對中位差、最小值、最大值、值域、峰度和平均值的標準誤差。

程式碼清單7-5 通過psych包中的describe( )計算描述性統計量

> library(psych)
> describe(mtcars[vars])
    vars  n   mean    sd median trimmed   mad   min    max  range skew
mpg    1 32  20.09  6.03  19.20   19.70  5.41 10.40  33.90  23.50 0.61
hp     2 32 146.69 68.56 123.00  141.19 77.10 52.00 335.00 283.00 0.73
wt     3 32   3.22  0.98   3.33    3.15  0.77  1.51   5.42   3.91 0.42
    kurtosis    se
mpg    -0.37  1.07
hp     -0.14 12.12
wt     -0.02  0.17

psych包和Hmisc包均提供了名為describe( )的函式。R在呼叫時優先呼叫最後載入的程式包中的。如果強制呼叫,可以使用Hmisc::describe( )

7.1.2 分組計算描述型統計量

在比較多組個體或觀測時,關注的焦點經常是各組的描述性統計資訊。在R中完整這個任務有託幹中方法,以變速箱型別各水平的描述性統計量開始。
在第5章中,介紹了整合資料的方法。可以使用aggregate( )函式(5.6.2節)來分組獲取描述性統計量。

程式碼清單7-6 使用aggregate( )分組獲取描述性統計量

> aggregate(mtcars[vars], 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(mtcars[vars], by=list(am=mtcars$am), sd)
  am      mpg       hp        wt
1  0 3.833966 53.90820 0.7774001
2  1 6.166504 84.06232 0.6169816

注意list(am=mtcars$am)</code>的使用。如果使用的是<code>list(mtcars$am),則am列被標註為Group.1而不是am。通過這個賦值指定了一個更有幫助的列標籤。如果有多個分組變數,可以使用by=list(name1=groupvar1, name2=groupvar2, ... , groupvarN)這樣的語句。
但是,aggregate( )僅允許在每次呼叫中使用平均數、標準差這樣的單返回值函式。無法一次返回若干個統計量。此時可以使用by( )函式,格式為:
by(data, INDICES, FUN)
其中data是一個數據框或矩陣,INDICES是一個因子或因子組成的列表,定義了分組,FUN是任意函式。

程式碼清單7-7 使用by( )分組計算描述性統計量

> 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( )函式的使用格式為:
summaryBy(formula, data=dataframe, FUN=function)
其中的formula接受以下的格式:
var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 + ... + groupvarN
~左側的變數是需要分析的數值型變數,而右側的變數是類別型的分組變數。function可為任何內建或使用者自編的R函式。

程式碼清單7-8 使用doBy包中的summaryBu( )分組計算概述統計量

> 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
1  0    19 17.14737  3.833966 0.01395038   -0.8031783   19 160.2632 53.90820
2  1    13 24.39231  6.166504 0.05256118   -1.4553520   13 126.8462 84.06232
      hp.skew hp.kurtosis wt.n  wt.mean  wt.stdev   wt.skew wt.kurtosis
1 -0.01422519  -1.2096973   19 3.768895 0.7774001 0.9759294   0.1415676
2  1.35988586   0.5634635   13 2.411000 0.6169816 0.2103128  -1.1737358

psych包中的describeBy( )函式可計算和describe相同的描述性統計量,只是按照一個或多個分組變數分層。

程式碼清單7-9 使用psych包中的describe.by( )`分組計算概述統計量

> library(psych)
> myvars <- c("mpg", "hp", "wt")
> describe.by(mtcars[myvars], mtcars$am)
> #describe.by棄用,請使用descriBy( )
$`0`
    vars  n   mean    sd median trimmed   mad   min    max  range  skew
mpg    1 19  17.15  3.83  17.30   17.12  3.11 10.40  24.40  14.00  0.01
hp     2 19 160.26 53.91 175.00  161.06 77.10 62.00 245.00 183.00 -0.01
wt     3 19   3.77  0.78   3.52    3.75  0.45  2.46   5.42   2.96  0.98
    kurtosis    se
mpg    -0.80  0.88
hp     -1.21 12.37
wt      0.14  0.18

$`1`
    vars  n   mean    sd median trimmed   mad   min    max  range skew
mpg    1 13  24.39  6.17  22.80   24.38  6.67 15.00  33.90  18.90 0.05
hp     2 13 126.85 84.06 109.00  114.73 63.75 52.00 335.00 283.00 1.36
wt     3 13   2.41  0.62   2.32    2.39  0.68  1.51   3.57   2.06 0.21
    kurtosis    se
mpg    -1.46  1.71
hp      0.56 23.31
wt     -1.17  0.17

describeBy( )不允許指定任意函式,適應性較差。若存在多個分組變數,可使用list(groupvar1, groupvar2, ... , groupvarN)表示它們。但這僅在分組變數交叉口不出現空白單元時有效。

最後,使用5.6.3節中描述的reshape包靈活地按組匯出描述性統計量。首先,使用:
dfm <- melt(dataframe, mesure.vars=y, id.vars=g)
融合資料框。其中的dataframe包含這資料,y是一個向量,指明瞭要進行概述的數值型變數(預設使用所有的變數),而g是由一個或多個分組變數組成的向量。然後使用:
cast(dfm, groupvar1 + groupvar2 + ... + variable ~ . , FUN)
重鑄資料。分組變數以+號分隔,這裡的variable只取其字面含義,而FUN是一個任意函式。

程式碼清單7-10 通過reshape包分組計算概述統計量

> library(reshape)
> dstats <- function(x) (c(n=length(x), mean=mean(x), sd=sd(x)))
> #定義函式dstats,求均值,標準差
> dfm <- melt(mtcars, measure.vars=c("mpg", "hp", "wt"), id.vars=c("am", "cyl"))
> #資料融合,按照am,cyl對mpg、hp、wt融合統計
> cast(dfm, am + cyl + variable ~ ., dstats)
> 按照am,cyl分組,對其變數求均值、標準差
   am cyl variable  n       mean         sd
1   0   4      mpg  3  22.900000  1.4525839
2   0   4       hp  3  84.666667 19.6553640
3   0   4       wt  3   2.935000  0.4075230
4   0   6      mpg  4  19.125000  1.6317169
5   0   6       hp  4 115.250000  9.1787799
6   0   6       wt  4   3.388750  0.1162164
7   0   8      mpg 12  15.050000  2.7743959
8   0   8       hp 12 194.166667 33.3598379
9   0   8       wt 12   4.104083  0.7683069
10  1   4      mpg  8  28.075000  4.4838599
11  1   4       hp  8  81.875000 22.6554156
12  1   4       wt  8   2.042250  0.4093485
13  1   6      mpg  3  20.566667  0.7505553
14  1   6       hp  3 131.666667 37.5277675
15  1   6       wt  3   2.755000  0.1281601
16  1   8      mpg  2  15.400000  0.5656854
17  1   8       hp  2 299.500000 50.2045815
18  1   8       wt  2   3.370000  0.2828427
> dfm
   am cyl variable   value
1   1   6      mpg  21.000
2   1   6      mpg  21.000
3   1   4      mpg  22.800
4   0   6      mpg  21.400

7.1.3 結果的視覺化

分佈特徵的數值刻畫很重要,但不能代替視覺呈現。

7.2 頻數表和列聯表

類別型變數的頻數表和列聯表,以及相應的獨立性檢驗、相關性的度量、圖形化展示結果的方法。
本節資料來自vcd包中的Arthritis資料集。表示一項風溼性關節炎新療法的雙盲臨床實驗結果。
資料展示:

> #資料展示
> library(vcd)
載入需要的程輯包:grid
> head(Arthritis)
  ID Treatment  Sex Age Improved
1 57   Treated Male  27     Some
2 46   Treated Male  29     None
3 77   Treated Male  30     None
4 17   Treated Male  32   Marked
5 36   Treated Male  46   Marked
6 23   Treated Male  58   Marked

7.2.1 生成頻數表

R中提供了用於建立頻數表和列聯表的若干方法。部分列於表中:

表7-1 頻數表和列聯表函式.png

1. 一維列聯表

可以使用table( )函式生成簡單的頻數統計表。示例如下:

> mytable <- with(Arthritis, table(Improved))
> #with限定資料框
> mytable
Improved
  None   Some Marked 
    42     14     28 
> prop.table(mytable)
> #表中條目轉化為比例值
Improved
     None      Some    Marked 
0.5000000 0.1666667 0.3333333 
> prop.table(mytable)*100
> #百分比
Improved
    None     Some   Marked 
50.00000 16.66667 33.33333 

2. 二維列聯表

對於二維列聯表,table( )函式的使用格式為:
mytable <- table(A, B)
其中A是行變數,B是列變數。
xtabs( )函式還可使用公式風格的輸入建立列聯表,格式為:
mytable <- xtabs(~ A + B, data=mydata)
其中mydata是一個矩陣或資料框。要進行交叉分類的變數應出現在公式的右側(即~符號右方),以+作為分隔符。若某個變數解除安裝公式左側,則其為一個頻數向量(在資料已經被表格化時很有用)。

> 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

這裡的下標2指代table( )語句中的第二個變數。

各單元格比例可用如下語句獲取:

> prop.table(mytable)
         Improved
Treatment       None       Some     Marked
  Placebo 0.34523810 0.08333333 0.08333333
  Treated 0.15476190 0.08333333 0.25000000

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( )時,預設為表中所有變數建立邊際和
> 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
> #新增列和

注意:table( )函式預設忽略缺失值(NA)。要在頻數統計中將NA視為一個有效的類別,請設定引數useNA="ifany"

使用gmodels包中的CrossTable( )函式是建立二維列聯表的第三種方法。CrossTable( )函式仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二維列聯表。

程式碼清單7-11 使用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 |           | 
--------------------|-----------|-----------|-----------|-----------|

CrossTable( )函式有很多選項,可以計算(行、列、單元格)的百分比,指定小數位數;進行卡方、Fisher和McNemar獨立性檢驗;計算期望和(皮爾遜、標準化、調整的標準化)殘差;將缺失值作為一種有效值;進行行和列標題的標註;生成SAS或SPSS風格的輸出。相見help(CrossTable)

3、多維列聯表

table( )xlab( )都可以基於三個或更多的類別型變數生成多維列聯表。margin.table( )prop.table( )addmargins( )函式可以自然地推廣到高於二維的情況。另外,ftable( )函式可以緊湊地輸出多維列聯表。
程式碼清單7.12 多維列聯表

> library(vcd)
載入需要的程輯包:grid
> 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) #ftable輸出緊湊表格
                 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
> 

比例將被新增到不在prop.table( )呼叫中的下標上(上例中是第三個下標,或Improve)。在最後一個表格中看出,因為那裡為第三個下標添加了邊際和。

列聯表可以展示組成表格的各種變數組合的頻數或比例。

7.2.2 獨立性檢驗

R提供了多種檢驗類別型變數獨立性的方法。

1、卡方獨立性檢驗

chisq.test( )函式可對二維表的行變數進行卡方獨立性檢驗。
程式碼清單7-13 卡方獨立性檢驗

> library(vcd)
> mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
> chisq.test(mytable)

        Pearson's Chi-squared test

data:  mytable
X-squared = 13.055, df = 2, p-value = 0.001463

> mytable <- xtabs(~ Improved + Sex, data = Arthritis)
> chisq.test(mytable)

        Pearson's Chi-squared test

data:  mytable
X-squared = 4.8407, df = 2, p-value = 0.08889

Warning message:
In chisq.test(mytable) : Chi-squared近似演算法有可能不準

第一個計算結果:p<0.01,患者接收的治療和改善的水平可能存在某種關係;
第二個計算結果:p>0.05,患者性別和改善情況之間互相獨立。

2、Fisher精確檢驗

fisher.test( )函式可進行FIsher精確檢驗。Fisher精確檢驗的厡假設是:邊界固定的列聯表中行和列是互相獨立的。呼叫格式為:
fisher.test(mytable)
其中mytable為一個二維列聯表。

> mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
> fisher.test(mytable)

        Fisher's Exact Test for Count Data

data:  mytable
p-value = 0.001393
alternative hypothesis: two.sided

這裡的fisher.test( )函式可以在任意行列數大於等於2的二維列聯表上使用,但不能用於2×2列聯表。

3、Cochran-Mantel-Haenszel檢驗

mentelhaen.test( )函式可用來進行Cochran-Mantel-Haenszel卡方檢驗,其原假設是,兩個名義變數在第三個變數的每一層中都是條件獨立的。
下列程式碼檢驗治療情況和改善情況在性別的每一水平下是否獨立。此檢驗假設不存在三階互動作用(治療情況×改善情況×性別)。

> mytable <- xtabs(~ Treatment + Improved + Sex, data = Arthritis)
> mantelhaen.test(mytable)

        Cochran-Mantel-Haenszel test

data:  mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647

> 

結果表明,患者接收的治療與得到的改善在性別的每一水平下並不獨立(即,分性別來看,用藥治療的患者較接收安慰劑的患者有了更多的改善)。

7.2.3 相關性度量

顯著性檢驗評估了是否存在充分的證據以拒絕變數見相互獨立的厡假設。如果可以拒絕原假設,問題轉向用以衡量相關性強弱的相關性度量。
vcd包中的assocstats( )函式可以用來計算二維列聯表的phi係數、列聯絡數和Cramer’s V係數。
程式碼清單7-14 二維列聯表的相關性度量

> 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 

總體來說,較大的值意味著較強的相關性。vcd包也提供了一個kappa( )函式,可以計算混淆矩陣的Cohen’s kappa值以及加權的kappa值。(混淆矩陣可以表示兩位評判者對於一系列物件進行分類所得結果的一致程度。)

7.2.4 結果的視覺化

vcd包中擁有優秀的、用於視覺化多維資料集中類別型變數間關係的函式,可以繪製馬賽克圖和關聯圖(參見11.4節)。ca包中的對應分析函式允許使用多種幾何表示(Nenadic & Greenacre,2007)可視地探索列聯表中行和列之間的關係。

7.2.5 將錶轉換為扁平格式

擁有一個列聯表如何轉化為原始資料,提供一個示例:
程式碼清單7-15 通過table2flat將錶轉換為扁平格式

> table2flat <- function(mytable) {
+   df <- as.data.frame(mytable)
+   rows <- dim(df)[1]
+   cols <- dim(df)[2]
+   x <- NULL
+   for (i in 1:rows) {
+     for (j in 1:df$Freq[i]) {
+       row <- df[i, c(1:(cols-1))]
+       x <- rbind(x, row)
+     }
+   }
+   row.names(x) <- c(1:dim(x)[1])
+   return(x)
+ }

此函式可以接收一個R中的表格(行列數任意)並返回一個扁平格式的資料框。也可以使用這個函式輸入別處的表格。
以下表為例進行轉換輸入

7-16 原始表格.png

程式碼清單7-16 使用table2flat( )`函式轉換已發表的資料

> treatment <- rep(c("Placebo", "Treated"), times=3)
> improved <- rep(c("None", "Some", "Marked"), each=2)
> Freq <- c(29, 13, 7, 17, 7, 21)
> #縱向資料
> mytable <- as.data.frame(cbind(treatment, improved, Freq))
> mydata <- table2flat(mytable)
> head(mydata)
  treatment improved
1   Placebo     None
2   Placebo     None
3   Placebo     None
4   Placebo     None
5   Treated     None
6   Placebo     Some

7.3 相關

相關係數可以用來描述定量變數之間的關係。相關係數的符號(±)表明關係的方向(正相關或負相關),其值表示關係的強弱程度。
本節使用R基礎安裝中的state.x77資料集,它提供了美國50個州在1977年的人口、收入、文盲率等資料。此外使用psych包和ggm包。

7.3.1 相關型別

R可以計算多種相關係數,包括Pearson相關係數、Spearman相關係數、Kendall相關係數、偏相關係數、多分格(polychoric)相關係數和多系列(polyserial)相關係數。

1、Prarson、Spearman和Kendall相關係數

Pearson積差相關係數衡量了兩個定量變數之間的線性相關程度。Spearman等級相關係數則衡量分級定序變數之間的相關程度。Kendall’sTau相關係數也是一種非引數的等級相關度量。
cor( )函式可以計算這三種相關關係,而cov( )函式可用來計算協方差。兩個函式的引數有很多,其中與相關係數的計算有關的引數可以簡化為:
cor(x, use= , method= )

  • x是矩陣或資料框;
  • use指定缺失資料的處理方式(all.obs假設無缺失資料,遇到將報錯;everything缺失資料相關係數計算為missing;complete.obs行刪除;pairwise.complete.obs成對刪除);
  • method指定相關係數的型別,pearson、spearman或kendall。
  • 預設的引數為use="everything"method="pearson"

程式碼清單7-17 協方差和相關係數

> states <- state.x77[, 1:6]
> cov(states)
              Population      Income   Illiteracy     Life Exp      Murder
Population 19931683.7588 571229.7796  292.8679592 -407.8424612 5663.523714
Income       571229.7796 377573.3061 -163.7020408  280.6631837 -521.894286
Illiteracy      292.8680   -163.7020    0.3715306   -0.4815122    1.581776
Life Exp       -407.8425    280.6632   -0.4815122    1.8020204   -3.869480
Murder         5663.5237   -521.8943    1.5817755   -3.8694804   13.627465
HS Grad       -3551.5096   3076.7690   -3.2354694    6.3126849  -14.549616
                HS Grad
Population -3551.509551
Income      3076.768980
Illiteracy    -3.235469
Life Exp       6.312685
Murder       -14.549616
HS Grad       65.237894
> cor(states)
            Population     Income Illiteracy    Life Exp     Murder
Population  1.00000000  0.2082276  0.1076224 -0.06805195  0.3436428
Income      0.20822756  1.0000000 -0.4370752  0.34025534 -0.2300776
Illiteracy  0.10762237 -0.4370752  1.0000000 -0.58847793  0.7029752
Life Exp   -0.06805195  0.3402553 -0.5884779  
            
           

相關推薦

R語言實戰基本統計分析

第7章 基本統計分析 在資料被組織成合適的形式後,可以使用圖形探索資料,接下來是使用數值描述每個變數的分佈,然後則是兩兩探索所選擇變數之間的關係。 本章將評述用於生成基本的描述性統計量和推斷統計量的R函式。 7.1 描述性統計分析 本節介紹分析連續型變

R語言進行詞雲統計分析

R語言進行詞雲統計分析 本文章從爬蟲、詞頻統計、視覺化三個方面講述了R語言的具體應用,歡迎大家共同談論學習 1、使用 rvest 進行資料的爬取 #如果沒有,先安裝rvest包 install.packages("rvest") library(rvest) url <- "http://www.

R語言學習基本語法

ggplot2圖形之基本語法: ggplot2的核心理念是將繪圖與資料分離,資料相關的繪圖與資料無關的繪圖分離 ggplot2是按圖層作圖 ggplot2保有命令式作圖的調整函式,使其更具靈活性 ggplot2將常見的統計變換融入到了繪圖中。 ggplot的繪圖有以下幾個特

R語言學習聚類分析

1.動態聚類:k—means 基本思想: (1)選擇K個點作為質心 (2)將每個點指派到最近的質心,形成K個類 (3)重新計算每個類的質心 (4)重複2—3知道質心不發生變化 例項: 優缺點: (1)有效率且不易受初始值的影響 (2)不能處理非球形簇 (3)不能處理不同尺

R語言實戰 - 基本統計分析(1)- 描述性統計分析

4.3 summary eas 方法 func -- 4.4 1.0 6.5 > vars <- c("mpg", "hp", "wt") > head(mtcars[vars]) mpg hp wt Maz

R語言- 基本統計分析

kruskal 最大的 turn clas 技術 ria 大於 stat pair 目的:   1.描述性統計分析   2.頻數表和;列連表   3.相關系數和協方差   4.t檢驗   5.非參數統計 在上一節中使用了圖形來探索數據,下一步就是給出具體的數據來描述每個變量

大數據分析學習使用R語言實戰機器學習視頻課程

https aid 通過 原理 har fsg follow mdf 學習 大數據分析學習之使用R語言實戰機器學習網盤地址:https://pan.baidu.com/s/1Yi9H6s8Eypg_jJJlQmdFSg 密碼:0jz3備用地址(騰訊微雲):https://s

R語言進行基本統計分析

1. 描述性統計分析 使用自帶的summary()函式 > myvars <- c("mpg","hp","wt") > summary(mtcars[myvars]) mpg hp wt Mi

R語言基本統計分析

描述性統計分析 #利用(mtcars)資料集,我們提取出英里數(mpg),馬力(hp),車重(wt) > myvars <- c("mpg","hp","wt") > head(

統計統計分析R語言方面的圖書

統計之都的成員編著、翻譯了大量關於統計分析和R語言方面的圖書。 已出版 讀者可以點選下面每本書的連結進入該書的的頁面,下載隨書程式碼,我們還會不定期釋出圖書的勘誤,也歡迎讀者留言提問。 免費電子書 即將出版 陳逸波翻譯:The R Bookfrom: http://cos.name/books/

R語言實戰 - 基本數據管理(3)

cat taf str logs 合並 exc number country lob 8. 數據排序 > leadership$age [1] 32 45 25 39 NA > newdata <- leadership[order(leadership

[讀書筆記] R語言實戰 (四) 基本數據管理

mean 圖片 數值 函數 nbsp 一個 img order 分享 1. 創建新的變量 mydata<-data.frame(x1=c(2,2,6,4),x2=c(3,4,2,8)) #方法一 mydata$sumx<-mydata$x1+mydat

R語言實戰 R的使用

works 結束 runif 摘要 區分大小寫 目錄 r語 work text 註意 R區分大小寫 R的對象可以是任何東西(數據、函數、圖形、分析結果等) 賦值:<-(常用),=(少用); #註釋 可使用上下方向鍵查看已輸入命令的歷史記錄 工作空間 工作空間(w

R語言實戰 創建數據集(第二章,各種數據結構)

前三 實戰 創建 data div mask 結果 搜索 全局 數據集 2.1數據集概念 概念:通常是由數據構成的矩形數據 不同行業對數據集的行和列叫法不同 行業人 行 列 統計學家 觀測(observation) 變量(variable) 數據庫分析師 記錄

R語言實戰 圖形初階(第三章)-- 初識

space wid spa 開啟 display tps ping microsoft 目標 圖形初級 3.1 使用圖形 在交互式會話中,通過組條輸入語句構建圖形,直至得到想要的效果 attach(mtcars)

(數據科學學習手劄19)R基本統計分析技巧總結

misc 總結 4.6 內部 red margin adjust 條件 置信區間 在獲取數據,並且完成數據的清洗之後,首要的事就是對整個數據集進行探索性的研究,這個過程中會利用到各種描述性統計量和推斷性統計量來初探變量間和變量內部的基本關系,本篇筆者便基於R,對一些常用的數

R-基本統計分析--描述性統計分析

及其 pre dice 數據集 returns length 平均值 sun 52.0 描述性統計分析主要包括 基本信息:樣本數、總和 集中趨勢:均值、中位數、眾數 離散趨勢:方差(標準差)、變異系數、全距(最小值、最大值)、內四分位距(25%分位數、75%分位數) 分布

R-基本統計分析--獨立、相關性及其檢驗

獨立性檢驗-(卡方、Fisher) 獨立性檢驗_百度百科 https://baike.baidu.com/item/%E7%8B%AC%E7%AB%8B%E6%80%A7%E6%A3%80%E9%AA%8C/4031921?fr=aladdin R之獨立性檢驗 -

R語言開發協方差分析瞭解下

我們通常使用迴歸分析來建立描述預測變數變數對響應變數的影響的模型,有時,如果我們有類似於是/否或男/女等值的分類變數,簡單迴歸分析為分類變數的每個值提供多個結果。在這種情況下,我們可以通過使用分類變數和預測變數來研究分類變數的影響,並比較分類變數的每個級別的迴歸線。 這樣的分

R筆記20181103基本統計分析

一、描述性統計分析 1、峰度與偏度 網上找了一些相關資料 2、aggregate() aggregate(x, by, FUN, ..., simplify = TRUE, drop = TRUE) 將x按by分組後將FUN作用於每一個子集上,其中by是一個l