R语言_基本统计分析

#基本统计分析#整体描述性统计分析，针对数值变量attach(mtcars)opar = par(no.readnoly=TRUE)d = mtcars[c("mpg","hp","wt")]head(d)#summary#较标准正态分布呈现正偏，且较平。（偏度为正，峰度为负）summary(d)plot(density(mpg))#describe#多了峰度，偏度等数据library(psych)describe(d)#分组描述统计，针对数值变量#aggregate,fun = c(mean,sd)aggregate(d,by=list(am=am),mean)#by, func = self-designeddstats <- function(x){    c(mean=sapply(x,mean),sum=sapply(x,sum))}tapply(d$mpg,factor(am),dstats)by(d,factor(am),dstats)by(d,factor(am),summary)#summaryBylibrary(doBy)dstats <- function(x){    c(mean=mean(x),sum=sum(x))}summaryBy(mpg+hp+wt~am,data=mtcars,FUN=dstats)#describe.by #多了峰度，偏度等数据,不允许指定函数，适应度较差library(psych)describeBy(d,am)#melt and casrlibrary(reshape)dstats <- function(x){    c(mean=mean(x),sum=sum(x),length=length(x))}dfm = melt(mtcars,measure.vars = c("mpg","hp","wt"),           id.vars=c("am","cyl"))cast(dfm,am+cyl+variable~.,dstats)#频数表和列联表,针对类别变量#函数总概table(var1,var2)xtabs(formula,data) #根据一个公式和一个矩阵或者数据框创建n维列联表prop.table(table,margins) #根据margins定义的边际列表将表中条目表示为分数形式margin.table(table,margin) #依据margin定义的边界计算和addmargins(table,margins) #将margin（默认求和结果）放入表中ftable(table) #创建一个紧凑的平铺式的列联表#一维列联表#table默认忽略缺失值，若不则useNA="ifany"table(cyl)prop.table(table(cyl))#二维列联表table(cyl,am)mytable = xtabs(~am+cyl,data=mtcars)margin.table(mytable,margin = 1) #行prop.table(mytable,margin = 1) #行prop.table(mytable,margin = 2) #列prop.table(mytable) #行列所占比例#添加边际和的二维列联表addmargins(mytable)addmargins(prop.table(mytable))addmargins(prop.table(mytable,1),2)addmargins(prop.table(mytable,2),1)#CrossTable生成二维列联表install.packages("gmodels")library(gmodels)CrossTable(am,cyl)#多维列联表mytable = xtabs(~am+cyl+gear,data=mtcars)mytableftable(mytable)margin.table(mytable,1)margin.table(mytable,3)margin.table(mytable,c(1,3))ftable(prop.table(mytable,c(1,2))) ftable(addmargins(prop.table(mytable,c(1,2)),3) )#将表转化为扁平格式table2flat <- function(mytable) {    df = as.data.frame(mytable)    rows = dim(df)[1]    cols = dim(df)[2]    x = NULL    counts=0    for (i in 1:rows) {        for (j in 1:df$freq[i]) {            counts = counts+1            row = df[i,c(1:(cols-1))]            x = rbind(x,row)        }        print(c(i,counts))    }    row.names(x) <- c(1:dim(x)[1])    return (list(x,counts))}treatment = rep(c("Placebo","Treated"),times=3)improved  = rep(c("None","Some","Marked"),each=2)freq = c(29,13,7,17,7,21)mytable = as.data.frame(cbind(treatment,improved,freq))mytable$freq = as.numeric(as.character(mytable$freq))mydata = table2flat(mytable)#独立性检验，描述类别变量独立性#卡方独立性检验#卡方备注：#p值表示从总体中抽取样本行变量与列变量相互独立的概率，# p<0.01，概率非常小，所以拒绝相互独立的原假设# p>0.05,概率不够小，没有足够理由说明原来的两个变量是不独立的#产生警告的原因，是6个单元格（男性，一定程度改善）有一个小于5，可能使卡方无效library(vcd)mytable = xtabs(~Treatment+Improved,data=Arthritis)chisq.test(mytable)  #治疗和改善不独立 p<0.01mytable = xtabs(~Sex+Improved,data=Arthritis)chisq.test(mytable)  #性别和改善独立 p>0.05#Fisher精确检验#原假设是：边界固定的列联表中行和列是相互独立的mytable = xtabs(~Treatment+Improved,data=Arthritis)fisher.test(mytable)#Cochran-Mantel-Haenszel检验#原假设是：两个名义变量在第三个变量的每一层中都是条件独立的#下面检验治疗情况和改善情况在性别的每一个水平下是否独立，检验不存在三阶交互作用#结果表明：患者接受的治疗与得到的改善在性别的每一个水平下并不独立#即：分性别看，用药治疗的患者较接受安慰剂的患者得到更多改善mytable = xtabs(~Treatment+Improved+Sex,data=Arthritis)mantelhaen.test(mytable)#相关性的度量#上面的显著性评估目的是判断变量间是否相互独立#若不，则接着衡量相关性的强弱#共得到了phi,列联，Cramer‘s V系数，较大意味着相关性越强library(vcd)mytable = xtabs(~Treatment+Improved,data=Arthritis)chisq.test(mytable)assocstats(mytable)#相关性#上述的独立性检验主要描述类别变量的独立性#针对定量变量，使用相关性去描述#Pearson 衡量两个定量变量之间的线性相关程度#Spearman 衡量分级定序变量之间的相关程度#Kendall 非参数的等级相关度量#cor(x,use=,method=) use制定缺失数据的处理方式，method制定相关系数的类型，默认use=everything,method=person#但是结果并不指明相关系数是否显著不为0，即根据样本数据是否有足够的证据得出整体相关系数不为0states=state.x77[,1:6]cov(states)cor(states)cor(states,method="spearman")x = states[,c("Population","Income","HS Grad")]y = states[,c("Life Exp","Murder")]cor(x,y)#偏相关#指控制一个或多个定量变量时，另外两个定量变量之间的相互关系install.packages("ggm")library(ggm)pcor(c(1,5,2,3,6),cov(states))  #控制2,3,6的影响，判断1,5的相关系数#相关性的显著检验#原假设：变量不相关，相关系数为0#cor.test(x,y,alternative=,method=)cor.test(states[,3],states[,5])#计算相关矩阵并进行显著性检验library(psych)corr.test(states,use="complete")#t检验#关注连续型变量的组间比较，类别型变量参考上文独立性检验部分#例子：新药治疗的患者相比旧药是否有更大程度改善；新工艺是否比旧工艺制造的不合格产品更少#独立样本的t检验#假设：两个总体的均值相等，并且从正态总体中取得#下面进行假设方差不等的双侧检验，比较南方和非南方的监禁概率#可以拒绝相同监禁概率的假设，因为p<0.01library(MASS)t.test(Prob~So,data=UScrime)#非独立样本的t检验#假设：组件的差异呈现正态分布#P值反映了如果实际均值相等，那么获得一个差异如此大的样本的概率小于2.2e-16library(MASS)sapply(UScrime[c("U1","U2")],function(x) (c(mean=mean(x),sd=sd(x))) )with(UScrime, t.test(U1,U2,paired=TRUE))#多于两组的情况#假设数据从正态总体中独立抽样而得ANOVA分析#组件差异的非参数检验#如果数据无法满足t检验或者anova的参数假设，一般采用非参数方法#例如：结果变量在本质上就严重偏斜或呈现有序关系#两组的比较#若两组数据独立，可以使用Wolcoxon秩和检验（Mann-Whitney U检验）。来评估观测是否是从相同概率分布中抽的#即：在一个总体中获得更高得分的概率是否比另一个总体更大#评价：是非独立样本t检验的一种非参数替代方法。适用于两组成对数据和无法保证正态性假设的情景。#1with(UScrime,by(Prob,So,median))wilcox.test(Prob~So,data=UScrime)#2#在本例中，含参的t检验和其作用相同的非参数检验得到了相同的结论。#当t检验的假设合理时，参数检验的功效更强（更容易发现存在的差异）。非参数检验在假设非常不合理（等级有序数）更有效sapply(UScrime[c("U1","U2")],median)with(UScrime, wilcox.test(U1,U2,paired = TRUE))#非参数多于两组的比较class = state.regionvar = state.x77[,c("Illiteracy")]mydata = as.data.frame(cbind(class,var))rm(class,var)install.packages("npmc")library(npmc)summary(npmc(mydata),type="BF")#组件差异的可视化