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