power analysis
power= 1 - P(Type II error) = probability of finding an effect that is there
power↑ w/ sample size↑
SDとSEMの使い分け
SD: standard deviation
SEM: standard error of means
SEM=SD/√nが成立。
こちらから抜粋
Although the SD and the SEM are related (SEM=SD/√n), they give two very different types of information. Whereas the SD estimates the variability in the study sample, the SEM estimates the precision and uncertainty of how the study sample represents the underlying population. In other words, the SD tells us the distribution of individual data points around the mean, and the SEM informs us how precise our estimate of the mean is. It is therefore inappropriate and incorrect to present data only as the mean (SEM).
安定マッチング
医師臨床研修マッチングについて考える。
採用されているのは、Gale-Shapleyの安定マッチングシステム。
・安定なマッチングとは
今ペアになっていない 2 人で「今のペアを維持するよりもその 2 人が結婚した方が互いに幸せである」という組合せがないことが条件である。
・安定なマッチングだからといって、フルマッチする保証は無い
(絶望の定理: ある安定マッチングにおいてペアを作られなかった人は、どんな安定マッチングを用いてもペアを作ることができない)
・プロポーズをする側に圧倒的に有利なシステムである
・こちらにRソースを見つけたが、久しぶりにRに触れるこの状況ではかなり理解に時間を要しそう…
家系データシュミレーション
アレル頻度とは
ハプロタイプは、アレルの数珠つなぎ(tandem: 縦につないだもの)
N############### ### 1 ############### # あるマーカー # アレルの数 Na<-4 # アレル頻度 Pa<-c(0.3,0.2,0.1,0.4) # Pa<-c(3,2,1,4)のようにして、あとで頻度の和が１になるようにしてもよい # 念のため、アレル頻度の和が１になるように補正する Pa<-Pa/sum(Pa) ############### ### 2 ############### # 一様乱数からの発生 Pa<-runif(Na) # アレル頻度が1/Na付近に固まるので… # そうならないように道具をウェブ上から導入する install.packages("MCMCpack") library(MCMCpack) # Na個の足して１になる乱数を１セット得る Pa<-rdirichlet(1,rep(1,Na)) # できたPaの中身を見る Pa # 和が1担っていることを確認する sum(Pa) ############### ### 3 ############### # 同一染色体上にある多型の数Ns Ns<-4 # Ns個のマーカーのそれぞれのアレルの数をNaS NaS<-c(3,4,5,2) # アレル頻度をPaS PaS<-list() for(i in 1:Ns){ PaS[[i]]<-c(rdirichlet(1,rep(1,NaS[i]))) } NaS PaS hap1<-hap2<-rep(0,Ns) #父親由来、母親由来のハプロタイプ２つ用意する for(i in 1:Ns){ hap1[i]<-sample(1:NaS[i],1,prob=PaS[[i]]) #1〜NaS[i]番目までのアレルの中から１つ選ぶ。選び出される確率はPaS[i] hap2[i]<-sample(1:NaS[i],1,prob=PaS[[i]]) } hap1 hap2 ############### ### 4 ############### #関数にする(毎回上を使うのはめんどくさい) # 入力はマーカーのアレルNaS、アレル頻度PaSと作成するハプロタイプの本数N MakeHaplotype<-function(NaS,PaS,N){ # haplotypeはN x length(NaS)の行列に納める H<-matrix(0,N,length(NaS)) for(i in 1:length(NaS)){ H[,i]<-sample(1:NaS[i],N,prob=PaS[[i]]) } H } # 使ってみる N<-2 # 一人分の２本のハプロタイプ Haps<-MakeHaplotype(NaS,PaS,N) Haps ############### ### 5 ############### # 組換えをシミュレーションするために、２本のハプロタイプを作る # Ns個のマーカーの(Ns-1)個の「間」について、モルガン単位で遺伝的距離を定める Dg<-runif(Ns-1) # 組換えを起こして２つの組換え体を作ろう # マーカー間ごとに、交叉がおきるかどうかを確率的に決める rands.crossing<-runif(length(Dg)) # 交差がおきた箇所を取り出す Cross<-which((Dg-rands.crossing)>0) Cross # 交差が起きたかどうかの情報を用いて組換え体を作る hap1<-Haps[1,] hap2<-Haps[2,] rec.hap1<-hap1 rec.hap2<-hap2 if(length(Cross)>0){ for(i in 1:length(Cross)){ # 交差箇所より尾部について、rec.hap1,rec.hap2で交換する tmp.tail1<-rec.hap1[(Cross[i]+1):Ns] tmp.tail2<-rec.hap2[(Cross[i]+1):Ns] rec.hap1[(Cross[i]+1):Ns]<-tmp.tail2 rec.hap2[(Cross[i]+1):Ns]<-tmp.tail1 } } Cross hap1 hap2 rec.hap1 rec.hap2 ############### ### 6 ############### #やっぱり関数化 # 2本のハプロタイプと、多型間遺伝的距離から、組換え体２本を返す Rec.Haplotype<-function(hap1,hap2,d){ Ns<-length(hap1) rands.crossing<-runif(length(d)) # 交差がおきた箇所を取り出す Cross<-which((d-rands.crossing)>0) # 交差が起きたかどうかの情報を用いて組換え体を作る rec.hap1<-hap1 rec.hap2<-hap2 if(length(Cross)>0){ for(i in 1:length(Cross)){ # 交差箇所より尾部について、rec.hap1,rec.hap2で交換する tmp.tail1<-rec.hap1[(Cross[i]+1):Ns] tmp.tail2<-rec.hap2[(Cross[i]+1):Ns] rec.hap1[(Cross[i]+1):Ns]<-tmp.tail2 rec.hap2[(Cross[i]+1):Ns]<-tmp.tail1 } } matrix(c(rec.hap1,rec.hap2),byrow=TRUE,nrow=2) } # 使ってみる rec.haps<-Rec.Haplotype(Haps[1,],Haps[2,],Dg) rec.haps ############### ### 7 ############### # 父・母のジェノタイプをランダムに作る father<-MakeHaplotype(NaS,PaS,2) mother<-MakeHaplotype(NaS,PaS,2) father mother # 配偶子を作る father.gamates<-Rec.Haplotype(father[1,],father[2,],Dg) mother.gamates<-Rec.Haplotype(mother[1,],mother[2,],Dg) # 父由来・母由来の配偶子をランダムに選んで子のジェノタイプを決める from.father<-father.gamates[sample(1:2,1),] from.mother<-mother.gamates[sample(1:2,1),] offspring<-matrix(c(from.father,from.mother),byrow=TRUE,nrow=2) offspring father mother offspring ############### ### 8 ############### #関数化 MakeOffspring<-function(father,mother,d){ # 配偶子を作る father.gamates<-Rec.Haplotype(father[1,],father[2,],d) mother.gamates<-Rec.Haplotype(mother[1,],mother[2,],d) # 父由来・母由来の配偶子をランダムに選んで子のジェノタイプを決める from.father<-father.gamates[sample(1:2,1),] from.mother<-mother.gamates[sample(1:2,1),] offspring<-matrix(c(from.father,from.mother),byrow=TRUE,nrow=2) offspring } # 使ってみる offspring<-MakeOffspring(father,mother,Dg) father mother offspring ############### ### 9 ############### #「子は親よりも大きいID」をつけることに責任を持つ(よほど複雑な家系図でない限り)←プログラムを書きやすくするために # 両親と一人っ子 M.Oyako.Trio<-matrix(c(0,0,0,0,1,2),byrow=TRUE,ncol=2) M.Oyako.Trio # １人とその両親とそのさらに両親 M.3gen<-matrix(c(0,0,0,0,0,0,0,0,1,2,3,4,5,6),byrow=TRUE,ncol=2) M.3gen # 一組のいとこ # いとこペアの１人、その両親、その片方の両親(共有される祖父母) # 共有された祖父母のもう一人の子、その配偶者、その子(いとこペアのもう一人) M.cousins<-matrix(c(0,0,0,0,1,2,0,0,1,2,0,0,3,4,5,6),byrow=TRUE,ncol=2) M.cousins ############### ### 10 ############### # M.cousins[,1]はM.cousings[,1]の第１列のこと n.rows<-length(M.cousins[,1]) # nrows<-nrow(M.cousins)でも同じこと # 人数=n.rows 分のジェノタイプは マーカー数 x 2 行列が人数分あるので、それを納めるリストを作る sim.G<-list() # 第1行から第n.rows行まで、順番に処理するループ for(i in 1:n.rows){ # ここで処理される人はi行目の人 # その人の両親は tmp.father<-M.cousins[i,1] tmp.mother<-M.cousins[i,2] # tmp.father,tmp.motherのジェノタイプは if(tmp.father==0){ father.genotype<-MakeHaplotype(NaS,PaS,2) }else{ # 親のジェノタイプはこのループ処理ですでに決まっている father.genotype<-sim.G[[tmp.father]] } if(tmp.mother==0){ mother.genotype<-MakeHaplotype(NaS,PaS,2) }else{ mother.genotype<-sim.G[[tmp.mother]] } # 両親のジェノタイプが決まったので、この人のジェノタイプを決める sim.G[[i]]<-MakeOffspring(father.genotype,mother.genotype,Dg) } sim.G M.cousins ############### ### 11 ############### # 引数はM.cousins,NaS,PaS,Dg Simulate.Genotype<-function(M,NaS,PaS,Dg){ # M.cousins[,1]はM.cousings[,1]の第１列のこと n.rows<-length(M.cousins[,1]) # nrows<-nrow(M.cousins)でも同じこと # 人数=n.rows 分のジェノタイプは マーカー数 x 2 行列が人数分あるので、それを納めるリストを作る sim.G<-list() # 第1行から第n.rows行まで、順番に処理するループ for(i in 1:n.rows){ # ここで処理される人はi行目の人 # その人の両親は tmp.father<-M.cousins[i,1] tmp.mother<-M.cousins[i,2] # tmp.father,tmp.motherのジェノタイプは if(tmp.father==0){ father.genotype<-MakeHaplotype(NaS,PaS,2) }else{ # 親のジェノタイプはこのループ処理ですでに決まっている father.genotype<-sim.G[[tmp.father]] } if(tmp.mother==0){ mother.genotype<-MakeHaplotype(NaS,PaS,2) }else{ mother.genotype<-sim.G[[tmp.mother]] } # 両親のジェノタイプが決まったので、この人のジェノタイプを決める sim.G[[i]]<-MakeOffspring(father.genotype,mother.genotype,Dg) } sim.G } # 使ってみる sim.G<-Simulate.Genotype(M.cousins,NaS,PaS,Dg) sim.G ############### ### 12 ############### # 染色体ごとにNaS,PaS,Dgを作る # ここでは簡単のために、すべて同じものとする num.Chromosomes<-22 NaS.list<-list() PaS.list<-list() Dg.list<-list() for(i in 1:num.Chromosomes){ NaS.list[[i]]<-NaS PaS.list[[i]]<-PaS Dg.list[[i]]<-Dg } sim.Gs<-list() for(i in 1:num.Chromosomes){ sim.Gs[[i]]<-Simulate.Genotype(M.cousins,NaS.list[[i]],PaS.list[[i]],Dg.list[[i]]) } ############### ### 13 ############### # 入力値はM.cousins,NaS.list,PaS.list,Dg.list Simulate.Genotype.Multi<-function(M,NaS.list,PaS.list,Dg.list){ num.Chromosomes<-length(NaS.list) sim.Gs<-list() for(i in 1:num.Chromosomes){ sim.Gs[[i]]<-Simulate.Genotype(M.cousins,NaS.list[[i]],PaS.list[[i]],Dg.list[[i]]) } sim.Gs } # 使ってみる sim.Gs<-Simulate.Genotype.Multi(M.cousins,NaS.list,PaS.list,Dg.list) ############### ### 14 ############### # 繰り返し回数 Niter<-10 sim.Gs.list<-list() for(i in 1:Niter){ sim.Gs.list[[i]]<-Simulate.Genotype.Multi(M.cousins,NaS.list,PaS.list,Dg.list) }
まとめ
●必要な変数
No. markers→genetic distance btw markerも必要
No. allels
Allele.freq
familytree
●simulationする
●あとはfamilyの数を増やしたり、減らしたり
数学いらずの医科統計学
数学いらずの医科統計学 より
PART4 連続変数
p.61
■散らばりや分布を示すためのデータのグラフ化
散布図
Tm<- mean(c(37, 36, 37.1, 36.2, 37.3, 36.8, 37.0, 36.3, 36.9, 36.7, 36.8)) #観測したデータの平均をTmとする t1<- rnorm(130, Tm, 0.8) #t1はTmの集団から、130個ランダムに抽出したもの T1<- cbind(rep(1,length(t1)), t1) t2<- rnorm(12, Tm, 0.8) #t2はTmの集団から、12個ランダムに抽出したもの T2<- cbind(rep(2,length(t2)), t2) T<- rbind(T1,T2) plot(T, axes=FALSE, xlab=' ', ylab='体温(℃)') axis(1, c(0:3), c(" ", "全データ", "一部のデータ", " ")) axis(2, seq(0, 40, by=5))
箱ヒゲ図
par(mfrow=c(1,2)) Tm<- mean(c(37, 36, 37.1, 36.2, 37.3, 36.8, 37.0, 36.3, 36.9, 36.7, 36.8)) #観測したデータの平均をTmとする t1<- rnorm(130, Tm, 0.8) #t1はTmの集団から、130個ランダムに抽出したもの t2<- rnorm(12, Tm, 0.8) #t2はTmの集団から、12個ランダムに抽出したもの boxplot(t1, t2, ylim=c(30,40), axes=FALSE, xlab=' ', ylab='体温(℃)') axis(1, c(0:3), c(" ", "全データ", "一部のデータ", " ")) axis(2, seq(0, 40, by=5)) boxplot(t1, t2, range=0, ylim=c(30,40), axes=FALSE, xlab=' ', ylab='体温(℃)') # 髭の範囲をデータの最大値から最小値までにする axis(1, c(0:3), c(" ", "全データ", "一部のデータ", " ")) axis(2, seq(0, 40, by=5))
頻度分布ヒストグラム
Tm<- mean(c(37, 36, 37.1, 36.2, 37.3, 36.8, 37.0, 36.3, 36.9, 36.7, 36.8)) #観測したデータの平均をTmとする t1<- rnorm(130, Tm, 0.8) #t1はTmの集団から、130個ランダムに抽出したもの par(mfrow=c(2,2)) hist(t1,main="defalt") hist(t1,breaks=seq(34,40,1), xlab='体温(℃)', main="breaks=seq(34,40,1)") #binの範囲をどんどん小さく… hist(t1,breaks=seq(34,40,0.5), xlab='体温(℃)', main="breaks=seq(34,40,0.5)") hist(t1,breaks=seq(34,40,0.1), xlab='体温(℃)', main="breaks=seq(34,40,0.1)")
