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-designed
dstats <- 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)
#summaryBy
library(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 casr
library(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)
mytable
ftable(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.01
mytable = 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，即根据样本数据是否有足够的证据得出整体相关系数不为0
states=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.01
library(MASS)
t.test(Prob~So,data=UScrime)
#非独立样本的t检验
#假设：组件的差异呈现正态分布
#P值反映了如果实际均值相等，那么获得一个差异如此大的样本的概率小于2.2e-16
library(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检验的一种非参数替代方法。适用于两组成对数据和无法保证正态性假设的情景。
#1
with(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.region
var = state.x77[,c("Illiteracy")]
mydata = as.data.frame(cbind(class,var))
rm(class,var)
install.packages("npmc")
library(npmc)
summary(npmc(mydata),type="BF")
#组件差异的可视化