リスク分布のプロット。

  • リスク
    • 遺伝要因と環境要因の和とする。
    • 遺伝要因と環境要因は独立とする。
    • 遺伝要因も環境要因も平均0で分散vg,veの正規分布とする。
#Nは集団の人数
N<-100000
#gvfrac 遺伝率
gvfrac<-0.8
#遺伝要因を正規分布で仮定する
g<-rnorm(N,mean=0,sd=sqrt(gvfrac))
#環境要因を正規分布で仮定する。分散は1-遺伝率
e<-rnorm(N,mean=0,sd=sqrt(1-gvfrac))
#遺伝要因と環境要因の和が個人のリスク
m<-g+e
#疾患の有病率を与える
prevalence<-0.01
#リスクが大きいと発病するとし、上位から有病率分が発病するとしてそれをリスク閾値とする
threshold<-1-prevalence
#発病閾値のリスク値を取り出す
quant<-quantile(m,threshold)
#発病閾値以上の個体の番号リスト
case<-which(m>quant)
#集合演算を用いて、それ以外の個体の番号リスト作る
cont<-setdiff(c(1:length(m)),case)
#ケースの遺伝要因の分布がよく見えるようにしたいので、まず、ケースのみでその遺伝要因のヒストグラムを描いて、その描出条件を保管する
hcase<-hist(g[case])
#全個体の遺伝要因の分布をする場合の描出条件も保管する
hall<-hist(g)
#ケースの描出条件のうち、ヒストグラムの各棒の幅を用いて、全個体のヒストグラムを描きたいので、棒の割り振り条件を作る
br<-seq(from=min(hall$breaks),to=max(hall$breaks),by=hcase$breaks[2]-hcase$breaks[1])
#棒の幅をケースのみ描出条件に変更して全個体のヒストグラムを描いたとしたときの描出条件をほぞ関する
h<-hist(g,breaks=br)
#コントロールを青の斜線で描く。ケースの棒の最大長までで描き、それより多い部分ははしょる条件
hist(g[cont],col="blue",breaks=br,ylim=c(0,max(hcase$counts)),density=20,main="MAIN TITLE",xlab="Xlabel",ylab="Ylabel")
#重ね描き
par(new=T)
#ケースは赤の斜線。
hist(g[case],col="red",breaks=br,ylim=c(0,max(hcase$counts)),density=20,main="",xlab="",ylab="")
    • 全員を描く場合には、最後の3行を次のように変える
#コントロールを青の斜線で描く。ケースの棒の最大長までで描き、それより多い部分ははしょる条件
hist(g[cont],col="blue",breaks=br,ylim=c(0,max(h$counts)),density=20,main="MAIN TITLE",xlab="Xlabel",ylab="Ylabel")
#重ね描き
par(new=T)
#ケースは赤の斜線。
hist(g[case],col="red",breaks=br,ylim=c(0,max(h$counts)),density=20,main="",xlab="",ylab="")