ラベル R の投稿を表示しています。 すべての投稿を表示
ラベル R の投稿を表示しています。 すべての投稿を表示

2014年1月9日木曜日

フィッシャーイグザクトメソッド

fisher exact method(fisher検定)はノンパラメトリック検定の1つであり、少数例の時でも背景因子の偏りを確かめることが出来る。

正規性を仮定できる時のみに用いることが出来る上位互換としてカイ二乗検定がある。

研究に使う場合はp値に関しては保守的に考えるべきであると考えているし、ミススペシファイしてしまうリスクを負う必要はないと思っており、正直カイ二乗検定は捨ててしまってもいいと考えている。

さて、fisher検定をRで用いるときはbaseの中にfisher.testという関数がある。

行列等を指定し、大きな行列となる時は作業につかう領域をworkspaceで指定する。

例えば

fisher.test(matrix(2:5,2,2),workspace=100000)

という具合である。多くの場合はp値が欲しいことが多いので

fisher.test(matrix(2:5,2,2),workspace=100000)$p.value

とすれば得ることが出来る。

たまにオッズが欲しくなるが、その場合は

fisher.test(matrix(2:5,2,2),workspace=100000)$estimate

オッズの信頼区間が欲しい場合は(2かける2の時)

fisher.test(matrix(2:5,2,2),workspace=100000)$conf.int


tableのそれぞれのセルの数が知りたければfisher.testをtableに変更すればよい。

もちろんfor文でも使えるので、幾つもの検定を一気に行うことも可能。

しかし、p値が小さすぎるときはデータに戻ることが大切。

そもそもマーカー陰性がいなかったり罹患した人がいなかったりするとオッズが変な値になる(計算出来ない)

計算結果はすぐに出てくるが、解釈をちゃんと行うこと。

それが一番大事。



多分いちばん多い質問。FFFFFAQに認定。



2013年12月24日火曜日

Rでのデータの読み込み

Rでデータを読み込む方法はいろいろあるが個人的には2つしか使わない。

read.table("clipboard",header=T)
か
read.csv(file.choose(),header=T)


上はExcelとかからコピーしてそれをそのまま乗っけたい場合
下はcsvファイルを読み込みたいとき。

区切りの認識はオプションで

sep="\t"   ","  " "

等でつけることができる。もっとあるけど実務上はこのくらいで十分。

前向き試験ののデータでは比較的欠損等のデータを扱うことが多いが後ろ向きの場合はかなり乱暴なデータベースで有ることも覚悟の上。

まあそもそも前向き試験の解析をRで行うこと自体があまりないのだけれども。

だから欠損とかをどう扱うかの方がレクチャーの満足度をあげるかもしれないな。

多くの教科書にはその辺りのことは書いていない、書いているけど、その時はよく検討して欲しい。というような内容。そもそもいろんな場合が考えられすぎて一般化しにくいんだと思うけど。

SPSS、JMPともにこの辺りはかなりうまくやってくれる(丸文字等もきっちりよんでくれる)のであまり重要な話ではないけれど。

Rを使っている先生からは必ずと言っていいほど読み込めません。相談は一度はあるので。

QAに認定。ッて感じで。

2013年11月23日土曜日

There is always light behind the clouds.

全然ブログかけてなかった。明日は久しぶりに何もない休み。めざましかけずにとりあえず寝る。


11日 
統計ソフトの検討を行い情報収集。主にJMPとSPSSについて。夜はよく相談を受ける診療科の予演会。2週間後本番ということで、診療科としても一気に発表ムード?いえ、私のほうがそういう質問が増えるので、時間を管理して、すべての質問に応えることが出来るように。

12日
私自身上で述べた学会で研究グループの統計解析を担当することになっていたので、それ関連について調べる。アポイントありの統計相談が2件、突発的相談が4件ほど。

13日
N大の講義。本日はP2関連。かなり内容が盛りだくさんで面白かった。理論の整理をする時間が欲しい。合間を見つけて行おう。

14日
某協会の方が来院されるということで、準備等。建設的な議論となったと思います。(自分が総評価するのもおこがましいですが・・・)
また新たな前向き試験の話が1件。なかなか興味深く歯ごたえのある相談。エンドポイントの設定が肝。
あと治験のキックオフが1件。ヒアリングで質問した内容がいろいろと検討されていた。自分の臨床の勉強だけでなく、少しだけ貢献できているかなとか思って嬉しくなる。いや、油断は禁物。油断?
アポあり相談2件。あまり研究をしたことがない診療科の先生からの相談。(人数の関係でバイアスあり。)興味深い。また勉強することが増えた。今日の診療とかup to dateとか読み漁る。

15日
アポあり相談3件。1件は少し前にいなくなった先生から。非常勤で週1で来られているらしい。今日の相談の先生は3人が3人共熱意がすごくある。比較するものでもないかもしれないけど、すごいパッションを感じる。いや、みな有るのだと思うけど、私には特にそう感じる。
ヒアリングが1件。内容は言えない。
博士論文の発表を人に聞いてもらう。もうスグだ。

16日17日
博士論文の発表の準備。胃が痛い。

18日
今週と迫った学会の予演会。統計の部分、体裁等を少しコメントさせていただく。多分指摘はされないんだろうけど、でも、やっぱりああ統計的な曖昧さがないなあといろんな人に思って欲しい。
自分の発表の方もデータ解析前大詰め。しかし明日は博士論文の審査。終わったらやろう。

19日
審査当日。準備はしたつもりだったが、あまり良くなかった。半期伸ばすことにする。申し訳ない。悔しい。ただ現実を受け止め前に進む。正直まだあまりうまく整理できていないため、多くをかけない。
帰宅後データの解析のソースづくり。寝たのは・・・。

20日
気持ちの整理がつかなくても、次の発表は来る。データ解析の一日。ちょこちょこ突発的な相談や雑務に追われながらも、行う。解析は私一人のため何度かおこないたいと考えていたため、夕方までに一度完成させ、夜にもう一度。予防接種も行った。心はあまりでも体はそこまで悪い状態ではないようだ(あまり寝る時間がないけど)

21日
朝4時に起きて最終確認。5時に研究者の皆さんに送付。
その後少し統計相談を受けてから、夕方前に学会会場へ。
会場に向かう間にメールを頂き、追加解析。
なんとかコアメンバーの先生に助けていただきながら発表スライドを完成。
発表も私の担当の部分は無難に終わった。その後臨床医の先生のスライドの解釈を聞いてもっと勉強しないとダメだと思った。
その後、ご飯をご一緒させていただき、その後別グループに合流させて頂く。
明日は学会を楽しもう。

22日
解析から開放されしばし学術的のみをもって学会を堪能。
途中院内からの連絡に対応しながらも、ほぼ予定通り聞くことが出来た。
院外で研究をご一緒させて頂いている先生とも多く対面でお話させていただいたし、新しいネタがあるからという話も聞くことが出来た。数いる統計家の中から私に声をかけてくださるのは嬉しい限りだ。
実は明日も東京なのだが、さすがにいろんな疲れがドバドバ出てきているので、一度帰宅することにする。

23日
築地でセミナー。久しぶりに座学だったのとこの2週間の怒涛の様なあれでちょっと・・・・。でも、新たな知識も入ってきたし、明日の業務に活かしていこうと思います。というより少し知らない言葉とかフォロー出来てない通知とかがあったので、確認かな。うん。まだまだすべてをフォロー出来てない。勉強。で、夜は後輩と飯食って今これを書いている。

ということでお休み。いい夢は見たくない。私は安眠を望む。






2012年11月12日月曜日

並べ替える時はsortをいきなり使わないように。

データをソートするときは、sortを用いるのではなく、一度orderで順位を確認したあと、それを用いてソートする癖をつけよう。

そうでないと、2変量以上の時データがぐちゃぐちゃになってよくわからなくなる。

分割表で0のセルがあるとRcmdrでは解析できない。

Rコマンダでは背景因子の分析をするときに2*2で0のセルが一つ以上あれば解析できない。

多分内部ではchisq.testを使っているからだと思うが、医学のデータではかなり少数の検討をしなければならないこともあり、0のセルがあってもやるしかないこともある。

その時はコマンドfisherを使ってやるべきなのだが、僕がいつも使っているfisherというものは群馬大学の青木先生という方が作って公開しているもので、このページにかいている。
http://aoki2.si.gunma-u.ac.jp/R/fisher.html

ここの


source("http://aoki2.si.gunma-u.ac.jp/R/src/fisher.R", encoding="euc-jp")

という1分をRの画面に貼り付ければ出来る。(RコマンダではなくRコマンダではなく本体)

そして

x <- matrix(c(50,30,20,40), byrow=TRUE, ncol=2)
fisher(x)

とすれば計算してくれる。

xは行列、ここでは2*2の分割表を指定する。

byrowってのは行を優先してください。

ncol=2というのは列は2までにしてください。ってこと。


数字をc()で囲んでいるのはおまじないとして覚えてください。

あと、もしかするとfisher(x)とするとエラーが出てしまう可能性があります。

大概の場合そのエラーの文章をコピーしてgoogle先生に見てもらうと同じエラーを起こしている人がかならずいます。

確認してみてください。

2012年11月7日水曜日

file.choose()

file.choose()
と打てば、Office系の開くと同じように対話的にファイルを選ぶことが出来るようになる。
ただ、このコマンドのみでファイルを選んでもそのファイルの絶対パスを読みこんでくれるだけなので、本来は、

read.table(file.choose(),header=T)

などとして使うことが多い。もちろんread.csvでも可能。

dat1<-read.table(file.choose(),header=T)

などとして

head(dat1)
summary(dat1)

などを確認していくのがよい。

2012年11月5日月曜日

カプランマイヤー

データは何が起こるかわからないので、エクセルで作った場合でも最終的に値で貼り付け。

そして、日本語は取り除く方がよい。

変数に名前をつける癖は作ってるほうが良い。

膨大であればあるほど。

また出来ればいつも同じディレクトリでやるようにした方がいいが、それを意識するのが煩わしいと感じる人もいる。(そもそも自分がどこのディレクトリで作業するかはwindowsでは気にしないことが多い)

そのため、クリップボードからコピーするとき以外は、

read.csv(file.choose(),header=T)

を使うべきである。

そのままカプランマイヤーを書く場合は

library(survival)

で関数を呼び出し、(ないと言われればネットにつなげた状態でinstall.packages(survival)をやる。)

plot(survfit(Surv(time,event)~group,data=data),main=c(""),xlab=...)

という感じで。(Survの引数この順番だったかな・・・)

カプランマイヤーの拡張はまた質問された時にでもまとめます。

Bayesian continuous reassessment method (CRM) designs for Phase I dose escalation trials.

Rに既に実装済み。

bcrmという関数なので。
教科書通り、ネットにつながっている状態で

install.packages("bcrm")
library(bcrm)

で完了。
色々オプションがあるが、わからない人は使わないほうがいいので、これ以上はまとめない。

どうしても使いたければ統計家に。

2012年11月1日木曜日

Rコマンダーについて 追記有り さらに追記あり

職場でRコマンダーでの解析の講座のサポートをした。

そこで出てきたよくわからないエラー。昨日1日調べたが、原因不明だった。

・エクセルがそもそも読み込めない。
→エクセルを読みに行くためのパッケージがない?そのパッケージ名がわからず。

・エクセルでカテゴリカル・データを読み込んだが、層別因子として認識されない。
→エクセルの情報をクリップボードに貼り付けてやれば可能であった。そのため、エクセルの読み込み方の問題?
・rJavaというパッケージが入っていない。
→インストールで解決

・Macの場合Rcmdrがインストール出来無い?
→Macもってないのでわかんない。使いふるしでいいから余ってないかな。ただ、できてた人もいた模様。


また実習があるようなので、それまでに解決したい。やっぱ英語で使うべきか-?


追記12・11・01
X11をダウンロードして
ttp://cran.r-project.org/bin/macosx/tools/
から
ttp://cran.r-project.org/bin/macosx/tools/tcltk-8.5.5-x11.dmg
をダウンロード→
Rを起動→上のイラストからX11を起動→
library(Rcmdr)
で行けるみたい。


X11とか学部の時の授業で出てきたなーと思ったくらいだった気がする。

とりあえずこれで一件落着かな?


追記12/11/02
Mountain Lionだとうまくいかないとの報告有り。
MacOSが近くにないので、確認できないが・・・。

かうか?そこまでする必要はないと思うけど、Mac使いの先生多いんですよね・・・。

2012年10月18日木曜日

強制投入した場合の推定値の信頼区間

強制投入で解析したら推定値の信頼区間でないじゃないかという要望から。

confint(モデル式)

で出てきます。

2012年10月17日水曜日

青木先生の関数fisher exact testを使う時の注意点

http://aoki2.si.gunma-u.ac.jp/R/fisher.html
このページから。


source("http://aoki2.si.gunma-u.ac.jp/R/src/fisher.R", encoding="euc-jp")

をコピペしたらfisher exact testが出来る。

xに適当なデータを入れてやってみると

> fisher(x)
 以下にエラー asymptotic() :  関数 "printf" を見つけることができませんでした 

と出てくる事がある。

この場合はこのページ


の真ん中あたりで回答があり、

printf <- function(fmt, ...)
{
cat(sprintf(fmt, ...))
}

と打ち込めば使えるようになります。

何度か間違ってしまったのでまとめ。

all.R入れろって話ですかね。うーん。

あと、レトロな解析な場合は基本t検定とカイ自乗検定は使わずにマンホイットニー(ウィルコクソン)とフィッシャーでいったほうがいいと思います。


2012年10月5日金曜日

強制投入法以外で変数候補を探索する。Rで。


まずRの画面で

library(MASS)

と入力し実行。(ダウンロードは必要ないと思います。)

この宣言でステップワイズ法が使えるようになります。

デフォルトでstepという関数があるけど個人的にちょっと結果がみにくいと思っている。

計算自体は

stepAIC(coxph(Surv(time,outcome)~factor1+factor2+...,data=解析したいデータセットの名前),direction="backward")

で行えます。(環境によっては改行が入るかもしれませんが1行です)


coxph()の部分は考えるモデルで強制投入の時に作ったものと同じ。もちろんコックスハザード以外も使えますが、医学論文のための手法としてはこれかロジスティックでしょう。

あと

cox.analysis<-coxph(Surv(time,outcome)~factor1+factor2+...,data=data)


などとしていただいて、


stepAIC(cox.analysis,direction="backward")


としていただいても構いません。

後ろのdirectionの部分はbackward,forward,bothの3種類を指定できます。

減らしていく方法、増やしていく方法、どちらもやっていく方法です。


医学論文では基本的にいわゆるp値でしかみないけど、AICって指標を用いることでどの要因が効いているのかを評価できる方法もあります。

AICってのは日本人の赤池先生が作って世界中で使われている指標です。