||
这是我最近研究工作中遇到的一个有趣现象:有一组含对照与控制样本的多维数据在主成分分析里没有看到差异,但进行差异分析时发现几乎所有维度上都有差异。这就会形成解释数据上的悖论:如果我用主成分分析来进行探索式数据分析,会认为对照组与控制组很难区分;但如果我用单一维度的差异分析去对比,即使经过错误发现率控制,在绝大多数维度上都能看到差异。这里的问题在于,因为你每个维度上都显示了差异,按说整体降维后可视化差异应该很明显才对。那么,这两组数据究竟算有差异还是没差异?
这里我用模拟仿真来重现下这个问题。首先,我们模拟两组数据,对照组与控制组均有100维,都存在1.2倍均值差异,每组十万个点:
library(genefilter)# 100个维度np <- rnorm(100,100,100) z <- c()for(i in 1:100){ case <- rnorm(100000,mean = np[i],sd=abs(np[i])*0.1+1) control <- rnorm(100000,mean=np[i]*1.2,sd=abs(np[i])*0.1+1) zt <- c(case,control) z <- cbind(z,zt) }
我们来看下主成分分析结果:
pca <- prcomp(z)# 主成分贡献summary(pca)$importance[,c(1:10)]
## PC1 PC2 PC3 PC4 PC5 PC6 ## Standard deviation 127.9706 34.74543 30.08979 28.33013 25.78264 25.46286 ## Proportion of Variance 0.4817 0.03551 0.02663 0.02361 0.01955 0.01907 ## Cumulative Proportion 0.4817 0.51719 0.54382 0.56743 0.58698 0.60605 ## PC7 PC8 PC9 PC10 ## Standard deviation 24.53831 23.80237 23.11194 22.56009 ## Proportion of Variance 0.01771 0.01666 0.01571 0.01497 ## Cumulative Proportion 0.62376 0.64043 0.65614 0.67111
plot(pca$x[,1],pca$x[,2],pch=19,col=c(rep(1,100000),rep(2,100000)))
colnames(z) <- 1:ncol(z) re <- colttests(z,factor(c(rep(1,100000),rep(2,100000)))) re$adj <- p.adjust(re$p.value,'BH') sum(re$adj<0.05)
## [1] 100
很明显差异。下面我们要做些变化:
z <- c()for(i in 1:100){ case <- rnorm(100000,mean = np[i],sd=abs(np[i])*0.5+1) control <- rnorm(100000,mean=np[i]*1.2,sd=abs(np[i])*0.5+1) zt <- c(case,control) z <- cbind(z,zt) } pca <- prcomp(z)# 主成分贡献summary(pca)$importance[,c(1:10)]
## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 184.0775 158.44812 144.37309 134.64138 125.6758 ## Proportion of Variance 0.0798 0.05913 0.04909 0.04269 0.0372 ## Cumulative Proportion 0.0798 0.13893 0.18802 0.23071 0.2679 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 122.32096 117.99105 114.59217 111.63828 108.53031 ## Proportion of Variance 0.03524 0.03279 0.03093 0.02935 0.02774 ## Cumulative Proportion 0.30315 0.33594 0.36687 0.39622 0.42396
plot(pca$x[,1],pca$x[,2],pch=19,col=c(rep(1,100000),rep(2,100000)))
colnames(z) <- 1:ncol(z) re <- colttests(z,factor(c(rep(1,100000),rep(2,100000)))) re$adj <- p.adjust(re$p.value,'BH') sum(re$adj<0.05)
## [1] 100
目前差异已经开始互相融合了,好,进一步加大力度。
z <- c()for(i in 1:100){ case <- rnorm(100000,mean = np[i],sd=abs(np[i])*1.2+1) control <- rnorm(100000,mean=np[i]*1.2,sd=abs(np[i])*1.2+1) zt <- c(case,control) z <- cbind(z,zt) } pca <- prcomp(z)# 主成分贡献summary(pca)$importance[,c(1:10)]
## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 418.51939 355.72299 338.84637 303.01048 298.57129 ## Proportion of Variance 0.07504 0.05421 0.04919 0.03933 0.03819 ## Cumulative Proportion 0.07504 0.12925 0.17843 0.21777 0.25596 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 286.12208 278.05519 270.55114 261.78349 258.70506 ## Proportion of Variance 0.03507 0.03312 0.03136 0.02936 0.02867 ## Cumulative Proportion 0.29103 0.32415 0.35551 0.38487 0.41354
plot(pca$x[,1],pca$x[,2],pch=19,col=c(rep(1,100000),rep(2,100000)))
colnames(z) <- 1:ncol(z) re <- colttests(z,factor(c(rep(1,100000),rep(2,100000)))) re$adj <- p.adjust(re$p.value,'BH') sum(re$adj<0.05)
## [1] 100
目前这两组差异数据在前两个主成分上已经看不出来了。
其实这也算不上悖论,主要线索就是主成分的方差贡献,而我做的改动也在相对标准偏差上,第一个是10%,第二个是50%,第三个是120%。也就是说,当单一维度上方差已经大到均值水平时,前几个主成分展示的就不会是本来存在的均值差异。也就是说组内方差已经大于组间方差了,而主成分分析的前几个主成分展示的主要方差的方向,此时就很难看到差异了。
另外一个问题就在为什么每个维度上做t检验都能看到差异?是不是这个检验太灵敏了?这里我们可以换下非参检验试一下:
re <- apply(z, 2, function(x) wilcox.test(x~factor(c(rep(1,100000),rep(2,100000))))) pv <- sapply(re, function(x) x$p.value) table(pv)
## pv ## 0 3.91700237397975e-307 1.51138626356535e-306 ## 4 1 1 ## 1.69033357540812e-305 3.69816700489022e-305 7.41041359986221e-304 ## 1 1 1 ## 1.19907386314392e-302 1.6979081670706e-302 2.11441218478598e-300 ## 1 1 1 ## 1.02923123985138e-299 2.84715946155774e-299 1.13983910431486e-298 ## 1 1 1 ## 2.84511424476962e-298 1.36481875863161e-297 2.4138926412321e-297 ## 1 1 1 ## 5.01777165508798e-297 6.3740319190573e-296 1.06512806333407e-295 ## 1 1 1 ## 2.9684345038554e-295 1.61778936275495e-294 2.60338020017055e-294 ## 1 1 1 ## 1.91558533686995e-293 2.10282840972079e-292 8.9660345549283e-291 ## 1 1 1 ## 1.23746174866693e-290 9.7772064637835e-290 5.54917392263547e-289 ## 1 1 1 ## 1.49892579574228e-288 9.35563548715912e-288 1.08649454476636e-287 ## 1 1 1 ## 3.16786218854795e-287 4.69279905849196e-287 1.118282060232e-286 ## 1 1 1 ## 1.1962582225918e-286 7.05752716298586e-286 2.42440953972578e-285 ## 1 1 1 ## 2.13314214586901e-284 2.47347628881802e-284 5.56306148900296e-284 ## 1 1 1 ## 6.48930682676215e-284 1.03069189166261e-283 2.06944343269889e-283 ## 1 1 1 ## 2.21015902959185e-283 9.04094542944513e-283 6.8633065287707e-282 ## 1 1 1 ## 1.42377796179449e-281 3.12369877273024e-281 4.66807820660733e-281 ## 1 1 1 ## 1.09540654209284e-280 5.98881828424942e-280 1.77128558057562e-279 ## 1 1 1 ## 1.19605377826797e-278 1.38296886675744e-278 2.38617223132791e-278 ## 1 1 1 ## 7.84788712325474e-278 9.53098284619264e-278 1.1170495994422e-277 ## 1 1 1 ## 5.91455469655551e-277 6.52415199674643e-276 1.07167093679512e-274 ## 1 1 1 ## 1.31256681444049e-274 1.50209913056876e-274 5.65627008544608e-274 ## 1 1 1 ## 1.07915989034337e-273 1.41257365167248e-273 6.74115421225619e-273 ## 1 1 1 ## 7.37969169980536e-272 7.50240786629847e-272 7.98499354414677e-272 ## 1 1 1 ## 4.41356905563368e-271 6.75885482781668e-271 7.56920082447655e-271 ## 1 1 1 ## 1.63370604905236e-270 6.11794474873258e-270 1.33732619109681e-269 ## 1 1 1 ## 1.04857332193612e-268 2.65251222934155e-268 1.48710002018571e-267 ## 1 1 1 ## 9.78882899780169e-267 1.74052952159569e-266 4.73167179342876e-266 ## 1 1 1 ## 5.16985694132384e-266 7.13131540530313e-264 1.98256549564857e-263 ## 1 1 1 ## 7.05483245806848e-263 1.09357212567035e-262 1.33842846539623e-259 ## 1 1 1 ## 2.57013184589882e-254 5.98676375371462e-250 3.13570702579997e-239 ## 1 1 1 ## 5.72466365564095e-236 2.46892905766707e-226 5.75155691148484e-223 ## 1 1 1 ## 4.95859083943006e-214 2.04379312687636e-160 4.87202652382768e-157 ## 1 1 1 ## 1.08969775204526e-129 ## 1
结果很直观,依然全是差异。那么问题在哪里?其实就是差异本身的性质,这个差异太小,同时样品量又太大,结果这个客观存在的差异就会在大样本的情况下被检验出来。那么皮球就踢回来了,这个微弱但存在的差异有没有物理意义或科学意义?
如果这个差异出现在无关紧要的地方,那么即使有差异可能也没多少意义。打比方我们对比了两个国家的人口拇指指甲宽度,发现其中一个国家的人口拇指指甲平均长度比另一个国家宽0.5毫米。从统计意义上存在显著差异,但实际意义上几乎为零。不过如果这个客观存在的差异有实际意义,那么想用探索式数据分析找到或者发现就不容易了。
什么时候会出现这个问题呢?就是样本量特别大的时候,或者说现在。我们现在把样品量降下来看看:
# 不同样品 n <- c(50,100,1000,5000,10000,50000) par(mfrow=c(2,3))for(t in n){ z <- c() for(i in 1:100){ case <- rnorm(t,mean = np[i],sd=abs(np[i])*1.2+1) control <- rnorm(t,mean=np[i]*1.2,sd=abs(np[i])*1.2+1) zt <- c(case,control) z <- cbind(z,zt) } pca <- prcomp(z) # 主成分贡献 print(summary(pca)$importance[,c(1:10)]) plot(pca$x[,1],pca$x[,2],pch=19,col=c(rep(1,t),rep(2,t)),main=paste(t,'samples')) colnames(z) <- 1:ncol(z) re <- colttests(z,factor(c(rep(1,t),rep(2,t)))) re$adj <- p.adjust(re$p.value,'BH') print(sum(re$adj<0.05)) }
## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 456.07143 424.5170 391.09387 355.43286 336.11953 ## Proportion of Variance 0.08715 0.0755 0.06408 0.05293 0.04733 ## Cumulative Proportion 0.08715 0.1626 0.22674 0.27966 0.32700 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 311.41814 307.23097 291.02704 282.45216 279.76759 ## Proportion of Variance 0.04063 0.03955 0.03549 0.03343 0.03279 ## Cumulative Proportion 0.36763 0.40718 0.44266 0.47609 0.50888
## [1] 0 ## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 442.27870 393.34622 370.6548 348.04532 320.82168 ## Proportion of Variance 0.08216 0.06499 0.0577 0.05088 0.04323 ## Cumulative Proportion 0.08216 0.14715 0.2049 0.25573 0.29896 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 309.48976 305.33719 295.9924 277.85593 267.75739 ## Proportion of Variance 0.04023 0.03916 0.0368 0.03243 0.03011 ## Cumulative Proportion 0.33919 0.37835 0.4152 0.44758 0.47769
## [1] 4 ## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 422.01184 366.44795 341.67684 304.49015 298.23717 ## Proportion of Variance 0.07616 0.05743 0.04992 0.03965 0.03804 ## Cumulative Proportion 0.07616 0.13359 0.18351 0.22316 0.26120 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 284.93725 283.75241 273.16411 266.97608 257.2485 ## Proportion of Variance 0.03472 0.03443 0.03191 0.03048 0.0283 ## Cumulative Proportion 0.29592 0.33035 0.36226 0.39274 0.4210
## [1] 96 ## PC1 PC2 PC3 PC4 PC5 PC6 ## Standard deviation 418.54777 357.8062 338.015 300.82293 298.73619 289.19724 ## Proportion of Variance 0.07512 0.0549 0.049 0.03881 0.03827 0.03587 ## Cumulative Proportion 0.07512 0.1300 0.179 0.21783 0.25610 0.29197 ## PC7 PC8 PC9 PC10 ## Standard deviation 279.78262 272.06504 262.07130 256.52711 ## Proportion of Variance 0.03357 0.03174 0.02945 0.02822 ## Cumulative Proportion 0.32553 0.35728 0.38673 0.41495
## [1] 100 ## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 415.28187 355.71525 339.53774 300.95094 298.23751 ## Proportion of Variance 0.07398 0.05428 0.04945 0.03885 0.03815 ## Cumulative Proportion 0.07398 0.12825 0.17771 0.21656 0.25471 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 286.94753 275.76851 271.76577 263.1141 259.25091 ## Proportion of Variance 0.03532 0.03262 0.03168 0.0297 0.02883 ## Cumulative Proportion 0.29003 0.32265 0.35433 0.3840 0.41286
## [1] 100 ## PC1 PC2 PC3 PC4 PC5 ## Standard deviation 416.94971 354.95062 337.60963 303.41448 299.91093 ## Proportion of Variance 0.07452 0.05401 0.04886 0.03946 0.03856 ## Cumulative Proportion 0.07452 0.12853 0.17739 0.21685 0.25541 ## PC6 PC7 PC8 PC9 PC10 ## Standard deviation 286.65155 277.67511 270.11921 261.86270 258.63352 ## Proportion of Variance 0.03522 0.03305 0.03128 0.02939 0.02867 ## Cumulative Proportion 0.29063 0.32369 0.35496 0.38436 0.41303
## [1] 100
这里我们可以看到,在样品单一维度有差异但方差比较大的情况下,当样品数超过维度后,差异分析跟主成分分析就已经出现灵敏度差异了。也就是说,当样品数增加后,类似主成分分析这种探索式数据分析已经看不出预设差异了。同时我们又注意到,在样品量少于维度时,差异分析功效又是不足的,也就是根本测不到预设差异。
这里有人可能觉得是不是主成分分析只考虑了维度间线性组合,如果我换一种非线性方法会不会还能看到差异?没问题,我们继续模拟:
library(umap) n <- c(50,100,1000,5000) par(mfrow=c(2,2))for(t in n){ z <- c() for(i in 1:100){ case <- rnorm(t,mean = np[i],sd=abs(np[i])*1.2+1) control <- rnorm(t,mean=np[i]*1.2,sd=abs(np[i])*1.2+1) zt <- c(case,control) z <- cbind(z,zt) } umap <- umap(z) plot(umap$layout,col=c(rep(1,t),rep(2,t)),pch=19,main=paste(t,'samples')) }
上面是用流形学习里的umap来进行的可视化,可以看出情况一致,毕竟我们模拟时差异并非是非线性的,所以流形学习也没法捕捉这种差异。
此处我想说的是,如果我们过分依赖差异分析,在大数据量背景下,大概率会检测到一些客观存在但微弱的差异,这些差异只在均值水平才能发现,具体到个体差异里很难看出来,属于统计意义可能大于实际意义的范畴。
这种情况不大可能发生在样品数小于维度数的场景下,那个场景里最大的问题是差异分析功效不足,增大样品量总是有益的。但却可以发生在维度远低于样品数,此时大数定律反而成了负面作用,因为我们总能发现差异或者说规律,但却缺少用专业知识来筛选差异实际意义的手段。
很显然,此时主成分分析或聚类分析这些方法都没啥用,这些自上而下的方法根本就看不出数据中的异质性。不过,这倒可能帮我们过滤掉哪些微弱差异,只能暂时认为这些差异就是无关紧要的。但此时会遇到的问题就是我开头说的悖论,明明差异分析发现了大量差异,但整体就好像糊成一片,如果数据分析经验不足可能就不知道咋解释了。
这种微弱的差异本质就是 type M 型错误,也就是如果效应比较弱,我们测到了也没实际意义。不仅仅是数据分析,在我们日常生活中这种错误也很常见,例如财经新闻在每天收市后的报道经常就是“受某某影响,今天大盘下降零点几个百分点”。零点几个百分点其实是属于随机过程,并不受某某影响,但即使是财经专家也要对着一堆噪音去找解释。更搞笑的是,这些专家的追随者在听了专家点拨后,马上也看到了潜在的规律性,然后带着优越感去跟其他人说你们能力不行,想不到这么远之类的。但问题是本来也没规律性啊,或者说这个所谓的规律本身的不确定性是大于确定性的。也就是当样本间方差足够大时,样本间均值差异的实际意义就很微弱了。好比两国人均GDP有显著差异,但两个国家内部贫富差异悬殊,此时与其比对人均GDP,不如把两国穷人划分成一类,富人划分成一类来进行对比,此时发现的问题可能更有实际意义。
我一般不讨论社会科学规律,很大原因在于这些学科没有完成科学化改造,其规律性更像是一种自我实现的过程。也就是我首先构建一套逻辑自洽的理论,然后用理论指导实践,怎么做怎么对,用实际行动证明理论可行性。但如果同一个学科存在两种逻辑都自洽且都经过事实检验的理论,那你很难说哪种就是客观规律,或者说不存在最优路线,这就跟科学规律不一样了。自然科学规律是相对唯一的,做不到既同意理论A又同意理论B,最后一定会被更一般性的理论C统一起来。社会科学里的很多所谓规律其实是很主观的学术观点或认知体系,好一点的可以通过观察实验来证明,极端的就单纯讲逻辑自洽,在自己理论圈子里满地打滚抱团取暖,而他们眼中的金科玉律也许就像是上面我模拟的那样,客观存在但实际没意义,或者压根分类就不合理。
要警惕那些“大”数据带来的发现。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-8 09:35
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社