環境與生態統計||統計假設

統計思維本質上講是歸納性的，它不像數學是演繹性的要藉助於幾條定理來推理，所以統計假設是統計分析和推斷的基礎。這個基礎的意思是說，我們所做的統計是在滿足這些條件之後的，而這些基礎，往往在我們的教科書上草草提過之後，我們就忘了，然後以後的模型就是拿來就用。

其實這是不合理的，因為所有的模型都是有使用條件的，條件不滿足，模型是不能用的。本章就介紹三種統計學上的基本假設：

• 正態假設
• 獨立性假設
• 方差齊性假設

正態假設

所謂正態假設就是在統計分析之前假設變數是符合正態分佈的？為什麼是正態分佈？

• 簡單：兩個變數（均值、標準差）即可表述
• 廣泛存在：中心極限定理
• 生態環境多可近似為對數正態分佈

什麼是正態分佈？

在隨機收集來自獨立來源的資料中，通常觀察到資料的分佈是正常的。 這意味著，在繪製水平軸上的變數的值和垂直軸中的值的計數時，我們得到一個鐘形曲線。 曲線的中心代表資料集的平均值。 在圖中，百分之五十的值位於平均值的左側，另外五十分之一位於圖的右側。 統稱為正態分佈。

在R中，生成服從正態分佈的隨機變數序列可通過rnorm()函式實現，服從正態分佈的隨機變數在某點的和取值分別用dnorm()和pnorm()實現。

若隨機變數X服從均值為μ，標準差為σ的概率分佈，記為：

則其概率密度函式為：

正態分佈的概率密度函式曲線呈鐘形，對於單變數的正態分佈而言，約68.3%數值分佈在距離平均值有1個標準差之內的範圍，約95.4%數值分佈在距離平均值有2個標準差之內的範圍，以及約99.7%數值分佈在距離平均值有3個標準差之內的範圍。該理論也稱為“68-95-99.7法則”或“經驗法則”。

1.dnorm()函式
該函式給出給定平均值和標準偏差在每個點的概率分佈的高度(密度值)。

# Create a sequence of numbers between -100 and 100 incrementing by 0.1.
x <- seq(-100, 100, by = .1)

# Choose the mean as 2.5 and standard deviation as 0.5.
y <- dnorm(x, mean = 0, sd = 15)
plot(x,y)
abline(h=0,lty=3)
abline(h=0,lty=3)

2.pnorm()函式
該函式給出正態分佈隨機數小於給定數值的概率。它也被稱為“累積分佈函式”。

y <- pnorm(x, mean = 0, sd = 15)
plot(x,y)

3.qnorm()函式
該函式採用概率值，並給出其累積值與概率值匹配的數字值。

# Create a sequence of probability values incrementing by 0.02.
x <- seq(0, 1, by = 0.02)

# Choose the mean as 2 and standard deviation as 3.
y <- qnorm(x, mean = 2, sd = 1)
plot(x,y)

4.rnorm()函式
該函式用於生成分佈正常的隨機數，它將樣本大小作為輸入，並生成許多隨機數。我們繪製直方圖以顯示生成數字的分佈。

# Create a sample of 50 numbers which are normally distributed.
y <- rnorm(50)
hist(y, main = "正態分佈")

R中還提供了許多函式用於檢驗某隨機變數是否服從正態分佈，如Shapiro-Wilk normality檢驗、直方圖或者QQ圖，分別對應R中shapiro.test()、hist()和qqnorm()函式。QQ圖的軸是用資料估計出的分位數，軸是用標準正態分佈的分位數。兩者擬合的越好，資料越接近正態分佈。

#生成服從標準正態分佈的隨機變數x，樣本數為1000
x <- rnorm(1000, mean = 0, sd = 1)
#SW檢驗:樣本數量需控制在3~5000
shapiro.test(x) 

Shapiro-Wilk normality test

data:  x
W = 0.99844, p-value = 0.5178

par(mfrow=c(2,2))
#繪製QQ圖和直方圖
qqnorm(x)
hist(x)
#隨機生成長度為100的服從均勻分佈的隨機變數
y <- runif(100, min = 0, max = 1)
#SW檢驗
shapiro.test(y)


    Shapiro-Wilk normality test

data:  y
W = 0.95249, p-value = 0.001215

#繪製QQ圖和直方圖
qqnorm(y)
hist(y)

獨立性假設

什麼是資料的獨立性？我的理解就是在抽樣時，取得的每一個抽樣值均不受其它抽樣值的影響。從這個觀點看，在沒有特殊因素影響的條件下，有放回抽樣就是獨立的，而無放回抽樣就是不獨立的。

在生態學中由於生態系統往往有生態梯度，生態資料也容易引起不獨立的情況。為什麼要求資料是獨立的？這跟各種統計分析方法的前提(假設)條件有關。回想一下數理統計書裡，在引入定理時，第一句就是設X獨立同分布,這就是前提，沒有這個假設，後面的推導就是錯誤的。

如果發現資料不是獨立的，建議僅能使用執行圖展示資料採集結果的狀況，可以使用點估計值來估計引數，但不要計算置信區間上下限，不要進行假設檢驗。在這種情況下，最重要的任務不是討論如何進行後面的分析，而是首先搞清為什麼會出現不獨立的情況。

怎樣來驗證資料的獨立性呢？兩種方法，遊程檢驗和自相關係數。

# 獨立性檢驗
# R提供了多種檢驗類別型變數獨立性的方法。本節中描述的三種檢驗分別為卡方獨立性檢驗、
# Fisher精確檢驗和Cochran-Mantel–Haenszel檢驗。
 
# 1. 卡方獨立性檢驗
# 你可以使用chisq.test()函式對二維表的行變數和列變數進行卡方獨立性檢驗
 
library(vcd)
mytable<-xtabs(~Treatment+Improved,data = Arthritis)
mytable
# 
# > mytable
# Improved
# Treatment None Some Marked
# Placebo   29    7      7
# Treated   13    7     21
 
chisq.test(mytable)
 
# > chisq.test(mytable)
# 
# Pearson's Chi-squared test
# 
# data:  mytable
# X-squared = 13.055, df = 2, p-value = 0.001463
 
mytable<-xtabs(~Improved+Sex,data = Arthritis)
 
chisq.test(mytable)
 
# > chisq.test(mytable)
# 
# Pearson's Chi-squared test
# 
# data:  mytable
# X-squared = 4.8407, df = 2, p-value = 0.08889
 
# 這裡的p值表示從總體中抽取的樣本行變數與列變數是相互獨立的概率。
# 在結果1中，患者接受的治療和改善的水平看上去存在著某種關係（p < 0.01）。2 而患者性別
# 和改善情況之間卻不存在關係（p > 0.05）
# 
# 由於1的概率值很小，所以你拒絕了治療型別和治療結果相互獨立的原假設
# 
# 由於2的概率不夠小，故沒有足夠的理由說明治療結果和性別之間是不獨立的
 
# 2. Fisher精確檢驗
# 可以使用fisher.test()函式進行Fisher精確檢驗。Fisher精確檢驗的原假設是：邊界固定
# 的列聯表中行和列是相互獨立的。其呼叫格式為fisher.test(mytable)，其中的mytable是一個二維列聯表。
 
mytable<-xtabs(~Treatment+Improved,data = Arthritis)
fisher.test(mytable)
# > fisher.test(mytable)
# 
# Fisher's Exact Test for Count Data
# 
# data:  mytable
# p-value = 0.001393
# alternative hypothesis: two.sided
 
# 
# 3. Cochran-Mantel—Haenszel檢驗
# mantelhaen.test()函式可用來進行Cochran—Mantel—Haenszel卡方檢驗，其原假設是，兩
# 個名義變數在第三個變數的每一層中都是條件獨立的。下列程式碼可以檢驗治療情況和改善情況在
# 性別的每一水平下是否獨立。此檢驗假設不存在三階互動作用（治療情況×改善情況×性別）
 
mytable<-xtabs(~Treatment+Improved+Sex,data=Arthritis) mantelhaen.test(mytable)

# Cochran-Mantel-Haenszel test
# 
# data:  mytable
# Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647
# 結果表明，患者接受的治療與得到的改善在性別的每一水平下並不獨立

方差齊性假設

總體之間的方差相等時比較總體均值時的必要假設。如果標準差不同，比較均值的意義就沒那麼大了。方差齊性是做方差分析的前提（在最小二乘法的框架下）。

檢驗資料方差齊性的方法有很多：

• Bartlett檢驗 - 如果我們的資料服從正態分佈，那麼這種方法將是最為適用的。對於正態分佈的資料，這種檢驗極為靈敏；而當資料為非正態分佈時，使用該方法則很容易導致假陽性誤判。

> require(graphics)
> plot(count ~ spray, data = InsectSprays)
> bartlett.test(InsectSprays$count, InsectSprays$spray)

    Bartlett test of homogeneity of variances

data:  InsectSprays$count and InsectSprays$spray
Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05

> bartlett.test(count ~ spray, data = InsectSprays)

    Bartlett test of homogeneity of variances

data:  count by spray
Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05

• Levene檢驗 - 相較於Bartlett檢驗，這一方法更為穩健。這一方法被封裝於car程式包中。

library(car)
> with(Moore, leveneTest(conformity, fcategory))
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2   0.046 0.9551
      42               
> with(Moore, leveneTest(conformity, interaction(fcategory, partner.status)))
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  1.4694 0.2219
      39               
> leveneTest(conformity ~ fcategory*partner.status, data=Moore)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  1.4694 0.2219
      39               
> leveneTest(lm(conformity ~ fcategory*partner.status, data=Moore))
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  1.4694 0.2219
      39               
> leveneTest(conformity ~ fcategory*partner.status, data=Moore, center=mean)
Levene's Test for Homogeneity of Variance (center = mean)
      Df F value Pr(>F)
group  5  1.7915 0.1373
      39               
> leveneTest(conformity ~ fcategory*partner.status, data=Moore, center=mean, trim=0.1)
Levene's Test for Homogeneity of Variance (center = mean: 0.1)
      Df F value Pr(>F)
group  5  1.7962 0.1363
      39               

• Fligner-Killeen檢驗 - 這是一個非引數的檢驗方法，完全不依賴於對分佈的假設。

> require(graphics)
> 
> plot(count ~ spray, data = InsectSprays)
> fligner.test(InsectSprays$count, InsectSprays$spray)

    Fligner-Killeen test of homogeneity of variances

data:  InsectSprays$count and InsectSprays$spray
Fligner-Killeen:med chi-squared = 14.483, df = 5, p-value = 0.01282

> fligner.test(count ~ spray, data = InsectSprays)

    Fligner-Killeen test of homogeneity of variances

data:  count by spray
Fligner-Killeen:med chi-squared = 14.483, df = 5, p-value = 0.01282

對於上述所有的檢驗，我們的原假設都為“變數的總體方差全部相同”；備擇假設則為“至少有兩個變數的總體方差時不同的”。

require(graphics)
require(mvabund)
## Load the tikus dataset:
data(tikus)
tikusdat <- mvabund(tikus$abund)
year <- tikus$x[,1]

par(mfrow=c(2,2))
## Plot mean-variance plot for a mvabund object with a log scale (default):
meanvar.plot(tikusdat)  

## Again but without log-transformation of axes:
meanvar.plot(tikusdat,log="")   

## A mean-variance plot, data organised by year, 
## for 1981 and 1983 only, as in Figure 7a of Warton (2008):
is81or83 <- year==81 | year==83
meanvar.plot(tikusdat~year, subset=is81or83, col=c(1,10))

參考：

為什麼要假設變數為正態分佈？
統計基礎篇之四：為什麼樣本的獨立性如此重要？
獨立性檢驗的基本思想和初步應用
方差分析為何要進行方差齊次檢驗？如果方差不同會有什麼影響？
為什麼要進行方差齊性檢驗？
統計學檢驗——正態性檢驗和方差齊性檢驗等
R語言正態分佈
R語言——正態分佈
R之獨立性檢驗