2014年度マイコースプログラム

4回生K.K. @統計遺伝学分野

10月9日(木)

傾向スコア分析 傾向スコアについてRでやっているサイトがあったので、それをほぼもろパクリして自分でやってみた。 install.packages("Matching") library(Matching) data(lalonde) logi <- glm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u…

10月8日(水)

検出力について 数理統計の教科書を読んでいたら、検定における「検出力」に関する内容があって、すぐにはピンと来なかったので、Rで検出力を求める関数をつくってみた。 (ここへ来てより道しすぎな感は否めないが…。)正規母集団 N(mu,1)からn個の無作為…

10月3日(金)

plot3dについて 先生のアドバイスを頂いて、plot3dについてようやく分かった。あるサイトに library(rgl) z <- matrix(1:20, ncol = 4, nrow = 5) x <- c(1, 2, 3, 4) # z 行列の各行の座標 y <- c(1, 2, 3, 4, 5) # z 行列の各列の座標 plot3d(x, y, z) こ…

10月2日(木)

回帰モデルによる交絡の制御 2 ただの入力ミスだった。 k0 <- 0.2 k1 <- 0.3 k2 <- 2 N <- 10000 n <- 100 out. <- function(x){ k0+k1*x[1]+k2*x[2]+rnorm(1) } a1.d <- sample(20:80,N,replace=T) a2.d <- sample(0:1,N,replace=T) a3.d <- apply(cbind(a1…

10月1日(水)

ロジスティック回帰モデル 3Dプロットは僕の理解の範疇を超えているようである。意味が分からなさすぎて発狂しそうになったので、諦めてperspでグラフをつくってみた。 age <- c(18,21,22,25,26,28,33,34,35,37,42,47,55,56,58,61,65,68,75,77) out <- c(0…

9月30日(火)

ロジスティックモデル 「ロスマンの疫学」p.303のデータにおいて、モデル式をリスク=とし、最尤法による係数推定を行った。 age <- c(18,21,22,25,26,28,33,34,35,37,42,47,55,56,58,61,65,68,75,77) out <- c(0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,1,1,1,1,1) d …

9月29日(月)

今日で、「ロスマンの疫学」を最後まで読み切った。 とりあえず、Rでやってみたい題材には困らなさそうである。ただ、統計の教科書で回帰分析について読んでから始めたいと思ったので、あとは統計の教科書を読んでいた。徐々に残りの期間も少なくなってきて…

9月22日(月)

今日は、「ロスマンの疫学」の第10章の残りと、第11章(交互作用についての章)を読んだ。標準化についてRでやってみようと思っていたが、基準を決めて重み付き平均を出すだけで、そんなにややこしい話ではないみたいだったので、スルーした。交互作用は興…

9月19日(金)

データの層化 2 kを0から1までの実数として、曝露率をe=0.5+k*xとし、k(すなわち交絡と曝露率の相関の強さ)によって層化の影響がどのように変わるか調べてみた。 N <- 1000 e <- 0.5 b <- 0.05 c.max <- 0.5 alpha <- 0.1 n <- 7 myfunc <- function(k){ d…

9月18日(木)

データの層化 交絡因子と曝露率に関係があったほうが層化の威力が分かりやすいと思ったので、e=c+0.25として計算してみた。 N <- 1000 e <- 0.5 b <- 0.05 c.max <- 0.5 alpha <- 0.1 n <- 5 myfunc <- function(){ deta <- matrix(nrow=N,ncol=3) deta[,1] …

9月17日(水)

昨日は、「ロスマンの疫学」の第7章から第10章の途中まで読み進めた。第9章で、この前のケース・コントロール研究のシミュレーションの結果にも関連しそうな話題も出てきたが、統計学的な理解が追い付かなかったので、一旦保留にした。とりあえず今日は、第1…

9月12日(金)午後

リード・フロストモデル 拡張版 いまちょうどデング熱が流行っているということで、蚊によって媒介される病気のモデルをつくってシミュレーションしてみた。 (完全に「やってみた」だけなので、だからどうということはないのだが…。)仮定として、(1)集団…

9月12日(金)午前

リード・フロストモデルによるシミュレーション 一応ネットでも調べておくか…と思って「リード・フロストモデル」と検索したら、自分の昨日の記事が4番目くらいに出てきて唖然とした。 ネット検索において、価値のある情報が上位にくるわけではないというの…

9月11日(木)

累積症例対照研究のシミュレーション まとめ 2 結論から言うと、どうして分散がこのように変化するのかの説明は出来なかった。ある多項分布にしたがう4つの確率変数のかけ算・割り算をしたときに、計算結果はどういった分布を示すか、という問題に帰着できて…

9月10日(水)午後

累積症例対照研究のシミュレーション まとめ まず、「ロスマンの疫学」に累積症例対照研究は、発症する人が多くなるとオッズ比を過大評価としてしまうと書いてあったので、中央値を見てその確認をしてみる。p/q=2のとき(グラフのx軸はpとなっているが、qの…

9月10日(水)午前

累積症例対照研究のシミュレーション 3 先生からの助言を受けて、先生のおっしゃる通り1~4の4つの数字で集団を表すというのはもっとも効率のよいアイデアだと思ったので、それで箱ひげ図を描く関数を作り直してみた。下は、同じ処理を今までのものにさせた…

9月9日(火)午後

密度症例対照研究のシミュレーションまとめ 2 各パラメータにおける分散についてグラフにしてみた。◎e1についてe1が0.5に近い値をとっているときのほうが分散の値が小さいことは確かなようである。◎p/qについてp/qの値が大きいほうが分散の値は大きくなった…

9月9日(火)午前

累積症例対照研究のシミュレーション 2 関数のバグはとれた。 mat<-function(N,e,p,q){ a<-N*10000 b<-e*10000 r<-runif(a,0,1) r[1:b][r[1:b]>=1-p]<-2 r[(b+1):a][r[(b+1):a]>=1-q]<-2 s<-rep(c(1,0),c(b,a-b)) return(rbind(r,s)) } box<-function(N,e,p,…

9月8日(月)午後

密度症例対照研究のシミュレーション 2 それぞれの結果に対する理由の説明を試みたが、思ったよりも難しくて、ちょっといい説明が分からなった。 ミーティングでの話題にするとして、一旦保留にすることにした。結果を見て思ったのは、思ったよりもバラつき…

9月8日(月)午前

密度対症症例研究のシミュレーション 結局土日は研究室には来なかった。 けど、土日はある程度休んでもいいんじゃないかと思い始めてきたので、進行度がやばい!ってなったらまた考えることにする。一応土日の間ちょくちょくケース・コントロール研究のモデ…

9月5日(金)午後

ガリガリ計算して確率の対数をとって、もっともらしいpとqの値をそれぞれ0.01~0.05まで0.0050刻みで選べる関数をつくった。 p<-runif(1,0.01,0.05) q<-runif(1,0.01,0.05) Surv.time<-function(p,q,l){ cnt<-1 repeat{ if (runif(1)< p) {x<- c(cnt,1);brea…

9月5日(金)午前

生存曲線のシミュレーション 4 昨日のコード、ところどころ修正したミスを書き換えてなかった部分があったので、ちゃんとしたものを記録として残しておく。 Surv.time<-function(p,q,l){ cnt<-1 repeat{ if (cnt==l){x<- c(cnt,0);break} else if (runif(1)<…

9月4日(木)午後

生存曲線のシミュレーション 2 について、 ri<-function(d,t){ if (length(d[,1][d[,1]>=t])==0)return(1) else return(length(d[,1][d[,1]>=t]) }としておかないと、length(d[,1][d[,1]>=t])==0のときにS.hatがNaNとなってしまうことに気が付いた。(「1」…

9月4日(木)午前

生存曲線のシミュレーション 時刻を0,1,2…とし、時間1の間の死亡率を、脱落率をとするモデルを考えると、 Surv.time<-function(p,q,l){ cnt<-1 repeat{ if (cnt==l){x<- c(cnt,0);break} else if (runif(1)< p) {x<- c(cnt,1);break} else if (runif(1)< q) …

9月3日(水)午後

「ロスマンの疫学」を読んで、そこに出てきた話題を調べたりしているうちに、午後は終わってしまった。 マイコースの目標について 「ロスマンの疫学」では、”疫学とは通常「疾病頻度の分布と決定要因に関する学問」、あるいはもっと簡単に「病気の発生に関す…

9月3日(水)午前

「ロスマンの疫学」を読み進めているところである。先生とのミーティングで、「ロスマンの疫学」を使って最終的にどうしたいのか、というビジョンがあったほうがいいという話があったので、もう少し考えておくことにする。また、昨日の生存曲線について、は…

9月2日(火)午後

生存曲線 3 まず、カプラン・マイヤー法について簡単にまとめる。カプラン・マイヤー法とは、時刻までの生存率をと推定するものである。(ちなみに上の式書くのかなり苦労した。)ここで、式の中の は、時点 における死亡数、は時点 の直前まで死亡しうる可…

9月2日(火)午前

mychisq.test()の検証 (やり直し) 先生のご指摘をいただいて、同じデータセットを用いてそれぞれのp値を出した方がいいことに気づく。どうしてこんなことに気が付かなかったのかと思う。というわけで、やり直し。 s <- function(){ sample(1:3,20,replace=…

9月1日(月)午後

mychisq.test()の検証 2 母分散が同じがどうかの検定をしてみる。 s <- function(){ sample(1:3,20,replace=TRUE) } s.t <- function(){ table(s())}#テーブル作成 x<-c() for(i in 1:100)x<-c(x,mychisq.test(s.t())) y<-c() for(i in 1:100)y<-c(y,chisq.t…

9月1日(月)午前

いよいよ正規のマイコースプログラム期間開始である。気を引き締めてかかりたいところである。 mychisq.test()の自作 とりあえず求めるものは出来上がった。 mat1<-function(n,N){ if (N==0) return (matrix(0,ncol=n)) else if(n==1)return(matrix(N)) else…