最近有朋友问到这个问题,说是计算某个性状的遗传力,发现这个性状不符合正态分布,问怎么转化才可以分析???
我的经验是一般不做转化,因为转化后结果可能与你想要的不一样,包括遗传力,育种值排名,稍微的偏态也是可以分析的。数据量较少的严重的偏态,转化后也是无能为力的。
在模型假定中,数据需要符合正态分布,在动物模型分析中,如果对数据进行转化,对育种值和方差组分变化的情况,进行测试。
library(asreml)
library(learnasreml)
data("animalmodel.dat")
data("animalmodel.ped")
dat = animalmodel.dat
ped = animalmodel.ped
dat[dat==0] = NA
# 原始数据分布
hist(dat$BWT)
方差组分:
# 原始数据,计算动物模型方差组分
head(dat)
ainv = ainverse(ped)
m1 = asreml(BWT ~ BYEAR + SEX, random = ~ vm(ANIMAL,ainv),data=dat)
summary(m1)$varcomp
育种值blup:
# 计算原始数据育种值blup
blup = tiqu_blup(coef(m1)$random)
head(blup)
结论:
# 原始数据减去一个常数:10
# 方差组分不变,育种值不变
dat$a = dat$BWT - 10
hist(dat$a)
m2 = asreml(a ~ BYEAR + SEX, random = ~ vm(ANIMAL,ainv),data=dat)
summary(m2)$varcomp
blup1 = tiqu_blup(coef(m2)$random)
head(blup1)
cor(blup$effect,blup1$effect)
表型数据分布:
方差组分:
结论:
# 原始数据,log10(x)转化
# 方差组分变化,育种值与原来育种值相关性系数:0.9
dat$a = log10(dat$BWT)
hist(dat$a) # 偏离正态分布
m2 = asreml(a ~ BYEAR + SEX, random = ~ vm(ANIMAL,ainv),data=dat)
summary(m2)$varcomp
blup1 = tiqu_blup(coef(m2)$random)
cor(blup$effect,blup1$effect)
plot(blup$effect,blup1$effect)
表型数据分布:
方差组分:
和原育种值相关性分析:
结论:
# 原始数据,log10(x +10)转化
# 方差组分变化,育种值与原来育种值相关性系数:0.99
dat$a = log10(dat$BWT+10)
hist(dat$a) # 符合正态分布
m2 = asreml(a ~ BYEAR + SEX, random = ~ vm(ANIMAL,ainv),data=dat)
summary(m2)$varcomp
blup1 = tiqu_blup(coef(m2)$random)
cor(blup$effect,blup1$effect)
plot(blup$effect,blup1$effect)
表型数据分布:
方差组分:
育种值相关性分析:
方差组分变化:
遗传力变化:
所以,不要轻易做数据转化,多多考虑加大数据量,或者提高数据质量吧。统计方法不能拯救失败的数据,只能告诉你数据分析为什么会失败,这是我说的。
下面是Fisher说的:
科学试验做完后再找统计学家分析数据,如同病人死了找医生进行尸体解剖,医生会告诉你病人死的原因。同样,统计学家会告诉你试验失败的原因。 ------Ronald Fisher
All models are wrong, but some models are useful | 模型皆有误,或尤建奇功。 George E. P. Box 在终极的分析中,一切知识都是历史;在抽象的意义下,一切科学都是数学;在理性的基础上,所有的判断都是统计学。 C.R.Rao 统计学是对令人困惑费解的问题做出数字设想的艺术。 David Freedma 统计具有非凡的能力处理各种复杂的问题,它需要非常精细的方法和小心翼翼的解释。当人类科学探索者在问题的丛林中遇到难以逾越的障碍时,唯有统计工具可为其开辟一条前进的通道。 Francis Galton 对统计学的一知半解常常会造成一些不必要的上当受骗,对统计学的一概排斥往往会造成某些不必要的愚昧无知。 C.R.Rao
所以那些星座迷,以星座不配甩了你的妹子千万不能要。
欢迎关注我的公众号:育种数据分析之放飞自我