R语言实战-基本方法
一.基本图形
1.条形图
1.1 简单条形图 使用barplot(height,....),其中height 是向量或者矩阵，默认垂直条形图，horiz = TRUE（水平条形图），可以使用par() 函数对R的默认图形做出大量修改
barplot(height, width = 1, space = NULL,
names.arg = NULL, legend.text = NULL, beside = FALSE,
horiz = FALSE, density = NULL, angle = 45,
col = NULL, border = par("fg"),
main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
xlim = NULL, ylim = NULL,...)
1.2堆砌条形图和分组条形图
height 是矩阵时，barplot 绘制的是堆砌或者分组条形图，beside = FASLE是默认值，每一列生成一个条形，每个值是堆砌的“子条”高度。设置beside = TRUE,那么每一列的各个值相互并排形成一组
1.3棘状图
棘状图对堆砌图进行了重缩放，每个图形的高度都为1，每一段的高度表示这个值在这一列的比例，可以使用vcd包中的spine()绘制
> install.packages("vcd")
> library(vcd)
载入需要的程辑包：grid
> attach(Arthritis)
> counts <- table(Treatment,Improved)
> counts
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
> spine(counts,main ="Spinogram Example")
> detach(Arthritis)
2.饼图
使用相对较少，人们对长度的判读比面积更加精确，可以使用 pie(x,labels)创建
使用plotrix包中的pie3D()可以创建3维饼图
install.packages("plotrix")
> library(plotrix)
> data <- c(1,2,3,4)
> tag <- c("a","b","c","d")
> pie3D(data,labels = tag,explode = 0.1,main="3d chart")
3.核密度图和箱线图
一种用于观察连续型变量分布的有效方法，使用plot(x) ,使用lines()函数可以在已存在的图形中叠加一条曲线
使用sm包中的sm.density.compare()函数可以向图形叠加两组或者更多的核密度图，但用于可视化分布和组间差异的更好图例是还是箱线图
箱线图:通过绘制连续变量的五数总括（最大值、上四分位数（第75%位数）、中位数（第50%位数）、下四分位数（第25%分位数）、最小值），箱线图可以直观的看出离群点（范围+-（1.5）* IQR(上四分位减下四分位)）
使用 boxplot(formula,data=dataframe)
其中,formula是一个公式， 一个示例：y~ A,这将为类别型变量A的每个值并列地生成数值型变量y的箱线图，公式 y~ A*B 表示将类别型变量A和B所有水平的两两组合生成数值型变量y的箱线图 ，dataframe 是数据集
下面来以mtcars数据集为例
mtcars
如果想要直观的显示出自动挡和手动挡在分别在4、6、8缸发动机的耗油量数据，如何用箱线图表示？
#创建汽车缸数因子
> mtcars$cyl.f <- factor(mtcars$cyl,levels = c(4,6,8),labels = c("4","6","8"))
#创建变速箱类型因子
> mtcars$am.f <- factor(mtcars$am,levels = c(0,1),labels = c("auto","standard"))
#绘图
> boxplot(mpg ~am.f*cyl.f,data = mtcars,varwidth= TRUE,col=c("gold","darkgreen"),main="MPG distribution by auto type",xlab="auto type",ylab="Miles per gallon")
效果如下
不同变速箱类型和汽车缸数对油耗的影响
4.点图
使用dotchart(x,labels=)
x 是一个数值向量，labels 是每个点的标签组成的向量
二.基本统计分析
R语言有许多用于描述生成基本描述性和推断统计量
2.1.描述性统计量
2.1.1 summary()
summary()函数提供了最小值、四分位数和均值
> testCars <- c("mpg","wt","hp")
> summary(mtcars[testCars])
mpg wt hp
Min. :10.40 Min. :1.513 Min. : 52.0
1st Qu.:15.43 1st Qu.:2.581 1st Qu.: 96.5
Median :19.20 Median :3.325 Median :123.0
Mean :20.09 Mean :3.217 Mean :146.7
3rd Qu.:22.80 3rd Qu.:3.610 3rd Qu.:180.0
Max. :33.90 Max. :5.424 Max. :335.0
2.1.2 sapply()
sapply(x,FUN,options) 可以用于计算所选择的任意描述性统计量，x:数据框(矩阵)，FUN为任意函数，指定options,它们将传递给FUN,这里的FUN 函数可以是 mean()，sd()，var()，min()，max()，median()等等
可以使用sapply函数生成基础安装不带的描述性变量 eg skew :偏度 kurtosis ：峰度
myfuction <- function(x,na.omit = FALSE){
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))
}
myvars <- c("mpg","hp","wt")
sapply(mtcars[myvars],myfuction)
2.1.3 许多开源包中也提供了很多计算描述性统计量的函数
2.1.3.1 Hmisc 包中的describe() :
函数可返回 平均值、分位数、及5个最大最小值
2.1.3.2 pastecs包中的 stat.desc() ：
提供了种类繁多的描述性统计量，使用格式为：stat.desc(x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95)
install.packages("pastecs")
> library(pastecs)
> mydata <- c("mpg","hp","wt")
> stat.desc(mtcars[mydata])
mpg hp wt
nbr.val 32.0000000 32.0000000 32.0000000
nbr.null 0.0000000 0.0000000 0.0000000
nbr.na 0.0000000 0.0000000 0.0000000
min 10.4000000 52.0000000 1.5130000
max 33.9000000 335.0000000 5.4240000
range 23.5000000 283.0000000 3.9110000
sum 642.9000000 4694.0000000 102.9520000
median 19.2000000 123.0000000 3.3250000
mean 20.0906250 146.6875000 3.2172500
SE.mean 1.0654240 12.1203173 0.1729685
CI.mean.0.95 2.1729465 24.7195501 0.3527715
var 36.3241028 4700.8669355 0.9573790
std.dev 6.0269481 68.5628685 0.9784574
coef.var 0.2999881 0.4674077 0.3041285
2.1.4 分组计算描述性统计量
使用psych包中的describeBy()
install.packages("psych")
library(psych)
myvars <- c("mpg","hp","wt")
describeBy(mtcars[myvars],list(am=mtcars$am))
2.2 类别型变量的独立性检验
R中提供了许多独立性检测的方法，这里介绍三种： 1. 卡方独立性检验，2Fisher精确检验 , 3. Cochran-Mantel-haenszel 检验
首先将数据集生成频数表，训练数据集来自vcd 包中的Arthritis 数据集
library(vcd)
head(Arthritis)
ID Treatment Sex Age Improved
1 57 Treated Male 27 Some
2 46 Treated Male 29 None
3 77 Treated Male 30 None
4 17 Treated Male 32 Marked
5 36 Treated Male 46 Marked
6 23 Treated Male 58 Marked
使用table(A,B)生成二维列联表，A是行变量,B是列变量，或者使用xtabs()函数生成eg：
#mytable <- xtabs(~Treatment+Improved,data=Arthritis)
mytable <- xtabs(~ Arthritis$Treatment+Arthritis$Improved,data=Arthritis)
> mytable
Arthritis$Improved
Arthritis$Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
如果要生成三维列联表
调用 mytable <- xtabs (~ A+B+C，data = mydata)
2.2.1.卡方独立性检测
chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
mytable <- xtabs(~ Arthritis$Treatment+Arthritis$Sex,data=Arthritis)
> chisq.test(mytable)
Pearson's Chi-squared test with Yates' continuity correction
data: mytable
X-squared = 0.38378, df = 1, p-value = 0.5356
结果1中p-value=0.001463 <0.01 说明接收治疗和改善水平中存在某种关系，而患者接收治疗和性别中p-value = 0.5356 >0.05 说明他们没有关联性
2.2.2 Fisher精确检验
mytable <- xtabs(~ Arthritis$Treatment+Arthritis$Improved,data=Arthritis)
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided
注意：fisher.test()可以在任一行列数大于等于2的二维列联表上使用，但不能用用2x2的列联表
2.2.3 Cochran-Mantel-Haenszel 检验
注意，此条件检测两名义变量在第三变量的每一层都是条件独立的，不存在三阶交互作用
threetable <-xtabs(~ Arthritis$Treatment +Arthritis$Improved+Arthritis$Sex,data = Arthritis)
> mantelhaen.test(threetable)
Cochran-Mantel-Haenszel test
data: threetable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647
2.3 相关性检验
1.相关性的显著性检验
2.3.1使用col.test(x,y,alternative,method=)
x,y 为需要检验相关性的变量，alternative表示指定双侧检验还是单侧检验，col.test ()每次只能检验一种相关关系，使用psych包中的corr.test ()可以一次检测多个变量间的相关性
2.3.2 t检验
独立的检验
使用t.test(y~x,data) ,y是一个数值型变量，x 是一个二分量(只有两个取值)
#数据集是MASS包中的UScrime(刑法制度对犯罪率的影响)，几个重要参数
#Prob:监禁的概率 ；
#So :是否位于南方
#U1 :(14-25岁年龄段城市男性失业率)
#U2:(35~39岁年龄段城市男性失业率)
library(MASS)
t.test(UScrime$Prob~UScrime$So,data=UScrime)
Welch Two Sample t-test
data: UScrime$Prob by UScrime$So
t = -3.8954, df = 24.925, p-value = 0.0006506
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.03852569 -0.01187439
sample estimates:
mean in group 0 mean in group 1
0.03851265 0.06371269
非独立样本的t检验
判断较年轻男性U1是否比较年长男性U2更高
sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),sd=sd(x))))
U1 U2
mean 95.46809 33.97872
sd 18.02878 8.44545
> with(UScrime,t.test(U1,U2,paired = TRUE))
Paired t-test
data: U1 and U2
t = 32.407, df = 46, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
57.67003 65.30870
sample estimates:
mean of the differences
61.48936
我们可以发现，两个变量的均值的差异为61.5 足够大，可以保证拒绝年长和年轻男性的平均失业率相同的假设，年轻男性的失业率更高，要获得两个变量总体均值相同的概率微乎其微，概率为2.2e-16
