R语言之统计专题

作者: 天秤座的机器狗 | 来源:发表于2018-09-16 15:39 被阅读36次

图片资源来自B站之阿雷边学边教

1.常用的两种抽样方法

image.png

1.随机抽样(sample函数)

data <- read.table("ehbio.DESeq2.all.DE.entrez",sep="\t") #导入差异基因
data <- data[order(data$V2),]#分类排序
index <- sample(1:nrow(data),80) #根据差异基因文件的行数，取这个范围内0.1比例的随机数
data_sample <- data[index,] #根据随机数代表的行数，从data中把数据提取出来

data

image.png

index

image.png

data_sample

image.png

2.分层抽样


install.packages("sampling") #安装sampling包
library(sampling) #加载包

data <- read.table("ehbio.DESeq2.all.DE.entrez",sep="\t") #导入差异基因
data <- data[order(data$V2),]#分类排序
colnames(data)=c("counts","sample") #更改列名
summary(data)  #查看data统计数据
# counts                             sample   
# Min.   :       12   untrt._higherThan_.trt:383  
# 1st Qu.:     5232   untrt._lowerThan_.trt :423  
# Median :    23030                               
# Mean   :  1203176                               
# 3rd Qu.:    84509                               
# Max.   :105369247      

a <- strata(data,stratanames = "sample",size = c(38,42),method = "srswor") #使用sampling中的strata函数
                                                          #stratanames就是分层的名字，size就是每层取得个数
  
#，“srswor”就是不放回抽样   

data_strata <- data[a$ID_unit,] #根据行号提取数据

a的结果

image.png

data_strata的结果

image.png

2.描述性统计

image.png

1.summary函数
（1）数据源单列情况

a <- summary(mtcars$cyl)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 4.000   4.000   6.000   6.188   8.000   8.000

最小值、第一4分位数、中位数、平均值、第三4分位数、最大值

比如，如果要提取min

a[[1]]
#[1] 4

（2）数据源多列情况

b <- mtcars[c("cyl","wt","hp")]
d <- summary(b)
# cyl              wt              hp       
# Min.   :4.000   Min.   :1.513   Min.   : 52.0  
# 1st Qu.:4.000   1st Qu.:2.581   1st Qu.: 96.5  
# Median :6.000   Median :3.325   Median :123.0  
# Mean   :6.188   Mean   :3.217   Mean   :146.7  
# 3rd Qu.:8.000   3rd Qu.:3.610   3rd Qu.:180.0  
# Max.   :8.000   Max.   :5.424   Max.   :335.0

比如，如果要提取cyl中的min值

e <- d[,1][1]
# "Min.   :4.000  " 
f <- strsplit(e,":")[[1]][2]
as.numeric(f)
#[1] 4

2.pastecs包的stat.decs函数来做描述性统计

install.packages("pastecs")
library("pastecs")
stat.desc(b)
 # cyl          wt           hp
# 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            4.0000000   1.5130000   52.0000000
# max            8.0000000   5.4240000  335.0000000
# range          4.0000000   3.9110000  283.0000000
# sum          198.0000000 102.9520000 4694.0000000
# median         6.0000000   3.3250000  123.0000000
# mean           6.1875000   3.2172500  146.6875000
# SE.mean        0.3157093   0.1729685   12.1203173
# CI.mean.0.95   0.6438934   0.3527715   24.7195501
# var            3.1895161   0.9573790 4700.8669355
# std.dev        1.7859216   0.9784574   68.5628685
# coef.var       0.2886338   0.3041285    0.4674077

比如，要提取cyl的mean值

g[,1][4]
#[1] 4

3.table计数
先看一下数据集

> mtcars$cyl
 [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4

table一下

> table(mtcars$cyl)

 4  6  8 
11  7 14

结果就是4出现11次
6出现7次
8出现14次

3.正态分布和抽样分布

image.png

1.什么是正态分布

image.png

在R中模拟

library(ggplot2)
#模拟生成一个正态分布，均值为7000，标准差为2000，数目为10000的男性毕业生收入
data1 <- rnorm(10000,mean=7000,sd=2000)
data1 <- as.data.frame(data1)
data1[,2]="male"
names(data1) <- c("income","gender")
p1 <-ggplot(data1,aes(income))+geom_histogram(stat="density")


#模拟生成一个正态分布。均值为5000，标准差为2000，数目为10000的女性毕业生收入
data2 <- rnorm(10000,mean=5000,sd=2000)
data2 <- as.data.frame(data2)
data2[,2]="female"
names(data2) <- c("income","gender")

newdata <- rbind (data1,data2)



p2 <- ggplot(newdata,aes(income))+geom_histogram(stat="density")+facet_grid(gender~.)

第一个图

image.png

第二个图

image.png

可见μ值越大，图像越向右偏移

如果将男性女性均值设定成一样，男性标准差为2000，女性为3000
则

library(ggplot2)
#模拟生成一个正态分布，均值为7000，标准差为2000，数目为10000的男性毕业生收入
data1 <- rnorm(10000,mean=7000,sd=2000)
data1 <- as.data.frame(data1)
data1[,2]="male"
names(data1) <- c("income","gender")



#模拟生成一个正态分布。均值为7000，标准差为3000，数目为10000的女性毕业生收入
data2 <- rnorm(10000,mean=7000,sd=3000)
data2 <- as.data.frame(data2)
data2[,2]="female"
names(data2) <- c("income","gender")

newdata <- rbind (data1,data2)



p<- ggplot(newdata,aes(income))+geom_histogram(stat="density")+facet_grid(gender~.)

image.png

可见，标准差越大，曲线越平坦

2.检验正态分布的三种方法

image.png

（1）qq图

qqnorm(data1[,1])

image.png

qqline(data1[,1],col="red")

image.png

如果所有数据的点都落在这条直线上，就说明这个数据符合正态分布
直观，但是不精确，毕竟是肉眼看

(2).夏皮罗检验（shapiro.test）
当w接近1，p>0.05时，说明数据符合正态分布

> shapiro.test(data1[,1])
Error in shapiro.test(data1[, 1]) : 样本大小必需在3和5000之间

所以，必须要对数据抽样，再使用夏皮罗检验

shapiro.test(sample(data1[,1],5000))
# Shapiro-Wilk normality test
# 
# data:  sample(data1[, 1], 5000)
# W = 0.99956, p-value = 0.3166

因为涉及到抽样
所以，希望能多抽取几次，看看情况
这就需要写一个循环，并且将p值保存（如果只想看看p值的话）
注意shapiro检验得到的结果是列表形式的


save_data <- c()
for(i in 1:1000){
  result_test <- shapiro.test(sample(data1[,1],5000))
  save_data <- append(save_data,result_test[[2]])
}
length(save_data[save_data<=0.05])
#[1] 20

1000次才有20次不符合，还是可以认为是符合正态分布的

（3）.ks检验，D接近0，且p>0.05时，数据符合正态分布
ks.test()除了要验证的数据之外，还需要生成一个与待验证数据的均值和标准差以及数量一致的符合正态分布的标准数据集

 ks.test(data1[,1],rnorm(10000,mean=mean(data1[,1]),sd=sd(data1[,1])))

# Two-sample Kolmogorov-Smirnov test
# 
# data:  data1[, 1] and rnorm(10000, mean = mean(data1[, 1]), sd = sd(data1[, 1]))
# D = 0.0106, p-value = 0.628
# alternative hypothesis: two-sided

3.关于正态分布的4个函数及其应用

image.png

已知某大学男性毕业生收入的平均值为7000，标准差为2000
求：
（1）若甲同学的收入大于80%的人，那么他的收入是多少？（由概率求值的问题）

> qnorm(0.8,mean=7000,sd=2000)
[1] 8683.242

（2）乙同学的收入为8500左右的概率是多少？（点概率）

> dnorm(8500,mean=7000,sd=2000)
[1] 0.0001505687

（3）已知丙同学的收入为9000，他的收入会比百分之多少的人高？（区间概率问题）

> pnorm(9000,mean=7000,sd=2000)
[1] 0.8413447

这么看来qnorm和pnorm是一对

概率密度函数是指某个点对应一个概率，而累计概率密度函数（分布函数）是指区间的概率累加。为什么要定义分布函数，是因为在很多情况下，我们并不想知道在某样东西在某个特定的值的概率，顶多想知道在某个范围的概率，于是，就有了分布函数的概念。

4.标准正态分布

image.png

其实R里面可以直接实现，不用查表
这里体现的是区间概率，所以用的是pnorm

> pnorm(0.95,mean=0,sd=1)
[1] 0.8289439

image.png

其实就是将X标准化

zdata=(data1[,1]-mean(data1[,1]))/sd(data1[,1])
zdata=as.data.frame(zdata)
ggplot(zdata,aes(zdata))+geom_histogram(stat="density")

image.png

这样就将data1数据转变成的标准正态分布

5.抽样分布

image.png

6.中心极限定理

image.png

在R中验证一下：
（1）当总体的分布服从正态分布时，任意样本，无论样本容量多大，样本均值的抽样分布都服从正态分布

c_data <- c()
for (i in 1:1000){
  sample_data <- sample(data1[,1],30)
  c_data <-append(c_data,mean(sample_data))
}
c_data <-as.data.frame(c_data)
ggplot(c_data,aes(c_data))+geom_histogram(stat="density")

#因为夏皮洛检验需要的是向量，所以将c_data转换成向量
c_data <- as.matrxi(c_data)
c_data <- as.vector(c_data)


# shapiro.test(c_data)
# Shapiro-Wilk normality test
# 
# data:  c_data
# W = 0.99798, p-value = 0.2753

image.png

需要注意的是，向量转化成数据框，可以直接就是as.data.frame()
但是！
数据框转成向量，需要想转换成矩阵，再转换成向量
as.matrix()--------->as.vector()
(2).当总体的分布不服从正态分布时，只要样本的容量足够大，样本均值的抽样分布就会服从近似服从正态分布

data <- read.table("ehbio.DESeq2.all.DE.entrez",sep="\t")
library(ggplot2)
ggplot(data,aes(data$V1))+geom_histogram(stat="density")


save_data <- c()

for(i in 1:10000){
  sample_data <- sample(data[,1],30)
  save_data <- append(save_data,mean(sample_data))
}


ggplot(data=NULL,aes(save_data))+geom_histogram(stat="density")

image.png

但是如果将样本量改为20，则

image.png

4.参数估计和假设检验

image.png

（1）参数和统计量

image.png

所以参数肯定就是对总体的描述
统计量肯定就是对样本的描述

（2）参数估计

image.png

#生成一个-3到3的序列
x <- seq(-3,3,length.out = 50)
y <- dnorm(x,0,1)

#绘图
ggplot(data=NULL,aes(x,y))+geom_line(col="red")

#置信水平
c_level <- 0.95

#置信水平下的占比区间
a1 <- qnorm((1-c_level)/2)
# > a1
# [1] -1.959964

ggplot(data=NULL,aes(x,y))+geom_line(col="red")+geom_vline(xintercept = 0)+geom_vline(xintercept = a1)+
  geom_vline(xintercept = -a1)+
  annotate("text",-2.23,0.01,label=(1-c_level)/2)+ annotate("text",2.23,0.01,label=(1-c_level)/2)



#求出置信区间
paste("置信区间:[",a1,",",-a1,"]",sep="")
#[1] "置信区间:[-1.95996398454005,1.95996398454005]"

这样就声成了一个显著性水平为0.05下的置信区间

（3）各类均值的差别

image.png

（4）用样本方差来估计总体方差

image.png

（5）假设检验

image.png

在R中来实现

#p值计算
#z统计量为正数时
z <-0.8

#对应的单向累计概率
1-pnorm(z,0,1)
#[1] 0.2118554

#对应的双向累计概率，即p值
2*(1-pnorm(z,0,1))
#[1] 0.4237108


#z统计量为负数的时候
z <- -2.2
#对应的单向累计概率
pnorm(z,0,1)
#[1] 0.01390345
#对应的双向累计概率，即p值
2*pnorm(z,0,1)
#[1] 0.0278069



#假设检验
data <- read.csv("lunzhou.csv")
x <- mean(data[,1])
#[1] 40.33867
z <- (x-40)/(2/sqrt(30))
#[1] 0.9274769
p <- (1-pnorm(z))*2
#[1] 0.353679

5.t检验

image.png

1.疑惑与问题

image.png

2.总体均值的t检验

image.png

（1）单样本t检验

data <- read.csv("单样本t检验.csv")
#原假设H0： μ=5000 备择假设H1： μ!=5000
t.test(data$income,alternative = "two.sided",mu=5000)
# One Sample t-test
# 
# data:  data$income
# t = 0.11681, df = 27, p-value = 0.9079
# alternative hypothesis: true mean is not equal to 5000
# 95 percent confidence interval:
#   4544.461 5510.539
# sample estimates:
#   mean of x 
# 5027.5 



#结论p>0.05,接受原假设，也就是说明在0.05的显著性水平下，有证据证明
#该校该专业的毕业生的平均收入和该校毕业生的平均收入标准5000元之间没有显著差异

（2）两独立样本t检验

#两独立样本t检验
#原假设H0：μ1=μ2   备择假设H1：μ1!=μ2
#做法1
data2 <- read.csv("独立样本t检验.csv")
t.test(income~gender,data2,var.equal=TRUE)

# Two Sample t-test
# 
# data:  income by gender
# t = 0.7564, df = 54, p-value = 0.4527
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   -393.0071  869.2214
# sample estimates:
#   mean in group 男 mean in group 女 
# 5027.500         4789.393 


#做法2

data2 <- split(data2,data2$gender)
male <- data2$男
female <- data2$女
t.test(male$income,female$income,var.equal = TRUE)
# Two Sample t-test
# 
# data:  male$income and female$income
# t = 0.7564, df = 54, p-value = 0.4527
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   -393.0071  869.2214
# sample estimates:
#   mean of x mean of y 
# 5027.500  4789.393 


#结论p>0.05,接受原假设，也就是说明在0.05的显著性水平下，有证据证明
#该校该专业的女生毕业生的平均收入和该校男生毕业生的平均收入之间没有显著差异

（3）配对样本t检验

#配对样本t检验
#原假设H0： μ1=μ2   备择假设H1：μ1!=μ2

#做法1
> data3 <- read.csv("配对样本T检验.csv")
> View(data3)
> t.test(data3$培训前销售额,data3$培训后销售额,paired=TRUE)

# Paired t-test
# 
# data:  data3$培训前销售额 and data3$培训后销售额
# t = -8.1428, df = 27, p-value = 9.553e-09
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -1673.942 -1000.130
# sample estimates:
# mean of the differences 
# -1337.036 


#做法2
t.test(data3$培训前销售额-data3$培训后销售额,mu=0)
#   One Sample t-test
# 
# # data:  data3$培训前销售额 - data3$培训后销售额
# # t = -8.1428, df = 27, p-value = 9.553e-09
# # alternative hypothesis: true mean is not equal to 0
# # 95 percent confidence interval:
# # -1673.942 -1000.130
# # sample estimates:
# # mean of x 
# # -1337.036 


#结论p<0.05,拒绝原假设，也就是说明在0.05的显著性水平下，有证据证明
#培训前的销售额和培训后的销售额有显著差异

3.回首t检验

image.png

6.方差分析

image.png

1.学习目标

image.png

2.什么是方差分析

image.png

3.为什么学习方差分析

image.png

多组之间用t检验，犯错概率大，因此不适合

4.方差分析相关概念阐述

image.png

5.单因素方差分析

image.png

#单因素方差分析
#导入数据
data <- read.csv("data16.csv")

#步骤1 对数据进行初步了解
#查看因素的水平及其频数分布
table(data$培训方案)
# > table(data$培训方案)
# 
# A      B 无培训 
# 15     15     15 

#看变量数和观测数
str(data)
# > str(data)
# 'data.frame': 45 obs. of  3 variables:
#   $ 销售员  : int  1 2 3 4 5 6 7 8 9 10 ...
# $ 培训方案: Factor w/ 3 levels "A","B","无培训": 1 1 1 1 1 1 1 1 1 1 ...
# $ 销售业绩: int  5628 5202 5500 5692 5100 5446 5640 5618 5965 5990 ...



#查看各组均值和标准差
aggregate(销售业绩~培训方案,data,mean)
aggregate(销售业绩~培训方案,data,sd)


# > aggregate(销售业绩~培训方案,data,mean)
# 培训方案 销售业绩
# 1        A 5553.600
# 2        B 7528.933
# 3   无培训 4523.267
# > aggregate(销售业绩~培训方案,data,sd)
# 培训方案 销售业绩
# 1        A 299.2481
# 2        B 315.3797
# 3   无培训 324.2286

#也可以看看数据的图形分布
library(ggplot2)
ggplot(data,aes(销售员,销售业绩,color=培训方案))+geom_point()
       
ggplot(data,aes(销售员,销售业绩,color=培训方案))+geom_boxplot()




#步骤2：方差齐性检验
#方差齐性检验
bartlett.test(销售业绩~培训方案,data)
#左边写因变量，右边写自变量

# > bartlett.test(销售业绩~培训方案,data)
# 
# Bartlett test of homogeneity of variances
# 
# data:  销售业绩 by 培训方案
# Bartlett's K-squared = 0.089252, df = 2, p-value = 0.9564
#说明方差没有显著性差异


#步骤3：方差分析
data_aov <- aov(销售业绩~培训方案,data)
summary(data_aov)
# > summary(data_aov)
# Df   Sum Sq  Mean Sq F value Pr(>F)    
# 培训方案     2 69987803 34993902   356.9 <2e-16 ***
#   Residuals   42  4117931    98046                   
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#p值小于0.05，说明培训方案是影响销售业绩的主要因素




#多重比较
TukeyHSD(data_aov)
par(las=1)
par(mar=c(3,5,3,4))
plot(TukeyHSD(data_aov))

#F分布
x <- seq(0,5,length.out = 1000)
y <- df(x,2,42)
ggplot(data=NULL,aes(x,y))+geom_line()+geom_vline(xintercept = qf(0.95,2,42))+
  annotate("text",x=3.8,y=0.1,label=qf(0.95,2,42))

image.png

6.双因素方差分析

image.png

#双因素方差分析
#数据源
data2 <- ToothGrowth

#步骤1 对数据进行初步了解
attach(data2)
#查看因素的水平及其频数分布
table(supp,dose)
# dose
# supp 0.5  1  2
# OJ  10 10 10
# VC  10 10 10




#看变量数和观测数
str(data2)
# 'data.frame': 60 obs. of  3 variables:
#   $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
# $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
# $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...


#查看各组的均值和标准差
aggregate(len~supp+dose,data2,mean)
aggregate(len~supp+dose,data2,sd)
# > aggregate(len~supp+dose,data2,mean)
# supp dose   len
# 1   OJ  0.5 13.23
# 2   VC  0.5  7.98
# 3   OJ  1.0 22.70
# 4   VC  1.0 16.77
# 5   OJ  2.0 26.06
# 6   VC  2.0 26.14
# > aggregate(len~supp+dose,data2,sd)
# supp dose      len
# 1   OJ  0.5 4.459709
# 2   VC  0.5 2.746634
# 3   OJ  1.0 3.910953
# 4   VC  1.0 2.515309
# 5   OJ  2.0 2.655058
# 6   VC  2.0 4.797731



#步骤2 方差齐性检验
#方差齐性检验（interaction可以将多个自变量折叠成一个单一的变量用于表示不同变量之间的因素组合）
bartlett.test(len~interaction(supp,dose),data2)
# > bartlett.test(len~interaction(supp,dose),data2)
# 
# Bartlett test of homogeneity of variances
# 
# data:  len by interaction(supp, dose)
# Bartlett's K-squared = 6.9273, df = 5, p-value = 0.2261

#步骤3 方差分析
dose <- factor(dose)
data2_aov <- aov(len~supp+dose+supp:dose)
summary(data2_aov)

# Df Sum Sq Mean Sq F value   Pr(>F)    
# supp         1  205.4   205.4  15.572 0.000231 ***
#   dose         2 2426.4  1213.2  92.000  < 2e-16 ***
#   supp:dose    2  108.3    54.2   4.107 0.021860 *  
#   Residuals   54  712.1    13.2                     
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

detach(data2)



#结论：豚鼠的牙齿是受喂食方法和坏血栓含量这两个因素影响的，而且非常显著

7.回顾学习目标

image.png

R语言之统计专题

1.常用的两种抽样方法

2.描述性统计

3.正态分布和抽样分布

4.参数估计和假设检验

5.t检验

6.方差分析

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

R语言

生信

R语言可视化