R语言实战-第七章基本统计分析

#第7章 基本统计分析
#描述性统计分析
myvars <- c("mpg","hp","wt")
head(mtcars[myvars])
#方法1 基础函数summary
summary(mtcars[myvars])
#利用sapply或者apply可以对描述性统计量进行计算sapply(x,FUN,options)
sapply(mtcars[myvars],fivenum)#分别是最小 下四分位 中位数 上四分位 最大值
#FUN可以是mean length range var max quantile...

#方法2 packages
#如Hmisc::describe 
#pastecs::stat.desc 
#  psych::describe 

library(Hmisc)
#Hmisc::describe可以返回变量 观测数量 缺失值 唯一值的数目
#均值 分位数 5个最大和最小值
Hmisc::describe(mtcars[myvars])

library(pastecs)
#pastecs::stat.desc
#若basic=TRUE(默认值)，则计算所有值的 空值 缺失值的数量，以及
#最小值 最大值 值域 总和
#若desc=TRUE(默认值)，则计算中位数 平均值 平均数的标准误 置信度为95%的置信区间
#方差 标准差 变异系数
#若norm=TRUE(非默认)，则返回正态分布统计量，
#包括偏度skewness和峰度kurtosis和shapiro-wilk正态检验结果
pastecs::stat.desc(mtcars[myvars],norm=TRUE)

#stat.desc()原始的输出用科学计数法来显示，可以用round()调整输出数字
#正态分布的偏度系数和峰度系数均为0.
#偏度系数为正表示分布中低值较多;偏度系数为负表示分布中高值较多。
#峰度系数为正表示分布峰的形状较正态分布更尖、分布形状相对“陡峭”;峰度系数为负表示分布峰的形状较正态分布更平、分布形状相对“矮胖”。
#峰度系数和偏度系数离0越远,表明数据越不太可能符合正态分布。
#stat.desc()的运行结果中有skew.2SE和kurt.2SE,
#分别对应偏度系数除以2倍标准误和峰度系数除以2倍标准误,
#即转换后的Z值再除以2（这里将1.96近似为2处理了）,
#因此skew.2SE绝对值大于1，就意味着偏度系数为0的差异性检验具有统计学意义,
#kurt.2SE同理。
qplot(wt,data = mtcars)#快速绘图hitogram可以看看分布情况

library(psych)
#非缺失值的数量 平均值 标准差 中位数 截尾均值 绝对中位差 
#最小值 最大值 域值 偏度 峰度 标准误
psych::describe(mtcars[myvars])

#分组计算描述性统计量
#方法1 aggregate(data,by=list(name1=groupvar1,nameN=groupvarN...),mean/sd...)
aggregate(mtcars[myvars],by=list(gear=mtcars$gear),mean)
#aggregate仅允许在每次调用中使用平均数、标准差等这样的单返回值函数
#无法一次返回若干个统计量,用by函数可以

#方法2 by(data,INDICES,FUN)
mystate <- 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^2)/n
  kurt <- sum((x-m)^4/s^4)/n-3
  return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))
}

dstats <- function(x)sapply(x, mystate)
by(mtcars[myvars],mtcars$am,dstats)
data <- data.frame(age=rep(c(21:24),25),
                   sex=rep(c("M","F"),50),
                   level=round(runif(100,min=10,max=50),digits = 2),
                   code=c(1:100))

by(data["level"],data$sex,dstats)#注意by后第一个参数是[]而不是$

myvar1 <- c("age","level","code")
by(data[myvar1],data$sex,dstats)
data[1,1] <- NA #可见此函数不会忽略NA
by(data[myvar1],data$sex,dstats)

#分组计算的扩展
#doBy::summaryBy(formula,data=dataframe,FUN=function)
#formula接受这样的格式var1+var2...需要分析的变量~groupvar1+groupvar2...分组变量
library(doBy)
doBy::summaryBy(hp+mpg+wt~am,data=mtcars,FUN=mystate)

#psych::describeBy 此函数不允许指定函数，因此普适性较低
psych::describeBy(mtcars[myvars],list(am=mtcars$am))

#频数表和列联表
#用到的包为vcd中的Arthritis数据集
library(vcd)
head(Arthritis)
str(Arthritis)

#生成频数表
#table() xtabs() prop.table() margin.table() addmargins() ftable()

#一维列联表
#table()
mytable <- with(Arthritis,table(Improved))
mytable
#prop.table()可以将频数转化成比例值
prop.table(mytable)
#prop.table()*100可以将频数转化成百分比
prop.table(mytable)*100

#二维列联表
#table(A,B)AB分别是行和列变量，默认忽略NA值,若要将其视为一个有效的类别,则useNA="ifany"
mytable2 <- with(Arthritis,table(Improved,Sex))
mytable3 <- table(Arthritis$Improved,Arthritis$Sex)
#xtabs(~A+B+...,data=mydata)
#prop.table() margin.table()生成边际频数和比例
mytable4 <- xtabs(~Improved+Sex,data=Arthritis)#~左边若是某个变量，则为一个频数向量
margin.table(mytable4)
prop.table(mytable4)
margin.table(mytable4,2)
prop.table(mytable4,2)
#1 2分别代表table语句中的第1/2个变量

#addmargin为表格添加边际和
addmargins(mytable4)
addmargins(margin.table(mytable4,2))
addmargins(prop.table(mytable4,1))#默认为表中所有变量仅创造各行/列的边际和(和为1)
addmargins(prop.table(mytable4,1),1)

#二维列联表
library(vcd)
head(Arthritis)

library(gmodels)
gmodels::CrossTable(Arthritis$Treatment,Arthritis$Improved,digits=2,prop.r = TRUE)
help(CrossTable)
#可以计算(行 列 单元格)的百分比，指定小数位数；
#进行卡方 Fisher McMemar独立性检验；
#计算期望和(皮尔逊 标准化 调整的标准化)残差；
#将缺失值作为一种有效值；进行行和列标题的标注；生成SAS或SPSS风格

#三维列联表
#table() xtabs() prop.table() margin.table() addmargins() ftable()都可以
mytable <- xtabs(~Treatment+Sex+Improved,data = Arthritis )#频数
ftable(mytable)
margin.table(mytable,c(1,3)) #边际频数
ftable(prop.table(mytable,c(1,2)))
ftable(addmargins(prop.table(mytable,c(1,2)),3))

#独立性检验
#chisq.test() fisher.test() mantelhaen.test()
mytable <-  xtabs(~Treatment+Improved,data = Arthritis)
mytable2 <-  xtabs(~Improved+Sex,data = Arthritis)
chisq.test(mytable)#p<0.01 拒绝原假设，表明不独立，说明治疗情况和改善情况不独立
chisq.test(mytable2)#治疗独立情况与性别独立

fisher.test(mytable)#不能用于2*2的列联表，可以用于任意大于等于2的二维列联表

mytable3 <-  xtabs(~Treatment+Improved+Sex,data = Arthritis)
mantelhaen.test(mytable3)#原假设是两名义变量在第三个变量的每一层中都是条件独立的
#p<0.01分性别看 治疗情况和改善情况不独立，即会相互影响

#相关性的度量(独立性检验后如果结果是不独立，则可以进行相关性的度量)
#vcd包中的assocstats函数 计算phi系数 列联系数 Cramer's V 系数
library(vcd)
vcd::assocstats(mytable)
#Phi系数：phi =1，表明变量完全相关；phi = 0，表明变量相互独立；
#列联系数：适用于大于2×2列联表，当表中两变量相互独立时，系数为0。
#其最大值依赖于行数和列数，但是不可能大于1。
#Cramer’s V系数：两个变量相互独立时，等于0；两个变量完全相关时，等于1。
#总体来说，较大的值意味着较强的相关性。

#相关
#方形相关矩阵
states <- state.x77[,1:6]
cov(states)#方差和协方差
#cor(x,use,method)
#use=all.obs/everything默认/complete.obs/pairwise.complete.obs 
#method=pearson默认/spearman/kendall
cor(states)
cor(states,method = "spearman")

#非方形相关矩阵
x <- states[,c("Population","Income","Illiteracy","HS Grad")]
y <- states[,c("Life Exp","Murder")]
cor(x,y)

#偏相关
library(ggm)
#ggm::pcor(u,S) u是数量向量，前两个数值表示要计算的相关系数的变量下标
#其余数量为条件变量(即需要排除的变量)，S为变量的协方差阵
colnames(states)
ggm::pcor(c(1,5,2,3,6),cov(states))
#在控制了收入、文盲率、高中毕业率的影响时，人口和谋杀率之间的相关系数为0.346

#其它类型相关
#polycor::hetcor()

#相关性的显著性检验，判断样本相关系数对总体相关程度的代表性
#原假设为变量之间不相关(即总体的相关系数为0)

#cor.test(x,y,alternative=two.side/less/greater, method= pearson/spearman/kendall)只能检验一种相关关系
cor(states[,3],states[,5])
cor.test(states[,3],states[,5],alternative = "greater",method = "pearson")
#p越小 越拒绝原假设，接受备择假设，表明变量间相关

#psych::corr.test
library(psych)
psych::corr.test(states,use = "complete")
#其它显著性检验
ggm::pcor.test

#t检验(两样本检验)
#独立样本t检验
#调用格式t.test(y~x,data)或t.test(y1,y2)
library(MASS)#UScrime数据集
t.test(Prob~So,data=UScrime)

#非独立样本t检验
with(UScrime,t.test(U1,U2,paired=TRUE))
sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),sd=sd(x))))

#方差分析(多于两组的情况,第9章细讲)

#非参数检验
#两组比较
#wilcoxon轶和检验
#调用格式wilcox.test(y~x,data)或wilcox.test(y1,y2)
#独立
wilcox.test(Prob~So,UScrime)
with(UScrime,by(Prob,So,median))

#非独立
with(UScrime,wilcox.test(U1,U2,paired=TRUE,exact=FALSE))
#数据中有相同的数字，成为结（tie）,结中的秩为按照升幂排列后位置的平均值，用exact=FALSE处理
sapply(UScrime[c("U1","U2")],median)

#多于两组的比较
#调用格式kruskal.test(y~A,data) 
#friedman.test(y~A|B,data)
states <- data.frame(state.region,state.x77)
kruskal.test(Illiteracy~state.region,data=states)