《R语言实战》学习笔记seventh

时间：2014-05-28 09:36:48 阅读：365 评论：0 收藏：0 [点我收藏+]

标签：des c class blog code a

由于在准备软考中级数据库系统工程师外加巩固SQL Server 2012，所以拖了好久一直没继续学R 下去

所以今天重开R 的战事

这次是关于基本统计分析的内容，即关于用于生成基本的描述性统计量和推断统计量的R 函数

首先，将着眼于定量变量的位置和尺度的衡量方式

然后将是生成类别型变量的频数表和列联表的方法（以及连带的卡方检验）

接下来将考察连续型和有序型变量相关系数的多种形式

最后转而通过参数检验（t检验）和非参数检验（Mann-Whitney U检验、Kruskal-Wallis检验）方法研究组间差异

描述性统计分析

在这一节里将关注分析连续型变量的中心趋势、变化性和分布形状的方法

> 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

使用的数据集是Motor Trend的mtcars数据集

首先查看所有32种车型的描述性统计量，然后按照变速箱类型am和汽缸数cyl考差描述性统计量；变速箱类型是一个以0表示自动档、1表示手动档来编码的二分变量。而汽缸数可为4、5和6

方法云集

> 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( )函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计

> 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))
+                   }
> 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

对于sapply( )函数，其使用格式：sapply(x, FUN, options)

其中x是数据库或矩阵，FUN为一个任意的函数，如果指定了options，将被传递给FUN；插入的函数可以用mean, sd, var, min, max, median, length, range, quantile; 函数fivenum( )可返回图基五数总括（Tukey‘s five-number summary, 即最小值、下四分数、中位数、上四分位数、最大值）；R 的基础安装中没有提供偏度和峰度的计算，但是可以自行添加

如果希望单纯地忽略缺失值，那么应当使用sapply(mtcars[vars], mystats, na.omit=TRUE)

然后再来一些扩展把

除了基础安装之外，R 中许多包也提供了计算描述性统计量的函数，其中包括Hmisc, pastecs, psych

> library(Hmisc)
载入需要的程辑包：grid
载入需要的程辑包：lattice
载入需要的程辑包：survival
载入需要的程辑包：splines
载入需要的程辑包：Formula
 
载入程辑包：‘Hmisc’
 
下列对象被屏蔽了from 
‘package:plyr’:
 
    is.discrete, summarize
 
下列对象被屏蔽了from 
‘package:base’:
 
    format.pval, round.POSIXt, trunc.POSIXt, units
 
> describe(mtcars[vars])
mtcars[vars] 
 
 3  Variables      32  Observations
-----------------------------------------------------------------------------
mpg 
      n missing  unique    Mean     .05     .10     .25     .50     .75 
     32       0      25   20.09   12.00   14.34   15.43   19.20   22.80 
    .90     .95 
  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  unique    Mean     .05     .10     .25     .50     .75 
     32       0      22   146.7   63.65   66.00   96.50  123.00  180.00 
    .90     .95 
 243.50  253.55 
 
lowest :  52  62  65  66  91, highest: 215 230 245 264 335 
-----------------------------------------------------------------------------
wt 
      n missing  unique    Mean     .05     .10     .25     .50     .75 
     32       0      29   3.217   1.736   1.956   2.581   3.325   3.610 
    .90     .95 
  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正态检验

> 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( )函数，它可以计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度、平均值的标准差

> library(psych)
 
载入程辑包：‘psych’
 
下列对象被屏蔽了from 
‘package:boot’:
 
    logit
 
下列对象被屏蔽了from 
‘package:Hmisc’:
 
    describe
 
下列对象被屏蔽了from 
‘package:ggplot2’:
 
    %+%
 
> 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

再是分组计算描述性统计量

> 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

这里用到了aggregate( )函数，注意list(am=mtcars$am)，如果使用list(mtcars$am)，则am列将被标注为Group.1而不是am；如果想标注更多列，by=list(name1=groupvar1, name2=groupvar2, ... , groupvarN)这样的语句

遗憾的是，aggregate( )仅允许每次调用中使用平均数、标准差这样的单返回值函数，它无法一次返回若干统计量，要完成这项任务，可以使用by( )函数，格式为：by(data, INDICES, FUN)

其中data是一个数据框或矩阵，INDICES是一个因子或因子组成的列表，定义了分组，FUN是任意函数

（书中的例子在R 上不知道怎么运行也都是报错，所以忽略）

还是说一下扩展把

doBy包和psych包提供了分组计算描述性统计量的函数

doBy包中summaryBy( )函数使用格式：summaryBy(formula, data=dataframe, FUN=function)

其中formula接受以下格式：var1+var2+var3+ ... +varN ~ groupvar1+groupvar2+ ... + groupvarN

~左侧的变量是需要分析的数值型函数，而右侧的变量是类别型的分组变量

> library(doBy)
载入需要的程辑包：MASS
> 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包describe.by( )函数可计算和describe相同的描述性统计量，只是按照一个或多个分组变量分层

> library(psych)
> describe.by(mtcars[vars], mtcars$am)
group: 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
--------------------------------------------------------- 
group: 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
警告信息：
describe.by 
is deprecated.  Please use the describeBy function 

reshape包也可以灵活地按组导出描述性统计量

> library(reshape)
 
载入程辑包：‘reshape’
 
下列对象被屏蔽了from 
‘package:plyr’:
 
    rename, round_any
 
下列对象被屏蔽了from 
‘package:reshape2’:
 
    colsplit, melt, recast
> dstats <- function(x)(c(n=length(x), mean=mean(x), sd=sd(x)))
> dfm <- melt(mtcars, measure.vars=c("mpg", "hp", "wt"),
+        id.vars=c("am", "cyl"))
> cast(dfm, am+cyl+variable ~ ., dstats)
   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

当然也许这种方式最为简洁动人

再是结构的可视化

对于定量变量，有直方图、密度图、箱线图、点图等

然后到了频数表和列联表

> library(vcd)
> 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

先是生成频数表

用于创建和处理列联表的函数有以下：

table(var1, var2, ... , varN) -- 使用N个类别型变量（因子）创建一个N维列联表

xtabs(formula, data) -- 根据一个公式和一个矩阵或数据框创建一个N维列联表

prop.table(table, margins) -- 依margins定义的边际列表将表中条目表示为分数形式

margin.table(table, margins) -- 依margins定义的边际列表计算表中条目和

addmargins(table, margins) -- 将概述为margins(默认是求和结果)放入表中

ftable(table) -- 创建一个紧凑的“平铺”式列联表

1、一维列联表

> mytable <- with(Arthritis, table(Improved))
> 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)

xtabs( )函数格式为：mytable <- xtabs(~ A + B, data=mydata)

其中A为行变量，B为列变量

> mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
> mytable
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

> 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

> 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

> 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(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

gmodels包中的CrossTable( )函数是创建二维列联表的第三种方法，CrossTable( )函数仿照SAS中的PROC FREQ或SPSS中的CROSSTABS的形式生成二维列联表

> 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( )函数功能强大，详见help(CrossTable)

接着是多位列联表了啊

> 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、卡方独立性检验

> 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
 
警告信息：
In chisq.test(mytable) : Chi-squared近似算法有可能不准

2、Fisher精确检验

> 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

3、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.6323, 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   : 0.394 
Contingency Coeff.: 0.367 
Cramer‘s V        : 0.394 

vcd包中的assocstats( )函数可以计算二维列联表的phi系数、列联系数和Cramer‘s V系数

另外vcd包中的kappa( )函数可以计算混淆矩阵的Cohen‘s kappa值以及加权kappa值（混淆矩阵可以表示两位评判者对于一系列对象进行分类所得到结构的一致程度）

再是结构的可视化

比如R 中的vcd包、ca包中的函数

然后是将表转化为扁平格式

（但是我觉得这部分有点偏，所以省略啊）

相关

1、Pearson、Spearman、Kendall相关

Pearson积差相关系数衡量了两个定量变量之间的线性相关程度；Spearman等级相关系数则衡量分级定序变量之间的相关程度；Kendall‘s Tau相关系数也是一种非参数的等级相关度量

cor( )函数可以计算这三种相关系数，cov( )函数可用来计算协方差；两个函数的参数很多，简化为cor(x, use= , method= )

其中x为数据框或矩阵；use为指定缺失数据的处理方式，可选的方式为all.obs(假设不存在缺失数据，遇到缺失值报错)、 everything(遇到缺失值时，相关系数的计算结构将被设为missing)、complete.obs(行删除)以及pairwise.complete.obs(成对删除，pairwise deletion)

默认参数为use="everything", method="pearson"

> 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  1.00000000 -0.7808458
Murder      0.34364275 -0.2300776  0.7029752 -0.78084575  1.0000000
HS Grad    -0.09848975  0.6199323 -0.6571886  0.58221620 -0.4879710
               HS Grad
Population -0.09848975
Income      0.61993232
Illiteracy -0.65718861
Life Exp    0.58221620
Murder     -0.48797102
HS Grad     1.00000000
> 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

> x <- states[,c("Population", "Income", "Illiteracy", "HS Grad")]
> y <- states[,c("Life Exp", "Murder")]
> cor(x,y)
              Life Exp     Murder
Population -0.06805195  0.3436428
Income      0.34025534 -0.2300776
Illiteracy -0.58847793  0.7029752
HS Grad     0.58221620 -0.4879710

2、偏相关

偏相关是指控制一个或多个定量变量时，另外两个定量变量之间的相互关系；可以用到ggm包中的pcor( )函数计算偏相关系数，函数调用格式为pcor(u, s)

> library(ggm)
载入需要的程辑包：igraph
 
载入程辑包：‘ggm’
 
下列对象被屏蔽了from 
‘package:Hmisc’:
 
    rcorr
 
> #在控制了收入、文盲率和高中毕业率时
> #人口和谋杀率的偏相关系数
> pcor(c(1,5,2,3,6),cov(states))
[1] 0.3462724

3、其他类型的相关

polycor包中的hetcor( )函数可以计算一种混合的相关矩阵，详见help(polycor)

再是相关性的显著性检验

可以使用cor.test( )函数对Pearson、Spearman、Kendall相关系数进行检验，简化后的使用格式为：cor.test(x, y, alternative= , method= )

其中x和y为要检验相关性的变量，alternative则用来指定进行双侧检验或单侧检验（取值为"two.side","less","greater"），而method用以指定要计算的相关类型（"pearson","kendall","spearman"）；当研究的假设为总体的相关系数小于0时，请使用alternative="less", 大于0，="greater", ="two.side"为总体相关系数不等于0

> cor.test(states[,3], states[,5])
 
        Pearson‘s product-moment correlation
 
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 

cor.test每次只能检验一种相关关系，但是psych包中提供的corr.test( )函数可以一次做更多的事情；corr.test( )函数可以为Pearson、Spearman、Kendall相关计算相关矩阵和显著性水平

> library(psych)
> corr.test(states, use="complete")
Call:corr.test(x = 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
 
 To see confidence intervals of the correlations, print with the short=FALSE option

参数use=的取值可为‘pairwise"或"complete"(分别表示对缺失值执行成对删除或行删除)；method=取值可为"pearson"(默认值)、"spearman"或"kendall"

其他显著性检验

在多元正态性的假设下，psych包中的pcor.test( )函数可以用来检验在控制一个或多个额外变量时两个变量之间的条件独立性，使用格式为：pcor.test(r, q, n)

其中r是由pcor( )函数计算得到的偏相关系数，q为要控制的变量数（以数值表示位置），n为样本大小

此外，psych包中的r.test( )函数提供了多种实用的显著性检验方法，此函数可用来检验：

某种相关系数的显著性；两个独立相关系数的差异是否显著；两个基于一个共享变量得到的非独立相关系数的差异是否显著；两个基于完全不同的变量得到的非独立相关系数的差异是否显著

最后是关于相关系数的可视化

主要是散点图、散点图矩阵、相关图等

t检验

一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设，这里假设两组数据是独立的，并且是从正态总体中抽得，检验的调用格式为：t.test(y ~ x, data)

其中y是一个数值变量，x是一个二分变量，调用格式为：t.test(y1, y2)

其中y1和y2为数值型向量（即各组的结果变量），可选参数data的取值为一个包含了这些变量的矩阵或数据框

在R 中t检验默认假定方差不相等，并使用Welsh的修正自由度（其实我不懂这句话）；但是可以通过添加一个参数var.equal=TRUE以假定方差相等，并使用合并方差估计，默认的备择假设是双侧的（即均值不等，但大小的方向不确定）；但是可以添加一个参数alternative="less‘/"greater"来进行有方向的检验

> 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 

非独立样本的检验

非独立样本的t检验假定组间的差异呈正态分布，检验调用的格式为：t.test(y1, y2, paired=TRUE)

其中y1和y2都是两个非独立组的数值向量

> library(MASS)
> sapply(UScrime[c("U1","U2")], function(x)(c(mean=mean(x), sd=sd(x))))
           U1       U2
mean 95.46809 33.97872
sd   18.02878  8.44545
> with(UScrime, t.test(U1, U2, paired=TRUE))
 
        Paired t-test
 
data:  U1 and U2
t = 32.4066, df = 46, p-value < 2.2e-16
alternative hypothesis: true 
difference in 
means is 
not equal to 0
95 percent confidence interval:
 57.67003 65.30870
sample estimates:
mean of the differences 
               61.48936 

对于多于两组的情况

如果能够假设数据是从正太总体中独立抽样而得的，那么可以使用方差分析（ANOVA）

组间差异的非参数检验

如果数据无法满足t检验或ANOVA检验的参数假设的话，可以转而使用非参数方法；举例来说，若结果变量在本质上就严重偏倚或呈现有序关系，就可以用了啊

两组的比较

若两组数据独立，可以使用Wilcoxon秩和检验（即Mann-Whitney U检验）来评估观测是否从相同的概率分布中抽得（即，在一个总体中获得更高得分的概率是否比另一个总体要大），调用格式为：wilcox.test(y ~ x, data)

其中y是数值型变量，x是一个二分变量，调用格式为：wilcox.test(y1, y2)

其中y1和y2为各组的结果变量，默认进行一个双侧检验，当然可以添加参数exact来进行精确检验，指定alternative="less"/"greater"进行有方向的检验

>  with(UScrime, by(Prob, So, median))
So: 0
[1] 0.038201
--------------------------------------------------------- 
So: 1
[1] 0.055552
> wilcox.test(Prob ~ So, data=UScrime)
 
        Wilcoxon rank sum test
 
data:  Prob by 
So
W = 81, p-value = 8.488e-05
alternative hypothesis: true 
location shift is 
not equal to 0

Wilcoxon符号秩检验是非独立样本t检验的一种非参数替代方法，适用于两组成对数据和无法保证正态性假设的情境，还可以添加参数paired=TRUE

> sapply(UScrime[c("U1","U2")],median)
U1 U2 
92 34 
> with(UScrime, wilcox.test(U1, U2, paired=TRUE))
 
        Wilcoxon signed rank test with continuity correction
 
data:  U1 and U2
V = 1128, p-value = 2.464e-09
alternative hypothesis: true 
location shift is 
not equal to 0
 
警告信息：
In wilcox.test.default(U1, U2, paired = TRUE) :
  cannot compute exact p-value with ties

多于两组的比较

如果无法满足ANOVA设计的假设，那么可以使用非参数方法来评估组间的差异，如果各组独立，则Kruskal-Wallis检验将是一种实用的方法，如果各组不独立（如重复测量设计或随机区组设计），那么Friedman检验会更合适

Kruskal-Wallis检验的调用格式为：kruskal.test(y ~ A, data)

其中y是一个数值型结果变量，A是一个拥有两个或更多水平的分组变量（若有两个水平，则它与Mann-Whitney U检验等价）

Friedman检验的调用格式为：friedman.test(y ~ A | B, data)

其中y是数值型结果变量，A是一个分组变量，而B是一个用以认定匹配观测的区组变量

此外，以上两检验中data为可选

> states <- as.data.frame(cbind(state.region, state.x77))
> kruskal.test(Illiteracy ~ state.region, data=states)
 
        Kruskal-Wallis rank sum test
 
data:  Illiteracy by 
state.region
Kruskal-Wallis chi-squared = 22.6723, df = 3, p-value = 4.726e-05

npmc包提供了所需要的非参数多组比较程序（一种是在控制犯第一类错误的概率的前提下，执行可以同步进行的多组比较，这样可以直接完成所有组之间的成对比较）

（但是很遗憾的是，我的电脑R 3.03和R 2.14都无法安装npmc包）

组间差异的可视化

这包括比如箱线图、核密度图啊

最后的最后，感谢作者Rober I. Kabacoff等，译者高涛、肖楠、陈钢

《R语言实战》学习笔记seventh,布布扣,bubuko.com

《R语言实战》学习笔记seventh

标签：des c class blog code a

原文地址：http://www.cnblogs.com/SSSSQD/p/3753841.html

踩

(0)

评论一句话评论（0）

分享档案

更多>

2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)

周排行