R語言統計分析篇

阿新 • • 發佈：2019-01-19

1.描述性統計分析

（1）方法雲集

通過summary，sapply（）計算描述性統計量

vars<-c("mpg","hp","wt")
vars
head(mtcars[vars])
#通過summary（）函式來獲取描述性統計量
summary(mtcars[vars])
#sapply（）函式，格式：sapply（x，FUN，options）x:資料框或矩陣;FUN：為任意的函式；options：若被指定，則將被傳遞給FUN
#通過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)/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)

通過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

library(pastecs)
stat.desc(mtcars[vars])

通過psych包中的describe（）計算描述性統計量

可計算非缺失值的數量、平均數、標準差、中位數、截尾均值、絕對中位差、最小值、最大值、值域、偏度、峰度和平均值的標準誤

library(psych)
describe(mtcars[vars])

（2）分組計算描述性統計量

使用aggregate（）分組獲取描述性統計量

注意：aggregate（）僅允許在每次呼叫中使用平均數、標準差這樣的單返回值函式

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

使用by（）分組計算描述性統計量

by（data,INDICES,FUN）

dstats<-function(x)(c(mean=mean(x),sd=sd(x)))
by(mtcars[vars],mtcars$am,dstats)

使用doBy包中的summaryBy（）分組計算概述統計量

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

使用psych包中的describe.by（）分組計算概述統計量

library(psych)
describeBy(mtcars[vars],mtcars$am)

使用reshape包分組計算概述統計量

library(reshape)
dstats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x)))
dfm<-melt(mtcars,measure.vasrs=c("mpg","hp","wt"),id.vars=c("am","cyl"))
cast(dfm,am+cyl+variable~.,dstats)

2.頻數表和列聯表

（1）生成頻數表

用於建立和處理列聯表的函式

函式	描述
table（var1，var2,…，varN）	使用N個類別變數（因子）建立一個N維列聯表
xtabs（forula，data）	根據一個公式和一個矩陣或資料框建立一個N維列聯表
prop.table（table，margins）	依margins定義的邊際列表中條目表示為分數形式
margin.table（table，margins）	依margins定義的邊際列表計算表中條目的和
addmargins（table，margins）	將概述邊margins放入表中
ftable（table）	建立一個緊湊的”平鋪“式列聯表

（2）一維列聯表

install.packages("vcd")
library(vcd)
mytable<-with(Arthritis,table(Improved))
mytable
#使用prop.table（）將頻數轉化為比例值：
prop.table(mytable)
#使用prop.table()*100轉化為百分比
prop.table(mytable)*100

（3）二維列聯表

table（）

格式： table（A，B）#A是行變數，B是列變數 xtabs（~A+B，data=mydata）#其中mydata是一個矩陣或資料框；要進行交叉分類的變數應出現在公式的右側（即~符號的右方），以+作為分隔符 eg：

mytable<-xtabs(~Treatment+Improved,data=Arthritis)
mytable
#可以使用margin.table（）和prop.table（）函式分別生成邊際頻數和比例
margin.table(mytable,1) #下標1指代table（）語句中的第一個變數
prop.table(mytable,1)
#列和列比例可以這樣計算
margin.table(mytable,2)#下標2指代table（）語句中的第二個變數
prop.table(mytable,2)
#各單元所佔比例可用如下語句獲取
prop.table(mytable)
#可使用addmargins（）函式為表格新增邊際和；使用addmargins（）時，預設行為是為表中所有的變數建立邊際和
addmargins(mytable)
addmargins(prop.table(mytable))
#僅新增各行的和
addmargins(prop.table(mytable,1),2)
#僅新增各列的和
addmargins(prop.table(mytable,2),1)

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

使用CrossTable生成二維列聯表

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

（4）多維列聯表

table（）

xtabs（）

margin.table（）

prop.table（）

addmargins（）

ftable（）

eg：

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

ftable(mytable)

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))
ftable(addmargins(prop.table(mytable,c(1,2)),3))*100 #得到百分比而不是比例，可將結果乘以100

（5）獨立性檢驗

名稱	原理	函式
卡方獨立性檢驗	對二維表的行變數和列變數進行卡方獨立性檢驗	chisq.test（）
Fisher精確檢驗	邊界固定的列聯表中行和列是相互獨立的注意：與許多統計軟體不同的是，這裡fisher.test（）函式可以在任意行列數大於等於2的二維旬聯表上使用，但不能用於2*2的列表	fisher.test（）
Cochran-Mantel-Haenszel檢驗	兩個名義變數在第三個變數的每一層中都是條件獨立的	mantelhaen.test（）

library(vcd)
mytable<-xtabs(~Treatment+Improved,data=Arthritis) #Treatment和Improved是否獨立
chisq.test(mytable)
mytable<-xtabs(~Improved+Sex,data=Arthritis)
chisq.test(mytable)

</pre><pre code_snippet_id="374527" snippet_file_name="blog_20140604_18_9851680" name="code" class="plain">mytable<-xtabs(~Treatment+Improved,data=Arthritis) 
fisher.test(mytable)

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

（6）相關性的度量

若變數間不獨立，衡量相關性強弱的相關性度量

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

較大的值意味著較強的相關性

vcd包中也提供了一個kappa（）函式，可以計算混淆矩陣的Cohen's kappa值以及加權的kappa值（混淆矩陣可表示兩位遮天蓋地者對於一系列物件進行分類所得結果的一致程度）

（7）結果視覺化

將錶轉換為扁平格式

通過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)
}
table2flat(mytable)

mytable<-xtabs(~Treatment+Sex+Improved,data=Arthritis)
mytable
mytable<-ftable(mytable)
df<-as.data.frame(mytable)
df
rows<-dim(df)[1]
rows
cols<-dim(df)[2]
cols
df$Freq[1]

3.相關

（1）相關型別

相關型別	定義	函式
Pearson積差相關係數	衡量兩個定量變數之間的線性相關程度	cor（x,use=,method=pearson）
Spearman等級相關係數	衡量分級定序變數之間的相關程度	cor（x,use=,method=spearman）
Kendall's Tau相關係數	是一種非引數的等級相關度量	cor（x,use=,kendall）

（2）協方差和相關係數

states<-state.x77[,1:6]
states
cov(states)
cor(states)
cor(states,method="spearman")

（3）偏相關

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

可使用ggm包中的pcor（）函式計算偏相關係數

函式呼叫格式為：

Pcor（u,S）

其中u是一個數值向量，前兩個數值表示要計算相關係數的變數下標，其餘的數值為條件變數（即要排除影響的變數）的下標。

S為變數的協方差陣

Eg：

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

（4）相關性的顯著性檢驗

格式：

Cor.test（x,y,alternative=,method=）

cor.test(states[,3],states[,5]) #cor.test每次只能檢驗一種相關關係，psych包中提供的corr.test（）可一次做更多的事

psych包中的pcor.test（）函式可以用來檢驗在控制一個或多個額外變數時兩個變數之間的條件獨立性，

格式為：pcor.test(r,q,n)

R是由pcor（）函式計算得到的偏相關係數

Q為要控制的變數數（以數值表示位置）

N為樣本的大小

4.t檢驗

（1）獨立樣本的t檢驗

檢驗的呼叫格式為：

t.test（y~x，data）

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

（2）非獨立樣本的t檢驗

在兩組的觀測之間相關時，獲得的是一個非獨立組設計（dependent groups design）

前-後測設計（pre-post design）

重複測量設計（repeated measures design）

呼叫格式：

t.test（y1,y2,paired=TRUE）

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

若想在多於兩個的組之間進行比較，若假設資料是從正態總體中獨立抽樣而得的，可使用方差分析（ANOVA）

5. 組間差異的非引數檢驗

（1）兩組的比較

若兩組資料獨立，可使用Wilcoxon秩和檢驗（Mann-WhitneyU檢驗）來評估觀測是否是從相同的概分佈中抽得的 呼叫格式為： wilcox.test（y~x，data）其中y是數值型變數，x是一個二分變數 呼叫格式為： wilcox.test（y1,y2）其中y1和y2為各組的結果變數

library(MASS)
with(UScrime,by(Prob,So,median))

wilcox.test(Prob~So,data=UScrime)

（2）多於兩組的比較

若無法滿足ANOVA設計的假設，可以使用非引數方法來評估組間的差異。

若各組獨立，Kruskal-Wallis檢驗是一種實用的方法

若各組不獨立（如重複測量設計或隨機區組設計），則Friedman檢驗更合適

Kruskal-Wallis呼叫格式為：kruskal.test（y~A，data）（y是數值型結果變數，A是一個分組變數，而B是一個用以認定匹配觀測的區組變數（blocking variable））

states<-as.data.frame(cbind(state.region,state.x77))
kruskal.test(Illiteracy~state.region,data=states)

非引數多組比較

install.packages("npmc")
library(npmc)
class<-state.region
var<-state.x77[,c("Illiteracy")]
mydata<-as.data.frame(cbind(class,var))
rm(class,var)
summary(npmc(mydata),type="BF")
aggregate(mydata,by=list(mydata$class),median)