2011-07-04から1日間の記事一覧

(12) 塩基置換

R

#(U,C,A,G)=(1,2,3,4) b<-1:4 #b<-c("T","C","A","G") codons<-as.matrix(expand.grid(b,b,b)) aas<-c( "F","F","L","L","S","S","S","S","Y","Y","X","X","C","C","X","W","L","L","L","L","P","P","P","P","H","H","Q","Q","R","R","R","R","I","I","I","M…

(11) 有害事象

R

# 有害事象 N<-500 # 単位期間にDという有害事象が起きる確率 pD<-0.1 D<-runif(N)/pD # 単位期間にVという有害事象が起きる確率 pV<-1 V<-runif(N,min=0,max=1)/pV t<-seq(from=0,to=1,length=100) xlim<-ylim<-c(0,1) plot(rep(0,length(t)),t,xlim=xlim,y…

(10) ポアソン過程(ランダムに起きること)

R

一定間隔で起きること、まったくでたらめに起きること(どの時点も同じ確率で起きうる〜ポアソン過程、生起間隔が指数分布)、起きるとなったら集中しておきがちなこと # ぽつりぽつりと起きる現象 m<-0.5 Ne<-20 proportionRunif<-0.3 proportionRnorm<-0.4 e…

(9) 2次元の的当て(2次元正規分布)

R

対立仮説が1点 # 2次元の的当て Npt<-50 xlim<-c(-1,1) ylim<-xlim ori<-c(0,0) mx<-0.2 my<-0 sdx<-1/5 sdy<-1/5 r<-1 # 2 は正規分布 rateParam<-0.1 XYs<-matrix(10,Npt,2) for(i in 1:Npt){ plot(ori[1],ori[2],cex=60,pch=19,col=1,xlim=xlim,ylim=yl…

(7) 対立仮説からの量的形質の分布、等分散・異分散(t-test)

R

等分散からのサンプリング # 対立仮説 N1<-5 N2<-5 m1<-500 m2<-510 v1<-10 v2<-10 # 異分散 Niter<-1000 ps1<-ps2<-rep(0,Niter) for(i in 1:Niter){ d1<-rnorm(N1,m1,sqrt(v1)) d2<-rnorm(N2,m2,sqrt(v2)) t.test.out1<-t.test(d1,d2,var.equal=TRUE) t.te…

(6) 帰無仮説からの量的形質の分布、等分散・異分散(t-test)

R

等分散からのサンプリング # 帰無仮説 N1<-5 N2<-5 m1<-500 m2<-500 v1<-10 v2<-10 # 異分散 Niter<-1000 ps1<-ps2<-rep(0,Niter) for(i in 1:Niter){ d1<-rnorm(N1,m1,sqrt(v1)) d2<-rnorm(N2,m2,sqrt(v2)) t.test.out1<-t.test(d1,d2,var.equal=TRUE) t.te…

(5) 量的形質の分布(t-test)

R

# 平均を変える N1<-8 N2<-8 m1<-50 m2<-50 v1<-5 v2<-5 d1<-rnorm(N1,m1,sqrt(v1)) d2<-rnorm(N2,m2,sqrt(v2)) mdeltas<-seq(from=-5,to=5,by=0.5) par(ask=TRUE) d1<-rnorm(N1,m1,sqrt(v1)) d2<-rnorm(N2,m2,sqrt(v2)) ylim<-range(c(d1,d2+min(mdeltas),d…

(4) ケース・コントロール比が大きいとき(コホート疫学)

R

# コホート # 2x2分割表 NOp<-FALSE N1<-10000 N2<-20 a<-1000 A<-c(rep(1,a),rep(2,N1-a)) A<-sample(A) xlim<-ylim<-c(0,1) bs<-seq(from=1,to=N2-1,by=1) t<-seq(from=0,to=1,length=100)*2*pi ra<-0.2 rb<-ra max<-0.25 may<-0.5 mbx<-0.75 mby<-0.5 …

(3) 対立仮説からの2x2表ランダムサンプリング

R

N1<-100 N2<-150 p1<-0.3 p2<-0.4 Niter<-10000 p.out<-rep(0,Niter) v1<-c(rep(1,N1),rep(2,N2)) for(i in 1:Niter){ v2<-c(sample(c(1,2),N1,replace=TRUE,prob=c(p1,1-p1)),sample(c(1,2),N2,replace=TRUE,prob=c(p2,1-p2))) #print(table(v1,v2)) #p.out…

(2) 帰無仮説からの2x2表ランダムサンプリング

R

N1<-100 N2<-150 p1<-p2<-0.3 Niter<-10000 p.out<-rep(0,Niter) v1<-c(rep(1,N1),rep(2,N2)) for(i in 1:Niter){ v2<-sample(c(1,2),N1+N2,replace=TRUE,prob=c(p1,1-p1)) #print(table(v1,v2)) #p.out[i]<-chisq.test(c(rep(1,N1),rep(2,N2)),c(v1,v2))$p.…

第6章 Rでシミュレーション

R

遺伝統計学集中講義のための資料の全体の目次 Rでシミュレーション (1) 2x2分割表(ケース・コントロール スタディ) (2) 帰無仮説からの2x2表ランダムサンプリング (3) 対立仮説からの2x2表ランダムサンプリング (4) ケース・コントロール比が大きい…

(1) 2x2分割表(ケース・コントロール スタディ)

R

# 2x2分割表 NOp<-FALSE N1<-100 N2<-100 a<-50 A<-c(rep(1,a),rep(2,N1-a)) xlim<-ylim<-c(0,1) bs<-seq(from=1,to=N2-1,by=1) t<-seq(from=0,to=1,length=100)*2*pi ra<-0.2 rb<-ra max<-0.25 may<-0.5 mbx<-0.75 mby<-0.5 p.out<-rep(0,length(bs)) fo…

(8) 量的形質の増減(paired t-test)

R

# paired test N<-50 par(ask=TRUE) Pre<-rnorm(N) d<-runif(N)^2 pairedTps<-rep(0,N+1) rankSumps<-rep(0,N+1) for(i in 0:N){ Post<-Pre+d*c(rep(-1,i),rep(1,N-i)) plot(rep(1,length(Pre)),Pre,pch=18,col=3,xlim=c(0,3),ylim=range(c(Pre+max(d),Pre-m…