二氧化碳吸收速率统计分析

俱往矣`

已于 2022-04-12 20:41:45 修改

阅读量6.3k

点赞数 5

分类专栏： R语言数据分析与挖掘文章标签：可视化 r语言数据分析统计学

于 2020-09-25 23:52:21 首次发布

本文链接：https://blog.csdn.net/weixin_43180762/article/details/108805893

版权

R语言数据分析与挖掘专栏收录该内容

33 篇文章 44 订阅

订阅专栏

二氧化碳吸收速率统计分析

在 R 自带的 datasets 包中包含一个 CO2（注意是大写字母）数据集（可通过代码“head(CO2)”查看数据的前 5 行，通过“?CO2”查看每个变量的具体意义），该数据是在某项实验中植物对二氧化碳的吸收情况的记录，共有 84 行 5 列。针对该数据集，请进行以下问题的分析：
（1）二氧化碳吸收速率（uptake）的数据是正态分布吗？如果不是，那是什么分布？
先绘制基础的分布图像：

head(CO2)
co2=CO2
#uptake指的是吸收率
#conc环境二氧化碳浓度的数值矢量（毫升/升）
plot(co2$uptake,type = "b")

在这里插入图片描述

hist(co2$uptake, prob = TRUE, col = "lightblue",main="Normal Distribution")

在这里插入图片描述
再进行精确的 Kolmogorov-Smirnov 检验：

library(nortest) 
lillie.test(co2$uptake)#精确的Kolmogorov-Smirnov检验
shapiro.test(co2$uptake)#夏皮罗-威尔克检验

结果如下：

> lillie.test(co2$uptake)#精确的Kolmogorov-Smirnov检验

	Lilliefors (Kolmogorov-Smirnov) normality test

data:  co2$uptake
D = 0.1077, p-value = 0.01745


> shapiro.test(co2$uptake)#夏皮罗-威尔克检验

	Shapiro-Wilk normality test

data:  co2$uptake
W = 0.94105, p-value = 0.0007908

p-value<0.05，说明数据不服从标准正态分布。

可绘制QQ图检验分布情况：

qqnorm(co2$uptake);qqline(co2$uptake)
par(mai=c(0.9, 0.9, 0.6, 0.2)) #图形边界参数
hist(co2$uptake,breaks = c(0,7.5, 15, 22.5, 30,37.5,45,52.5,60),freq=FALSE,density=10,angle=60)

在这里插入图片描述

经过检验可得知，并非是正态分布，而是混合正态分布，是一种双峰分布。

（2）针对上面检验出的二氧化碳吸收速率的数据分布形式，对其参数进行估计，计算出相应分布下的参数；

library(MASS)
attach(geyser)
#定义log-likelihood函数
LL<-function(params,data) 
{#参数"params"是一个向量，依次包含了五个参数：p,mu1,sigma1,mu2,sigma2.
  #参数"data"，是观测数据。
  t1<-dnorm(data,params[2],params[3])
  t2<-dnorm(data,params[4],params[5])
  #这里的dnorm()函数是用来生成正态密度函数的。
  f<-params[1]*t1+(1-params[1])*t2
  #混合密度函数
  ll<-sum(log(f))
  #log-likelihood函数
  return(-ll)
  #nlminb()函数是最小化一个函数的值，但我们是要最大化log-likeilhood函数，所以需要在“ll”前加个“-”号。
}
#用hist函数找出初始值
hist(co2$uptake,freq=F)
lines(density(co2$uptake))
#拟合函数####optim####
geyser.res<-nlminb(c(0.5,15,15,35,15),LL,data=co2$uptake,
                   
                   lower=c(0.0001,-Inf,0.0001,-Inf,-Inf,0.0001),
                   
                   upper=c(0.9999,Inf,Inf,Inf,Inf))
#初始值为p=0.5,mu1=15,sigma1=15,mu2=15,sigma2=15

#LL是被最小化的函数。

#data是拟合用的数据

#lower和upper分别指定参数的上界和下界。
#查看拟合的参数
geyser.res$par
#拟合的效果
X<-seq(0,50,length=100)
#读出估计的参数
p<-geyser.res$par[1]
mu1<-geyser.res$par[2]
sig1<-geyser.res$par[3]
mu2<-geyser.res$par[4]
sig2<-geyser.res$par[5]
#将估计的参数函数代入原密度函数。
f<-p*dnorm(X,mu1,sig1)+(1-p)*dnorm(X,mu2,sig2)
#作出数据的直方图
hist(co2$uptake,probability=T,col=0,ylab="Density",  
     ylim=c(0,0.04),xlab="Eruption waiting times")
#画出拟合的曲线
lines(X,f)

结果如下：

> geyser.res$par
[1]  0.4014459 15.7219998  4.1541111 34.9200993  5.7905600
> p
[1] 0.4014459
> mu1
[1] 15.722
> sig1
[1] 4.154111
> mu2
[1] 34.9201
> sig2
[1] 5.79056

参数估计前的分布图形：
在这里插入图片描述

参数估计后的图像： 在这里插入图片描述
（3）针对 4 个因子变量和二氧化碳吸收速率，分别进行 t 检验和单因素方差分析，对比二氧化碳的吸收速率在不同因素下的差异性；


co2$value<-gl(4,21)
head(co2)
#Treatment之间是否有显著性差异
t.test(uptake~co2$Treatment, data = co2, alternative = "less", var.equal = TRUE)
#Type之间是否有显著性差异
t.test(uptake~co2$Type, data = co2, alternative = "less", var.equal = TRUE)
##############观察两两之间是否有显著性差异###########
c1<-co2[which(co2$value==1|co2$value==2),]#1和2类
t.test(uptake~value, data = c1, alternative = "less", var.equal = TRUE)
c2<-co2[which(co2$value==1|co2$value==3),]#1和3类
t.test(uptake~value, data = c2, alternative = "less", var.equal = TRUE)
c3<-co2[which(co2$value==1|co2$value==4),]#1和4类
t.test(uptake~value, data = c3, alternative = "less", var.equal = TRUE)
c4<-co2[which(co2$value==2|co2$value==3),]#2和3类
t.test(uptake~value, data = c4, alternative = "less", var.equal = TRUE)
c5<-co2[which(co2$value==2|co2$value==4),]#2和4类
t.test(uptake~value, data = c5, alternative = "less", var.equal = TRUE)
c6<-co2[which(co2$value==3|co2$value==4),]#3和4类
t.test(uptake~value, data = c6, alternative = "less", var.equal = TRUE)

结果如下：

> #Treatment之间是否有显著性差异
> t.test(uptake~co2$Treatment, data = co2, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by co2$Treatment
t = 3.0485, df = 82, p-value = 0.9985
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
   -Inf 10.603
sample estimates:
mean in group nonchilled    mean in group chilled 
                30.64286                 23.78333 

> #Type之间是否有显著性差异
> t.test(uptake~co2$Type, data = co2, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by co2$Type
t = 6.5969, df = 82, p-value = 1
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 15.85208
sample estimates:
     mean in group Quebec mean in group Mississippi 
                 33.54286                  20.88333 

> 
> ##############观察两两之间是否有显著性差异###########
> c1<-co2[which(co2$value==1|co2$value==2),]#1和2类
> t.test(uptake~value, data = c1, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by value
t = 1.2061, df = 40, p-value = 0.8826
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 8.580289
sample estimates:
mean in group 1 mean in group 2 
       35.33333        31.75238 

> c2<-co2[which(co2$value==1|co2$value==3),]#1和3类
> t.test(uptake~value, data = c2, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by value
t = 3.5471, df = 40, p-value = 0.9995
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 13.83421
sample estimates:
mean in group 1 mean in group 3 
       35.33333        25.95238 

> c3<-co2[which(co2$value==1|co2$value==4),]#1和4类
> t.test(uptake~value, data = c3, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by value
t = 8.5846, df = 40, p-value = 1
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 23.34765
sample estimates:
mean in group 1 mean in group 4 
       35.33333        15.81429 

> c4<-co2[which(co2$value==2|co2$value==3),]#2和3类
> t.test(uptake~value, data = c4, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by value
t = 2.1861, df = 40, p-value = 0.9826
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 10.26737
sample estimates:
mean in group 2 mean in group 3 
       31.75238        25.95238 

> c5<-co2[which(co2$value==2|co2$value==4),]#2和4类
> t.test(uptake~value, data = c5, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by value
t = 6.9798, df = 40, p-value = 1
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
    -Inf 19.7831
sample estimates:
mean in group 2 mean in group 4 
       31.75238        15.81429 

> c6<-co2[which(co2$value==3|co2$value==4),]#3和4类
> t.test(uptake~value, data = c6, alternative = "less", var.equal = TRUE)

	Two Sample t-test

data:  uptake by value
t = 5.5033, df = 40, p-value = 1
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 13.24007
sample estimates:
mean in group 3 mean in group 4 
       25.95238        15.81429

经过检验，两两之间都有显著性差异。

试着对其可视化：

library(gplots)
par(family="STKaiti")
plotmeans(uptake~value,co2,col="red",main="")

在这里插入图片描述
由上图可发现，各种组合的四种水平下，差异较为明显，尤其是 1-1 与与 4-4相对比，相差过半。
单因素方差分析：

lamp.aov<-aov( uptake~value , data=co2) #单因素方差分析 
summary(lamp.aov) #提取方差分析表 • 
##使用lm()和anova()函数组合, 可得到相同的计算结果 
lamp.lm<-lm(uptake~value, data=co2) 
anova(lamp.aov)
##################多重比较并可视化#################
tky <- TukeyHSD(lamp.aov)
tky = as.data.frame(tky$value)
tky$pair<-rownames(tky)
library(ggplot2)
# Plot pairwise TukeyHSD comparisons and color by significance level
ggplot(tky, aes(colour=cut(`p adj`, c(0, 0.01, 0.05, 1), 
                           label=c("p<0.01","p<0.05","Non-Sig")))) +
  theme_bw(base_size = 16)+
  geom_hline(yintercept=0, lty="11", colour="grey30",size = 1) +
  geom_errorbar(aes(pair, ymin=lwr, ymax=upr), width=0.2,size = 1) +
  geom_point(aes(pair, diff),size = 2) +
  labs(colour="")+
  theme(axis.text.x = element_text(size = 14))

结果如下：

> summary(lamp.aov) #提取方差分析表 • 
            Df Sum Sq Mean Sq F value   Pr(>F)    
value        3   4579  1526.5   23.82 4.11e-11 ***
Residuals   80   5128    64.1                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


> anova(lamp.aov)
Analysis of Variance Table

Response: uptake
          Df Sum Sq Mean Sq F value    Pr(>F)    
value      3 4579.4 1526.46  23.816 4.106e-11 ***
Residuals 80 5127.6   64.09                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

方差分析可看出，各因子水平下，有显著性差异。
接下来可分析哪些种类之间平均值相差较大，可使用 TukeyHSD（）函数得到方差分析的结果。
在这里插入图片描述
可发现：2-1，3-2不太显著外，其余都有显著性差异。

（4）将 Treatment、conc、uptake 三个变量进行双因素方差分析，并对输出的结果进行解释说明。

#无交互作用
aov1<-aov(uptake~conc+Treatment,data=co2)
summary(aov1)
#验证是否有交互作用
aov2<-aov(uptake~conc*Treatment,data=co2)
summary(aov2)

结果如下：

> #无交互作用
> aov1<-aov(uptake~conc+Treatment,data=co2)
> summary(aov1)
            Df Sum Sq Mean Sq F value   Pr(>F)    
conc         1   2285  2285.0   28.77 7.55e-07 ***
Treatment    1    988   988.1   12.44 0.000695 ***
Residuals   81   6434    79.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> #验证是否有交互作用
> aov2<-aov(uptake~conc*Treatment,data=co2)
> summary(aov2)
               Df Sum Sq Mean Sq F value   Pr(>F)    
conc            1   2285  2285.0  28.554 8.38e-07 ***
Treatment       1    988   988.1  12.348  0.00073 ***
conc:Treatment  1     32    31.9   0.398  0.52979    
Residuals      80   6402    80.0                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1