R语言实战(第二版)第七章-基本统计分析

Peaches’Notepad

综合各路大神,讲了三遍可是把自己讲会了

代码们 开整


 


#7.1 描述性统计分析
#示例
myvars<-c("mpg","hp","wt")#生成向量
head(mtcars[myvars])#查看
dt<-mtcars[myvars]#赋值给dt

#7.1.1方法云集
#基础安装
summary(dt)
#图基五数总括
#   即最小值、下四分位数、中位数、上四分位数和最大值)针对其中一个变量
fivenum(dt$hp)

#apply(x,margin,FUN,...)将一个任意函数“应用”到矩阵、数组、数据框的任何维度上
apply(dt,2,mean)

#sapply(x,FUN)
      #自定义一个函数:
mystats<-function(x,na.omit=FALSE){
  if(na.omit)#如果忽略na
    x<-x[!is.na(x)]#判断 如果x是非缺失值,去掉了na
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(dt, mystats)

#7.1.2 更多方法


install.packages("Hmisc")
#Hmisc::describe(dt)
library(Hmisc)
describe(dt)

install.packages("pastecs")
# stat.desc(x, basic = TRUE, desc = TRUE, norm = FALSE, p = 0.95)
library(pastecs)
stat.desc(dt,norm = TRUE)

install.packages("psych")
library(psych)
describe(dt)

#7.1.3 分组计算描述性统计量
#如果使用的是list(mtcars$am),则am列将被标注为 Group.1而不是am.
aggregate(dt,by=list(mtcars$am),sd)
aggregate(dt,by=list(am=mtcars$am),sd)
#如果有多个分组变量,可以使用by=list(name1=groupvar1, name2=...
aggregate(dt,by=list(am=mtcars$am,cyl=mtcars$cyl),mean)

#aggregate只能使用单返回值函数
#使用by()分组计算描述性统计量 自定义函数

dstats<-function(x)sapply(x,mystats)
by(dt,mtcars$am,dstats)

#7.1.4分组计算的扩展
#1)doBy
install.packages("doBy")
library(doBy)
summaryBy(mpg+hp+wt~mtcars$am,data=mtcars,FUN = mystats)

#2)psych   
install.packages("psych")
library(psych)
describeBy(dt,list(am=mtcars$am))#与describe相同的统计量,不允许指定函数


#3)reshape
install.packages("reshape")
library(reshape)
dstats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x)))
#融合:dfm <- melt(dataframe, measure.vars = y, id.vars = g)
#y:数值型变量(默认使用所有变量),g:一个或多个分组变量组成的向量

#重铸:cast(dfm, groupvar1, groupvar2 + ... + variable - ., FUN)
#variable只取其字面含义(即仅表示重铸后数据框中的变量variable)
dfm<-melt(mtcars,measure.vasrs=c("mpg","hp","wt"),id.vars=c("am","cyl"))
cast(dfm,am+cyl+variable~.,dstats)

#7.2频数表和列联表
library(vcd)
head(Arthritis)

#7.2.1 生成频数表
#1)一维列联表
#table()函数
mytable1<-with(Arthritis,table(Improved))
mytable1
#prop()函数---转化为比例值
prop.table(mytable1)
prop.table(mytable1)*100

#2)二维列联表
#table(A,B)函数
mytable2<-table(Treatment=Arthritis$Treatment,Improved=Arthritis$Improved)
mytable2

#xtabs(~A+B,data=mydata)函数
mytable2<-xtabs(~Treatment+Improved,data=Arthritis)
mytable2

#margin.table()
margin.table(mytable2,1)#对行进行操作
margin.table(mytable2,2)

#prop.table()
prop.table(mytable2,1)
prop.table(mytable2,2)

##各单元格所占比例
prop.table(mytable2)

#概述边margins,为表格添加边际和
addmargins(mytable2)#(默认是为表中所有变量创建边际和)
addmargins(prop.table(mytable2))

addmargins(prop.table(mytable2,1),2)
addmargins(prop.table(mytable2,2),1)

#3)gmodels-CrossTable
#install.packages("gmodels")
library(gmodels)
CrossTable(Arthritis$Treatment,Arthritis$Improved)
#3)多维列联表
mytable3<-xtabs(~Treatment+Sex+Improved,data=Arthritis)
mytable3

#ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表
ftable(mytable3)
#边际频数
margin.table(mytable3,1)
margin.table(mytable3,2)
margin.table(mytable3,3)

margin.table(mytable3,c(1,3))##治疗情况x改善情况的分组边际频数

#治疗情况x性别的各类改善情况的比例,比例将被添加在不在prop.table调用的下标上
ftable(prop.table(mytable3,c(1,2)))
#治疗情况x性别的各类改善情况的边际和
##将被添加在不在prop.table调用的下标上
ftable(addmargins(prop.table(mytable3,c(1,2)),3))
ftable(addmargins(prop.table(mytable3,c(1,2)),3))*100

#7.2.2独立性检验:
#1、卡方独立性检验
#对数据框中的两列进行检验,是否独立
##原假设 相互独立

#治疗情况与改善情况
mytable2<-xtabs(~Treatment+Improved,data=Arthritis)
chisq.test(mytable2)
#治疗情况和改善情况不独立

#性别与改善情况
mytable3<-xtabs(~Improved+Sex,data=Arthritis)
mytable3
chisq.test(mytable3)
#性别和改善情况独立
#warning出现,因为表中单元格之一有一个小于5的值,可能会导致卡方近似无效

#Fisher精确检验
##不能用于2x2列联表
mytable2<-xtabs(~Treatment+Improved,data=Arthritis)
fisher.test(mytable2)


#Cochran-Mantel-Haenszel检验
#检验治疗情况和改善情况在性别的每一水平下是否独立(原假设独立)

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

#7.2.3相关性的度量

library(vcd)
assocstats(mytable2)



#7.3 相关

#7.3.1 相关的类型


# cor(x,use=,method= ) 
states<-state.x77[,1:6]
cor(states)#默认的method是pearson
cor(states,method="spearman")

cov(states)


#有目的地针对某两组变量之间的关系进行研究
#非方形的相关矩阵
x <- states[,c("Population", "Income", "Illiteracy", "HS Grad")]
y <- states[,c("Life Exp", "Murder")]
cor(x,y)

#偏相关:在控制一个变量时,另两个变量的相关性。
install.packages("ggm")
library(ggm)
colnames(states)
pcor(c(1,5,2,3,6),cov(states))

#7.3.2 相关性的显著性检验
#cor.test(x,y,alternative=less/greater/two side,method= )

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

#检验相关矩阵下的所有相关系数:corr.test()
library(psych)
corr.test(states,use="complete")
#若需要显示相关系数95%CI,以print()函数打印该检验结果并声明参数short=FALSE
print(corr.test(states,method = "pearson"),short = FALSE)

#偏相关系数的显著性检验
library(ggm)
#pcor.test(pcor()输出结果,q=控制变量序号,n=观测记录数)
pcor.test(r=pcor(c(1,5,2,3,6),cov(states))
,q=c(2,3,6),n=nrow(states))

#7.4 t检验
#7.4.1 独立样本t检验
#t.test(y~x,data) OR t.test(y1,y2)
library(MASS)
dt2<-UScrime
t.test(Prob~So,data=dt2)

#7.4.2 非独立样本的t检验
library(MASS)
sapply(dt2[c("U1","U2")], function(x)(c(mean=mean(x),sd=sd(x))))
#配对t检验
 with(dt2,t.test(U1,U2,paired=TRUE))

       
#7.5组间差异的非参数检验
#7.5.1 两组的比较
with(dt2, by(Prob, So, median))
#相互独立
wilcox.test(Prob~So,data=dt2)
#不相互独立
with(dt2, wilcox.test(U1, U2, paired = TRUE))

#7.5.2多于两组的比较
#各组独立
dts<-data.frame(state.region,state.x77)
kruskal.test(Illiteracy~state.region,data = dts)
#四个地区的文盲率各不相同

#各组不独立
dt3<-matrix(1:20,nrow = 4,ncol = 5)
dt3
friedman.test(dt3)


#作者拓展
source("http://www.statmethods.net/RiA/wmc.txt")              
states <- data.frame(state.region, state.x77)
wmc(Illiteracy ~ state.region, data=states, method="holm")




























  • 1
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值