R检验数据是否符合幂律分布

作者: 生信小鹏 | 来源:发表于2022-11-01 23:32 被阅读0次

R检验数据是否符合幂律分布
R检验数据是否符合正态分布
大数据时代的“小数据系列3 --Shapiro-Wilk检验
R语言统计检验
r语言S-W和K-S两种正态性检验
R语言-T检验、秩和检验、百分比检验、卡方检验
刘嘉《概率论》15
内容平台的“幂律分布”
python数据分析之卡方检验
Hypothesis test

在进行WGCNA分析中，确定软阈值并给基因之间的相关性，加上一个幂次后，使得这种基因之间的连接度符合幂次分布。那么如何确定，或者通过统计学的方式检验，相关的分布符合幂律分布呢？

以【WGCNA学习笔记】无尺度分布和软阈值中的计算为例，详细分解一下如何检验是否符合幂律分布。

原始计算

power <- 8
ADJ=abs(cor(datExpr,use = 'p'))^power

其中表达矩阵是行为样本名，列为基因名
这样，将相关性矩阵绝对值之后，都给了一个8次幂，结果如下

> dim(ADJ)
[1] 3600 3600
> ADJ[1:5,1:5]
             MMT00000044  MMT00000046  MMT00000051  MMT00000076  MMT00000080
MMT00000044 1.000000e+00 4.322615e-12 5.724787e-09 1.392642e-05 1.195759e-07
MMT00000046 4.322615e-12 1.000000e+00 1.022965e-02 5.504572e-13 2.455594e-10
MMT00000051 5.724787e-09 1.022965e-02 1.000000e+00 1.967975e-10 7.471142e-11
MMT00000076 1.392642e-05 5.504572e-13 1.967975e-10 1.000000e+00 8.623133e-14
MMT00000080 1.195759e-07 2.455594e-10 7.471142e-11 8.623133e-14 1.000000e+00

然后对每一列的基因进行求和，得到的就是这个基因跟其他基因相关性之和。减去1是排除了自身。

k=apply(ADJ,2,sum) -1
hist(k)

结果如下

> k[1:10]
MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080 MMT00000102 MMT00000149 
  0.0814313  15.6397241  10.3721409   0.5134981  11.9155349   3.1397842  16.5059843 
MMT00000159 MMT00000207 MMT00000212 
 15.6716271  22.4635281   1.0466731

横坐标是基因的连接度，纵坐标是拥有该连接点数的频率，大体看上去，是头大尾巴长的幂律分布。

使用R验证数据符合幂律分布

那么怎么验证一个数据分布符合幂律分布

image.png
图片来自无标度网络幂律分布

从这幅图可以看出，如果一个数据分布符合幂律分布，那么其自变量的对数和因变量的对数是呈现线性关系的。

这个时候，就将问题转换为，是否符合线性关系的假设检验

线性关系的假设检验

先把之前的数据进行整理

cut1=cut(k,20)

#计算每个区间的平均值
binned.k=tapply(k,cut1,mean)
binned.k

# 计算每个区间频率
freq1=tapply(k,cut1,length)/length(k)
freq1

这个时候就是要看区间平均值和对数和每个区间的频率的对数是否为线性关系

也就是， $log(k)$ 和 $log(freq k)$ 是否为线性关系

假设检验

$H_0: \beta_1 = 0$ （没有线性关系）
$H_0: \beta_1 ≠ 0$ （有线性关系）

myfit_log <- lm(log(freq1) ~ log(binned.k))
anova(myfit_log)
summary(myfit_log)

结果如下

> anova(myfit_log)
Analysis of Variance Table

Response: log(freq1)
              Df Sum Sq Mean Sq F value    Pr(>F)    
log(binned.k)  1 48.979  48.979  102.03 7.647e-09 ***
Residuals     18  8.640   0.480                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这是方差分析表，结果表明模型拟合数据较好（F value = 102.03， P < 0.001)

> summary(myfit_log)

Call:
lm(formula = log(freq1) ~ log(binned.k))

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7026 -0.4063  0.2258  0.5573  0.7013 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     1.2268     0.5885   2.084   0.0516 .  
log(binned.k)  -1.7110     0.1694 -10.101 7.65e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6928 on 18 degrees of freedom
Multiple R-squared:   0.85, Adjusted R-squared:  0.8417 
F-statistic:   102 on 1 and 18 DF,  p-value: 7.647e-09

这是参数估计结果，自变量 $log(k)$ 有显著意义（ P < 0.001），常数项（ P < 0.1），因变量 $log(freq k)$ 变异中85%的可以被 $log(k)$ 解释( $R^2 = 0.85$ )