
? 描述性統(tǒng)計(jì)分析
? 頻數(shù)表和列聯(lián)表
? 相關(guān)系數(shù)和協(xié)方差
? t檢驗(yàn)
? 非參數(shù)統(tǒng)計(jì)
7.1描述性統(tǒng)計(jì)分析
> myvars<-c("mpg","hp","wt")
> head(mtcars[myvars])#返回mtcars數(shù)據(jù)集中的三類
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
7.1.1
- 通過 summary() 計(jì)算描述性統(tǒng)計(jì)量。
> myvars<-c("mpg","hp","wt")
> summary(mtcars[myvars])#獲取描述下統(tǒng)計(jì)量:
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
> #包括最小值、最大值、四分位數(shù)和數(shù)值型變量的均值,
> #以及因子向量和邏輯型向量的頻數(shù)統(tǒng)計(jì)。
- 通過 sapply() 計(jì)算描述性統(tǒng)計(jì)量。
sapply(x, FUN, options)其中的 x 是你的數(shù)據(jù)框(或矩陣), FUN 為一個(gè)任意的函數(shù)。
> options(digits = 3)
> mystats<-function(x,na.omit=FALSE){
+ if(na.omit)
+ x<-x[!is.na(x)]#獲得除過缺失值的數(shù)值
+ 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.000 32.000 32.0000
mean 20.091 146.688 3.2172
stdev 6.027 68.563 0.9785
skew 0.611 0.726 0.4231
kurtosis -0.373 -0.136 -0.0227
7.1.2
- 通過 Hmisc 包中的 describe() 函數(shù)計(jì)算描述性統(tǒng)計(jì)量
> library(Hmisc)
> myvars<-c("mpg","hp","wt")
> describe(mtcars[myvars])#變量和觀測量、缺失值和唯一值得數(shù)目
#平均值、分位數(shù)和五個(gè)最大、五個(gè)最小的值
mtcars[myvars]
3 Variables #變量 32 Observations#觀測量
----------------------------------------------------
mpg
n missing distinct Info Mean
32 0 25 0.999 20.09
Gmd .05 .10 .25 .50
6.796 12.00 14.34 15.43 19.20
.75 .90 .95
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
32 0 22 0.997 146.7
Gmd .05 .10 .25 .50
77.04 63.65 66.00 96.50 123.00
.75 .90 .95
180.00 243.50 253.55
lowest : 52 62 65 66 91, highest: 215 230 245 264 335
----------------------------------------------------
wt
n missing distinct Info Mean
32 0 29 0.999 3.217
Gmd .05 .10 .25 .50
1.089 1.736 1.956 2.581 3.325
.75 .90 .95
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() 函數(shù)計(jì)算描述性統(tǒng)計(jì)量。
> options(digits = 3)
> library(pastecs)
> myvars<-c("mpg","hp","wt")
> stat.desc(mtcars[myvars])
mpg hp wt
nbr.val#所有值 32.00 32.000 32.000
nbr.null#空值 0.00 0.000 0.000
nbr.na#缺失值得數(shù)量 0.00 0.000 0.000
min#最小值 10.40 52.000 1.513
max#最大值 33.90 335.000 5.424
range#值域 23.50 283.000 3.911
sum#總和 642.90 4694.000 102.952
median#中位數(shù) 19.20 123.000 3.325
mean#平均數(shù) 20.09 146.688 3.217
SE.mean#平均數(shù)的標(biāo)準(zhǔn)誤1.07 12.120 0.173
CI.mean.0.95平均數(shù)置信度95%的置信區(qū)間
2.17 24.720 0.353
var#方差 36.32 4700.867 0.957
std.dev#標(biāo)準(zhǔn)差 6.03 68.563 0.978
coef.var#變異系數(shù) 0.30 0.467 0.304
- 通過 psych 包中的 describe() 計(jì)算描述性統(tǒng)計(jì)量
> library(psych)
> myvars<-c("mpg","hp","wt")
> describe(mtcars[myvars])
vars n mean sd median trimmed mad min
mpg 1 32 20.09 6.03 19.20 19.70 5.41 10.40
hp 2 32 146.69 68.56 123.00 141.19 77.10 52.00
wt 3 32 3.22 0.98 3.33 3.15 0.77 1.51
max range skew kurtosis se
mpg 33.90 23.50 0.61 -0.37 1.07
hp 335.00 283.00 0.73 -0.14 12.12
wt 5.42 3.91 0.42 -0.02 0.17
統(tǒng)計(jì)量結(jié)果包括:非缺失值的數(shù)量、平均數(shù)、標(biāo)準(zhǔn)差、中位數(shù)、截尾均值、絕對(duì)中位差、最小值、最大值、值域、偏度、峰度和平均值得標(biāo)準(zhǔn)誤。(忒豐富了,有點(diǎn)變態(tài)哈)
7.1.3 分組計(jì)算描述性統(tǒng)計(jì)量
- 使用 aggregate() 分組獲取描述性統(tǒng)計(jì)量
> myvars<-c("mpg","hp","wt")
> aggregate(mtcars[myvars],by=list(am=mtcars$am),mean)#返回平均值
am mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(am=mtcars$am),sd)#返回標(biāo)準(zhǔn)差
am mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
其中l(wèi)ist(am=mtcars$am)的作用是為第一列(am)指定一個(gè)標(biāo)簽,若是無指定,則是如下這樣:
> myvars<-c("mpg","hp","wt")
> aggregate(mtcars[myvars],by=list(mtcars$am),mean)
Group.1 mpg hp wt
1 0 17.1 160 3.77
2 1 24.4 127 2.41
> aggregate(mtcars[myvars],by=list(mtcars$am),sd)
Group.1 mpg hp wt
1 0 3.83 53.9 0.777
2 1 6.17 84.1 0.617
- 使用 by() 分組計(jì)算描述性統(tǒng)計(jì)量
by(data, INDICES, FUN),其中data是一個(gè)數(shù)據(jù)框或者矩陣,INDICSE是一個(gè)因子或因子組成的列表,定義了分組。FUN是任意函數(shù)。
> 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))
+ }
> 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.000 19.0000 19.000
mean 17.147 160.2632 3.769
stdev 3.834 53.9082 0.777
skew 0.014 -0.0142 0.976
kurtosis -0.803 -1.2097 0.142
---------------------------------------------------------------------
mtcars$am: 1
mpg hp wt
n 13.0000 13.000 13.000
mean 24.3923 126.846 2.411
stdev 6.1665 84.062 0.617
skew 0.0526 1.360 0.210
kurtosis -1.4554 0.563 -1.174
7.1.4 分組計(jì)算的擴(kuò)展
- 使用 doBy 包中的 summaryBy() 分組計(jì)算概述統(tǒng)計(jì)量
> 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.1 3.83 0.0140 -0.803 19 160 53.9 -0.0142
2 1 13 24.4 6.17 0.0526 -1.455 13 127 84.1 1.3599
hp.kurtosis wt.n wt.mean wt.stdev wt.skew wt.kurtosis
1 -1.210 19 3.77 0.777 0.976 0.142
2 0.563 13 2.41 0.617 0.210 -1.174
其中左側(cè)的變量是需要分析的數(shù)值型變量(mpg+hp+wt),右側(cè)(am)是類別型的分組變量。
- 使用 psych 包中的 describeBy() 分組計(jì)算概述統(tǒng)計(jì)量
> 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
其中describeBy() 函數(shù)不允許指定任意函數(shù)。
7.2 頻數(shù)表和列聯(lián)表(主要是類比型變量,而非描述性統(tǒng)計(jì)量)
> library(vcd)
載入需要的程輯包:grid
Warning message:
程輯包‘vcd’是用R版本3.3.3 來建造的
> 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
本節(jié)中的數(shù)據(jù)來自 vcd 包中的 Arthritis 數(shù)據(jù)集,上例為前幾個(gè)觀測。
7.2.1 生成頻數(shù)表
table(var1, var2, ..., varN)使用 N 個(gè)類別型變量(因子)創(chuàng)建一個(gè) N 維列聯(lián)表
xtabs(formula, data)根據(jù)一個(gè)公式和一個(gè)矩陣或數(shù)據(jù)框創(chuàng)建一個(gè) N 維列聯(lián)表
prop.table(table, margins)依 margins 定義的邊際列表將表中條目表示為分?jǐn)?shù)形式
margin.table(table, margins)依 margins 定義的邊際列表計(jì)算表中條目的和
addmargins(table, margins)將概述邊 margins (默認(rèn)是求和結(jié)果)放入表中
ftable(table)創(chuàng)建一個(gè)緊湊的“平鋪”式列聯(lián)表
1. 一維列聯(lián)表
> #一維列聯(lián)表
> mytable<-with(Arthritis,table(Improved))
> mytable
Improved
None Some Marked
42 14 28
轉(zhuǎn)化為比例值和百分比:
> prop.table(mytable)
Improved
None Some Marked
0.500 0.167 0.333
> prop.table(mytable)*100
Improved
None Some Marked
50.0 16.7 33.3
2.二維列聯(lián)表
對(duì)于二維列聯(lián)表,table()函數(shù)的使用格式為:mytable<-table(A,B),其中的 A 是行變量, B 是列變量。
使用mytable <- xtabs(~ A + B, data=mydata)時(shí),mydata是一個(gè)矩陣或數(shù)據(jù)框。
> mytable<-xtabs(~Treatment+Improved,data=Arthritis)
> mytable#對(duì)于Arthitis數(shù)據(jù)集,創(chuàng)建一個(gè)自己想要內(nèi)容的表,然后再利用其它函數(shù)進(jìn)行處理
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
生成邊際頻數(shù)和比例。
>margin.table(mytable,1)#下標(biāo)1指代table()語句中的第一個(gè)變量Treatment
Treatment
Placebo Treated
43 41
> prop.table(mytable,1)# 下標(biāo) 1 指代 table() 語句中的第一個(gè)變量Treatment
Improved
Treatment None Some Marked
Placebo 0.674 0.163 0.163
Treated 0.317 0.171 0.512
> margin.table(mytable,2)
Improved
None Some Marked
42 14 28
> prop.table(mytable,2)#下標(biāo)2指的是table中的Improved
Improved
Treatment None Some Marked
Placebo 0.69 0.50 0.25
Treated 0.31 0.50 0.75
各單元所占比例如下:
> prop.table(mytable)
Improved
Treatment None Some Marked
Placebo 0.3452 0.0833 0.0833
Treated 0.1548 0.0833 0.2500
使用 addmargins() 函數(shù)為這些表格添加邊際和。默認(rèn)情況是為表中所有變量創(chuàng)建邊際和。
> 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.3452 0.0833 0.0833 0.5119
Treated 0.1548 0.0833 0.2500 0.4881
Sum 0.5000 0.1667 0.3333 1.0000
僅添加各行/各列的和:
> addmargins(prop.table(mytable,1),2)
Improved
Treatment None Some Marked Sum
Placebo 0.674 0.163 0.163 1.000
Treated 0.317 0.171 0.512 1.000
> addmargins(prop.table(mytable,2),1)
Improved
Treatment None Some Marked
Placebo 0.69 0.50 0.25
Treated 0.31 0.50 0.75
Sum 1.00 1.00 1.00
使用 CrossTable 生成二維列聯(lián)表
> 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 | |
--------------------|-----------|-----------|-----------|-----------|
3.多維列聯(lián)表
table()和xtabs()都可以基于三個(gè)或更多的類別型變量生成多維列聯(lián)表。
margin.table()、prop.table()和addmargins()函數(shù)可以自然地推廣到高于二維的情況。
> mytable<-xtabs(~Treatment+Sex+Improved,data = Arthritis)
> mytable#生成各單元格的頻數(shù)
, , 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
邊際頻數(shù)。
> margin.table(mytable,1)#對(duì)使用公式 ~Treatement+Sex+Improve 創(chuàng)建的表,
#Treatment需要下標(biāo)1來引用
Treatment
Placebo Treated
43 41
> margin.table(mytable,2)#Sex需要下標(biāo)2來引用
Sex
Female Male
59 25
> margin.table(mytable,3)#Improved需要下標(biāo)3來引用
Improved
None Some Marked
42 14 28
治療情況( Treatment ) × 改善情況( Improved )的邊際頻數(shù)
> margin.table(mytable,c(1,3))
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
治療情況( Treatment ) × 性別( Sex )的各類改善情況比例
> ftable(prop.table(mytable,c(1,2)))
Improved None Some Marked
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875
Male 0.9091 0.0000 0.0909
Treated Female 0.2222 0.1852 0.5926
Male 0.5000 0.1429 0.3571
> ftable(addmargins(prop.table(mytable,c(1,2)),3))#以較為緊湊的方式輸出多維列聯(lián)表,
#并添加來邊際和。
Improved None Some Marked Sum
Treatment Sex
Placebo Female 0.5938 0.2188 0.1875 1.0000
Male 0.9091 0.0000 0.0909 1.0000
Treated Female 0.2222 0.1852 0.5926 1.0000
Male 0.5000 0.1429 0.3571 1.0000
7.2.2 獨(dú)立性檢驗(yàn)(對(duì)生成列聯(lián)表中的變量是否相關(guān)或獨(dú)立進(jìn)行檢驗(yàn))
** 1. 卡方獨(dú)立性檢驗(yàn)**(卡方檢驗(yàn)就是統(tǒng)計(jì)樣本的實(shí)際觀測值與理論推斷值之間的偏離程度,實(shí)際觀測值與理論推斷值之間的偏離程度就決定卡方值的大小,卡方值越大,越不符合;卡方值越小,偏差越小,越趨于符合,若兩個(gè)值完全相等時(shí),卡方值就為0,表明理論值完全符合。)
> #卡方獨(dú)立性檢驗(yàn)
> library(vcd)
> mytable<-xtabs(~Treatment+Improved,data = Arthritis)
> chisq.test(mytable)#對(duì)二維列表的行變量和列變量進(jìn)行檢驗(yàn)
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
> #本例中對(duì)治療情況和改善情況進(jìn)行檢驗(yàn),P<0.01說明存在一定的相關(guān)性。
> 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近似算法有可能不準(zhǔn)(信息量太少,且有一個(gè)小于5的值)
由于P=0.089>0.01,故說明治療結(jié)果和性別之間是不獨(dú)立的。
2. Fisher精確檢驗(yàn)(邊界固定的列聯(lián)表中行和列是相互獨(dú)立的)
> 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#替代假設(shè)
3. Cochran-Mantel-Haenszel檢驗(yàn)(兩個(gè)名義變量在第三個(gè)變量的每一層中都是條件獨(dú)立的)
> #兩個(gè)名義變量在第三個(gè)變量的每一層中都是條件獨(dú)立的
> 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 相關(guān)性的度量(數(shù)值越大表明相關(guān)性越強(qiáng))
二維列聯(lián)表的相關(guān)性度量(很是抽象)
> #二維列聯(lián)表的相關(guān)性度量
> 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
7.3 相關(guān)
7.3.1 相關(guān)的類型
1. Pearson、Spearman和Kendall相關(guān)
使用的函數(shù)為:cor(x, use= , method= ),以及cov()函數(shù)。
協(xié)方差和相關(guān)系數(shù)
> #協(xié)方差和相關(guān)系數(shù)
> options(digits = 3)
> states<-state.x77[,1:6]#使用R基礎(chǔ)安裝中的數(shù)據(jù)集,提供了大量的關(guān)于人口、收入、
#文盲率、預(yù)期壽命、
#謀殺率和高中畢業(yè)率數(shù)據(jù)。
> cov(states)#計(jì)算樣本協(xié)方差(所有變量之間兩兩計(jì)算相關(guān))
Population Income Illiteracy Life Exp Murder HS Grad
Population 19931684 571230 292.868 -407.842 5663.52 -3551.51
Income 571230 377573 -163.702 280.663 -521.89 3076.77
Illiteracy 293 -164 0.372 -0.482 1.58 -3.24
Life Exp -408 281 -0.482 1.802 -3.87 6.31
Murder 5664 -522 1.582 -3.869 13.63 -14.55
HS Grad -3552 3077 -3.235 6.313 -14.55 65.24
> cor(states)# 計(jì)算了Pearson積差相關(guān)系數(shù)
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.0000 0.208 0.108 -0.0681 0.344 -0.0985
Income 0.2082 1.000 -0.437 0.3403 -0.230 0.6199
Illiteracy 0.1076 -0.437 1.000 -0.5885 0.703 -0.6572
Life Exp -0.0681 0.340 -0.588 1.0000 -0.781 0.5822
Murder 0.3436 -0.230 0.703 -0.7808 1.000 -0.4880
HS Grad -0.0985 0.620 -0.657 0.5822 -0.488 1.0000
> cor(states,method = "spearman")# 計(jì)算了Spearman等級(jí)相關(guān)系數(shù)
Population Income Illiteracy Life Exp Murder HS Grad
Population 1.000 0.125 0.313 -0.104 0.346 -0.383
Income 0.125 1.000 -0.315 0.324 -0.217 0.510
Illiteracy 0.313 -0.315 1.000 -0.555 0.672 -0.655
Life Exp -0.104 0.324 -0.555 1.000 -0.780 0.524
Murder 0.346 -0.217 0.672 -0.780 1.000 -0.437
HS Grad -0.383 0.510 -0.655 0.524 -0.437 1.000
從表中可以看到收入和高中畢業(yè)率之間存在很強(qiáng)的正相關(guān)。Pearson相關(guān)關(guān)系為0.6199,Spearman相關(guān)關(guān)系為0.510,兩者均為極強(qiáng)的正相關(guān),但預(yù)期壽命和文盲率存在很強(qiáng)的負(fù)相關(guān)。Pearson相關(guān)關(guān)系為-0.588,Spearman相關(guān)關(guān)系為-0.555。
2. 偏相關(guān)(偏相關(guān)是指在控制一個(gè)或多個(gè)定量變量時(shí),另外兩個(gè)定量變量之間的相互關(guān)系)
使用pcor(u, S)函數(shù)。
> library(ggm)
> colnames(states)
[1] "Population" "Income" "Illiteracy" "Life Exp" "Murder" "HS Grad"
> pcor(c(1,5,2,3,6),cov(states))#前兩個(gè)1,5表示要計(jì)算相關(guān)系數(shù)
#的變量下標(biāo),其余數(shù)值為條件變量的下標(biāo)。
[1] 0.346#人口和謀殺率之間的相關(guān)系數(shù)0.346
7.3.2 相關(guān)性的顯著性檢驗(yàn)
使用cor.test(x, y, alternative = , method = )函數(shù),其中的 x 和 y 為要檢驗(yàn)相關(guān)性的變量, alternative 則用來指定進(jìn)行雙側(cè)檢驗(yàn)或單側(cè)檢驗(yàn)(取值為 "two.side" 、 "less" 或 "greater" ),而 method 用以指定要計(jì)算的相關(guān)類型( "pearson" 、"kendall" 或 "spearman" )。
- 檢驗(yàn)?zāi)撤N相關(guān)系數(shù)的顯著性
> #檢驗(yàn)?zāi)撤N相關(guān)系數(shù)的顯著性
> cor.test(states[,3],states[,5])
Pearson's product-moment correlation
data: states[, 3] and states[, 5]
t = 7, df = 50, p-value = 1e-08
alternative hypothesis: true correlation is not equal to 0#替代假設(shè),
95 percent confidence interval:#95%的置信區(qū)間
0.528 0.821
sample estimates:#樣本估計(jì)
cor
0.703
- 通過 corr.test 計(jì)算相關(guān)矩陣并進(jìn)行顯著性檢驗(yàn)
> 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