R基本统计分析

#一、基本统计分析之基本方法####
#R语言实战—第7章（描述性统计、频数表和列联表、相关及检验、t检验等）
--------------------------------------------------------------------
#1.1 描述性统计分析####

#例子：通过sapply()
#定义一个函数,其中自变量为x,这里选择不忽视缺失值
mystats<-function(x,na.omit=F) {
x<-x[!is.na(x)] #把x中所有不是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)) #返回值
}
mywars<-c('mpg','hp','wt')
sapply(mtcars[mywars],mystats)

#1.1.1 不分组计算描述统计量

##整体看描述统计量
mywars<-c('mpg','hp','wt')
head(mtcars[mywars])
summary(mtcars[mywars])

#通过Hmisc包的describe()
library('Hmisc')
mywars<-c('mpg','hp','wt')
describe(mtcars[mywars])

#通过psych包的describe()
library('psych')
describe(mtcars[mywars])

#通过pastecs包的stat.desc()
library('pastecs')
mywars<-c('mpg','hp','wt')
stat.desc(mtcars[mywars])

#1.1.2 分组计算描述统计量
##通过aggregate() 指定单一函数返回单一值
mywars<-c('mpg','hp','wt')
aggregate(mtcars[mywars],by=list(am=mtcars$am),mean)

#只分组,默认返回值 describeBy() 函数不允许指定任意函数
library('psych')
describeBy(mtcars[mywars],list(mtcars$am)) #增加分组变量

#通过分组,指定函数,一次返回若干统计值
#a 通过by分组
dstats<-function(x)sapply(x,mystats) #调用前面mystats的函数
myvars<-c('mpg','hp','wt')
by(mtcars[myvars],mtcars$am,dstats)

#b 通过summaryBy分组
#语法 summaryBy(formula=,data=,FUN=),formula左侧为数值型变量,右侧为类别型分组变量
library("doBy", lib.loc="D:/R/R-3.6.0/library")
summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)

-----------------------------------------------------------------------------
#1.2 频数表和列联表####
#使用table(var1,var2,...) #针对n个类别型变量(因子)创建n维列联表
#使用xtabs(formula,data) #根据一个公式和一个矩阵或数据框创建一个N维列联表

#考察类别型变量的频数表和列联表
library('vcd')
head(Arthritis)

#1.2.1 一维列联表
mytable<-table(Arthritis$Improved)
mytable<-with(Arthritis,table(Improved))

prop.table(mytable) #将频数转化为比例值(小数)
options(digits = 4)
prop.table(mytable)*100 #将比例转化为百分比

#1.2.2 二维列联表 p138
#两种方法
table(A,B) #A表示行变量,B表示列变量
mytable<-table(Arthritis$Improved,Arthritis$Treatment)

xtabs(~A+B,data=) #~之后的最后一个变量分类作为二维表的列
mytable<-xtabs(~Improved+Treatment,data=Arthritis) #此处Treatment中类别作为列

#1,2,3 根据二维表生成边际频数(结果同一维列联表)
#1.2.3.1 生成边际频数
mytable<-xtabs(~Improved+Treatment,data=Arthritis)
#指定单个变量
margin.table(mytable,1) #求和边际频数,下标1代表table语句中的第1个变量求和
#不指定变量
addmargins(mytable) #添加各行与各列和,addmargins()默认为所有变量创建边际和
addmargins(mytable,1) #仅添加各列的和
addmargins(mytable,2) #仅添加各列的和

#1.2.3.2 生成边际频数的比例
mytable<-xtabs(~Improved+Treatment,data=Arthritis)
prop.table(mytable) #各单元所占比例(小数)
prop.table(mytable)*100 #各单元所占百分比
#指定单个变量
prop.table(mytable,1) #根据二维表生成行与行的比例,1含义同上
#不指定变量
addmargins(prop.table(mytable)) #添加各单元边际百分比
addmargins(prop.table(mytable,1))*100#table语句中的第1个变量边际百分比和
addmargins(prop.table(mytable,1),2) #仅添加各行的和
addmargins(prop.table(mytable,1),1) #仅添加各列的和

#1.2.4 多维列联表 p141
#三维列联表
library('vcd')
mytable<-xtabs(~ Treatment+Sex+Improved,data=Arthritis)
ftable(mytable) #紧凑的方式输出多维列表

margin.table(mytable,1)
margin.table(mytable,c(1,2))
ftable(margin.table(mytable,c(1,2))) #注意和和上式结果显示区别

margin.table(prop.table(mytable),c(1,2))
ftable(margin.table(prop.table(mytable),c(1,2)))

addmargins(prop.table(mytable,c(1,2)),1)
addmargins(prop.table(mytable,c(1,2)),2)
addmargins(prop.table(mytable,c(1,2)),3) #注意三式区别

ftable(addmargins(prop.table(mytable,c(1,2)),1))
ftable(addmargins(prop.table(mytable,c(1,2)),2))
ftable(addmargins(prop.table(mytable,c(1,2)),3))
----------------------------------------------------------------------

#接下来检验列联表中变量是否相关或独立
#1.3 独立性检验####
#检验类别型变量独立性
#独立性检验原假设H0:相互独立;
#备择假设H1:相互不独立;可以进一步检验相关性 p143
#1.3.1 卡方独立性检验
mytable<-table(Arthritis$Improved,Arthritis$Treatment)
chisq.test(mytable) #检验二维表的行变量和列变量的独立性
mytable<-xtabs(~Improved+Sex,data=Arthritis)
chisq.test(mytable)

#1.3.2 fisher精确检验
#原假设H0:边界固定的列联表中的行和列是相互独立的
#fisher.test()可以在任意行列数大于等于2的二维表上使用，但不能用于2x2的列联表
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
fisher.test(mytable)

#1.3.3 Cochran-Mantel-Haenszel检验
#原假设:两个名义变量在第三个变量中的每一层中都是条件独立的
mytable<-xtabs(~Treatment+Improved+Sex,data=Arthritis)
ftable(mytable)
mantelhaen.test(mytable) #结果表明：患者接受的治疗与得到的改善在性别的每一个水平下并不独立

#1.3.4 相关性的度量
#拒绝独立性原假设(p<0.05)时,则转向衡量相关性强弱的相关性度量 p145
library('vcd')
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
assocstats(mytable) #总体来说,较大值意味着较强的相关性
---------------------------------------------------------------------

#1.4 相关性检验####
#1.4.1 相关系数p147:用来描述定量变量之间的线性相关程度,三种相关系数:
#pearson积差相关系数衡量了两个定量变量之间的线性相关程度;
#Spearman等级相关系数则衡量分级定序变量之间的相关程度;
#Kendall's Tau也是一种非参数的等级相关度量
#cor函数可以计算这三种相关系数
#cor(x,method='pearson') #method默认pearson
states<-state.x77[,1:6]
cor(states) #默认得到的是所有结果的方阵

#选择变量检测相关系数
x<-state.x77[,1:4]
y<-state.x77[,c('Life Exp','Murder')]
cor(x,y)

#1.4.2 相关性的显著性检验 p148
#前言:知道相关系数后,并不能说明总体相关系数显著不为0，
#还需要对相关系数进行显著性检验
#原假设Ho:变量间不相关(即总体相关系数为0)

#1 cor.test() 每次只能检验一种相关关系的显著性(两个变量间)
#用法:
cor.test(x,y, #x和y是要检验的相关性度量
alternative='two.side', #指定进行双侧检验或单侧检验,取值two.side/less/greater
method= ) #指定要计算的相关系数种类,同前

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

#2 corr.test() 一次检验多个变量间的相关关系的显著性
#corr.test()返回相关矩阵并进行显著性检验
library("psych", lib.loc="D:/R/R-3.6.0/library")
corr.test(x=states,
use='complete', # use取值'complete'或'pairwise'
#分别表示对缺失值执行成对删除或行删除
method='pearson') #默认,同前

#1.5 t检验####
#1.5.1独立样本的t检验
#关注的是结果变量y是连续型的变量的组间比较
#针对两组的独立样本t检可用于检验两个总体均值相等的假设,
#前提是假设两组数据是独立的,并且服从正态分布
library("MASS", lib.loc="D:/R/R-3.6.0/library")
t.test(y~x) #y因变量~x自变量
#H0:两个总体均值相等(u0=u1)
t.test(Prob~So,data=UScrime)
#t.test()中的t检验默认假定方差不相等，默认双侧检验
t.test(Prob~So,data=UScrime,
var.equal=TRUE #假定方差相等,使用合并方差估计
)