# 第七章 基本統計分析 # 本章內容 # 描述性統計分析 # 頻數表和列聯表 # 相關系數和協方差 # t檢驗 # 非參數統計# 7.1 描述性統計分析 # 本節中，我們將關注分析連續型變量的中心趨勢、變化性和分布形狀的方法。為了便于說明， 我們將使用第1章中Motor Trend 雜志的車輛路試(mtcars)數據集。我們的關注焦點是每加侖 汽油行駛英里數(mpg)、馬力(hp)和車重(wt)。 vars <- c("mpg", "hp", "wt") head(mtcars[vars]) # 我們將首先查看所有32種車型的描述性統計量，然后按照變速箱類型(am)和汽缸數(cyl)考察描述性統計量。變速箱類型是一個以0表示自動擋、1表示手動擋來編碼的二分變量，而汽缸 數可為4、5或6。# 7.1.1 方法云集 # 對于基礎安裝，你可以使用summary()函數來獲取描述性統計量。 # 代碼清單7-1 通過summary()計算描述性統計量 summary(mtcars[vars]) # summary()函數提供了最小值、最大值、四分位數和數值型變量的均值，以及因子向量和邏 輯型向量的頻數統計。 # 你可以使用第5章中的apply()函數或sapply()函數計算所選擇的任意描 述性統計量。對于sapply()函數，其使用格式為: # sapply(x, function, options) # 其中的x是你的數據框(或矩陣)，FUN為一個任意的函數。如果指定了options，它們將被傳遞 給FUN。你可以在這里插入的典型函數有mean、sd、var、min、max、median、length、range 和quantile。函數fivenum()可返回圖基五數總括(Tukey’s five-number summary，即最小值、 下四分位數、中位數、上四分位數和最大值)。 # 代碼清單7-2 通過sapply()計算描述性統計量 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) / nkurt <- sum((x - m) ^ 4 / s ^ 4) / n - 3return(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt)) } # 調用上面的函數 sapply(mtcars[vars], mystats) # 對于樣本中的車型，每加侖汽油行駛英里數的平均值為20.1，標準差為6.0。分布呈現右偏(偏 度+0.61)，并且較正態分布稍平(峰度-0.37)。# 擴展 # 若干用戶貢獻包都提供了計算描述性統計量的函數，其中包括Hmisc、pastecs和psych。 # Hmisc包中的describe()函數可返回變量和觀測的數量、缺失值和唯一值的數目、平均值、 分位數，以及五個最大的值和五個最小的值。 library(Hmisc) describe(mtcars[vars]) # 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)。 install.packages("pastecs") # 代碼清單7-4 通過pastecs包中的stat.desc()函數計算描述性統計量 library(pastecs) stat.desc(mtcars[vars]) # psych包也擁有一個名為describe()的函數，它可以計算非缺失值的數量、 平均數、標準差、中位數、截尾均值、絕對中位差、最小值、最大值、值域、偏度、峰度和平均 值的標準誤。 install.packages("psych") library(psych) describe(mtcars[vars]) # 在前面的示例中，psych包和Hmisc包均提供了名為describe()的函數。R如何知道該 使用哪個呢?簡言之，如代碼清單7-5所示，最后載入的程序包優先。在這里，psych在 Hmisc之后被載入，然后顯示了一條信息，提示Hmisc包中的describe()函數被psych 包中的同名函數所屏蔽(masked)。鍵入describe()后，R在搜索這個函數時將首先找 到psych包中的函數并執行它。如果你想改而使用Hmisc包中的版本，可以鍵入 7 Hmisc::describe(mt)。# 7.1.2 分組計算描述性統計量 # 在比較多組個體或觀測時，關注的焦點經常是各組的描述性統計信息，而不是樣本整體的描 述性統計信息。 # 代碼清單7-6 使用aggregate()分組獲取描述性統計量 aggregate(mtcars[vars], by = list(am = mtcars$am), mean) aggregate(mtcars[vars], by = list(am = mtcars$am), sd) # 注意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是任 意函數。 # 代碼清單7-7 使用by()分組計算描述性統計量 dstats <- function(x)(c(meann = mean(x), sd = sd(x))) # 以下語句會報錯： (list) object cannot be coerced to type 'double' by(mtcars[vars], mtcars$am, dstats) # 使用以下語句解決上面的異常 by(as.numeric(unlist(mtcars["mpg"])), mtcars$am, dstats) by(as.numeric(unlist(mtcars["hp"])), mtcars$am, dstats) by(as.numeric(unlist(mtcars["wt"])), mtcars$am, dstats)# 擴展 # doBy包和psych包也提供了分組計算描述性統計量的函數。同樣地，它們未隨基本安裝發布， 必須在首次使用前進行安裝。doBy包中summaryBy()函數的使用格式為: # summaryBy(formula, data = dataframe, FUN = function) # 其中的formula接受以下的格式: # var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 +... + groupvarN # 在~左側的變量是需要分析的數值型變量，而右側的變量是類別型的分組變量。function 可為任何內建或用戶自編的R函數。 install.packages("doBy") # 代碼清單7-8 使用doBy包中的summaryBy()分組計算概述統計量 library(doBy) summaryBy(mpg + hp + wt ~ am, data = mtcars, FUN = mystats) # 與前面的示例不同，describe.by()函數不允許指定任意函數，所以它的普適性較低。若 存在一個以上的分組變量，你可以使用list(groupvar1, groupvar2, ... , groupvarN) 來表示它們。但這僅在分組變量交叉后不出現空白單元時有效。 # 最后，可以使用5.6.3節中描述的reshape包靈活地按組導出描述性統計量。首先，使用: # dfm <- melt(dataframe, measure.vars = y, id.vars = g) # 融合數據框。其中的dataframe包含著數據，y是一個向量，指明了要進行概述的數值型變量(默 認使用所有變量)，而g是由一個或多個分組變量組成的向量。然后使用: # cast(dfm, groupvar1, groupvar2 + ... + variable - ., FUN) # 重鑄數據。分組變量以+號分隔，這里的variable只取其字面含義(即僅表示重鑄后數據框中的變量variable)，而FUN是一個任意函數。 library(reshape) dstats <- function(x) (c(n = length(x), mean = mean(x), sd = sd(x))) dfm <- melt(mtcars, measure.vars = c("mpg", "hp", "wt")) cast(dfm, am + cyl + variable ~ ., dstats)# 7.2 頻數表和列聯表 # 在本節中，我們將著眼于類別型變量的頻數表和列聯表，以及相應的獨立性檢驗、相關性的 度量、圖形化展示結果的方法。我們除了使用基礎安裝中的函數，還將連帶使用vcd包和gmodels 包中的函數。 # 本節中的數據來自vcd包中的Arthritis數據集。這份數據來自Kock & Edward(1988)，表 示了一項風濕性關節炎新療法的雙盲臨床實驗的結果。前幾個觀測是這樣的: library(vcd) head(Arthritis) # 治療情況(安慰劑治療、用藥治療)、性別(男性、女性)和改善情況(無改善、一定程度的改善、顯著改善)均為類別型因子1。下一節中，我們將使用此數據創建頻數表和列聯表(交叉的分類)。# 7.2.1 生成頻數表 # 表7-1 用于創建和處理列聯表的函數 # 函數 描述 # 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. 一維列聯表 # 可以使用table()函數生成簡單的頻數統計表。示例如下: mytable <- with(Arthritis, table(Improved)) mytable # 可以用prop.table()將這些頻數轉化為比例值: prop.table(mytable) # 或使用prop.table()*100轉化為百分比: prop.table(mytable) * 100# 2. 二維列聯表 # 對于二維列聯表，table()函數的使用格式為: # mytable <- table(A, B) # 其中的A是行變量，B是列變量。除此之外，xtabs()函數還可使用公式風格的輸入創建列聯表， 格式為: # mytable <- xtabs(~ A + B, data = mydata) # 其中的mydata是一個矩陣或數據框。總的來說，要進行交叉分類的變量應出現在公式的右側(即~ 符號的右方)，以+作為分隔符。若某個變量寫在公式的左側，則其為一個頻數向量(在數據已經 被表格化時很有用)。 # 對于Arthritis數據，有: mytable <- xtabs(~ Treatment + Improved , data = Arthritis) mytable # 你可以使用margin.table()和prop.table()函數分別生成邊際頻數和比例。行和與行比 例可以這樣計算: # 行和 margin.table(mytable, 1) # 列和 margin.table(mytable, 2) # 行比 prop.table(mytable, 1) # 列比 prop.table(mytable, 2) # 以上四個語句中的下標1指代table()語句中的第一個變量，下標2指代table()語句中的第二個變量。觀察表格可以發現，與接受安慰劑的個體中有顯著改善的16%相比，接受治療的個體中的51%的個體病情有了顯著的改善。 # 各單元格所占比例可用如下語句獲取: prop.table(mytable) # 你可以使用addmargins()函數為這些表格添加邊際和。例如，以下代碼添加了各行的和與各列的和: addmargins(mytable) addmargins(prop.table(mytable)) # 在使用addmargins()時，默認行為是為表中所有的變量創建邊際和。作為對照: addmargins(prop.table(mytable, 1), 2) # 僅添加了各行的和。類似地， addmargins(prop.table(mytable, 2), 1) # 添加了各列的和。在表中可以看到，有顯著改善患者中的25%是接受安慰劑治療的。 # table()函數默認忽略缺失值(NA)。要在頻數統計中將NA視為一個有效的類別，請設定參數useNA="ifany" # 使用gmodels包中的CrossTable()函數是創建二維列聯表的第三種方法。CrossTable() 函數仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二維列聯表。示例參閱代碼清單 7-11。 # 代碼清單7-11 使用CrossTable生成二維列聯表 install.packages("gmodels") library(gmodels) CrossTable(Arthritis$Treatment, Arthritis$Improved) # CrossTable()函數有很多選項，可以做許多事情:計算(行、列、單元格)的百分比;指 定小數位數;進行卡方、Fisher和McNemar獨立性檢驗;計算期望和(皮爾遜、標準化、調整的 標準化)殘差;將缺失值作為一種有效值;進行行和列標題的標注;生成SAS或SPSS風格的輸出。 參閱help(CrossTable)以了解詳情。# 如果有兩個以上的類別型變量，那么你就是在處理多維列聯表。我們將在下面考慮這種情況。 # 3. 多維列聯表 # table()和xtabs()都可以基于三個或更多的類別型變量生成多維列聯表。margin.table()、 # prop.table()和addmargins()函數可以自然地推廣到高于二維的情況。另外，ftable()函 數可以以一種緊湊而吸引人的方式輸出多維列聯表。 # 代碼清單7-12 三維列聯表 # 各單元格的頻數 mytable <- xtabs(~ Treatment + Sex + Improved, data = Arthritis) mytable ftable(mytable) # 邊際頻數 margin.table(mytable, 1) margin.table(mytable, 2) margin.table(mytable, 3) # 治療情況(Treatment) x 改善情況(Improved)的邊際頻數 margin.table(mytable, c(1, 3)) # 治療情況(Treatment) x 性別(Sex)的邊際頻數 ftable(prop.table(mytable, c(1, 2))) # 為第三個下標添加了邊際和 ftable(addmargins(prop.table(mytable, c(1, 2)), 3)) # 如果想得到百分比而不是比例，可以將結果表格乘以100。例如: ftable(addmargins(prop.table(mytable, c(1, 2)), 3)) * 100# 7.2.2 獨立性檢驗 # 1. 卡方獨立性檢驗 # 你可以使用chisq.test()函數對二維表的行變量和列變量進行卡方獨立性檢驗。 # 代碼清單7-13 卡方獨立性檢驗 library(vcd) # (1)治療情況和改善情況不獨立 mytable <- xtabs(~ Treatment + Improved, data = Arthritis) chisq.test(mytable) # (2)性別和改善情況獨立 mytable <- xtabs(~ Treatment + Sex, data = Arthritis) chisq.test(mytable) # 在結果(1)中，患者接受的治療和改善的水平看上去存在著某種關系(p < 0.01)。而患者性別和改善情況之間卻不存在關系(p > 0.05)(2) # 這里的p值表示從總體中抽取的樣本行變量與列變量是相互獨立的概率。由于(1)的概率值很小，所以你拒絕了治療類型和治療結果相互獨立的原假 設。由于(2)的概率不夠小，故沒有足夠的理由說明治療結果和性別之間是不獨立的。代碼清單7-13 中產生警告信息的原因是，表中的6個單元格之一(男性 - 一定程度上的改善)有一個小于5的值， 這可能會使卡方近似無效。# 2. Fisher精確檢驗 # 可以使用fisher.test()函數進行Fisher精確檢驗。Fisher精確檢驗的原假設是:邊界固定 的列聯表中行和列是相互獨立的。其調用格式為fisher.test(mytable)，其中的mytable是 一個二維列聯表。示例如下: mytable <- xtabs(~ Treatment + Improved, data = Arthritis) fisher.test(mytable) # 與許多統計軟件不同的是，這里的fisher.test()函數可以在任意行列數大于等于2的二維 列聯表上使用，但不能用于2×2的列聯表。# 3. Cochran-Mantel—Haenszel檢驗 # mantelhaen.test()函數可用來進行Cochran—Mantel—Haenszel卡方檢驗，其原假設是，兩 個名義變量在第三個變量的每一層中都是條件獨立的。下列代碼可以檢驗治療情況和改善情況在 性別的每一水平下是否獨立。此檢驗假設不存在三階交互作用(治療情況×改善情況×性別)。 mytable <- xtabs(~ Treatment + Improved + Sex, data = Arthritis) mantelhaen.test(mytable) # 結果表明，患者接受的治療與得到的改善在性別的每一水平下并不獨立(即，分性別來看， 用藥治療的患者較接受安慰劑的患者有了更多的改善)。# 7.2.3 相關性的度量 # 上一節中的顯著性檢驗評估了是否存在充分的證據以拒絕變量間相互獨立的原假設。如果可以拒絕原假設，那么你的興趣就會自然而然地轉向用以衡量相關性強弱的相關性度量。vcd包中的assocstats()函數可以用來計算二維列聯表的phi系數、列聯系數和Cramer’s V系數。 # 代碼清單7-14 二維列聯表的相關性度量 library(vcd) mytable <- xtabs(~ Treatment + Improved, data = Arthritis) assocstats(mytable) # 總體來說，較大的值意味著較強的相關性。vcd包也提供了一個kappa()函數，可以計算混淆矩陣的Cohen’s kappa值以及加權的kappa值。(舉例來說，混淆矩陣可以表示兩位評判者對于一 系列對象進行分類所得結果的一致程度。)# 7.2.4 結果的可視化 # R中擁有遠遠超出其他多數統計軟件的、可視地探索類別型變量間關系的方法。通常，我們 會使用條形圖進行一維頻數的可視化(參見6.1節)。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 <- NULLfor (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) } # 調用 table2flat(mytable) # 代碼清單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)# 7.3 相關 # 相關系數可以用來描述定量變量之間的關系。相關系數的符號(?)表明關系的方向(正相 關或負相關)，其值的大小表示關系的強弱程度(完全不相關時為0，完全相關時為1)。 # 本節中，我們將關注多種相關系數和相關性的顯著性檢驗。我們將使用R基礎安裝中的 state.x77數據集，它提供了美國50個州在1977年的人口、收入、文盲率、預期壽命、謀殺率和 高中畢業率數據。 # 數據集中還收錄了氣溫和土地面積數據，但為了節約空間，這里將其丟棄。你 可以使用help(state.x77)了解數據集的更多信息。除了基礎安裝以外，我們還將使用psych 和ggm包。 # 7.3.1 相關的類型 # 1. Pearson、Spearman和Kendall相關 # Pearson積差相關系數衡量了兩個定量變量之間的線性相關程度。Spearman等級相關系數則衡 量分級定序變量之間的相關程度。Kendall’s Tau相關系數也是一種非參數的等級相關度量。 # cor()函數可以計算這三種相關系數，而cov()函數可用來計算協方差。兩個函數的參數有 很多，其中與相關系數的計算有關的參數可以簡化為: # cor(x, use= , method = ) # 表7-3 cor和cov的參數 # 參數 描述 # x 矩陣或數據框 # use 指定缺失數據的處理方式。可選的方式為all.obs(假設不存在缺失數據——遇到缺失數據時將報 # 錯)、everything(遇到缺失數據時，相關系數的計算結果將被設為missing)、complete.obs # (行刪除)以及 pairwise.complete.obs(成對刪除，pairwise deletion) # method 指定相關系數的類型。可選類型為pearson、spearman或kendall # 默認參數為use="everything"和method="pearson"。 # 代碼清單7-17 協方差和相關系數 states <- state.x77[, 1:6] cov(states) cor(states) cor(states, method = "spearman") # 以上首個語句計算了方差和協方差，第二個語句則計算了Pearson積差相關系數，而第三個語句計算 了Spearman等級相關系數。舉例來說，我們可以看到收入和高中畢業率之間存在很強的正相關， 而文盲率和預期壽命之間存在很強的負相關。 # 請注意，在默認情況下得到的結果是一個方陣(所有變量之間兩兩計算相關)。你同樣可以 計算非方形的相關矩陣。觀察以下示例: x <- states[, c("Population", "Income", "Illiteracy", "HS Grad")] y <- states[, c("Life Exp", "Murder")] cor(x, y) # 當你對某一組變量與另外一組變量之間的關系感興趣時，cor()函數的這種用法是非常實用 的。注意，上述結果并未指明相關系數是否顯著不為0(即，根據樣本數據是否有足夠的證據得 出總體相關系數不為0的結論)。由于這個原因，你需要對相關系數進行顯著性檢驗(在7.3.2節中 闡述)。# 2. 偏相關 # 偏相關是指在控制一個或多個定量變量時，另外兩個定量變量之間的相互關系。你可以使用 ggm包中的pcor()函數計算偏相關系數。ggm包沒有被默認安裝，在第一次使用之前需要先進行 安裝。函數調用格式為: # pcor(u, s) # 其中的u是一個數值向量，前兩個數值表示要計算相關系數的變量下標，其余的數值為條件變量 (即要排除影響的變量)的下標。S為變量的協方差陣。這個示例有助于闡明用法: install.packages("ggm") library(ggm) # 在控制了收入、文盲率和高中畢業率時人口和謀殺率的偏相關系數 pcor(c(1, 5, 2, 3, 6), cov(states)) # 本例中，在控制了收入、文盲率和高中畢業率的影響時，人口和謀殺率之間的相關系數為 0.346。偏相關系數常用于社會科學的研究中。# 3. 其他類型的相關 # polycor包中的hetcor()函數可以計算一種混合的相關矩陣，其中包括數值型變量的 Pearson積差相關系數、數值型變量和有序變量之間的多系列相關系數、 # 有序變量之間的多分格相 關系數以及二分變量之間的四分相關系數。多系列、多分格和四分相關系數都假設有序變量或二 分變量由潛在的正態分布導出。請參考此程序包所附文檔以了解更多。# 7.3.2 相關性的顯著性檢驗 # 在計算好相關系數以后，如何對它們進行統計顯著性檢驗呢?常用的原假設為變量間不相關 (即總體的相關系數為0)。你可以使用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)。參考代碼清單7-18中的示例。 # 代碼清單7-18 檢驗某種相關系數的顯著性 cor.test(states[, 3], states[, 5]) # 這段代碼檢驗了預期壽命和謀殺率的Pearson相關系數為0的原假設。假設總體的相關度為0， 則預計在一千萬次中只會有少于一次的機會見到0.703這樣大的樣本相關度(即p = 1.258e?08)。 6 由于這種情況幾乎不可能發生，所以你可以拒絕原假設，從而支持了要研究的猜想，即預期壽命 和謀殺率之間的總體相關度不為0。 # 遺憾的是，cor.test每次只能檢驗一種相關關系。但幸運的是，psych包中提供的 corr.test()函數可以一次做更多事情。corr.test()函數可以為Pearson、Spearman或Kendall 相關計算相關矩陣和顯著性水平。 # 代碼清單7-19 通過corr.test計算相關矩陣并進行顯著性檢驗 library(psych) corr.test(states, use = "complete") # 參數use=的取值可為"pairwise"或"complete"(分別表示對缺失值執行成對刪除或行刪 除)。參數method=的取值可為"pearson"(默認值)、"spearman"或"kendall"。這里可以看到，人口數量和高中畢業率的相關系數(?0.10)并不顯著地不為0(p = 0.5)。# 其他顯著性檢驗 # 在7.4.1節中，我們關注了偏相關系數。在多元正態性的假設下，psych包中的pcor.test() 函數(這里可能是作者的筆誤，函數pcor.test事實上包含于ggm包中)可以用來檢驗在控制一個或多個額外變量時兩個變量之間的條件獨立性。使用格式為: # poor.test(r, q, n) # 其中的r是由pcor()函數計算得到的偏相關系數，q為要控制的變量數(以數值表示位置)，n為樣本大小。 # 在結束這個話題之前應當指出的是，psych包中的r.test()函數提供了多種實用的顯著性 檢驗方法。此函數可用來檢驗: # 某種相關系數的顯著性; # 兩個獨立相關系數的差異是否顯著; # 兩個基于一個共享變量得到的非獨立相關系數的差異是否顯著; # 兩個基于完全不同的變量得到的非獨立相關系數的差異是否顯著。 # 參閱help(r.test)以了解詳情。# 7.3.3 相關關系的可視化 # 以相關系數表示的二元關系可以通過散點圖和散點圖矩陣進行可視化，而相關圖 (correlogram)則為以一種有意義的方式比較大量的相關系數提供了一種獨特而強大的方法。# 7.4 t檢驗 # 7.4.1 獨立樣本的t檢驗 # 如果你在美國的南方犯罪，是否更有可能被判監禁?我們比較的對象是南方和非南方各州，因變量為監禁的概率。一個針對兩組的獨立樣本t檢驗可以用于檢驗兩個總體的均值相等的假設。這里假設兩組數據是獨立的，并且是從正態總體中抽得。檢驗的調用格式為: # t.test(y ~ x, data) # 其中的y是一個數值型變量，x是一個二分變量。調用格式或為: # t.test(y1, y2) # 其中的y1和y2為數值型向量(即各組的結果變量)。可選參數data的取值為一個包含了這些 變量的矩陣或數據框。與其他多數統計軟件不同的是，這里的t檢驗默認假定方差不相等，并使 用Welsh的修正自由度。 # 你可以添加一個參數var.equal=TRUE以假定方差相等，并使用合并方 差估計。默認的備擇假設是雙側的(即均值不相等，但大小的方向不確定)。你可以添加一個參 數alternative="less"或alternative="greater"來進行有方向的檢驗。 # 在下列代碼中，我們使用了一個假設方差不等的雙側檢驗，比較了南方(group 1)和非南 方(group 0)各州的監禁概率: library(MASS) t.test(Prob ~ So, data = UScrime) # 你可以拒絕南方各州和非南方各州擁有相同監禁概率的假設(p < .001)。 # 由于結果變量是一個比例值，你可以在執行t檢驗之前嘗試對其進行正態化變換。在本例 中，所有對結果變量合適的變換(Y/1-Y、log(Y/1-Y)、arcsin(Y)、arcsin(sqrt(Y)) 都會將檢驗引向相同的結論。數據變換詳述于第8章。# 7.4.2 非獨立樣本的t檢驗 # 非獨立樣本的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)))) with(UScrime, t.test(U1, U2, paired = TRUE)) # 差異的均值(61.5)足夠大，可以保證拒絕年長和年輕男性的平均失業率相同的假設。 年輕男性的失業率更高。事實上，若總體均值相等，獲取一個差異如此大的樣本的概率小于 0.000 000 000 000 000 22(即2.2e?16)。# 7.4.3 多于兩組的情況 # 如果想在多于兩個的組之間進行比較，應該怎么做?如果能夠假設數據是從正態總體中獨立 抽樣而得的，那么你可以使用方差分析(ANOVA)。ANOVA是一套覆蓋了許多實驗設計和準實驗 設計的綜合方法。就這一點來說，它的內容值得單列一章。你可以隨時離開本節轉而閱讀第9章。# 7.5 組間差異的非參數檢驗 # 如果數據無法滿足t檢驗或ANOVA的參數假設，可以轉而使用非參數方法。舉例來說，若結果變量在本質上就嚴重偏倚或呈現有序關系，那么你可能會希望使用本節中的方法。 # 7.5.1 兩組的比較 # 若兩組數據獨立，可以使用Wilcoxon秩和檢驗(更廣為人知的名字是Mann–Whitney U檢驗) 來評估觀測是否是從相同的概率分布中抽得的(即，在一個總體中獲得更高得分的概率是否比另 一個總體要大)。調用格式為: # wilcox.test(y ~ x, data) # 其中的y是數值型變量，而x是一個二分變量。調用格式或為: # wilcox.test(y1, y2) # 其中的y1和y2為各組的結果變量。可選參數data的取值為一個包含了這些變量的矩陣或數據框。默 認進行一個雙側檢驗。你可以添加參數exact來進行精確檢驗，指定alternative="less"或 alternative="greater"進行有方向的檢驗。 # 如果你使用Mann–Whitney U檢驗回答上一節中關于監禁率的問題，將得到這些結果: with(UScrime, by(Prob, So, median)) wilcox.test(Prob ~ So, data = UScrime) # 你可以再次拒絕南方各州和非南方各州監禁率相同的假設(p < 0.001)。 # Wilcoxon符號秩檢驗是非獨立樣本t檢驗的一種非參數替代方法。它適用于兩組成對數據和 無法保證正態性假設的情境。調用格式與Mann–Whitney U檢驗完全相同，不過還可以添加參數 paired=TRUE。讓我們用它解答上一節中的失業率問題: sapply(UScrime[c("U1", "U2")], median) with(UScrime, wilcox.test(U1, U2, paired = TRUE)) # 你再次得到了與配對t檢驗相同的結論。 # 在本例中，含參的t檢驗和與其作用相同的非參數檢驗得到了相同的結論。當t檢驗的假設合 理時，參數檢驗的功效更強(更容易發現存在的差異)。而非參數檢驗在假設非常不合理時(如對于等級有序數據)更適用。# 7.5.2 多于兩組的比較 # 如果無法滿足ANOVA設計的假設，那么可以使用非參數方法來評估組間的差異。如果各組獨立，則Kruskal—Wallis檢驗將是一種實用的方法。如果各組不獨立(如重復測量設計或隨機區 組設計)，那么Friedman檢驗會更合適。Kruskal–Wallis檢驗的調用格式為: kruskal.test(y ~ A, data) # 其中的y是一個數值型結果變量，A是一個擁有兩個或更多水平的分組變量(grouping variable)。 (若有兩個水平，則它與Mann–Whitney U檢驗等價。)而Friedman檢驗的調用格式為: # friedman.test(y ~ A | B, data) # 其中的y是數值型結果變量，A是一個分組變量，而B是一個用以認定匹配觀測的區組變量(blocking variable)。在以上兩例中，data皆為可選參數，它指定了包含這些變量的矩陣或數據框。 # 讓我們利用Kruskal–Wallis檢驗回答文盲率的問題。首先，你必須將地區的名稱添加到數據 集中。這些信息包含在隨R基礎安裝分發的state.region數據集中。 states <- as.data.frame(cbind(state.region, state.x77)) # 現在就可以進行檢驗了: kruskal.test(Illiteracy ~ state.region, data = states) # 顯著性檢驗的結果意味著美國四個地區的文盲率各不相同(p <0.001)。 # 雖然你可以拒絕不存在差異的原假設，但這個檢驗并沒有告訴你哪些地區顯著地與其他地區 不同。要回答這個問題，你可以使用Mann–Whitney U檢驗每次比較兩組數據。一種更為優雅的 方法是在控制犯第一類錯誤的概率(發現一個事實上并不存在的差異的概率)的前提下，執行可 以同步進行的多組比較，這樣可以直接完成所有組之間的成對比較。npmc包提供了所需要的非 參數多組比較程序。 # 代碼清單7-20 非參數多組比較 class <- state.region var <- state.x77[, c("Illiteracy")] mydata <- as.data.frame(cbind(class, var)) rm(class, var) # 安裝npmc包報錯： package ‘npmc’ is not available (for R version 3.5.2) install.packages("npmc") library(npmc) summary(npmc(mydata), type = "BF") # 調用了npmc的語句生成了六對統計比較結果(東北部對南部、東北部對中北部、東北部對。可以從雙側的p值(p.value.2s)看 處可以看到南部的文西部、南部對中北部、南部對西部，以及中北部對西部) 出南部與其他三個地區顯著不同，而其他三個地區之間并沒有什么不同。在 盲率中間值更高。注意，npmc在計算積分時使用了隨機數，所以每次計算的結果會有輕微的不同# 7.6 組間差異的可視化 # 在7.4節和7.5節中，我們關注了進行組間比較的統計方法。使用視覺直觀地檢查組間差異，同 樣是全面的數據分析策略中的一個重要組成部分。它允許你評估差異的量級、甄別出任何會影響結 果的分布特征(如偏倚、雙峰或離群點)并衡量檢驗假定的合理程度。R中提供了許多比較組間數 11 據的圖形方法，其中包括6.5節中講解的箱線圖(簡單箱線圖、含凹槽的箱線圖、小提琴圖)、6.4.1 節中疊加的核密度圖，以及在第9章中討論的評估檢驗假定的圖形方法。

總結

以上是生活随笔為你收集整理的《R语言实战》第7章的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。