JHS Pedals - Little Black Amp Box

これまでに掲載している記事とはかなり毛色が異なる，個人的な趣味に関する記事です。
今回はタイトルにあるように，JHS PedalsというメーカーのLittle Black Amp Box（LBAB）というエフェクターについての書きます。

JHS Pedals パッシブアッテネーター風ペダル Little Black Amp Box (国内正規品)

まずは，代理店であるLEP INTERNATIONALの製品ページより紹介文を引用いたします。*1
Little Black Amp Box - LEP INTERNATIONAL

Little Black Amp Boxは、エフェクトループを持つアンプのために設計されたパッシブのアッテネーター風ペダルです。

アンプにSend/Returnのエフェクトループがあるなら、このLittle Black Amp Boxをエフェクトループに設置するだけで、プリアンプで飽和させた音をそのまま音量だけを下げ、ベッドルームレベルでチューブの歪みを楽しむことができます！

紹介文にもあるように，アンプのセンドリターンに入れてマスターボリュームのように使える，電源不要のパッシブなエフェクターです。
新宿にある某ビシバシ楽器店で在庫処分セールをやっており，税込6000円で購入しました。*2

購入した理由は，現在E.W.SからSubtle Volume Control（SVC）というペダルを所有しており，エフェクターボード最後段にてマスターボリューム（ボリュームペダル）として使用しているのですが，これに似ているなぁミニサイズだし置き換えられないかなぁと考えたためです。

E.W.S. Subtle Volume Control

ということで，LBABについてのレビューを適宜SVCと比較しながら書いてみようかと思います。
このペダルについてはかなり地味な働きをするからか，日本語の情報はほとんどないようですし，一定の価値はあるでしょう。

まずは外観です。

見ての通りミニサイズで，ボリューム一発という非常にシンプルなペダルですね。

中身は次の写真のようになっており，外観の印象を裏切らない圧倒的シンプルさです。

スイッチクラフト製のように見えるモノラルジャックが2つに，Alpha electronics社製100KΩ・Aカーブのポットが1つ，それをつなぐワイヤーが数本といった感じでしょうか。
材料費を計算すると，6000円（または45$）という金額に対して思うところができてしまうので考えないことにします。

中身はどう違うのか

中身の違いについて確認するにあたって，今回おおいに参考にしたSVCの中身について言及なさっている日本語情報のブログをご紹介いたします。
（勝手に紹介しておりますので，問題がありましたらご一報ください）

My Things:E.W.S. Subtle Volume Control Wiring Inside (EWS) - livedoor Blog（ブログ）
以後このブログのことは便宜的にブログタイトルであるMy Thingsさんとお呼びします。

また，SVCの中身は次の写真のようになっています。
これについては，はじめブログに掲載する予定はなかったので手が映り込んでしまっている気の抜けた写真になっていますがご容赦ください。

私は電子部品には疎いですが，My Thingsさんの記事を拝見すると，SVCとLBABの違いは大きく次の2つだと言えそうです。

• 配線の違い

• 可変抵抗器（ポット）の違い

ひとつずつ確認してみます。

配線の違い

どちらもボリュームポットが1つにジャックが2つというシンプルな構成ですが，内部での配線のしかたに細かな違いがあるようです。
うーん，この違いが何をもたらすのか，うーん。

適当ではありますが，2つの中身を並べて画像にしてみました。

ジャック部分にあるTとSはそれぞれTipとSleeveを意味しています。
Tipはホット，つまり信号が通る部分となっており，Sleeveはアースに落ちる部分となっています。
また，可変抵抗の部分にある数字はそれぞれ端子の番号となっています。

これをふまえて，素人目で気づいたことを箇条書きしてみます。
この部分については間違いがあるかもしれないことをおことわりしておきます。

• LBABでは，input側のアースがポットの1番につながっておりミュートされるときはこの部分が機能して入ってくる信号がミュートされている？対してSVCで1番端子につながっているのはoutput側のtip端子であるため，出ていく信号がミュートされるかたちになる？

• ボリューム最大となる3番端子について，SVCではジャックから3番に行った信号が2番からアウトプットへ行くかたちをとり，LBABではジャックから2番を通って3番，アウトプットと流れる？

なんだか，素人目にはSVCの方が理にかなった配線をしているような気がしてきました……。
一般には2番からoutputへ通るらしいので，LBABで2番にinputのtipから直で信号が通っているように思われるのはよくわかりません。
この配線の仕方が，「アッテネーター風」を謳うこのペダルのキモなのでしょうか，JHS Pedalsが適当な仕事をするイメージはあまりないので謎は深まるばかりです。

ポットの違い

My Thingsさんでは，SVCのポットは25KΩ2W・Aカーブの特殊なポットであるとされていました。
見た目にも大きいポットで，ワット数も大きいようです。
一方LBABでは，100KΩ（おそらく0.06W？）・Aカーブのポットであり，抵抗値とワット数は大きく違うようです。
抵抗値とワット数の違いが計算上どのような作用をするのかはわかりませんが，（流れる電流の量が変わるのかなーくらいの印象です）次のセクションで述べる実用上の違いに影響するひとつの要因なのでしょう。

実用上どう違うのか

調べものをしながら中身についての違いについて書いてきましたが，先にも述べたように私自身そのあたりには疎いので地に足のついていないレビューになってしまった感が否めません。
ということで，プレイヤー目線で，実際につなげて使用し，違いについて検討してみました。
音源などはありません，すみません。

手持ちのアンプシミュレーターにあるFX LOOPにつないでみました。
アンプ―キャビネット―ディレイ―FX LOOP―フィルター，という流れでFX LOOPにSVC，LBABを単体で挟み込み，ギターを鳴らしたあとにSVCまたはLBABでボリュームを操作してみる，という手続きです。

同じAカーブのはずですが，聴覚上SVCの方が使い勝手のよいマスターボリュームとして機能しており，音量の調整はボリュームペダルを踏んでいるときの音量変化のようにスムーズです。
LBABでは音量が最大（5時くらい）から9時くらいまではあまり変化せず，8時くらいから一気に小さくなるような印象でした。
また，LBABをツマミ最大でつないでいると，アンプシミュレーターでのループ内のインプットレベルが上がり，クリップしてしまうという不思議な現象に直面しました。

結論として，エフェクターボードに組み込むボリュームペダルとして使おうと思うとSVCの方が圧倒的に使いやすいように思います。
LBABは私には中身も働き方も難しいエフェクターのようです。
が，ポットと配線をSVC風味に改造してしまえば，穴あけをする手間なく使いやすいボリュームペダルにメルヘンチェンジできるのではないのかしら，とぼんやり考えています。

アンプのセンドリターンに挟むという正しい使い方とはやや異なりますが，アンプシミュレーターのセンドリターンに入れてマスターボリュームとして使うというほぼ誤りのない使い方をしてみての感想でした。
購入失敗だったかな，という考えが一瞬頭をよぎりましたが，この記事を書くにあたって随分勉強になったのでよしとしましょう。
まだメルヘンチェンジコースも残っていますし。
今度時間のあるときに，スタジオに行って私の予定しているボード最後段での使用ではどうなのかを比較してみようと思います。
最後に，今回参考にしたサイトを備忘もかねて参考文献としてリンクしておきます。

追記（162014）

こちらのページに掲載されているLBABの中身と今回購入したものの中身とが，微妙に違う気がするんですけれども……。
2014年10月絶叫日記

ここではinputとoutputのSleeve同士がSVC同様が繋がっていますが，私の手元にあるものではそこは繋がっていないようで，何が何やら。
某SPICEのようにスルーボックスを作る場合にはSleeve同士は結線するようなのですが，これは大丈夫なのかしら。

参考文献

（再掲）SVCの中身を分析されていたブログ
My Things:E.W.S. Subtle Volume Control Wiring Inside (EWS) - livedoor Blog（ブログ）

ボリュームポットこと可変抵抗器の仕組みについて，アニメーションでわかりやすく説明されていた記事
『ポット（ポテンショメーター）』 もっと自由に大きな音で！OZIMAS Guitar Atelierのブログ

同じく可変抵抗について，回路図風味でわかりやすかったもの
可変抵抗を使ってみよう : こた電 - こたつぁーの電子工作

抵抗やワット数などについて
抵抗

*1:この紹介文はJHS Pedalsのメーカーサイトにある英語での紹介文を翻訳したものであるようです，1行省略してありますが。

*2:本国では45$らしいです。

van der Linden & Lewis（2015）．カンニングを事後オッズで

van der Linden, W.J., & Lewis, C（2015）．Bayesian checks on cheating tests. Psychometrika, 80, 689-706.

カンニングを統計的に発見する，という方法論があることを知らなかったのでなかなか新鮮でした。

背景

アメリカではNo Child Behind Left政策が打ち出されて以後，教員の評価においても生徒の成績が加味されるようになったそうで，生徒の成績を上げられない教員は最悪クビなんだとか。
そのため，日本でカンニングといったときにすぐ浮かんでくる，生徒間での見せ合いや盗み見といったものだけでなく，教員によるカンニングも見られるそうです。
教員のカンニングというと，たとえば回収した答案の改ざん，試験中に「ここ，もっとよく考えよう」といったアドバイス（これは日本でも目にする？），事前に入手したテスト情報を生徒に教えてしまう，などが挙げられます。
こうしたカンニングに対処するためにテスト会社ではいろいろとがんばってきたようで，今回の論文で取り上げる方法もその一環ということです。

カンニングの種類と対策

4種のカンニング

1. テスト受験者が他人の回答を写す（copying）
2. テスト完了後に行われる不正な回答の取り消し修正（fraudulent erasures）
3. あとから受ける受験者のために，テスト中に項目を覚えようとする
4. 未公開のテスト項目に関する事前知識

大まかな対策

個々のカンニングへの対処法はいろいろとあるようですが，詳しい情報はこの論文でレビューされていますので，どうぞ原典へ。
ここではほとんどどの方法でも共通している点を大まかな対策としてまとめます。
簡単にいえば，カンニングしたと思しき人（以後コピー者と呼びます）と，カンニングされただろう人（以後，ソースと呼びます）との間で一致した解答について，その数が偶然かそうでないかという統計的検定を行うという対策を取ります。

しかしながら，ソースの項目反応をどのように取り扱うかは昔からずっと課題として残ってきたようですね。
現在出されている3つの取り扱い方は，次のようなものです。

• ソースの不正解の項目のみに着目する

• ソースによる反応はすべて所与のものとする

• すべてランダムなものとする

このうち2番目，3番目を見るに，要はコピー者の反応におけるソースの反応を条件づけるか条件づけないかという問題になり，検定においても条件付けた確率と条件付けない確率のどちらを用いるかで恣意性が残ってしまうことが課題なのでしょう。

これをそもそも回避する方法として，この論文ではベイズアプローチによる事後オッズを用いることが提案されています。
ベイズでは選ぶまでもなくすべての確率が条件付きで表されるため，その点については恣意性が残らないというわけです。

回答コピーにおける事後オッズ

ここでは回答コピーに注目して，運用しやすいよう項目反応モデルとの接続を考えながら，カンニングをしてしまったとう事後オッズを考えますが，ここでは途中に必要となる仮定なり展開なりがあまり省略できず，煩雑になってしまうのでその記述的な性質を述べるの留めます。

提案されている事後オッズについては次の性質があります。

• コピー者とソースの2人の反応について，一致している部分のみを対象とすればよい（したがって一致していない部分を用いて，コピー者の，チートを含まない能力を推定できうる。

• ソースの能力には依存しないで推定ができる。

• 一致した項目のサブセットとして実際にコピーが行われた項目があるはずであり，本来実際コピーをしたサブセットは知ることができない。しかし，この方法では一致した項目における全組み合わせを検討して確率を足し上げるアプローチをとるので，実際にコピーした項目がわからずとも事後オッズを計算できる。

• しかし，この事後オッズは一致した項目の数の関数なので注意が必要

• 一致した項目の全組み合わせを考えるので，テスト項目が多くなると爆発的に計算時間がかかる。したがって，全テスト受験者についてルーティン的にチェックするという方法は推奨しない。カンニングの疑いがある生徒について，詳細にチェックするときの方法と言える。

面白いなーと思ったのが，実際にコピーした項目がわからなくても事後オッズを計算できるということでしょうか。
まだまだ知らない，統計と世界とのかかわり方があるんだなぁとしみじみ思いました。

橘田いずみさんが紹介する「絶対外さない！お取り寄せ餃子5選」

こんにち餃子，ということで2016年1月22日放送「白熱ライブ　ビビット」より。
絶対外さない！お取り寄せ餃子5選として，餃子評論家の橘田いずみさんが餃子を紹介されていました。
今回のポイントは，“おうちで焼いてもお店の味をそのまま再現できる”ことだそうです。
わたくし，個人的に橘田いずみさんを応援しています。

橘田いずみのザ・餃子

1．餃子激戦区　宇都宮のお取り寄せ餃子

曰く，見た目がびっくりで若い人におすすめ，彩花というお店ののキング餃子（5個入り450円）。
普通の餃子の約3倍らしく，一時期健啖家を自負していた自分としては非常に胃袋を刺激されました。
ところで橘田さん，宇都宮のお取り寄せ餃子のすべてを食べ比べたとのことで，さすがの貫禄です。

2．餃子激戦区「浜松」　10個＿＿＿円のお取り寄せ餃子

浜松餃子栄福というお店の餃子（10個300円）。
材料を手で粗めに切ることでしっかりとした食感が残る餡となり，隠し味の山芋も合間ってとてもおいしいそうです。
この値段と数はいっぱい食べたいおひとりさまからわいわい食べたい複数人まで，幅広く楽しめて素敵です。

3．世田谷区梅丘　羽根つき餃子

自宅でも簡単につくれる羽根つき餃子として紹介されていたのは，世田谷区梅丘にある羽根木餃子というお店の羽根木餃子（5個入り600円）で，箱（小130円，大150円）に入れてくれるのでプレゼント用にもよいとのこと。
このお店ではトマト餃子，パクチー餃子といったものも販売しており，お取り寄せは半年待ちだそうです（！）。

4．鹿児島　1個で3種類の食感が楽しめるご当地餃子

丸っこくてお菓子のような見た目の，かごしま黒豚餅餃子（14個入り1080円）。
火が通った状態で送られてくるので，オーブントースターで焦げ目がつくまで焼くだけ，つまり誰でも焼ける餃子なんだとか。
スープに入れて水餃子，揚げて揚げ餅と，いろんな楽しみ方が紹介されていました。

5．香川県　2年待ちのクロワッサン餃子

たれ屋というお店のクロワッサン餃子（生餃子）（50個入り2650円）。
14年かけて完成した，すべてに徹底的にこだわりぬいた餃子らしいです。
2年待ちとはいえ，食べてみたい……。

番外編，手軽に作れる絶品たれ

番組内で餃子の裏技として紹介されていた，たれのつくり方についても備忘をかねて以下にまとめます。

• お酢に胡椒のみのシンプルなたれ

• 砂糖・はちみつ・酢・白味噌を混ぜるだけのお手軽な味噌だれ

2016年の餃子は「もちもち」らしい

完全なる私得でしかない備忘録です。

スーパーJチャンネルの特集，達人おススメ！「進化する餃子」④より。
2016年の餃子は「カリカリ」から「もちもち」へ，とのことです。

餃子の達人として，餃子評論家こと橘田いずみさんが「もちもち」をコンセプトに紹介されていたお店の情報をまとめておきます。
橘田いずみさんは個人的に応援しています。

橘田いずみのザ・餃子

渋谷区恵比寿の花林糖餃子

醤油（540円），味噌（650円），パクチー（650円）と3種のスープにもちもちの皮をした餃子の入った「スープ餃子」，橘田さんは味噌がお好きだそうで。
もっちりとした皮だと，うどんのようにしっかりと食べ応えがありスープの味もしみる，そんな一品になりそうでおいしそうですね。
tabelog.com

新宿区飯田橋のPAIRON

台湾の水餃子をさらに進化させて，中にいろいろなものを入れた色鮮やかな餃子が目を惹くお店でした。
番組で橘田さんとレポーターの方が召し上がっていたのは白龍餃子，青龍餃子，黒龍餃子，レッドドラゴンという4種類で各390円。
個人的にはシナモンの入った白龍餃子が気になります。
tabelog.com

世田谷区成城学園前のギョーザブラザーズ

ミニ中華まんのようにも見える（放送では「白玉？」なんて言われていました）かわいい餃子，もちもちまるチーズ／キムチ餃子（421円）が紹介されていました。
白いチーズ入りの方の餃子は求肥のようだと，言っているのを聞いてうまい表現だなと感心した次第です。
定番の焼き餃子もカリカリ・ジューシーとのことで食べてみたいです。
tabelog.com

豊島区池袋の中国家庭料理楊2号店

四川料理のお店とのことで辛いものが食べられない私には縁がないかもしれない，と思っていましたが，水餃子が絶品らしく非常に気になります！
紹介されていたものは真っ赤なスープが見ているだけで痛い紅油水餃子（920円）でした。
食べログのレビューで知ったのですが，孤独のグルメにも登場したお店なんですね。
tabelog.com

もちもちとした皮の餃子，とても気になりますね。
ごはんと一緒に餃子を食べるのと同じよう感覚になる部分もあるのかもしれません。
餃子と言えばカリカリ派でしたが，年も明けたことですし新たな境地を開拓してみようかしら。

この記事を書くために深夜3時に特集を見直していてしみじみ思いましたが，夜中に見るものじゃないですね，本当に。

mvrnorm関数で「'Sigma' is not positive definite」というエラーが出たときの対処法

シミュレーションなどのために，多変量正規分布からのサンプリングを行う関数として{MASS}パッケージのmvrnorm関数があります。
この関数は次のようなプログラム例を通してサンプリングを行います。

library(MASS)

#3変量正規分布の例
#平均値として0が3つのベクトルを作る
mux=numeric(3)

#3*3でそれぞれ0.5の相関係数をとる相関行列
sigma<-matrix(0.5, 3, 3)
diag(sigma)<-1

#サンプルサイズ=100，3変数で多変量正規分布からのサンプリング
data<-mvrnorm(n=100, mu=mux, Sigma=sigma, empirical=FALSE)

empiricalというのは与えた相関行列を寸分たがわず再現するかしないか，というロジカルな引数で，TRUEだとcor(data)として相関係数を求めたときに与えた相関行列と完全に一致するため母相関行列のように使えるかと思います。
今回はFALSEとしているので，標本から推定した相関係数のように使うことができます。

さて，ここから本題ということで次のような状況を考えてみましょう。
x1-x10までの10個の変数があるときに，x1-x9までは相互に0.1で相関しており，x10だけは他の変数と0.5で相関しているという状況です。
これを前提として，10変量正規分布からのサンプリングを行うプログラムを次のように書きます。

library(MASS)

mux<-numeric(10)
sigma<-matrix(0.1, 10, 10)
sigma[10,]=0.5
sigma[,10]=0.5
diag(sigma)<-1

data<-mvrnorm(n=100, mu=mux, Sigma=sigma, empirical=FALSE)

そうすると，
Error in mvrnorm(n = 100, mu = mux, Sigma = sigma, empirical = FALSE) :
'Sigma' is not positive definite
というエラーメッセージが返されて，データが生成されません。

この原因の確認と対処法をいかに述べようと思います。

エラーの原因の特定

メッセージを読めばわかるので言うまでもないのですけれど，相関行列の設定が不適切ということになります。
ただ，それで終わっては芸がないので，もう少し細かく見ていきます。

まず，mvrnorm関数の中身は次のようになっていることを確認します。
見やすさを重視して，スクリーンショットにしてみました。

この中で今回関心があるのは，
stop("'Sigma' is not positive definite")
という条件式までの部分になります。

かつ相関行列の指定までですので，次の部分が問題となっていると考えてよいでしょう。

eS <- eigen(Sigma, symmetric = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")

これを，日本語に翻訳してコメントをつけてみます。

#与えた相関行列／共分散行列の固有値を求める
eS <- eigen(Sigma, symmetric = TRUE)

#固有値の計算結果（固有値，固有ベクトルが含まれる）の中から固有値だけを取り出す
ev <- eS$values

#もしすべての固有値が条件式を満たさなかった場合にはエラーメッセージを返す
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")

これらの内容から，今回のエラーは
if (!all(ev >= -tol * abs(ev[1L])))
という固有値に関する条件式を満たせなかったことが原因であることがわかります。

この条件式についてもうすこし細かく見てみます。
tolというのは上のスクリーンショットを見ていただければわかりますが，1e-06*1が割り当てられています。
abs関数は絶対値を返すものなので，abs(ev[1L])は1つ目の固有値の絶対値という意味になります。
また，!は否定演算子です。
したがって，この条件式は(-0.000001)*一つ目の固有値の絶対値を，得られた10個の固有値のすべてが越えなければFALSEを返すif条件式ということになります。

ためしに今回の状況における固有値の計算結果を見てみましょう。

> eS <- eigen(sigma, symmetric = TRUE)
> ev <- eS$values
> ev
[1] 2.9524175 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 0.9000000 -0.1524175

> tol=1e-06
> -tol * abs(ev[1L])
[1] -2.952417e-06

> ev[10]>=-tol * abs(ev[1L])
[1] FALSE

ひとつでもFALSEを返す場合にはエラーを返すようになっているので，かくして'Sigma' is not positive definiteというメッセージが表示されるわけです。

対策（力技）

さて，対策なのですが，要は与えた相関行列の組み合わせによっては固有値を計算する際に不適切なものが生じる恐れがある，ということなので，与える相関行列を最初から修正してしまえばいいじゃないかということになります。
そこで使うのが，{sfsmisc}パッケージのnearcor関数です。
nearcor {sfsmisc} | inside-R | A Community Site for R

これは不適切な相関行列を，cor関数でとりうる最も近い値にスムージングする関数とのことなので，素直に先ほどの相関行列をnearcor()の中に入れてからmvrnorm関数の相関行列として指定してあげれば，とりあえずエラーは消えるという力技的な解決法となっています。

お粗末様でした。

*1:=0.000001

加工肉は発がんリスクを上げるのか？

WHOが「加工肉には発がん性がある」といった旨の声明を出したことが世間で話題となっているようですね。www.bbc.com

日本人への影響は小さいだの，続報がどんどん出てきているようですが，個人的にこの声明自体どこまで本当なのかしらと疑問に思ったので元ネタを探してみようと思いました。
おそらくそうした声明を出すに至った研究結果があったのだろうと考えてちらっと探してみたところ，いくつかヒットしました。

その中でも，今回はCross, Leitzmann, Gail, Hollenbeck, Schatzkin, & Sinha (2007) という論文*1をつまみぐいしてみようかと思います。journals.plos.org

タイトルを眺めるに，赤身肉・加工肉の摂取と発がんリスクの関係に関する前向き調査といったところでしょうか。
簡単にですが概要を訳してみました。

論文の概要

Abstract

Background

赤身肉と加工肉はいくつかの解剖学的観点から発がんと関連付けられてきたが，しかしながら悪性の範囲と肉の摂取との関係についての前向き調査はなされてこなかった。
本論文では，赤身肉ないしは加工肉の摂取ががんのリスクを高めるかどうか，様々な観点から調査した。

Methods and Findings

• 1995-1996をベースラインとした，50～71歳の500,000人に対するコホート調査

• 食肉の摂取はベースラインで質問紙調査を行い推定

• Cox比例ハザード回帰で，赤身肉と加工肉を摂取している人の5分位におけるハザード比と95%信頼区間を推定

• 8.2年のフォローアップで，53,396件のがんが確認された

• 赤身肉摂取の最も高い20%と最も低い20%との間で，統計的に有意なリスク増加（20%～60%のレンジ）が食道・結直腸・肝臓・肺について見られた

• 加えて，加工肉を摂取する人の上位20%は，結直腸がんのリスクが20%，肺がんのリスクが16%高まることが示された

Conclusions

赤身肉と加工肉の両方の摂取が，結直腸と肺のがんとの間に正の関連が見られた。
加えて，赤身肉の摂取は食道と肝臓のがんのリスク増加に関連が見られた。

今回のお話し

事前のお勉強

まず，そもそも私の不勉強さの問題ではあるのですが，Cox比例ハザードモデルについてはあまり知りませんでした。
医療系で使うらしい手法だなーとか，そのくらいの印象です。
なので，これについて簡単に，手元にあった森田（2014）の『実証分析入門　データから「因果関係」を読み解く作法』の第15章で勉強してみようと思います。

実証分析入門 データから「因果関係」を読み解く作法

サバイバル分析

Cox比例ハザードモデルというのは，サバイバル分析の文脈で用いられるモデルだそうです。
サバイバル分析（Survival Analysis）*2について，これはある個体が死亡してしまった場合そこから先は観測されなくなってサンプルから落ちてしまう，という特徴を持ったデータを扱うときの手法の枠組みです。
サバイバル分析に関する詳細な説明はここでは割愛しますが，重要なのはこの枠組みではハザード率（Hazard Ratio; HR）というものを分析対象とすることだと言えます。
ハザード率は，森田（2014）では次のように定式化されています。

たとえば，a歳のある病気に罹患した患者を集めて30年間の継続調査を行ったとしましょう。
この病気が5年以内死亡率が非常に高いものだったとしたとき，分子や分母にあたるある時点での死亡者数や生存者数（確率）だけで考えると，25年目まで生きている人自体がそもそも非常に少ないためその時点での死亡率などを考えることに実質的な意味がなくなってしまいます。
そこで，ハザード比というものに変換する，要するに割合に直すことで刑事変化の問題を回避して本当の死亡率に焦点化するのがサバイバル分析で扱う対称の特徴と言えるかもしれません。

Cox比例ハザードモデル（Cox Proportion Hazard model）って？

さて，本題のCox比例ハザードモデルですが，このモデルは次のように定式化されます。

ここではベースラインハザードという，全個体で共通のものです。
また，右辺の指数の中には時間を表す添え字tが存在していないことが見てわかると思います。
つまり，ベースラインハザードという共通のものと，時点で不変である様々な変数=要因を含んだ部分の積としてハザード率が表現されているということになります。
比例ハザードという名前の意味ですが，今回のお勉強本の記述がわかりやすかったので（一部変更はありますが）そのまま引用します。

各個体のハザード率が共通要員であるベースラインの何倍か（倍率は説明変数で決まる），というベースラインハザードに比例する形で決まるので，比例ハザードと呼ばれる

このモデルは時間部分，すなわちベースラインハザードは推定せず，どの要因の影響が大きいかという観点に絞った分析を行うことが特徴的であると言えます。
ベースライン部分の情報として，存続期間を順序情報のみ用いたうえでパラメトリックな解析を行うことから，セミパラメトリックな手法であると考えることもできるようですね。
ハザード率の経時変化を考えずに要因の影響という関心のある部分のみをモデリングするという便利な手法ではありますが，欠点として次の3つがあげられていました。

1. 一部使わない情報がありもったいない
2. 比例ハザードという仮定自体が強いものである
3. 複数のリスク要因を独立とする仮定も強いものである

まとめると，ハザード率が全員共通のベースラインから，説明変数によって何倍になるかをセミパラメトリックな形でモデリングしたサバイバル分析の一手法ということですかね。

論文の図表を簡単に眺めてみる

Cox比例ハザードモデルをよく知らなかったのでいい勉強の機会だと思って記事を書き始めた，という事情がありここまででほとんど目的は達成してしまったのですが，せっかくなので少し今回取り上げた論文の結果を眺めてみようと思います。

たとえば，加工肉に関する分析を行ったTable 3ですが，基準となるQ1，すなわち加工肉を食べる頻度でソートした時の下位20%の人に対して，上位20%の人であるQ5の結直腸がんに関するハザード比は1.20でした。
表の脚注によると，説明変数としては年齢，性別，教育歴，結婚歴，がん経験，人種，BMI，喫煙，運動習慣，総摂取エネルギー量，飲酒，果物・野菜消費量を投入しているとのことです。

この論文ではがん発症をイベントとしているようですから，それまでにがんにかかっていなかった人がその時点でがんになる確率がハザード比としたとき，Q1の人に対してQ5の人は1.2倍（20%）多く発症リスクがあるということになるのでしょう。
ただ，信頼区間を見てみると1.07-1.30となっていますから，ほぼ変わらないかもしれないし，大きくて30%ほどなのかもしれません。
要旨で触れられていた肺がんについても，Q1に対してQ5のハザード比が1.16とのことですから16%とのことでしたが，信頼区間を見てみると1.06-1.26ということで個人的にはなんとも言い難いように思います。

検定結果は統計的に有意であったかもしれませんが，有意であることを以て真実であると主張できるわけではないですし，メタ分析などを通じて効果の大きさのより厳密な評価は必要かもしれませんね。
まだまだCox比例ハザードモデルへの理解が浅く，使い方も含めわかっていない部分が多いことから手続き自体に関する批判的検討ができないのが悔やまれますが，私個人の印象としてはそこまで気にしなくても大丈夫ではないかと感じました。
フォロー期間や説明変数を見ていて，赤身肉や加工肉の摂取が因果関係として発がんリスクを高めるという主張は，この結果からはまだ強すぎるようにも思われます。
もちろん論文では強い表現はつかっていないようですから，報道の仕方の問題もあるのかしら。

メタ分析をしている論文もあるようですから，後日気が向いたらabstractだけでも眺めてみようかと思います。journals.plos.org

*1:PLoS Med, 4

*2:イベントヒストリー分析（Event History Analysis），存続モデル（Duration model）とも

【R】任意の因子負荷から架空データを発生させる（1因子モデル編）

久しぶりの更新です。
更新したいネタはたくさんあるのですが，なかなか筆を執ることができません。

そんなどうでもいい個人的事情はおいておきまして，今回は備忘もかねて，統計ソフトRの{psych}パッケージを利用して任意の因子負荷行列から架空データを発生させる方法を書いておこうと思います。
{psych}パッケージのマニュアルはこちら。（※リンク先pdf注意）

大まかに次のような手順を踏みます。

1. 因子負荷行列を用意する
2. {psych}パッケージのsim関数を用いて相関行列に変換
3. {MASS}パッケージのmvrnorm関数で正規乱数を発生

とりあえずということで，まずは1因子モデルでの方法をこの記事にまとめます。
まとめるといっても手順は少ないので，ごく簡単に使う関数を載せて，簡単な具体例を付記するのみですが。
ちょっと手抜きをしており，各関数の引数は必要最小限のみ指定しています。

わかりやすさを重視して各セクションごとにコードを分割して記載，一番最後にコピーペーストや参照しやすさの観点からひとつにまとめたものを記載というスタイルといたします。

因子負荷の状況設定～相関行列への変換

まずはここからです。
せっかくなので，CiNiiでオープンアクセスできる心理尺度を参考に状況設定してみます。
今回は（自分の性格を鑑みて，ダメージ的な意味で）心に突き刺さるものがあったので教育心理学研究に投稿された次の論文を選んでみました。

小浜　駿（2010）．先延ばし意識特性尺度の作成と信頼性および妥当性の検証．教育心理学研究，58，325-337．ci.nii.ac.jp

この論文で作成された尺度の第1因子である「先延ばし中の否定的感情」を構成するとされる7項目について因子負荷を引用すると，.88, .82, .79, .77, .59, .57, .56という数字が得られます。
これを因子負荷ベクトルとして与えて相関行列を生成し，そこからデータを発生させるというわけです。*1

#psychパッケージの読み込み
library(psych)

#引っ張ってきた因子負荷をベクトルとして変数flの中に格納
fl<-c(0.88, 0.82, 0.79, 0.77, 0.59, 0.57, 0.56)

#psychパッケージのsim関数を使う
#今回は1因子モデルなのでsim.congeneric()を用いる
#生成した相関行列を変数cor.matの中に格納
cor.mat<-sim.congeneric(loads=fl)

任意の平均，相関行列から架空データ生成

ここまでで相関行列の生成ができましたので，このセクションではそこからデータを生成します。
{MASS}パッケージのmvrnorm関数で一発です，便利な時代になりました。

#MASSパッケージの読み込み
library(MASS)

#平均値のベクトルを用意する
#今回は7項目なので，各項目について計7つ平均があればよい
#今回は5件法とのことなので便宜的に3としました
mux<-c(3, 3, 3, 3, 3, 3, 3)

#架空データの発生
#作成したデータをデータフレームとするためにdata.frame()をかませています
#empiricalはTRUEならば指示通り（与えた通り）のデータを生成する模様
#nはサンプルサイズなので適当に与えてください
data<-data.frame(mvrnorm(n=100, mu=mux, Sigma=fl, empirical=TRUE))

これで，サンプルサイズが100の，7項目に関するデータが生成されているはずです。

コードのまとめ

2つに分割して記載しているので，ひとまとめにしたものも以下に用意しておきます。
説明のためのコメントが多いのは備忘録的性質の記事のためご愛敬ということで。

#パッケージの読み込み
library(MASS)
library(psych)

#引っ張ってきた因子負荷をベクトルとして変数flの中に格納
fl<-c(0.88, 0.82, 0.79, 0.77, 0.59, 0.57, 0.56)

#psychパッケージのsim関数を使う
#今回は1因子モデルなのでsim.congeneric()を用いる
#生成した相関行列を変数cor.matの中に格納
cor.mat<-sim.congeneric(loads=fl)

#平均値のベクトルを用意する
#今回は7項目なので，各項目について計7つ平均があればよい
#今回は5件法とのことなので便宜的に3としました
mux<-c(3, 3, 3, 3, 3, 3, 3)

#架空データの発生
#作成したデータをデータフレームとするためにdata.frame()をかませています
#empiricalはTRUEならば指示通り（与えた通り）のデータを生成する模様
#nはサンプルサイズなので適当に与えてください
data<-data.frame(mvrnorm(n=100, mu=mux, Sigma=fl, empirical=TRUE))

そのうち1～4因子まですべてまとめて因子負荷行列として与えたうえで架空データを生成できるようなものをまとめようという気持ちはあります。

論文で明記されているデータを再現するようなデータを発生させて，それを分析してみるなんていうのは勉強になりそうなのでいいですよね。

この記事では既存の関数だけを使ってやりくりしていますが，そのうち目的に応じて関数を組む必要も出てきそうです。

次の本，『みんなのR』と言いますが，バイブルクラスに出来がいいらしいので私も早く買います。
友人も即買って，べた褒めいていました。
なんだかんだ，Rもまだまだできない方なのが悩みの種です。

みんなのR -データ分析と統計解析の新しい教科書-

*1:本当は第4因子まであるので，4因子のモデルとして因子負荷行列を与えてデータ生成するほうが望ましいのだと思いますが，今回は単純化のために中から第1因子のみをピックアップして記事にしています。