sample関数 cbind関数

• 入力値
  • Nm:マーカー数
  • p:マーカーごとのアレル頻度。以下の例では、0-1の一様乱数でランダムに振ってある
  • f:マーカーごとのディプロタイプ頻度に関するHWEからの逸脱度
  • r:マーカーごとのリスクの強さ。以下の例では、平均0、標準偏差0.01で振ってある。すべてのマーカーは程度の差こそあれ、リスクに寄与するとみなしている。また、アレル別・ディプロタイプ別のリスクの強さは、個々のマーカーごとに集団全体でのリスクの平均値が（HWEならば）０になるように与えているので、実際には、このrで与えるのは、マーカーごとのリスクの集団全体での分散の多寡になる
  • Ns:サンプル人数

Nm<-1000
p<-runif(Nm)
f<-rnorm(Nm,mean=0.0001,sd=0.000)
r<-rnorm(Nm,mean=0,sd=0.01)
Ns<-1000

• 関数の使用

dataDiplotype<-SimDiplotypeData(p,f,Ns)
dataRisk<-SimRiskData(p,f,r,Ns)

mean(apply(dataDiplotype,1,sum))
var(apply(dataDiplotype,1,sum))
hist(apply(dataDiplotype,1,sum))

mean(apply(dataRisk,1,sum))
var(apply(dataRisk,1,sum))
hist(apply(dataRisk,1,sum))

• • SimDiplotypeData
    • dataDiplotypeは列数Nm、行数Nsの結果。３ディプロタイプが0,1,2（ホモ、ヘテロ、逆ホモ）で表されている
  • SimRiskData
    • dataDiplotypeは列数Nm、行数Nsの結果。マーカーごとにがディプロタイプ別のリスクの値が表されている
    • アレル頻度pの座位においては、アレルごとのリスクをr(1-p), -rpとおくと、Additive modelにおいて３ディプロタイプのリスクは2r(1-p),r(1-2p),-2rpとなり、リスクの平均値は０になり、リスクの分散は2r^2p(1-p)となるので、SimRiskDataでは、マーカーのリスクをアレル頻度ともう一つのパラメタでこのように与えることで、マーカーごとの集団全体でのリスクを０にしつつ、マーカーごとのリスクへの寄与分はrによって増減させている

SimDiplotypeData<-function(p,f,Ns){

Nm<-length(p)
ret<-rep(0,Ns)
for(i in 1:Nm){
dif<-p[i]*(1-p[i])*f[i]
pr<-c(p[i]^2+dif,2*p[i]*(1-p[i])-2*dif,(1-p[i])^2+dif)
if(i==1){
ret<-sample(c(0,1,2),Ns,replace=T,prob=pr)

}else{
ret<-cbind(ret,sample(c(0,1,2),Ns,replace=T,prob=pr))
}
}

return(ret)
}

SimRiskData<-function(p,f,r,Ns){

Nm<-length(p)
ret<-rep(0,Ns)
for(i in 1:Nm){
dif<-p[i]*(1-p[i])*f[i]
pr<-c(p[i]^2+dif,2*p[i]*(1-p[i])-2*dif,(1-p[i])^2+dif)
risk<-c(2*r[i]*(1-p[i]),(1-2*p[i])*r[i],-2*p[i]*r[i])
if(i==1){
ret<-sample(risk,Ns,replace=T,prob=pr)

}else{
ret<-cbind(ret,sample(risk,Ns,replace=T,prob=pr))
}
}
return(ret)
}