# R in Action: Chapter 7
#7.1 描述性统计
summary()     #获取描述性统计
apply()
sapply( x , fun , options)     #宽mean/sd/var/min/max/median/length/range/quantile
fivenum()
aggregate( x , by = list(name1=groupvar1,name2=groupvar2,.....) , fun)     
#一维用pastecs包里的stat.desc()直观点。
library(pastecs)
stat.desc(mtcars[vars], basic=T,desc=T, norm=T,p=0.95)   
#若basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和。若desc=TRUE（默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。最后，norm=TRUE（不是默认），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro–Wilk正态检验结果。这里使用了p值来计算平均数的置信区间（默认置信度为0.95）。
#reshape     分组进行统计
library(reshape)
dstats <- function(x) (c(n = length(x), mean = mean(x), sd = sd(x)))
dfm <- melt(mtcars, measure.vars = c("mpg", "hp", "wt"), id.vars = c("am", "cyl"))     #melt()融合为长数据型
#dfm1<-melt(mtcars, id = c("am", "cyl" ))
#dfm和dfm1的区别，measure.vars可以指明要进行改善的数值型变量（默认使用所有变量），dfm1是默认的。
cast(dfm, am + cyl + variable ~ ., c(length,mean,sd))


#7.2 二维频数表和列联表
table(x1, x2, .....)     #方便后续调用
xtable(~x1+x2+...)     #以公式风格形式
prop.tabel()     #对已求得的table结果进行频率转换


library(vcd)
options(digits = 2)
mytable <- with(Arthritis, table(Improved))     #频数统计表
prop.table(mytable)     #频率统计表，all
prop.table(mytable,1)     #频率统计表，行
prop.table(mytable,2)     #频率统计表，列
prop.table(mytable)*100     #百分比统计表


mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
margin.table(mytable, 1)     #生成边际频数,以行数值计算
prop.table(mytable, 1)     #生成比例，以行数值计算
margin.table(mytable, 2)     #生成边际频数,以列数值计算
prop.table(mytable, 2)     #生成比例，以列数值计算
prop.table(mytable)     #各单元所占比例


addmargins(mytable)     #添加边际和
addmargins(prop.table(mytable))
addmargins(prop.table(mytable, 1), 2)      #添加边际和,以行数值计算
addmargins(prop.table(mytable, 2, 1)     #添加边际和，以列数值计算
#table()默认忽略缺失值(NA)。要在频数统计中将NA视为一个有效的类别，设定useNA="ifany"。即计算时候是否添加空值


#二维的gmodels包中的crosstable结果直观点，format用sas
CrossTable(x, y, digits=3, max.width = 5, expected=FALSE, prop.r=TRUE, prop.c=TRUE,
           prop.t=TRUE, prop.chisq=TRUE, chisq = FALSE, fisher=FALSE, mcnemar=FALSE,
           resid=FALSE, sresid=FALSE, asresid=FALSE,
           missing.include=FALSE,
           format=c("SAS","SPSS"), dnn = NULL, ...)
library(gmodels)
with(Arthritis ,CrossTable(Treatment, Improved))


data(infert, package = "datasets")
CrossTable(infert$education, infert$induced, expected = TRUE)     #带期望值和皮尔森卡方检验
CrossTable(infert$education, infert$induced, expected = TRUE, format="SAS")  


#卡方独立性检验。对二维表
library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
chisq.test(mytable)
CrossTable(Arthritis$Treatment, Arthritis$Improved，expected = T)     #方法2，gmodels包
fisher.test(mytable)     #原假设为边界固定的列联表中行和列是独立的，用于任意行数列数大于等于2的二维列联表上，理论频数有小于5的。
#table()和xtabs()都可以基于3个或更多的类别型变量生成多维列联表
mytable <- xtabs(~Improved+Sex, data=Arthritis)
chisq.test(mytable)
CrossTable(Arthritis$Improved, Arthritis$Sex,expected = T ,format= "SAS")     
#H0相互独立；H1相互不独立。当p<0.05,相互不独立，存在关系；p>0.05,相互独立，不存在关系

• 如果个别字段的期望次数太低，会使概率分配无法近似于卡方分配。一般要求：自由度 {\displaystyle df>1}1" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/96282a7a0f25c2ae7ed911273a0301f80ae023d5" style="border: none; vertical-align: -0.671ex; display: inline-block; width: 6.797ex; height: 2.509ex;">时，期望次数小于5的字段不多于总字段的20%。
• 若自由度 {\displaystyle df=1}，且若期望次数 {\displaystyle <10} ，则近似于卡方分配的假设不可信。此时可以将每个观察值的离差减去 {\displaystyle 0.5} 之后再做平方，这便是叶氏连续性修正。
  

# Listing 7.14 - Measures of association for a two-way table
#相关性度量。值越大，相关性越强。phi系数、列联系数和Cramer's V系数
           library(vcd)
           mytable <- xtabs(~Treatment+Improved, data=Arthritis)
           assocstats(mytable)


#宽转长最原始列表


           table2flat <- function(mytable) {
             df <- as.data.frame(mytable)
             rows <- dim(df)[1]
             cols <- dim(df)[2]
             x <- NULL
             for (i in 1:rows) {
               for (j in 1:df$Freq) {
                 row <- df[i, c(1:(cols - 1))]
                 x <- rbind(x, row)
               }
             }
             row.names(x) <- c(1:dim(x)[1])
             return(x)
           }


  # Listing 7.16 - Using table2flat with published data
           treatment <- rep(c("Placebo", "Treated"), 3)
           improved <- rep(c("None", "Some", "Marked"), each = 2)
           Freq <- c(29, 13, 7, 7, 7, 21)
           mytable <- as.data.frame(cbind(treatment, improved, Freq))
           mydata <- table2flat(mytable)
           head(mydata)


#7.3 相关
求相关系数：cor()、pcor()
#pearson积差相关系数衡量了两个定量变量之间的线性相关程度。spearman等级相关系数则衡量分级定序变量之间的相关程度。kendall’s Tau相关系数是一种非参数的等级相关度量。
#cor( x , use = , method = )  
#指定缺失数据的处理方式。可选的方式为all.obs（假设不存在缺失数据——遇到缺失数据时将报错）、everything（遇到缺失数据时，相关系数的计算结果将被设为missing）、complete.obs（行删除）以及 pairwise.complete.obs（成对删除，pairwise deletion）
#指定相关系数的类型。可选类型为pearson、spearman或kendall


# Listing 7.17 - Covariances and correlations     
           states <- state.x77[, 1:6]
           cov(states)     #协方差
           cor(states)     #person积差相关系数
           cor(states, method="spearman")     #spearman登记相关系数
           cor(x,y)     #求两组变量之间的相关系数


#偏相关。指在控制一个或多个定量变量时，另外两个定量变量之间的相互关系。
#pcor(u , s)     #u是一个数值向量，前2个数值表示要计算相关系数的变量下标，其余的数值为条件变量（即要排除影响的变量）的下标，s为变量的协方差阵。
library(ggm)
pcor(c(1, 5, 2, 3, 6), cov(states))


#相关性显著性检验
检验相关系数：cor.test()、corr.test()、 pcor.test()、r.test()
cor.test(x,y, alternative= , method = )     #默认原假设为变量间不相关，0。alternative:<0用less，>0用greater，不等于0用two.side。method:person、spearman、kendall
#计算相关矩阵并进行显著性检验
library(psych)
corr.test(states, use = "complete")     #生成两个矩阵，第一个是相关系数，第二个是检验p值

#参数use=的取值可为"pairwise"或"complete"（分别表示对缺失值执行成对删除或行删

除）。参数method=的取值可为"pearson"（默认值）、"spearman"或"kendall"
pcor.test(r, q , n) #检验偏相关，r是pcor函数计算得到的偏相关系数，q为要控制的变量数（以数值表示位置），n为样本大小。
r.test()#help(r.test)

#7.4 t检验
假设两组数据是独立的，并且是从正态总体中抽得的，检验两总体均值相等的假设用t检验。
t.test( y ~ x , data) #y是数值变量，x是二分变量。
t.test( y1 , y2 )     #y1和y2为数值型向量。
t检验默认方差不相等，并使用Welsh的修正自由度，可以添加var.equel = T 假定方差相等。
alternative="less",alternative="greater"
#非独立样本t检验
前后测设计、重复测量设计等均为非独立。条件：假定组间的差异呈正态分布：
t.test(y1 , y2 ,paired = T)


#7.5 非参数检验
如果数据无法满足t检验或ANOVA的参数假设，可以转用非参数方法。
#7.5.1 两组比较（不独立，不正态），
两组独立，使用Wilcoxon秩和检验（U检验）评估观测是否是从相同的概率分布中抽得的
wilcox.test( y ~x , data);wilcox.test(y1,y2)
两组配对检验情况需增加paired=T
#7.5.2 多于两组的比较（不正态）
无法满足ANOVA设计的假设。
如果各组独立，用Kruskal-Wallis检验：kruskal.test( y~A , data)，A为分组变量
如果各组不独立，用Friedman检验：friedman.test( y ~A | B , data),A为分组变量，B为一个用以认定匹配观测的区组变量