R语言统计分析篇

最新推荐文章于 2024-05-17 08:30:00 发布

bestlinjiayin

最新推荐文章于 2024-05-17 08:30:00 发布

阅读量1.3k

点赞数

分类专栏： R

R 专栏收录该内容

3 篇文章 0 订阅

订阅专栏

1.描述性统计分析

（1）方法云集

通过summary，sapply（）计算描述性统计量

[plain] view plain copy

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（）函数计算描述性统计量

[plain] view plain copy

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

[plain] view plain copy

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

通过psych包中的describe（）计算描述性统计量

可计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误

[plain] view plain copy

library(psych)
describe(mtcars[vars])

（2）分组计算描述性统计量

使用aggregate（）分组获取描述性统计量

注意：aggregate（）仅允许在每次调用中使用平均数、标准差这样的单返回值函数

[plain] view plain copy

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

使用by（）分组计算描述性统计量

by（data,INDICES,FUN）

[plain] view plain copy

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

使用doBy包中的summaryBy（）分组计算概述统计量

[plain] view plain copy

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

使用psych包中的describe.by（）分组计算概述统计量

[plain] view plain copy

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

使用reshape包分组计算概述统计量

[plain] view plain copy

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）一维列联表

[plain] view plain copy

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：

[plain] view plain copy

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生成二维列联表

[plain] view plain copy

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

（4）多维列联表

table（）

xtabs（）

margin.table（）

prop.table（）

addmargins（）

ftable（）

eg：

[plain] view plain copy

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

[plain] view plain copy

ftable(mytable)

[plain] view plain copy

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（）

[plain] view plain copy

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

[plain] view plain copy

</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)

[plain] view plain copy

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

（6）相关性的度量

若变量间不独立，衡量相关性强弱的相关性度量

[plain] view plain copy

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

较大的值意味着较强的相关性

vcd包中也提供了一个kappa（）函数，可以计算混淆矩阵的Cohen's kappa值以及加权的kappa值（混淆矩阵可表示两位遮天盖地者对于一系列对象进行分类所得结果的一致程度）

（7）结果可视化

将表转换为扁平格式

通过table2flat将表转换为扁平格式

[plain] view plain copy

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）