part 2 基本方法

-------------第6章 基本图形------------------

#1.条形图
#一般是类别型（离散）变量
library(vcd)
help(Arthritis) #类风湿性关节炎新疗法研究结果
head(Arthritis)
count <- table(Arthritis$Improved)
barplot(count,main="simple bar plot",xlab = "improvement",ylab = "frequency")
barplot(count,horiz = T)
plot(Arthritis$Improved) #类别型因子可直接plot画

#堆积bar
counts <- table(Arthritis$Improved,Arthritis$Treatment) #列联表
counts
barplot(counts,col=c("red","yellow","green"),legend=rownames(counts))
#分组bar
barplot(counts,legend=rownames(counts),beside = T)
#棘状图（对堆积图重缩放，总高度为1，按比例显示）
spine(counts)

#微调参数：
par(mar=c(5,8,4,2)) #调边界
par(las=2) #旋转条形标签
barplot(count,
    cex.names=0.8, #缩小字体
    names.arg=c("A","b","c") #修改标签文本
    )


#2.饼图
x <- c(12,12,14,18,4)
y <- c("a","b","c","d","e")
pie(x,labels=y,col=rainbow(length(y)),main="pie") #即rainbow(5)

#3D饼图
library(plotrix)
pie3D(x,labels=y,col=rainbow(length(y)),main="3D pie",explode = 0.1)
#扇形饼图
fan.plot(x,labels=y)


#3.直方图：展示连续型变量
x <- rnorm(100)
hist(x)
hist(mtcars$mpg)
hist(mtcars$mpg,
     breaks = 12, #指定组数
     freq = F, #根据概率密度显示，而非频数
     col="blue",
     xlab = "x",
     main = "histogram"
     )
rug(jitter(mtcars$mpg)) #添加轴须图，jitter向每个数据点添加一个小的随机值，避免点重叠！
lines(density(mtcars$mpg),col="red",lwd=2) #添加概率密度曲线
box() #添加外框


#4.核密度图：观察连续型变量分布
#用于估计随机变量概率密度函数的一种非参数方法
plot(density(x))
plot(density(mtcars$mpg),main = "density test")
polygon(density(mtcars$mpg),
        col="red", #填充颜色
        border="blue") #边界颜色
rug(mtcars$mpg)

#多组核密度图比较
library(sm)
attach(mtcars)
sm.density.compare(mpg,cyl,xlab="gallon")
title(main="test")
legend(locator(1),levels(factor(cyl,levels=c(4,6,8))),fill = c(2:4)) 
#locator(1)表用鼠标点击交互式放置图例
detach(mtcars)

#5.箱线图
#连续型变量5数总括：最小值、下四分位数、中位数、上四分位数、最大值
#离群点：±1.5*IQR外的值，IQR即四分位距（上四分位数-下四分位数）
#正偏：上侧须比下侧须要长，反之负偏
boxplot.stats(mtcars$mpg) #输出构建图形的统计量
boxplot(mpg~cyl,mtcars, #两者顺序不能反
        notch=T, #含凹槽，若两个箱之间互不重叠，表明其中位数有显著差异。
        varwidth=T, #箱线宽度与其样本大小的平方根成正比
        col="red",
        main="Car mileage data",
        xlab="number of cylinders",
        ylab="miles per gallon",
        horizontal=T #反转坐标
        )

#两个交叉因子的箱线图
mtcars$a <- factor(mtcars$cyl,levels = c(4,6,8)) #不同气缸
mtcars$b <- factor(mtcars$am,levels = c(0,1)) #不同变速箱
boxplot(mpg~a*b,
        mtcars,
        varwidth=T,
        col=c("gold","darkgreen")
        )

#6.小提琴图
#小提琴图是箱线图和核密度图的结合
library(vioplot)
attach(mtcars)
vioplot(mpg[cyl==4],mpg[cyl==6],mpg[cyl==8], #将不同组分离到不同变量
        names = c("4cyl","6cyl","8cyl"),
        col = "gold"
        )
detach(mtcars)
#白点中位数，黑盒上四分位数到下四分位数，细黑线为须，外部形状即为核密度估计（镜像叠加）

#7.点图
dotchart(mtcars$mpg,labels = row.names(mtcars),cex = 0.7,main = "dotplot",xlab = "gallon")
#一般点图在经过排序、分组被不同符号和颜色区分开后最有用
x <- mtcars[order(mtcars$mpg),]
x$cyl <- factor(x$cyl)
x$color[x$cyl==4] <- "red"
x$color[x$cyl==6] <- "blue"
x$color[x$cyl==8] <- "darkgreen"
dotchart(x$mpg,
         labels = row.names(x),
         groups = x$cyl, #数据点分组
         gcolor = "black", #数字颜色
         color = x$color, #点和标签颜色
         cex=0.7,
         pch=19,
         man="dot\nplot test",
         xlab = "miles per gallon"
         )

-----------------第7章 基本统计分析-------------------------------

#参数检验：t检验
#非参检验：U检验，KW检验

#1.描述性统计分析
#方法一：
summary(mtcars[,c("mpg","hp","wt")])

#方法二：
#lapply/sapply结合FUN：mean/sd/var/min/max/median/length/range/quantile
mystats <- function(x,na.omit=F){
  if(na.omit)
    x <- x[!is.na(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(mtcars[,c("mpg","hp","wt")], mystats)


#方法三：
fivenum(mtcars$mpg) #不能同时指定多列

#其他方法：
library(Hmisc)
describe(mtcars[,c("mpg","hp","wt")]) #返回变量名、变量数目、缺失值、唯一值数目、平均值、分位数、最小和最大的五个数

library(pastecs)
stat.desc(mtcars[,c("mpg","hp","wt")]) #返回更多信息

library(psych)
describe(mtcars[,c("mpg","hp","wt")])
#不同包的同名函数，最后载入的包为先，如这里Hmisc包的describe被屏蔽（masked）
#最好用Hmisc::describe避免被屏蔽

#分组计算描述统计量:
#方法1：
aggregate(mtcars[,c("mpg","hp","wt")],
          by=list(am=mtcars$am), #若不指定列名，将默认为Group.1
          mean)
##只能用mean、sd等单返回值函数，一次无法返回多个统计量

#方法2：
by(mtcars[,c("mpg","hp","wt")],
   mtcars$am,
   function(x) sapply(x,mystats))
#一次返回多个统计量

#其他方法：
library(doBy)
summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)

psych::describeBy(mtcars[,c("mpg","hp","wt")],list(am=mtcars$am))

#2.频数表和列联表
library(vcd)
head(Arthritis)

#一维列联表
table(Arthritis$Improved)
prop.table(table(Arthritis$Improved)) #将table比例化
prop.table(table(Arthritis$Improved))*100

#二维列联表
table(Arthritis$Treatment,Arthritis$Improved) #默认忽略NA
mytable <- xtabs(~ Treatment+Improved,data = Arthritis);mytable
##边际频数
margin.table(mytable,1) #行
margin.table(mytable,2) #列
##比例
prop.table(mytable,1) #行
prop.table(mytable,2) #列
##添加边际之和
addmargins(mytable)
addmargins(prop.table(mytable))
addmargins(prop.table(mytable),1)
addmargins(prop.table(mytable),2)

#多维列联表
mtable <- xtabs(~ Treatment+Sex+Improved,data = Arthritis);mtable
ftable(mtable) #易读的方式
margin.table(mtable,1)
margin.table(mtable,2)
margin.table(mtable,3)
margin.table(mtable,c(1,3))
ftable(addmargins(prop.table(mtable,c(1,2)),3))*100

#独立性检验：原假设都是独立
##1）卡方检验
chisq.test(mytable) #p<0.05不独立（存在关系）
chisq.test(margin.table(mtable,c(1,2))) #p>0.05独立

##2）.Fisher精确检验
fisher.test(mytable) #p<0.05不独立，但不能用于2x2列联表

##3）.Cochran-Mantel-Haenszel检验（CMH卡方检验）
mantelhaen.test(mtable) #用于三维列联表

#相关性度量
#以上两个变量不独立（存在关系），度量相关性
assocstats(mytable) #越大相关性越强

#3.相关
head(state.x77)
help("state.x77")

##相关的类型：
#person：两个定量变量
#spearman：分级定序变量
#kendall：非参
cor(x,use ="everything" ,method ="pearson" ) #默认
#use指定缺失值处理方式：all.obs(不允许NA，否则报错),everything(missing), complete.obs(行删除), pairwise.complete.obs(成对删除)

cor(state.x77[,1:6])
cor(state.x77[,1:6],method = "spearman")
cor(state.x77[,c(1,2,3,6)],state.x77[,c(4,5)]) #两组变量
cov(state.x77[,1:6]) #协方差

#偏相关：控制一个或多个变量时，其他两个变量间的关系
library(ggm)
colnames(state.x77[,1:6])
pcor(c(1,5,2,3,6),cov(state.x77[,1:6])) #第一个参数为数值向量，前两个是待计算的变量下标，其他是条件变量下标

##相关性的显著性检验
cor.test(x,y,alternative = "two.side",method = "pearson") #默认
#alternative参数：less(研究假设总体相关系数小于0时)，greater(研究假设总体大于0时)，two.side
cor.test(state.x77[,3],state.x77[,5])
mycor <- cor.test(state.x77[,3],state.x77[,5])
str(mycor)
mycor$ estimate
mycor$p.value
#缺点：一次只能检验一种相关关系
cor.test(state.x77[,1:6]) #error

library(psych)
corr.test(state.x77[,1:6],use = "complete")
#一次检验多种变量之间的关系,返回两个cor和pvalue两个矩阵
corr.test(state.x77[,3],state.x77[,5],use = "complete")
#use参数：pairwise/complete, method参数：pearson(default)/spearman/kendall

#检验在控制一个或多个变量时，其他两个变量间的独立性
psych::pcor.test(r,q,n)


##4.t检验
library(MASS)
head(UScrime)
help("UScrime")

###1)独立样本
t.test(y~x,data) #y为数值型变量，x为二分变量
t.test(y1,y2) #y1,y2都为数值型变量
t.test(Prob~So,data=UScrime)

###2)非独立样本
#两个变量之间相关
t.test(y1,y2,paired = T)
sapply(UScrime[c("U1","U2")], function(x)(c(mean=mean(x),sd=sd(x))))
with(UScrime,t.test(U1,U2,paired = T))

#多组比较anova参考后续章节

##5.组间差异的非参检验
###1)两组
#独立样本：wilcoxon秩和检验（U检验）
wilcox.test(y ~ x, data) #y是数值，x是二分变量
with(UScrime,by(Prob,So,median))
wilcox.test(Prob~So,data = UScrime)

#非独立样本：wilcoxon符号秩检验(适用于成对数据和非正态情况)
sapply(UScrime[c("U1","U2")], median)
with(UScrime,wilcox.test(U1,U2,paired = T))

#当t检验的假设合理时，参数检验的功效更强，更容易发现差异；非参数检验在假设不合理时更适用，如等级有序数据。

###2)多组
#one-way design 单向设计，如美国四个地区的文盲率比较。
#各组独立：Kruskal-Wallis检验
kruskal.test(y~A,data) #y为数值变量，A为两个及其以上水平的分组变量（两个水平等于U检验）
states <- data.frame(state.region,state.x77)
kruskal.test(Illiteracy~state.region,data = states)
#p<0.01表明四个地区的文盲率各不相同，但无法知道具体哪些地区显著不同（只能每次比较两组）

#所有组之间成对比较，可以用本书作者写的wmc函数：
source("http://www.statmethods.net/RiA/wmc.txt") #函数下载
wmc(Illiteracy~state.region,data = states,method = "holm")
wmc #查看源码

#各组不独立(如重复测量设计或随机区间设计)：Friedman检验
friedman.test(y~A|B,data) #A为分组变量，B为匹配变量的区组变量