二値形質の狭義遺伝率 - ryamadaの遺伝学・遺伝統計学メモ

二値形質の場合には、上述の量的形質の場合に算出した、遺伝因子と環境因子との和であるphenoをそのまま表現型にするわけにはいかないので工夫が必要である
一つのモデルは、閾値を設定しそれより大きいと発病(値 1)、そうでなければ非発病(値0)とするものである
やってみる
上の量的形質と同じリスク計算をし、表現型決定だけ変えたもの

######## case/control
# No. loci
n.marker <- 100
# population size
N <- 10^5
# Genotype matrix
g <- matrix(0,N,n.marker)
# allele freq of markers
f <- 0.1 + runif(n.marker) * 0.4
# HWE and random genotype generation
for(i in 1:n.marker){
	g[,i] <- sample(0:2,N,replace=TRUE,prob=c(f[i]^2,2*f[i]*(1-f[i]),(1-f[i])^2))
}
# Risk of each marker's allele
r <- rnorm(n.marker)
# Total genetic risk when addtive model
R <- rep(0,N)
for(i in 1:n.marker){
	R[which(g[,i]==1)] <- R[which(g[,i]==1)] + r[i]
	R[which(g[,i]==2)] <- R[which(g[,i]==2)] + r[i]*2

}
# narrow-sense heritability
h <- 0.7
# environmental variance 
eV <- (1-h)/h * var(R)
# sum of genetic and environmental risks
R. <- R + rnorm(N,0,sqrt(eV))
# phenotype is based on threshold model
#pheno <- R.
rnk <- rank(R.)
pheno <- rep(0,N)
# 有病率
prevalence <- 0.01
pheno[which(rnk > N * (1-prevalence))] <- 1

# linear regression
lm.out <- lm(pheno~g)
# phenotype value that is expected by genetic factors only
pheno. <- lm.out$fitted
# narrow-sense heritability
h1 <- var(pheno.)/var(pheno)
h1

> h1
[1] 0.05043295

注意するべきは、狭義遺伝率がぐっと小さくなっていること
これは分散を用いた狭義遺伝率が線形回帰的な枠組みを使っているのに対し、0/1表現型ではロジスティック回帰の方がよいことに対応する(はず)
では、線形回帰ではなく、ロジスティック回帰でやってみよう

# logistic regression
glm.out <- glm(pheno~g,family=binomial)
# phenotype value that is expected by genetic factors only
pheno. <- glm.out$fitted
# narrow-sense heritability
h1 <- var(pheno.)/var(pheno)
h1

> h1
[1] 0.2889976

線形回帰の場合よりはずいぶん上がったけれど、0.7からは程遠い
二値形質での場合は、遺伝率100というのは、なかなかに難しいはず。なぜなら、閾値モデルがそもそも連続値を0/1に峻別するという方法であって、そのこと自体に無理があるから