Statistical Code for Clinical Research Papers in a High-Impact Specialist Medical Journal

インパクトの大きな医学研究で使用した統計ソフトのコードを公開するべきだという論評

Assel M, Vickers AJ. Statistical Code for Clinical Research Papers in a High-Impact Specialist Medical Journal. Ann Intern Med. 2018;168:832–833. doi: 10.7326/M17-2863

<原文は自前のデータを提示してはいますが、ほぼ、思うところを述べている論評です>


1.      定量的医学研究・分析の中核となるはずのコーディングを教える教育プログラムが少ない、としているが、あるのか日本には? (医者でちょっと深くやろうとすると、システムエンジニアや数理統計学者の真似みたいなことに向かおうとするような人もいる。自分のサイロに閉じこもって医学研究と言う当初の目的から離れていったりしていないかい?)

First, software practices and principles should become a core part of biostatistics curricula, regardless of the degree (under- or postgraduate) or subject (biostatistics, public health, or epidemiology). Given that students will have to write code when they perform analyses as practicing investigators, we question why so few degree programs in quantitative medical science teach good coding practice.


2.      教室内で、統計解析のコードをピアレビューすべきだというが、教室内にいるのかな自分以外でコードがわかる人が? 仮にいたとして、他人の書いたプログラムを読むのはかなり厳しい。それより、ダブルプログラミングして答え合わせする方が現実的。

Second, statistical code should undergo intramural peer review. Colleagues should routinely share code to receive constructive criticism just as they share drafts of scientific papers.


3.      自分自身が「エレガント」だと思うような綺麗なコードは、おそらく他人から見て意味不明で、何でこれで正しい答えが出るのが不思議なやつになっている(と思う)。若干冗長でも、何をやっているのか解りやすいほうが後日のデバックや学習用には向いている。

Finally, code associated with published research should be archived. Doing so would not only improve transparency and reproducibility but also help to ensure that investigators write better-quality code…. We believe well-written code has more than cosmetic value and that dirty code may lead to scientific errors.


We urge the medical research community to take immediate remedial action」:自分では他の人ができないようなことができるから、他人より一歩前に出ているという認識があって、皆ができるようになると「なぁーんだ」てな感じで埋もれてしまいそうだから、実は現状がちょうど良いというような気もしなくない。

# ちなみに、この論評で推進を論じているcode share、「格安航空会社のホームページでチケットを予約したのに、空港で待っていたのは国内大手の機材だった」感じです。

 

 

 

TRUMPET ENSEMBLE SEED

SEED – trumpet ensemble

2018/4/27(FRI) @ Casa Classica

ぱんだウインドオーケストラのトランペットパートメンバーによるアンサンブルSEEDのライブに行ってきました。出演されたのは、笠原日向さん・片野和泉さん・佐藤玲伊奈さん・重井吉彦さん・中村諒さんの5名。会場には若い方が多かった中、舞台目の前に陣取ったテーブルに私と、メンバーの片野さんのお母さんが世代が上で保護者席の様な感じになっていました。そこまで広くない会場で、トランペット5本が目の前で鳴ると音圧がすごい迫力でした。中村さんの相変わらずの安定した音は聴いていて安心できます。笠原さん、通常は耳に突き刺さるようなハイノートを楽しませてくださいますが、この日は「3万円ほどでヤフオクで落札した」というバストランペットにトロンボーンのマウスピースを付けてアンサンブルに入っておられました。みなさんおしゃべりは苦手な様ですが、前半おしゃべりが少なくどんどん進んでいくと少しバテから回復する時間がないのでは、という感じでした。チラシは重井さんがスマホアプリで作成されたものだそうです。

https://drive.google.com/file/d/0B23Rvhh_yGe1SkRMeUtUV2U5a19MdURVSmcyTTIxdnpaXzdz/view?usp=sharing

全体を録画したつもりなのですが、うまく取れていたのは最後の曲だけでした。(失礼)

なお、この楽曲はJASRAC DB上はPublic Domainなのですが、youtubeに掲載しましたところ「Hexacorp (music publishing)」というところが何らかの権利を有していると主張してメールを送ってきました。動画の掲載は継続していてよいが、広告が現れる場合がある(広告収入がHexacorpに入るものと思われます)と言う連絡を受けました。継続して掲載可能という事ですので特に異議を主張していません。将来条件が変わり削除を求められるようなことになりましたら、削除するかもしれませんのであしからず。

in-store-live of Sabão (ex Hysteric Blue)

Sabão (ex Hysteric Blue) インストアライブ

実は、このSabaoというグループ、寡聞にして2018年現在活動を続けておられることすら認識していなかったのですが、今回のライブの数日前にたまたま検索する気になって、このインストアライブを知ることになりました。目的は1999年紅白歌合戦でも披露した「春~spring~」 です。その紅白を見た当時、我が家は私の米国留学中で、数少ない日本からの放送(当時)で聴いた曲です。当時、米国から見て太平洋の向こう側のNHKホールのステージで歌っておられて、文字通り手の届かないような場所にいらっしゃったお二人が、本当に手の届きそうな目の前で曲を演奏しているのを見て感動です。そして、約20年と言う長い年月が経ってますます歌声に磨きがかかっていたのも素敵でした。特にファンという訳でもなかったので、大ヒットしたspring以外の他の曲はよく知らないわけですが、今回のライブで少し聴いてみようかなと言う気になりました。

リハーサル風景の動画

ライブ中からサイン会終了までは撮影・録音・録画禁止という事でアナウンスされていました。このページに貼ってある画像や動画は、リハーサルの音出しのと時にほんの少しだけ、録画・撮影しました。じつは音出しの時にはちゃんとバランスよく音が出来上がっていたのに、本番が始まるとボーカルの音がほとんど聞こえません。にもかかわらず、vocalのtamaさんはミキサーに向かって、音量を下げるようにと言う合図を何回も出します。別のスタッフは客席の後ろからミキサーのところまで駆け寄ってきて何か指示しています。2曲目まで終わったところで、解ったのですが、内音はvocalの声量が大きくてバランスが悪かったという事です。で、外音にはvocalが出ていなかったこともステージに伝わりました。結局DJのマイクをvocalが使用することで、残りの2曲は無事バランスよく楽しむことができました。4曲目が目的の「春~Spring~」でした。

 

 

PlayPlay

R のアンインストール MacOS

MacOS での統計ソフトRのアンインストール

MacにインストールしたRが古くなり、パッケージが動かない状況も見受けられるので一旦アンインストールして、その上で最新版をインストールしようと思いました。でも、RってApplestoreからインストールしたアプリでないのと、アンインストーラが付属していないので、どうやって削除するのか。まずは下調べです。

<https://cran.r-project.org/doc/manuals/r-release/R-admin.html#Uninstalling-under-macOS>

4.2 Uninstalling under macOS

R for macOS consists of two parts: the GUI (R.APP) and the R framework. The un-installation is as simple as removing those folders (e.g. by dragging them onto the Trash). The typical installation will install the GUI into the /Applications/R.app folder and the R framework into the /Library/Frameworks/R.framework folder. The links to R and Rscript in /usr/local/bin should also be removed.

If you want to get rid of R more completely using a Terminal, simply run:

sudo rm -rf /Library/Frameworks/R.framework /Applications/R.app \
/usr/local/bin/R /usr/local/bin/Rscript

The installation consists of up to four Apple packages:18 org.r-project.R.el-capitan.fw.pkg, org.r-project.R.el-capitan.GUI.pkg, org.r-project.x86_64.tcltk.x11 and org.r-project.x86_64.texinfo. You can use pkgutil –forget if you want the Apple Installer to forget about the package without deleting its files (useful for the R framework when installing multiple R versions in parallel), or after you have deleted the files.

Uninstalling the Tcl/Tk or Texinfo components (which are installed under /usr/local) is not as simple. You can list the files they installed in a Terminal by

pkgutil –files org.r-project.x86_64.tcltk.x11
pkgutil –files org.r-project.x86_64.texinfo

These are paths relative to /, the root of the file system.

上記を自動翻訳するとこんな感じです。

4.2 macOSでのアンインストール

R for macOSは、GUI(R.APP)とRフレームワークの2つの部分で構成されています。アンインストールは、それらのフォルダを削除するだけで簡単です(例えば、ゴミ箱にドラッグするなど)。標準的なインストールでは、GUIは/Applications/R.appフォルダに、Rフレームワークは/Library/Frameworks/R.frameworkフォルダにインストールされます。 / usr / local / binにあるRおよびRscriptへのリンクも削除する必要があります。
ターミナルを使用してRをより完全に削除したい場合は、次のコマンドを実行します。
sudo rm -rf /Library/Frameworks/R.framework /Applications/R.app \
/usr/local/bin/R /usr/local/bin/Rscript
インストールは、最大4つのAppleパッケージで構成されています:18 org.r-project.R.el-capitan.fw.pkg、org.r-project.R.el-capitan.GUI.pkg、org.r-project.x86_64 .tcltk.x11とorg.r-project.x86_64.texinfo。 ファイルを削除せずにパッケージをAppleインストーラに忘れさせるには(複数のRバージョンを並列にインストールするときにRフレームワークに役立つ)pkgutil –forgetを使用することができます。Tcl / TkやTexinfoコンポーネント(/ usr / localにインストールされている)をアンインストールするのは簡単ではありません。ターミナルにインストールしたファイルを

pkgutil –files org.r-project.x86_64.tcltk.x11
pkgutil –files org.r-project.x86_64.texinfo

これらは、/に関連するパスで、ファイルシステムのルートです

ん? 次のコマンドだけでTcl/Tk以外を完全に削除できるということかな?もしそうなら、これのコピペでいいじゃん。
> sudo rm -rf /Library/Frameworks/R.framework /Applications/R.app \
/usr/local/bin/R /usr/local/bin/Rscript
この rm -rf のコマンドはそのディレクトリ以下のファイルが警告・確認なしに削除されるため、必要なファイルが削除される危険があります。どうぞ、使用される際にはご注意ください。
なお、次のコマンドをLinux上で実行した際のスクリーンを記録した動画が公開されています。視聴者の反応からも、そのコマンドの危険さが感じられると思います。
rm -rf /

LearnBayes R パッケージによる臨床試験と市販後の副作用データの解析

RのLearnBayes パッケージで臨床試験と市販後の副作用データを解析してみました

ここに提示されるデータは、治験・市販後の実データを元にしていますが、実際のデータではありません。LearnBayesというRのパッケージはベイズ推計を行えるような関数を提供しています。説明書を読みながら、このパッケージの一部機能を使ってみてみました。


1.海外試験の結果が事前にあって、国内の副作用情報が得られた

海外の治験のデータがすでにあり、その後国内で試験を行って副作用データが新たに得られた。と言う場合を想定した解析をしてみた。

海外臨床試験では、被験者135例中50例に副作用が報告された。普通の頻度の解析によると副作用が報告される割合は、本剤使用者の37.0% (95%CI; 28.9 – 45.8)となる。

> binom.test(50, 135)

Exact binomial test

data: 50 and 135
number of successes = 50, number of trials = 135, p-value = 0.00328
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.2888945 0.4576672
sample estimates:
probability of success
0.3703704

この計算結果をもとに、β分布を推定する。95%信頼区間の上限(p=0.975)を与えるxの値は0.4576672であり、最尤値(p=0.5)を与えるxの値は0.3703704。この2点を与えれば事前確率のβ分布のパラメータが推定できる。

quantile2 <- list(p=.975, x=.458)
quantile1 <- list(p=0.5, x=.370)
beta_parm <- beta.select(quantile1, quantile2)
a <- beta_parm[1]
b <- beta_parm[2]

国内臨床試験では海外より被験者が少なく、23例が安全性解析対象症例となった。23例中11例で副作用が報告された。この結果と事前推定されているβ分布のパラメータを元に前後の分布を図示する。

# 事後分布
curve(dbeta(x,a+s,b+f), from=0, to=1, xlab=”p”,ylab=”Density”,lty=1,lwd=4)
# 実験結果から出てくるイベント発生頻度の尤度
curve(dbeta(x,s+1,f+1),add=TRUE,lty=2,lwd=4)
# 事前分布
curve(dbeta(x,a,b),add=TRUE,lty=3,lwd=4)
legend(.7,9,c(“Prior”,”Likelihood”,”Posterior”), lty=c(3,2,1),lwd=c(3,3,3))

国内試験のデータの分布(荒い破線)のデータが得られた後も、事後の分布(実線)は事前分布(細かい破線)と大きく変わらないようだ。

# 事後の分布による推定

> # 90%信頼区間
> qbeta(c(0.05, 0.95), a+s, b+f)
[1] 0.3222927 0.4550214
> # 点推定
> qbeta(0.5, a+s, b+f)
[1] 0.3872558

事後の確率は0.39 (90%CI; 0.32 – 0.46)で、国内のデータ0.48 (90%CI; 0.30 – 0.66)とずいぶんずれているように見える。

# 国内データのみによる推定

> binom.test(11, 23, conf.level=.9)

Exact binomial test

data: 11 and 23
number of successes = 11, number of trials = 23, p-value = 1
alternative hypothesis: true probability of success is not equal to 0.5
90 percent confidence interval:
0.2960934 0.6648524
sample estimates:
probability of success
0.4782609

2. 治験の結果をもとに製造販売が承認され、市販後に医薬品が使われて安全性データが集積した

上記の臨床試験の結果を元に承認申請がなされ、製造販売の承認を得た。市販後はこの医薬品は治験に参加した被験者の数よりはるかに多い患者さんに使用されて多くの安全性データを得た。治験の結果を事前分布、市販後のデータを加えて事後の分布を集計してみる。

上記の分布を事前分布の確率密度関数のパラメータ推計に用いる

# 90%信頼区間の上限 0.95のx値は0.455
quantile4 <- list(p=.95, x=.455)
quantile4$x <- qbeta(0.95, a+s, b+f)
# 事前確率(点推定; p=0.5)は0.387であった
quantile3 <- list(p=0.5, x=.387)
quantile3$x <- qbeta(0.5, a+s, b+f)
beta_parm <- beta.select(quantile3, quantile4)
a <- beta_parm[1]
b <- beta_parm[2]

国内の市販後の調査に登録され安全性解析対象とされた患者は2072例であり、そのうち964例に副作用が報告された。この結果と事前推定されているβ分布のパラメータを元に前後の分布を図示する。

s <- 964
f <- 2072-s
# 分布を見てみよう
# 事後分布
curve(dbeta(x,a+s,b+f), from=0, to=1, xlab=”p”,ylab=”Density”,lty=1,lwd=4)
# 実験結果から出てくるイベント発生頻度の尤度
curve(dbeta(x,s+1,f+1),add=TRUE,lty=2,lwd=4)
# 事前分布
curve(dbeta(x,a,b),add=TRUE,lty=3,lwd=4)
legend(.7,35,c(“Prior”,”Likelihood”,”Posterior”), lty=c(3,2,1),lwd=c(3,3,3))

 

市販後の調査のデータの分布(荒い破線)のデータが得られた後、事後の分布(実線)は事前分布(細かい破線)から大きく右へシフトした。

# 90%信頼区間
> qbeta(c(0.05, 0.95), a+s, b+f)
[1] 0.4427958 0.4776115
> # 点推定
> qbeta(0.5, a+s, b+f)
[1] 0.4601712

図から受ける印象の通り、事後の確率は0.46 (90%CI; 0.44 – 0.48)は、事前のデータ0.39 (90%CI; 0.32 – 0.46)とずれているように見える。実は、事後の確率は、国内治験のデータ0.48 (90%CI; 0.30 – 0.66)を精密にしたように見える。

3. まとめ

ここまでやってみた感想

  1. 事前の試験・調査と事後とで大きく規模(被験者数)が違う場合、被験者数の大きな方の試験・調査の結果が事後の分布に反映しているだけではないか
  2. それなら単に規模の大きな試験・調査の結果の分布を集計しても同じではないか、あるいは、単純に足したような併合解析でも結果は似たようなものになるのではないか
  3. この様な例では、国内外で副作用の報告頻度が異なるという仮説を導きたい
  4. 結局事後分布って何を示しているものなんだろうか?
  5. なぜ、β分布を想定するのだろうか?

疑問がより具体化したので良しとしよう。

 

Childhood leukemia and residential magnetic fields

高圧電線と小児白血病

1988年に公開されたSavitzらの論文(症例対象研究)1のデータを触っています。

Caseは白血病症例、Controlはそれ以外の住民。そして、Exposureは閾値を3milligauss (mG)とした場合の居住地での磁気フィールドへの曝露の有無。データは次のようになります。

これを基にβ分布で症例とコントロールの分布の1万例のシュミレーションを行うプログラム。

そして、結果です。

そのヒストグラムを書くと次のように分布します。

ま、これとは別に普通にFisherのテストをすると次のような結果になっていて、オッズ比は3.48、95%CI 0.51 – 18.89 と幅広い分布をすることが計算できます。

1.
Savitz D, Wachtel H, Barnes F, John E, Tvrdik J. Case-control study of childhood cancer and exposure to 60-Hz magnetic fields. Am J Epidemiol. 1988;128(1):21-38. [PubMed]

‘Rcpp’ という名前のパッケージはありません

RでインストールしたはずのRcppパッケージが認識されないという悩ましい状況

> library(Rcpp)
library(Rcpp) でエラー: ‘Rcpp’ という名前のパッケージはあり

インストールを3回試してみて、ダメでした。ダメだった時のログですが、有用なエラーメッセージはありません。コンソールを見ている限り何がダメなのかわかりません。

install.packages(‘Rcpp’)
パッケージを ‘C:/Users/Oshima/Documents/R/win-library/3.3’ 中にインストールします
(‘lib’ が指定されていないため)
— このセッションで使うために、CRAN のミラーサイトを選んでください —
URL ‘https://cloud.r-project.org/bin/windows/contrib/3.3/Rcpp_0.12.15.zip’ を試しています
Content type ‘application/zip’ length 4357370 bytes (4.2 MB)
downloaded 4.2 MB

パッケージ ‘Rcpp’ は無事に展開され、MD5 サムもチェックされました

ダウンロードされたパッケージは、以下にあります
C:\______\Local\Temp\RtmpUfqmSC\downloaded_packages

しばらくやった後、これは、先人も経験があるはずと調べたところ、それらしい記事に当たりました。リンク先の柴田啓文先生、ありがとうございます。

セキュリティソフトを10分間無効にしてインストールすると、わたしも成功しました。

 

Armenian Dance part1

アルメニアンダンス パート1

楽曲の元になったと思われるKomitas (Gomidas)の曲集を見つけたので、スクリーンショットをまとめてみました。

(スクリーンショットは表示されたりされなかったりがあるので、曲集のダウンロードリンクも作成しました。)


VOICES OF ARMENIA (Geghard Ensemble, Pahayan)

Translate »