Part1 相關統計量說明

峰度係數Coefficientof kurtosis

峰度係數（Kurtosis）用來度量資料在中心聚集程度。在正態分佈情況下，峰度係數值是3（但是SPSS等軟體中將正態分佈峰度值定為0，是因為已經減去3，這樣比較起來方便）。>3的峰度係數說明觀察量更集中，有比正態分佈更短的尾部；<3的峰度係數說明觀測量不那麼集中，有比正態分佈更長的尾部，類似於矩形的均勻分佈。

偏度係數skew

SEMean 是 Standard error ofthe mean的縮寫，

標準誤差平均值，也叫平均數標準誤差，是描述均數抽樣分佈的離散程度及衡量均數抽樣誤差大小的尺度。SE Mean的計算公式如下：

均方誤差MeanSquared Error, MSE

數理統計中均方誤差是指引數估計值與引數真值之差平方的期望值，記為MSE。MSE是衡量“平均誤差”的一種較方便的方法，MSE可以評價資料的變化程度，MSE的值越小，說明預測模型描述實驗資料具有更好的精確度。與此相對應的，還有均方根誤差RMSE、平均絕對百分誤差等等。

標準偏差StdDev,Standard Deviation

標準偏差反映數值相對於平均值(mean) 的離散程度。

變異係數（Coefficientof Variation）

統計百科參考

Part2 R中基本統計函式實現

R語言中除本身有的獲取統計量的方法summary()之外，能得到描述性統計量的包有Hmisc、pastecs和psych。這裡使用的資料是R中已有車輛路試（mtcars）資料集，挑取其中的幾個欄位，英里數(mpg)、馬力(hp)和車重(wt)來做後續的示例資料集。

vars <- c(“mpg”,”hp”,”wt”)
summary(mtcars[vars])
#統計結果有最小值、最大值、平均值、上四分位數、下四分位數
 
library(Hmisc)
describe(mtcars[vars])
#統計結果有總數、缺失數、唯一值、平均值、各個分位數、最大值最小值五個
 
library(pastecs)
stat.desc(mtcars[vars])
#統計結果有總數、null數、NA數、最小、最大、差值、和、平均值、0.95置信區間均值、方差、標準差、變異係數

其方法
desc=FALSE時，基本統計量總數、null數、NA數、最小、最大、差值、和norm=TRUE時，多六個個正態分佈統計量，包括偏度和峰度（以及它們的統計顯著程度）和Shapiro–Wilk正態檢驗結果。

library(psych)
describe(mtcars[vars])
#缺失值的數量、平均數、標準差、中位數、截尾均值、絕對中位差、最小值、最大值、值域、偏度、峰度和平均值的標準誤差
 

注意：Hmisc包中describe()函式被新載入的包psych中的describe函式遮蔽，要使用Hmisc中的函式時，使用下面程式碼Hmisc::describe(mt)

獲取分組統計量

aggregate()

aggregate(mtcars[vars], by = list(am =mtcars$am), mean)
aggregate(mtcars[vars], by = list(am =mtcars$am), sd)

缺點：只能獲取平均值、標準差這樣的單值統計量。

by()可以自定義函式

dstats <- function(x){c(max=max(x),min=min(x))}
by(mtcars[vars], mtcars$am, dstats)

doBy包中的summaryBy()可以自定義函式

library(doBy)
summaryBy(mpg + hp + wt ~ am, data = mtcars,FUN = mystats)

psych包中的describe.by()函式不能自定義函式，就是describe()方法中的那些統計量。若存在一個以上的分組變數，你可以使用list(groupvar1, groupvar2, ... , groupvarN)。

library(psych)
describe.by(mtcars[vars], mtcars$am)

利用reshape包中的cast方法來自定義分組統計量的一個例子如下：

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

頻數表和列聯表

資料來自vcd包中的Arthritis資料集。這份資料來自Kock & Edward （1988），表示了一項風溼性關節炎新療法的雙盲臨床實驗的結果。

install.packages("vcd")
library(vcd)
head(Arthritis)

R中提供的建立頻數表和列聯表的方法有下面幾種：

一維列聯表

mytable <- with(Arthritis,table(Improved))#某一列的頻數統計表
mytable
prop.table(mytable)#根據頻數統計表得到百分比表
prop.table(mytable)*100

二維列聯表

mytable <- xtabs(~ Treatment+Improved,data=Arthritis)
mytable

margin.table(mytable, 1)#按行計數
prop.table(mytable, 1) #按行求百分比

margin.table(mytable, 2)#按列計數
prop.table(mytable, 2)#按列求百分比

prop.table(mytable)#行列合計求百分比
addmargins(mytable)#對計數結果列求和
addmargins(prop.table(mytable))#對百分比結果列求和
addmargins(prop.table(mytable, 1), 2)#對百分比結果按行求和
addmargins(prop.table(mytable, 2), 1)#對百分比結果按列求和

使用gmodels包中的CrossTable()函式是建立二維列聯表的第三種方法。CrossTable()

函式仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二維列聯表。CrossTable()函式有很多選項，可以做許多事情：計算（行、列、單元格）的百分比；指定小數位數；進行卡方、Fisher和McNemar獨立性檢驗；計算期望和（皮爾遜、標準化、調整的標準化）殘差；將缺失值作為一種有效值；進行行和列標題的標註；生成SAS或SPSS風格的輸出。可參閱help(CrossTable)以瞭解詳情。

install.packages("gmodels")
library(gmodels)
CrossTable(Arthritis$Treatment,Arthritis$Improved)

多維列聯表

mytable <- xtabs(~Treatment+Sex+Improved, data=Arthritis)#生成三維表
mytable
ftable(mytable)#以一種更緊湊的方式生成三維表
margin.table(mytable, 1)# 邊際頻數表
margin.table(mytable, 2)
margin.table(mytable, 3)
margin.table(mytable, c(1,3))# 治療情況（Treatment）×改善情況（Improved）的邊際頻數
ftable(prop.table(mytable, c(1, 2)))# 治療情況（Treatment）× 性別（Sex）的各類改善情況比例
ftable(addmargins(prop.table(mytable, c(1,2)), 3))

獨立性檢驗

卡方獨立性檢驗

chisq.test()函式對二維表的行變數和列變數進行卡方獨立性檢驗。

library(vcd)
mytable <- xtabs(~Treatment+Improved,data=Arthritis)#治療情況和改善情況
chisq.test(mytable)
mytable <- xtabs(~Improved+Sex,data=Arthritis)#性別和改善情況
chisq.test(mytable)

Fisher精確檢驗

fisher.test()函式可進行Fisher精確檢驗，也是針對二維列聯表進行的檢驗。Fisher精確檢驗的原假設是：邊界固定的列聯表中行和列是相互獨立的。

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

Cochran-Mantel—Haenszel檢驗

mantelhaen.test()函式可用來進行Cochran—Mantel—Haenszel卡方檢驗，其原假設是，兩

個名義變數在第三個變數的每一層中都是條件獨立的。

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

上面p-value=0.0006647，可以知道，患者接受的治療與得到的改善在性別的每一水平下並不獨立，即不管性別是什麼，治療情況和改善情況都存在關係。

相關性的度量

R提供了多種檢驗類別型變數獨立性的方法，這裡講到三種檢驗分別為卡方獨立性檢驗、Fisher精確檢驗和Cochran-Mantel–Haenszel檢驗。

vcd包中的assocstats()函式可以用來計算二維列聯表的phi係數、列聯絡數和Cramer’s V係數。

library(vcd)
mytable <- xtabs(~Treatment+Improved,data=Arthritis)
assocstats(mytable)

這些相關性係數的介紹及其計算方式不詳細說明了，提供一個參考：

Phi係數

Φ相關係數的大小，表示兩因素之間的關聯程度。當Φ值小於0.3時，表示相關較弱；當Φ值大於0.6時，表示相關較強。適用於2×2表即四格表。

列聯絡數

C係數的大小，表示兩因素之間的關聯程度。C係數可能的最大值依賴於列聯表的行列數，並隨行列增大而增大，只有行列數相同的列聯表比較其列聯絡數的大小。

cramer’s V係數

cramer’s V係數是Phi係數的修正值，適用於四格表，也適用於大於2×2表格，當行列有一個為二維時，V係數就等於Phi係數。

多維聯錶轉換為扁平表格

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)
}
treatment <- rep(c("Placebo","Treated"), 3)
improved <- rep(c("None","Some", "Marked"), each = 2)
Freq <- c(29, 13, 7, 7, 7, 21)
mytable <- as.data.frame(cbind(treatment,improved, Freq))
mydata <- table2flat(mytable)#
head(mydata)

相關

相關係數可以用來描述定量變數之間的關係。相關係數的符號（+、-）表明關係的方向（正相關或負相關），其值的大小表示關係的強弱程度（完全不相關時為0，完全相關時為1）。這一部分介紹多種相關係數和相關性的顯著性檢驗方法。使用的資料是R中測試資料state.x77資料集，它提供了美國50個州在1977年的人口、收入、文盲率、預期壽命、謀殺率和高中畢業率資料。

R可以計算多種相關係數，包括Pearson相關係數、Spearman相關係數、Kendall相關係數、偏相關係數、多分格（polychoric）相關係數和多系列（polyserial）相關係數。接下來依次看這些資料。

Pearson、Spearman和Kendall相關

cor()函式可以計算這三種相關係數，而cov()函式可用來計算協方差。這兩個方法的引數如下：

states <- state.x77[, 1:6]
cov(states)#計算states資料集的協方差
cor(states)# 計算states資料集的相關係數，預設pearson
cor(states, method="spearman")# 計算states資料集的spearman係數

我們能看到收入和高中畢業率之間存在很強的正相關，文盲率和預期壽命之間存在很強的負相關，文盲率和謀殺率存在很強的正相關。還可以自己定義想要計算相關性的行列。

x <- states[, c("Population","Income", "Illiteracy", "HS Grad")]
y <- states[, c("Life Exp","Murder")]
cor(x, y)

偏相關

偏相關是指在控制一個或多個定量變數時，另外兩個定量變數之間的相互關係。

install.packages("ggm")
library(ggm)
pcor(c(1, 5, 2, 3, 6), cov(states))
pcor(c(3, 5, 1, 2, 6), cov(states))

上面兩個結果一是在控制了收入、文盲率和高中畢業率的影響時，人口和謀殺率之間的相關係數為0.3462724。二是在控制了人口、收入和高中畢業率的影響時，文盲率和謀殺率之間的相關係數為0.5989741，這個值已經很大了，說明這兩者之間的關係很大。

相關性的顯著性檢驗

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時，應使用alternative="greater"。預設情況下，假設為alternative="two.side"（總體相關係數不等於0）。

cor.test(states[, 3], states[, 5])

cor.test()方法有一個不足就是一次只能檢驗一種相關關係的，psych包中的corr.test()函式可以一次檢驗更多變數之間的相關性值和顯著水平。其使用格式如下：

corr.test(x, y, use=, method=)

引數use=的取值可為"pairwise"或"complete"（分別表示對缺失值執行成對刪除或行刪除）。引數method=的取值可為"pearson"（預設值）、"spearman"或"kendall"。

library(psych)
corr.test(states, use ="complete")

t檢驗

這裡使用MASS包中的UScrime資料集。它包含了1960年美國47個州的刑罰制度對犯罪率影響的資訊。我們感興趣的結果變數為Prob（監禁的概率）、U1（14~24歲年齡段城市男性失業率）和U2（35~39歲年齡段城市男性失業率）。類別型變數So（指示該州是否位於南方的指示變數）將作為分組變數使用。

獨立樣本的t檢驗

library(MASS)
t.test(Prob ~ So, data=UScrime)

美國的南方犯罪，是否更有可能被判監禁？物件是南方和非南方各州，因變數為監禁的概率。結果是拒絕南方各州和非南方各州擁有相同監禁概率的假設（p < .001）。

非獨立樣本的t檢驗

library(MASS)
sapply(UScrime[c("U1","U2")], function(x) (c(mean = mean(x), sd = sd(x))))
with(UScrime, t.test(U1, U2, paired =TRUE))

較年輕（14~24歲）男性的失業率是否比年長（35~39歲）男性的失業率更高？在這種情況下，這兩組資料並不獨立。差異的均值（61.5）足夠大，可以保證拒絕年長和年輕男性的平均失業率相同的假設。年輕男性的失業率更高。

參考：《R語言實戰》

有任何問題建議歡迎指正，轉載請註明來源，謝謝！