R語言進行基礎統計分析（一）

知白守黑。發表於2020-09-25

R語言進行基礎統計分析（一）

在 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

可發現，對於 conc 和 Treatment 都具有顯著性差異。
考慮交叉影響後：由上述結果可以看出，兩種因素水平下的樣本均值是有顯著性差異的，但它們之間相互作用檢驗p值為 0.52979 ，接受原假設，說明它們之間不存在互動作用。

社交網路分析的 R 基礎：（一）初探 R 語言
2022-02-05
R語言經典統計分析
2024-10-16
R語言
用R語言進行時間序列ARMA模型分析
2023-05-07
R語言模型
R語言-Survival analysis（生存分析）
2019-05-31
R語言
R語言資料質量分析
2024-03-21
R語言
使用R語言分析微信好友
2018-10-05
R語言
多元統計分析01：多元統計分析基礎
2021-10-19
Python文字統計與分析從基礎到進階
2024-05-06
Python
基於R語言的航空公司客戶價值分析
2020-05-14
R語言
r語言
2019-10-18
R語言
R語言入門與資料分析
2024-04-20
R語言
【R語言入門】R語言環境搭建
2021-09-09
R語言
dart語言基礎
2021-07-30
Dart
C語言基礎
2020-10-25
C語言
[Go]Go 語言基礎拾遺（一）
2022-01-18
Go
SQL語言基礎(資料控制語言)
2022-11-21
SQL
R 語言使用
2024-06-10
c語言的基礎
2024-03-17
C語言
【Go語言基礎】slice
2023-02-21
Go
e語言基礎01
2024-03-27
【01】C語言基礎
2020-09-30
C語言
Julia語言程式基礎
2021-04-16
python程式語言基礎
2019-01-24
Python
Go語言基礎-序言
2018-08-09
Go
c語言基礎的一些小技巧
2020-09-26
C語言
SQL語言基礎(SELECT語句)
2022-11-18
SQL
社交網路分析的 R 基礎：（四）迴圈與並行
2022-02-09
並行
OLAP引擎：基於Druid元件進行資料統計分析
2021-04-05
UI元件
用Python語言進行時間序列ARIMA模型分析
2023-05-08
Python模型
《R語言入門與資料分析》——向量索引
2020-10-02
R語言索引
R語言構建層次分析模型不看一下嗎~
2018-05-13
R語言模型
統計學基礎(一）
2018-10-22
零基礎轉行嵌入式——C語言
2024-05-28
C語言
學程式語言難嗎？零基礎想轉行進入IT行業，學哪個好呢？
2020-11-08
行業
Go語言基礎語法總結
2020-07-24
Go
C語言基礎函式
2024-07-24
C語言函式
C語言基礎-指標
2024-06-27
C語言指標
c語言基礎學習
2024-06-24
C語言

R語言進行基礎統計分析（一）

R語言進行基礎統計分析（一）

相關文章