环境与生态统计||统计假设

2018-10-29 本文已影响82人周运来就是我

统计思维本质上讲是归纳性的，它不像数学是演绎性的要借助于几条定理来推理，所以统计假设是统计分析和推断的基础。这个基础的意思是说，我们所做的统计是在满足这些条件之后的，而这些基础，往往在我们的教科书上草草提过之后，我们就忘了，然后以后的模型就是拿来就用。

其实这是不合理的，因为所有的模型都是有使用条件的，条件不满足，模型是不能用的。本章就介绍三种统计学上的基本假设：

正态假设
独立性假设
方差齐性假设

正态假设

所谓正态假设就是在统计分析之前假设变量是符合正态分布的？为什么是正态分布？

简单：两个变量（均值、标准差）即可表述
广泛存在：中心极限定理
生态环境多可近似为对数正态分布

什么是正态分布？

在随机收集来自独立来源的数据中，通常观察到数据的分布是正常的。这意味着，在绘制水平轴上的变量的值和垂直轴中的值的计数时，我们得到一个钟形曲线。曲线的中心代表数据集的平均值。在图中，百分之五十的值位于平均值的左侧，另外五十分之一位于图的右侧。统称为正态分布。

在R中，生成服从正态分布的随机变量序列可通过rnorm()函数实现，服从正态分布的随机变量在某点的 $pdf$ 和 $cdf$ 取值分别用dnorm()和pnorm()实现。

若随机变量X服从均值为μ，标准差为σ的概率分布，记为：
$X∼N(μ,σ^2)$
则其概率密度函数为：
$f(x)=\frac{1}{σ\sqrt{2\pi}}e^{\frac{-(x- \mu)^2}{2 σ^2}}$
正态分布的概率密度函数曲线呈钟形，对于单变量的正态分布而言，约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图的 $y$ 轴是用数据估计出的分位数， $x$ 轴是用标准正态分布的分位数。两者拟合的越好，数据越接近正态分布。

#生成服从标准正态分布的随机变量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之独立性检验

环境与生态统计||统计假设

正态假设

独立性假设

方差齐性假设

猜你喜欢

热点阅读