#描述性统计分析
myvars <- c("mpg","hp","wt")
head(mtcars[myvars])
# summary函数
# 数值型变量，提供min,max,四分位数（1st QU，Median,3rd Qu）,mean
# 因子向量，逻辑向量，频数统计
summary(mtcars[myvars])
#sapply函数
#自定义函数,对数值，剔除空值，计算平均数，方差，偏度，峰度，并返回
mystats <- function(x,na.omit=FALSE){
if(na.omit)
x <- x[!is.na(x)] #剔除x中非空值
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调用自定义函数
myvars <- c("mpg","hp","wt")
sapply(mtcars[myvars], mystats)
#Hmisc包,describe函数
#3 variables 3个变量(3列)，32个观测值（行）
#独立变量：mpg,hp,wt
#每个变量下：
#n: 表示总行数
#missiong: 缺失值个数
#unique:非重复值个数
#mean 均值
#百分位数：0.05,0.10，0.2,0.5,0.75,0.9.0.95
#最小的五个值,最大的五个值
library("lattice")
library("survival")
library("Formula")
library("ggplot2")
library("Hmisc") #首先完成关联依赖的加载
myvars <- c("mpg","hp","wt")
describe(mtcars[myvars])
#pastecs包,stat.desc函数
#stat.desc(x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95)
#当basic=True保持默认值：
#nbr.val 所有值的数量
#nbr.null 空值的数量
#nbr.na 缺失值的数量
#min,max 最小值，最大值
#range 值域
#sum 求和
#当desc=TRUE保持默认值：
#median 中位数
#mean 平均数
#SE.mean 平均数的标准误
#CI.mean.0.95 平均数置信度为95%的置信区间
#var 方差
#std.dev 标准差
#coef.var 变异系数
#当norm=TRUE不保持默认值，默认值为FALSE
#skewness 偏度
#skew.2SE 偏度显著程度
#kurtosis 峰度
#kurt.2SE 峰度显著程度
#normtest.W 正态检验结果
#normtest.p 正态检验结果
#p=0.95 置信度为0.95
#范例一
install.packages("pastecs")
library(pastecs)
myvars <- c("mpg","hp","wt")
stat.desc(mtcars[myvars],norm = TRUE)
#psych包，describe函数
library("psych")
#加载psych之后，提示Hmisc包中的describe函数被psych包中的describe函数覆盖；
#最后载入的程序包优先
#采用Hmisc::describe()
myvars <- c("mpg","hp","wt")
#返回参数介绍
# n 非缺失值的数量
# mean 平均数
# sd 标准差
# median 中位数
# trimmed 截尾均值
# mad 绝对中位差
# min 最小值
# max 最大值
# range 值域
# skew 偏度
# kurtosis 峰度
# se 平均值的标准误
describe(mtcars[myvars])
#分组计算描述性统计量
myvars <- c("mpg","hp","wt")
#基于mtcars$am分组，对mtcars[mpg,hp,wt]三列，分别求均值
#am 对应分组列，0,1均为组值
#am=mtcars$am,假如不赋值列名，将会被标注为Group.1
aggregate(mtcars[myvars],by=list(am=mtcars$am),mean)
#基于mtcars$am分组，对mtcars[mpg,hp,wt]三列，分别求标准差
#am 对应分组列，0,1均为组值
#aggregate() 仅允许在每次调用中使用平均数，标准差这样的单返回值函数，无法一次返回若干个统计量
aggregate(mtcars[myvars],by=list(am=mtcars$am),sd)
#使用by分组计算描述性统计量
dstats <- function(x)sapply(x,mystats) #嵌套函数定义
myvars <- c("mpg","hp","wt")
by(mtcars[myvars],mtcars$am,dstats) #规避aggregate仅能返回单个返回值问题
#doby包中的summaryBy()分组计算概率统计量
library(doBy)
#summaryBy(formula,data=dataframe,FUN=function)
#var1+var2+var3~groupvar1+groupvar2+...
#~左侧的变量是需要分析的数值型变量
#~右侧的变量是类别性的分组变量
summaryBy(mpg+hp+wt~am,data=mtcars,FUN = mystats)
#psych包中describeBy
library(psych)
myvars <- c("mpg","hp","wt")
#基于分组，分别统计各列的各种统计指标
#统计指标：
#n,mean,sd,median,trimmed,mad,min,max,range,skew,kurtosis,se
#不允许指定任意函数，普适性较低
describeBy(mtcars[myvars],list(am=mtcars$am))
#针对类别性变量的图表
#数据源
library("vcd")
head(Arthritis)
#一维列联表
#table函数默认忽略缺失值NA,如频数统计时需要视为一个有效的类别，请useNA="ifany"
mytable <- with(Arthritis,
table(Improved)
) #table 函数,频数统计
prop.table(mytable) #基于频数统计，转成成比例值,0-1间的小数形式表示
prop.table(mytable)*100 #基于频数统计，转成成比例值,1-100间的整数形式表示
#二维列联表
#mytable <- xtabs(~A+B,data=mydata)
#mydata 是一个矩阵或者数据框
#交叉分类的变量应出现在公式的右侧，即~符号的右方，以+作为分隔符
#A的值 行变量 第一个变量
#B的值 列变量 第二个变量
mytable <- xtabs(~Treatment+Improved,data=Arthritis) #生成二维频数表格
margin.table(mytable,1) #基于二维频数统计，二次加工，按照行，生成边际频数（按行求和）
margin.table(mytable,2) #基于二维频数统计，二次加工，按照列，生成边际频数（按列求和）
prop.table(mytable,1)#基于二维频数统计，二次加工，按照行，生成当前频数/边际频数（按行求和）
prop.table(mytable,2)#基于二维频数统计，二次加工，按照列，生成当前频数/边际频数（按列求和）
prop.table(mytable)#基于二维频数表格，不指定行维度，列维度时，所有元素得到的比例总和为1
addmargins(mytable) #基于二维频数表格，追加一列，保存基于行求和的值；追加一行，保存基于列求和的值
addmargins(prop.table(mytable))#基于百分比表格，追加一列，保存基于行求和的%值；追加一行，保存基于列求和的%值
addmargins(prop.table(mytable,1),2)#基于行百分比的表格，追加一列，保存行的和值
addmargins(prop.table(mytable,2),1)#基于列百分比的表格，追加一行，保存列的和值
#使用CrossTable生成二维列联表
library(gmodels)
#结果解释(不同隔断中，相同行位置，或者相同列位置值相加)
#第一行
#整数部分剥离
#相当于addmargins(mytable)),返回二维频数边际统计表
#行方向：29+7+7=43 13+7+21=41
#列方向：29+13=42,7+7=14,7+21=28
#42+14+28=43+41=84
#第二行
#Chi-square contribution
#第三行
#相当于addmargins(prob.table(mytable,1),2)
#行方向: 0.674+0.163+0.163=1,0.317+0.171+0.512=1
#第四行
#相当于addmargins(prob.table(mytable,2),1)
#列方向: 0.690+0.310=1,0.500+0.500=1,0.250+0.750=1
#第五行
#相当于addmargins(prob.table(mytable))
#行方向：0.345+0.083+0.083=0.512,0.155+0.083+0.250=0.488
#列方向: 0.345+0.155=0.500,0.083+0.083=0.167,0.083+0.250=0.333
#0.512+0.488=0.500+0.167+0.333=1
CrossTable(Arthritis$Treatment,Arthritis$Improved)
#多维列联表
#原生成二维列联表的函数，直接推广到多维
#生成三维数据源，频数表
#Treatment 可理解成原始行维度，满足块条件下的统计
#sex 可理解成原始列维度，满足块条件下的统计
#Improved 只能是第三维度，块标识
mytable <- xtabs(~Treatment+Sex+Improved,data=Arthritis)
ftable(mytable) #区别于xtabs的多维表格输出格式，这个更紧凑些
#margin.table 单维度
margin.table(mytable,1) #根据变量1(行变量)， 多个块计算边际值，横向求和
margin.table(mytable,2) #根据变量2(列变量)， 多个块计算边际值，纵向求和
margin.table(mytable,3) #根据变量3(块变量)， 单块计算边际值，块内求和
#margin.table 超过单维度
margin.table(mytable,c(1,3)) #相当于忽略第二个维度，仅根据第1，第3的维度的条件，统计频数
ftable(prop.table(mytable,c(1,2))) #第三个维度仍为块维度，比较难，没太懂
ftable(addmargins(prop.table(mytable,c(1,2)),3)) #结合上面一行理解
ftable(addmargins(prop.table(mytable,c(1,2)),3))*100 #结合上一行理解
#独立性检验
#显著性检验评估了是否存在充分的证据以拒绝变量间相互独立的原假设
#卡方独立性检验
library(vcd)
mytable <- xtabs(~Treatment+Improved,data=Arthritis)
#卡方假设检验
#用于检验原假设
#范例一
#p-value p值，p值落入小概率范围(p<0.01)则拒绝原假设
#df 自由度
chisq.test(mytable)
#范例二
#p-value p值，p值未落入小概率范围(p>0.01)则接受原假设，两者相互独立
mytable <- xtabs(~Treatment+Sex,data=Arthritis)
chisq.test(mytable)
#fisher 精确检验
mytable <- xtabs(~Treatment+Improved,data=Arthritis)
fisher.test(mytable) #可以用于任意行数大于等于2的二维列联表上使用，但是2行，2列不行
#Cochran-Mantel-Haenszel检验
#其原假设：两个名义变量在第三个变量的每一层中都是条件独立的
#范例一
mytable <- xtabs(~Treatment+Improved+Sex,data=Arthritis)
mantelhaen.test(mytable)
#相关性的度量
#如果拒绝了相互独立的原假设，开始转向衡量相关性强弱和相关性度量
#二维列联表的相关性度量
#phi-coefficient: fi相关系数
#Contigency Coeff: 列联相关系数
#Cramer's V : V相关系数
#较大值意味着较强的相关性
library(vcd)
mytable <- xtabs (~Treatment+Improved,data=Arthritis)
assocstats(mytable)
#结果的可视化
#相关
#相关系数可以用来描述定量变量之间的关系。相关系数的符号表明关系的方向（正相关或负相关），
#其值的大小表示关系的强弱程度(完全不相关时0，完全相关时为1)
#Pearson，积差相关系数，衡量了两个定量变量之间的线性相关程度
#spearman，等级相关系数则衡量了分级定序变量之间的相关程度
#Kendall's Tau相关系数也是一种非参数的等级相关度量
#协方差
states <- state.x77[,1:6]
cov(states)
cor(states) #协方差，默认pearson系数
cor(states,method="spearman") #spearman相关系数
#相关系数矩阵默认为一个方阵，所有变量之间两两计算相关，也可以计算非方形的相关矩阵
x <- states[,c("Population","Income","Illiteracy","HS Grad")]
y <- states[,c("Life Exp","Murder")]
cor(x,y)
#偏相关
#控制一个或多个定量变量时，另外两个定量变量之间的相互关系
library(ggm)
colnames(states) #获取列名
#pcor(u,S),偏相关系数常用于社会科学的研究中
#u是一个数值向量，前两个数值表示要计算相关系数的变量下标，其余的数值为
#条件变量(即要排除影响的变量)，S为变量的协方差阵
pcor(c(1,5,2,3,6),cov(states)) #控制其他变量，求得Poputation与Murder之间的相关关系
#返回0.346
#相关性的显著性检验
#cor.test(x,y,alternative=,method=),每次只能检验一种相关关系
#x,y 要检验相关性的变量
#alternative 则用来指定进行双侧检验或者单侧检验，
#greater,当研究的总体的相关系数大于0时
#less，当研究的总体的相关系数小于0时
#two.side 总体相关系数不等于0
#原假设：检验了预期寿命和磨砂率的Pearson相关系数为0的原假设
# p-value = 1.258e-08,拒绝原假设，即预期寿命和磨砂率之间的总体相关度不为0
cor.test(states[,3],states[,5])
#psych包 corr.test函数
library(psych)
#Correlation matrix 相关矩阵，对称阵
#Probability values 显著性检验
corr.test(states,use="complete")
#独立样本的t检验
#一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设
t.test(y~x,data)
#范例一
library(MASS)
#原假设：南方各州和非南方各州拥有相同监禁概率的假设
#p-value = 0.0006506，p<0.01,落入拒绝区域
#正态化变换，Y/1-Y,log(Y/1-Y),arcsin(Y),arcsin(sqrt(Y))
t.test(Prob~So,data=UScrime)
#非独立样本的t检验
t.test(y1,y2,paired = TRUE)
library(MASS)
sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),sd=sd(x))))
#p-value < 2.2e-16，落入拒绝域
#可以拒绝年长和年轻男性的平均失业率相同的假设
with(UScrime,
t.test(U1,U2,paired=TRUE)
)
#如果想在多于两个的组之间进行比较，假设数据是从正态总体中独立抽样而得的，那你
#可以使用方差分析(ANOVA),如果数据无法满足t检验或者ANOVA的参数检验，可以转而使用
#非参数方法
#wilcoxon秩和检验，是非独立样本t检验的一种非参数替代方法
#试用于两组成对数据和无法保证正态性假设的情景
#范例一
with(UScrime,
by(Prob,So,median)
)
wilcox.test(Prob~So,data=UScrime)
#多于两组的比较
#kruskal-wallis检验
#kruskal.test(y~A,data) y是一个数值型结果变量，A是一个拥有两个或更多水平的分组变量
#范例一
states <- data.frame(state.region,state.x77)
#p-value =4.726e-05,落入拒绝域，拒绝文盲率相同的假设
kruskal.test(Illiteracy~state.region,data=states)
#范例二 InternetOpenUrl failed: '无法与服务器建立连接'
source("http://www.statmethods.net/RiA/wmc.txt")
#执行失败， InternetOpenUrl failed: '无法与服务器建立连接'
7497
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
