• タグ別アーカイブ: COVID-19
  • 先を読む話

    まだまだ暑い日が続いていますが,いかがお過ごしでしょうか.いつになったら涼しくなるかと東京の10日間天気を見れば,まだしばらくは最高気温が30度以上で最低気温も25度以上の日が続くようでがっかりですね.実は,今年の夏は冷夏になるのではないかと,東北の古老が言っていたのです.7月に台風が一つも発生しなかったりして異常とは言われていましたが,過去の経験では,こういう年は不作に見舞われるのだそうです.気象庁の長期予報では猛暑であるとなってましたから,どうなることかと期待してました.古の知識が現代の科学技術を凌ぐというようなドラマチックな展開もあり得ただけに,科学的な予測に軍配が上がって内心ほっとしています.

    この暑さのせいか飴玉が売れないんだそうですね.暑さというよりも,コロナで通勤や通学の時間が減ったからなんだとかとも聞きます.確かに自宅で飴を舐めることはあまりないですね.どんな世の中になっても,食料品だけは不況に強いと思っていたのですが,ビールとかハムなども売れなくなっています.ビールなんかは自宅での消費は伸びているに違いないので,外食がシェアの多くを占めているんでしょうか.その他にも,一時期ほどではないにしても,化粧品が売れなくなったのはご存知だと思います.マスクと相性の悪い口紅だけは未だに売れないのだそうです.真夏にもマスクをする日が来るなどとは全く予想だにしていませんでした.先が見通せない時代になったんだということを強く感じています.

    先を読むことの難しさを示すために,実験計画セミナーの冒頭で枕にしているエピソードがあります.このブログで,以前斎藤一の話をしましたが,同じ新撰組の近藤勇に言及するところから始めます.司馬遼太郎の『新撰組血風録』で有名な近藤勇ですが,講談では「今宵の虎徹は地に飢えている。」という台詞が有名です.因みに彼が持っていた虎徹は偽物だったそうですが,この虎徹というのは、刀鍛冶の名前です.戦国時代から江戸時代初期,ちょうど宮本武蔵と時代が重なります.この虎徹ですが,元々は甲冑鍛冶だったそうで,材料に中古の甲冑を混ぜるという独特の製法が特徴だったとあります.それゆえ「古鉄」と呼ばれ,虎徹の名に転じたとか.この製法故に,虎徹の鎧はとても頑丈だったのです.おそらく鉄の含有炭素量などが一般的な甲冑の材料と異なっていたのではないかと思われます.


    甲冑鍛治として戦国の世に名を馳せた虎徹ですが,関ヶ原の戦いの後に豊臣勢も絶えたことから,50歳にして甲冑鍛治としての名声を捨て,刀鍛冶として出直すために江戸に移ります.天下太平の世の中では,武士の役割は戦士から政治へ変わっていくであろう,もはや甲冑は時代遅れのものとなり,これからは武士の精神的存在として刀のみが生き残っていくはず,そう考えたからです.凄い読みですね.確かに江戸時代を通じてそのようになって行ったのですから.今でこそ50歳の平均余命は30年以上ありますが,当時はおそらく20年くらいでしょうか.このコラムによれば,江戸時代の60歳の日本人の平均余命がほぼ14歳とあります. そのタイミングで心機一転して新しい仕事に挑んだという精神力には感服します.もちろんそれまでのキャリアは無駄にはなりませんでした.頑丈な甲冑を作るノウハウを刀に持ち込むことで,斬り合っても刃こぼれのない頑丈な刀であるとの名声は後の世に残っています.近藤勇が虎徹を欲しがったのも頷けます.


    もう一つ,先を読む話を続けますと,コロナによるリモートワークの普及も後押しして,アップルの株価が,ついに時価総額で2兆ドル突破したというニュースがありました.

    アップルの株価の変遷を知る面白いサイトがあります.What if I had bought Apple stock instead?という,もしもあの時Macを買わないで,その代わりにアップルの株を買っていたら今いくらになっているかがここに示されています.残念ながら,2012年4月からアップデートされていませんが,現在の株価は2012年の約5倍と考えていいでしょう.試しに,私が購入したPowerBook G3 300 (PDQ – Late 1998)は1998-09-01のリリースで当時$4999.0しました.(冒頭に当時の写真を載せましたが,まだ捨てずにしまってあります.)このサイトによれば,あのときそのお金でアップル株を買っていれば,今では$351365×5ですから,なんと1.9億円弱になります.計算間違いかもと思ったのですが,ほんとにそうなります.今新しいPowerBookが欲しいんですが,そのお金で600台近く買えます.このサイトを眺めるに,つくづく自分には先を読む能力がないとがっかりします.


    コロナ時代にも,虎徹のように先を読む能力を持つものが頭一つ突き抜けるはずです.とはいえ,凡人には真似できるものでもありません.唯一の手段が,天気予報のようにデータをもとに科学的に予測することです.そのためのデータを取得する方法が実験計画なんだよ,というお話でセミナーを始めると,気のせいかもしれませんが,より熱心に話を聞いてくれるように感じてます.


    なんだか雑なブログ記事になっていましましたが,暑い日が続いてやる気が失せている状況ですので大目に見てやってください.まだ暑い日が続きそうですので,皆様も体調管理にお気をつけくださいね.


    それでは.


  • アオサとかキャベツとか

     コロナ流行の初期に,感染防止に有効とされる食品がメディアに取り上げられました.覚えてますか?私の記憶によれば,納豆とかアオサとかがありました.アオサは某大学が公式に発表したものの,その研究のスポンサーがアオサのサプリを販売している会社だったことも問題視され,大学側が取り下げたはずです.一番問題だったのは,インフルエンザウィルスを使った実験結果でヒトコロナウィルスに対する増殖抑制効果を議論したことでしょうか.私のような素人でも乱暴と思いましたが,大学側の言い訳では,「ヒトコロナウィルスを含む各種ウィルスで効果があったから,新型コロナウィルスにも効くんじゃね.」としか言ってないのに,勝手に「新型コロナウィルスに効果があった」と曲解されてしまったとのことです.苦しい言い訳です.


    捏造ではないと思いますし,騙すつもりもなかったのでしょう.ただ,人の弱みにつけ込むことになってしまったという結果だけを見れば,真性の詐欺と変わりないように思います.そもそも,あのタイミングで発表すれば,この研究結果を参照する形で良い広告が打てます.出張で宿泊したときなど,早朝のTVショッピング番組を見ることがあるのですが.あの枠のほとんどが健康食品か美容関連の商品です.これらの番組では,ナントカ大学の研究によれば,この**が**に効くことが証明されています,などというカットが必ず入ってます.おそらく,最初からそれに使うことを目的とした研究もあるのではないでしょうか.

    フードファディズムという特定の食品に健康との関連性があると信じ込んでしまう認知バイアスがあります.面白いのは健康に良いというバイアスと良くないというバイアスとがあることです.後者の有名どころでは「グルテン・フリー ダイエット」でしょう.いわゆるベジタリアンもこの範疇かと.前者の例は,枚挙にいとまがありません.コロナ関連だけでも,アオサや納豆の他にもいろいろありましたね.当時は中国から帰国したインド人に感染者がいなかったことから,カレーが効果あるのではないかとか,岩手県に感染者がいなかったので,わんこそばがいいのではないかとか冗談も,そうとわかっていれば楽しかったものです.普段はフードファディズムなど寄せ付けないリテラシーを備えた人でも,何らかのパニック状態では簡単に騙されてしまいます.最近は食品だけに留まらず,この手の怪しい情報に振り回される大衆心理をインフォデミックなどと呼んでいるようですが,特定の用語ができてしまうほど,世界中で見られる一般的な現象なのでしょう.もちろんその背後には,従来は一部の人のものだったインターネットによる情報の拡散が誰でも気軽にできるようになったことにあります.

    コロナだけでなく311のときにもこの現象はやはり食品で多く観察できました.緊急事態宣言の頃はそもそも経済活動が停止してしまったので.下火になっていましたが,最近またコロナ予防に効く食品の記事を目にすることが多くなったように感じます.そのほとんどが,何らかの研究結果に基づくものではあるのですが,その研究を検証する姿勢がそもそもありません.メタノールがコロナ感染予防に効くとのデマで800人もの人が亡くなっているという事実を考えると危険なものを感じます.食品だから実害はないだろうという甘えもあるのだとは思いますが,白インゲンの食中毒事件も記憶に新しいところです.それに例のバナナダイエットのときのように特定の食品の購入が困難になるのは困ったことです.

    白インゲンならスーパーから消えても個人的には困りません.きゅうりはなくてもいいかなと言うところですが,キャベツは困ります.それというのも「コロナ感染症、キャベツときゅうりで死亡リスク低下との指摘」という記事が出ていたからです.下記に引用します.

    キャベツを食べる量が1日1g増えると死亡リスクが13.6%下がった。きゅうりを食べた場合は15.7%も低下する可能性があるという。

    オリジナルソースの著者は,WHO傘下の世界呼吸器疾病連盟元会長であるとWHOと関連付けて,仏モンペリエ大学医学部名誉教授という肩書とともにちゃんと権威付けもされて紹介されています.論文が参照されているわりには,リンクが貼られていないのは対象読者から仕方ないのかもしれませんが,そもそもライターは原著を読んでいないのではないかと疑いたくなります.

    原著はこちらです.Association between consumption of vegetables and COVID-19 mortality at a country level in Europe

    定番のmedRxivですね.いつも楽しく拝見させていただいています.記事ではmedRxivを世界的な医学論文公開サイトと紹介していますが,一番肝心なプレプレスということが抜け落ちています.査読前の論文は,査読に落ちれば単なる自由研究です.いや,それどころではなく間違っている可能性がある研究とも言えます.(査読に通ったからといって間違っていないとは言えませんが.)その段階の研究を一般読者を対象にして記事を書くのは大変危険です.それはどんなことでも記事にできてしまうからです.例えば,Moon Hoaxと呼ばれる有名なアポロ計画が捏造だったと主張する陰謀論があります.それを統計的に証明する研究結果をプレプレスすることは可能です.従って,それを記事にすることもできるわけです.プレプレスの目的は,研究成果の先行性を示すことにあるわけですから,最新の情報という意味では報道サイドにとって貴重なソースであることは間違いありません.中に正しい研究が埋まっているのも事実です.ですが,少なくとも一般人向けの報道ではプレプレスをソースにすべきでないと思います.それというのも,数理統計学の論文と違って,この類の疫学的な論文では,一般人にも結論だけは理解できてしまうからです.しかも,この記事のように間違った解釈,おそらく意図的に間違っているのでしょう,で読者を誘導することはいかがなものでしょう.

    この記事の解釈が間違っていることは,原著を読めばはっきりわかります.まず第一に,研究の手法が,欧州各国を対象に調査をしたところ,としか書かれていませんが,アブストラクトにははっきり書かれています.それによれば,欧州食品安全機関(EFSA)のデータベースの各国の野菜の消費量(アブラナ科野菜,ホウレンソウ,キュウリ,ズッキーニ,レタス,トマトなど)に加えてEuroStatの国内総生産(GDP),人口密度,64歳以上の人々の割合などの各種データを合わせて,各国のCOVID-19による死亡率を説明することを試みたようです.要するに実験したわけではなく,データベースからダウンロードしたデータをモデリングした研究だったわけです.

    この手の論文のお決まりとして,Discussionのlimitationのところにこうあります.

    (本研究で示された)関連性は因果関係を意味するものではなく,食物消費パターンは他の未知の要因の代理なのかもしれない.

    要するに擬似相関かもよということです.それなら何でこんなことしたのかというと,意訳するとこんな感じでしょうか.

    パンデミック分布の指標が更に研究されるまで,各国の地理的な(コロナ死亡率)の変動は説明できないけど,この研究は,国の単位でコロナ死亡率と食物消費を関連付る最初の試みとしては面白いよね?

    著者としては,安全弁を設けたのでしょうが,そのような意歳を一切無視して,「きゅうりを食べるとコロナでの死亡率が15.7%低下する」などと報道されては著者もたまったものではないですね.しかも,某クリニックの院長先生とか栄養管理士のエンドースメントまでついています.最も,このことはこの方達には責任はないのかもしれません.この論文をどう思いますか?なぜこんな結果になるのでしょうか?などと聞かれれば真面目な方だったら理屈を考えますよね.その返事をそのまま,このようなこの研究を肯定するコメントとして報道側が掲載した可能性もあります.もしかしたら,複数の人にコメントを求めて,そのような記事にしやすい回答だけを載せたのかもしれません.

    もちろん,このようなわかりやすい例であれば,直感的におかしいと気づく人が多いと思います.ですが,そうではない人が世の中の多数なのだということが問題です.例のイソジン騒動でも明らかになりました.プレプレスの論文は一般向けには報道すべきではないのはもちろんですが,それ以上に問題なのは,報道側のバイアスです.報道の引き出し効果とでも言いましょうか.「都道府県別新型コロナウイルス感染者数マップ」を公開してくださているジャッグジャパン株式会社の大濱﨑さんのTweetによれば,某社から緊急事態宣言の解除が第二波の原因であることを解説してほしいと依頼されたそうです.そこで「そうではない」ことをデータをつけてレポートを送ったところ,何の返事も来なくなったとのこと.報道以前に礼儀にも欠けるとこのような話は他にもいろいろと聞きます.

    さて,この論文では,過剰分散を考慮して死亡率をモデル化するために,死亡率をquasi-Poisson回帰モデルで分析しています.計数データのモデリングは,JMPでは一般化線形モデルで『分布』を「Poisson」にすることで実施できますが,データが過分散である場合は,『リンク関数』の下にある「過分散に基づく検定と信頼区間」にチェックを入れるとモデルに過分散パラメータを含めるとができます.こうすることで「モデル全体の検定」レポートにある「適合度統計量」にそのカイ二乗値を自由度で割った値の「過分散」という列が出てきます.この値が過分散パラメータで,残差逸脱度と残差自由度との関係性の指標で,大きければ過分散であることを意味します.

    Poisson分布と二項分布では過分散データに注意することが教科書には書いてあります.over-dispersion検定をすることを推奨している文献もあります.その結果として,データが実際に過分散だった場合に,その計数データをいかにモデリングするか?という問題に直面することになります.この場合によく使われるのが,擬似尤度を導入して,分布に分散の調整のためのパラメータを加えた,負の二項分布やquasi-Poisson分布を指定することです.後者は日本では擬似ポアソン分布とも呼ばれますが,勉強不足で自分にはあまり使い分けも違いもわかってないです.確か久保先生の青本に書かれていたと記憶しているので,詳細が知りたい方は参照されてください.自分の解釈としては,要するにデータのばらつき具合の推定値による誤差の調整のようなもんでしょうか.このモデリングはJMP Proならば『分布』のところに「Negative Binomial」というのがあって,それで実施可能と覚えているのですが,うろ覚えです.もしかしたら違うソフトかも.

    そもそもポアソン分布にとって,過分散は本質的なのかもしれません.母数はλだけで,期待値と分散がともにλだということを思い出してください.平均値とともに分散が大きくなるので等分散性を満たさないということは,何かの回数を記録したカウントデータににとって自然なのでしょう.因みに,指定した分布のもとで線形化するための変換関数が「リンク関数」で,ポアソン分布の場合は対数変換とするのが通常です.

    過分散のモデルに与える影響は,誤差が小さく推定されるので,有意差が出やすくなり,反対にその補正の結果,有意差が見逃されやすくなると考えればいいと思います.この論文では安全側に倒して擬似ポアソン分布回帰モデリングを実施したのでしょう.本当はこのデータを元に結果をJMPで再現したかったのですが.あいにくすべてのデータが公開されていないようですし,JMP Proが手元にないので...そのかわりと言っては何なんですが,この論文の Table3をグラフビルダーで素早くグラフにする練習をしてみました.その結果が冒頭のグラフです.確かに,きゅうりとキャベツの効果が目立ちますが,カリフラワーはどうなんだよ,とかいろいろ気づくことがあります.

    このグラフの作成にはPDF の読み込みから初めて,3分かかりませんでした.JMP15のPDFからテーブルを読み込む機能は英語だとうまくいくことが多いですね.medRxivにはデータが添付されていて,統計分析までフォローできる論文も多いので.JMP初心者には良い練習材料になります.もちろん,その研究を批判的に読んで.統計リテラシーを磨くことも忘れないようにしたいですね.

    それでは時間となりましたので,これにて.


  • コロナ時代のメディアリテラシー

    最近の報道ではコロナ対策についての主張が二極化しています.以前から言われていた,コロナなんてただの風邪という意見と,現状の危機を認識しPCR検査を積極的に実施して,感染者を炙り出すことで,部分的にでもコミュニティをロックアウトせよという意見です.それに誘導されてでしょうか,私たちの間にもある意味対立が生じてしまっているのは残念なことです.二週間後には日本はどこどこの国のようになる・・・という従来からのマスコミの基調に乗せられている人々もいる一方で,渋谷ではクラスターフェスとかいうノーマスク,ノー自粛を標榜する団体の催しが開催されたりしています.そもそも,医療や政治あるいは経済の観点だけから最適解を考えても意味がないのは確実です.自動車の運転テクニックにヒール・アンド・トーというのがありますが,まさにアクセルとブレーキを一本の足で踏み分ける高度なテクニックが要される状況と言えましょうか.いずれにしても,何らかの妥協を前提とした準解を探さなければなりません.当然,立場の違う人々の間の交渉が必要となるので,そのお膳立てをするのが政治の一つの役割のはずです.そのような方向でのリーダーシップをぜひ期待したいところです.

    報道にはこの両者の主張が織り交ざっているのでやっかいです.両者は互いに思うところが違うわけなのですが,純粋に情報を求めているものにとって,正義や思惑はいい迷惑です.とはいえ,情報発信者の思惑を窺い知るという統計リテラシーを鍛える良い練習になります.例えば,ほぼすべての報道はコロナ感染者数しかも今までの累計をトップに掲げるので,およそ45500人ということを知っていると思いますが,そのうち既に31000人程度は退院していることをご存知でしょうか.重症者数が増加傾向にあると報道されています.確かにその傾向は見えていますが,現在のECMO治療を受けている患者数を知っている人は少ないのではないでしょうか?本日の朝見たところでは,12人です.因みに日本にECMOは1412台ありますから1%以下です.時間差のため追従性がよくないので,感染状況の把握にはむいているわけではありませんが,おそらく,ほとんどの人がこのデータを知ると安心するのではないでしょうか.だから報道されないのです.メディアが稼ぐためにはより刺激的な情報で,できれば視聴者に恐怖を感じさせなければなりません.その思惑が分かっていれば,昨今の報道を色々と楽しむことができます. 

    統計リテラシーというよりはメディアリテラシーと言った方が正確かもしれません.でも,根っこは同じです. 百聞は一件にしかずと言いますが,報道では視覚的情報が重要視されます.例えば.TVに棒グラフや円グラフが多いのもそれが理由だとか.そう言えばTVでは散布図はあまり見ないですよね.1秒で誰にでも分かってもらえるようにという配慮からと言えば聞こえはいいですが,そのことを逆手にとって,誰でも騙しやすいように棒グラフを使っている例が散見されます.『統計的問題解決入門』にも書きましたけど,3D円グラフとかは目の錯覚を起こしやすいグラフの典型です.円グラフの三次元化は厳禁であると,生徒さんには教えてますが,それはグラフ要素が面積のグラフでは,斜めから見ると観測値が比例しなくなってしまうからです. 

     視覚情報として最も効果が高いのが写真です.写真は加工できてしまうので,データとしての価値は低いのですが,例えば,1ヶ月以上前になりますが,この記事を見てください.
    「コロナ前の喧騒、程遠く 人気居酒屋、閉店も検討―宣言解除25日で1カ月・アメ横」 この写真のトリックがお分かりでしょうか.既に同様な手法が色々と暴かれていますよね.この上下2枚の写真は異なるレンズで撮影されています.専門用語で言うと,上の人が少ない写真は焦点深度が深いレンズで撮影されています.その証拠がこちら.アメ横商店街のStreet Viewですが,こちらが解除前の5月6日の写真の撮影場所です.一方,こちらが解除後に6月23日に人通りが戻ってきたと言う写真の撮影場所です.そんなバカなと思うかもしれませんが,ある点に注目するとそれとわかります.一つだけヒントを.上の写真で,中央の赤いカラオケ屋の看板を見つけてみてください.場所はここです.この赤い看板を下の写真でも探してみてください.わかりますか?アメ横のアーチの「横」と言う字に上が隠れているのがそうです.まったく違って見えますよね?しかも,このことをごまかすために下の写真では撮影地点をかなり後方にずらしてアメ横のアーチが同じよう見える撮影ポイントを選んでいるのです.

    因みに,この場所を知っている人ならお分かりのように,この先で上野中通商店街と合流します.どうやら,解除前もその合流地点より先はそれなりの人通りがあったようにも見えます. お店が閉まってれば商店街を通る人が少ないのは当たり前なので,人の流れが変わっただけなのかもしれません.真実はこの報道からはわかりませんが,少なくとも写真撮影に細工をしているのは間違いなく,問題はなぜこのような細工をしたのか?ということです.単にニュースになるから,絵になるからと言う理由だけかもしれませんが,こんなことをしてまでも私たちに何かを思い込ませたいのではないかと勘ぐってしまうのです.それはなにかはこの場では書かないでおきます.

    今回のパンデミックが世界史に残る事件となったことは間違いありません.そして,この事件では,報道関係者だけでなく,医療関係者や研究者の姿勢も問われています.医療器具や各種検査装置をはじめ世の中を支える様々な技術に携わる技術者も同様です.自分としては,広い意味で技術倫理を考える良いきっかけになってます.あるいは自分自身に問うてみていると言ってもいいかもしれません.「Do Artifacts Have Politics」という社者科学者のウィナー先生のエッセーがあって,Langdon Winner(2000)『鯨と原子炉』紀伊国屋書店にも所収されているのですが,そこでの考察はとても重要です.それは技術者がこの世の中に創出する製品には政治があるということです.

    現在進行形ではありますが,この事件をいかに糧にしていくかが私たちに問われていると思います.ニューノーマルと呼ばれる生活様式もその一つです.日本では普及していない高水温が使える温水洗濯機がヨーロッパでは珍しくないのは,ペストの感染防止にお湯で寝具を洗濯した名残だとも言われています.仕事でも,オンライン会議が増えて,出張が減っていくことでしょう.このような目に見えることだけにとどまらず,ポストコロナの考え方というのも重要と思っています. そのうち,それらをまとめられたらと考えています.

    少し脱線してしまいた.それでは,これにて本日は失礼します.


  • また抗体検査の話

    一段と世の中が騒がしくなってきましたね.そんな中で菅,今日もブログだけは平常運転でいってみたいと思います.そうはいっても,特に書くこともないので,先週に引き続き,時事ネタを強引にJMPと紐付けてるという恒例のチャレンジをしてみます.先週は,横須賀市のコロナ抗体検査の結果の解釈に,統計的でもないのに「統計的に」という枕詞をつけるなと文句を言いました.そもそも「統計的に」という言葉をつける言説には,騙されるなという警告の意味があると解釈しています.さて,最近は,様々な自治体で同じ試みが実施されているようです.数日前のニュースでは宇都宮市でも無症状者を対象に抗体検査を実施したようです.

    こちらがその記事です.これによれば,742人中3人の陽性者が出たとのことです.この場合も,2290人を無作為抽出してはいるものの,応じたのは742人ですからそのうちの1/3です.やはり,結果としては無作為ではなく,なんらかのバイアスがかかってしまっていると疑えますが,そのことは今日は脇に置いておきます.引っかかったのは,記事のこの部分です.


    以下引用
    この時点の感染率は0.40%だが、統計学的に精度の高い数字に補正すると1.23%になるとしている。6月1日時点の同市人口(51万8610人)に換算すると感染者数は推定6378人。「第1波」とされる時期に市が把握していた陽性者23人の277倍に当たるという。
    引用ここまで

    「統計学的に精度の高い数字に補正すると」という部分に注目してください.これどういう意味でしょうか.そこで先週と同じように,「一変量の分布」で信頼区間を求めてみます.(普段は英語のUIを使うことが多く,裏で分析途中のテーブルが大量に開いている状態で,日本語UIに切り替えるのが面倒なのでそのままですいません.)
    上限は1.18%ですから,精度の高いという1.23%と異なってます.この信頼区間はそこにも書いてあるように,スコア信頼区間,日本ではWilsonのスコア信頼区間と呼ばれる二項分布から計算される信頼区間です.

    おそらく,記事にある「統計学的に精度の高い数字」は信頼区間とは関係ないところからきているのか,それとも異なる統計量なのか?気になったので,他にもある信頼区間を求めてみることにしました.残念ながら,JMPでは「一変量の分布」プラットフォームで出力される信頼区間は限定的なので,Pythonを使ってみます.といってもプログラムを書く必要はなく,
    from statsmodels.stats.proportion import proportion_confint
    で簡単に実装できます.
    proportion_confint(count, nobs, alpha=0.05, method=’normal’)の引数のmethodのパラメータはデフォルトではいわゆるWald法による信頼区間を与える「normal」ですが,その他に次の種類が用意されています.
    • normal : asymptotic normal approximation
    • agresti_coull : Agresti-Coull interval
    • beta : Clopper-Pearson interval based on Beta distribution
    • wilson : Wilson Score interval
    • jeffreys : Jeffreys Bayesian Interval
    • binom_test : experimental, inversion of binom_test


    検証の結果から言いますと,おそらく1.23%という上限値はAgresti-Coullの信頼区間ではないかと予想できました.Agresti-Coullの信頼区間はWaldの信頼区間を修正したものですけど,サンプルサイズが大きいときは大差ないと認識していました.興味ある方は,どうぞガチの統計論文ですが,こちらをお読みください. 確かに値は異なりますが,統計学的に精度が高いというのは,全くのミスリーディングですね.統計学学的により精度の高い信頼区間と言うことはできるかもしれませんが,そもそも95%の信頼区間の上限値ですから,1.23%という値に意味を持たせてはいけません.

    因みに,日本では,Agresti-Coullと英語表記のまま書かれることが多いのですが,Agrestiの方はカテゴリカルデータ分析の大家であるフロリダ大学のアラン・アグレスティ先生で,Coullの方は先生のもとで当時Dr.だった現在はハーバード大学で生物統計学の教授のブレント・クール先生です.Wikiなんかではコウルと訳されていますし,有名なイギリスのCoulle Quartet は日本ではカウル弦楽四重奏団と呼んでいたりと混乱していますが,少なくとも人名ではクールと読むのが普通です.ですから,アグレスティ-コール信頼区間と呼ぶべきだと思うのですが...

    閑話休題.JMPではAgresti-Coull信頼区間は,私の知る限りでは直接出力できませんが,Agresti-Coull検定として,実験計画メニューの「計画の診断>標本サイズ/検出力」で『1標本割合』を選択すると登場します.デフォルトでは「方法」が『近似Agresti-Coull検定の正確検出力計算』となっているはずです.この手法で,今回の宇都宮市の結果を評価してみます.少なくとも0.4%と1%の違いが検出できなければ結論が無意味ですから,割合の仮説値を0.01として,検出力と標本サイズの関係を示したのが冒頭に示したグラフです.(結局,日本語UIに変更しました.)これによれば,0.8を得るための標本サイズは1700人程度でしょうか.全然少ないですね.逆に742人では検出力は0.25程度です.統計的に信頼できる結果ではないと言うことです.お金と時間の無駄だったとまでは言いませんが,やるのであれば政府主導でもう少し大規模に実施することが必要かと思います.個人的には,学校などである程度強制的に検査できなければバイアスの影響が気にはなるところです.

    それでは,本日はこれにて失礼します


  • 統計リテラシーと抗体検査

    今日読んだニュースによれば(今日のニュースではありませんけど),インテルが7nm半導体の技術開発が予定より6カ月遅れていることを明らかにしたことで,株価が1割近く下がりました.Mac搭載のCPUとして毎日お世話になっていますけど,つい先日もAppleがCPUをインテル製から内製へと移行していく方針が発表されたばかりです.自社の製品開発のロードマップを他社の開発の遅れに妨げられないようにすることが目的と言われています.おそらく,インテルも今後は外部への生産委託を拡大する方向へ舵を取っていくことになると思いますが,同じ業界に長いこといて,かつての同僚のアメリカ人なども働いているので他人事のようには思えません.



    7nmといえば,コロナウィルスのサイズが100nmですから,いかに微細なパターンということかお分かりと思います.このような高度なプロセス技術をもはや一社で抱えるのは難しいのかもしれません.それは,現代では設計者と大工とが分業するのがあたりまえの,家作りにも似ています.いつだったか「アテ」の話をブログに書いた記憶があるのですが,昔は大工の棟梁が山に行って材料の選定からしていたわけです.もちろん,そういう家は特別です.大量製品を前提とした工業分野では,理想なのだということもわかってます.でもですね,中国やアジア諸国と互角に渡り合うには,そこに巧の精神というか技が必須だと思うのですよ.大量生産の工業製品であっても,少なくともそれらをより強くタフに,そしてより高性能にしていくための鍵が上流から下流までを貫く巧の技です.

    その巧の技を阻害するものの一つが分業化です.ですから,同じ会社組織に臆していてもそこに十分なコミュニケーションが図られていなければ意味がありません.インテルでも設計と製造とが分業化してしまい,同じ会社であってもデータの受け渡し程度の関係になっているのかもしれません.このことは,他の会社組織を見ても想像に難くないです.もちろん,分業化には経営効率の観点からは多大なメリットもあるので,いかに部門間のコミュニケーションをとっていくかが巧の技にとって大切になります.そのベースになるのが統計学であり,それ以前の統計リテラシーと考えています.これについは,そのうち書きます.

    というわけで,近頃は統計リテラシーの講座なんかをいろいろ開いているのですが,今回のコロナ騒動で多くのネタが拾えるので,その整理に追われている毎日です.例えば,治療薬として期待されている「アビガン」の臨床研究についての藤田医科大学の報告が報道されました.明確な有効性は確認できなかった,即ち有意性なしとの結果について,サンプルサイズを大きくすれば有効性が出てくるかもしれないとのことです.まだ,実際の報告を読んでいないので,メディアの報道だけから判断するのは危険ですが,『JMPではじめるデータサイエンス』の読者の方であれば,事前に検定力をもとにサンプルサイズを決めてるべきであることお分かりと思います.望み通りの結果が得られないからといってサンプルサイズを増やすのは,特に今回のように人命に関わる重要な研究では口に出すのさえ憚られることです.間違った検定で間違った判断を下すことの影響ははかり知れません.研究チームには統計の専門家もいるはずなので,おそらくメディアの誤解によるものと思ってますが,正直驚きました.この件については,後日検証してみたいと思っています.


    他の報道でも,首を捻ることがありました.横須賀市が無作為に抽出した市民の抗体検査をしたのですが,その結果を受けて,市長が「人口に比例させれば、4000人が感染したことがあるということになり、非常に高い数字だと受け止めている。自分も感染している可能性があるという意識を持って行動をとってもらいたい」と発言しています.自分も感染している可能性があるという意識を持って行動するのはもちろん,その通りですが,統計リテラシーでこの結果を捉えると,いくつかの問題が見えてきます.

    1.まず第一に無作為抽出したのは2000人でそのうち964人が応じた結果だということです.このため,本人しかわからない不安な履歴がある人が応じている可能性があります.そこにバイアスがかかり陽性率が上がっていることが考えられます.本当の無作為ではないということです.
    2.別の報道によれば,こんなことも書かれています.「厚生労働省の抗体検査では2種類の試薬で両方陽性だった人を抗体保有者とみなしているのに対し、横須賀市は1種類だけの陽性判定でも保有者にカウントしている。」ご存知のように,抗体検査ではIgMという現在感染していることを示す抗体と,IgGという過去に感染した可能性を示す抗体を二値判定します.横須賀市では,おそらくIgGのみでも陽性であればカウントしたのでしょう.発症後二週間程度経たないとIgGは上昇しないのでIgM陽性だけをカウントするのなら理解できますが,IgGのみ陽性者をカウントするのは大いに疑問です.無症状のままに治癒した人も網にかけたいという気持ちはわかりますが,単なる風邪でも陽性になってしまうことがあるからです.もともと精度が低い抗体検査は発症後の治癒の状況を確認する目的で使うべきです.二つの抗体の判定結果を踏まえ,臨床所見や症状から最終的に判定する必要があります.検査手法に大いに疑義ありといったところでしょうか.ダメなデータをもとに正い判断は降せません.こんなデータならない方がマシです,というのは言い過ぎでしょうか.
    3.報道によれば,「上地克明市長は記者会見で「統計上は市民4000人が抗体を保有していることになり、非常に多い」と強調した。」とありますが,統計上はという言葉は安易に使うべきではないと思います.冒頭のJMPのレポートを見てください.今回の横須賀市の結果をもとに,「信頼区間」を示してみました.この結果によれば,「横須賀市の人口39万人のうち抗体保有者は,信頼度を95%とするとおよそ2200人から7400人の間にあると推定できる.」と言うべきであることがわかります.確かに21日時点での,感染者累計69人よりは多いです.ですが,69人は有症状の感染者数です.比較すべきなのは無症状の感染者数なのですが,もちろんそれは不明です.


    因みに,この信頼度を99%にすれば信頼区間は1833人から8869人となりますし,更に,10人のうち半分が偽陽性だとすると,下側は677人になります.非常に多いと判断するのではなく,こんなもんだくらいに受け止めればいいのでしょうが,このために440万円かける価値があるのかは検討すべきですね.


    この分析をして気づいたのですが,レポートの下に,「注:スコア信頼区間を使って計算」とありますよね.こういうところがJMPらしいですね.せっかくなので少し解説しておきますと,統計学の教科書では母比率の信頼区間を次式で習います.ここでnをサンプルサイズとして,pとπはそれぞれ標本比率と母比率で,zは所望のアルファに対応したz値です.95%なら1.96ですね.覚えてますか?
    ですが,この式はサンプルサイズが十分大きくて正規分布が二項分布の良い近似となる場合でないと正確ではありません.今回のようにnやpが小さいときは全く使えません.そこで2標本のカイ二乗検定いわゆるピアソンの検定から,求めた信頼区間がスコア信頼区間です.式をここに書こうかと思ったのですが,結構複雑なので,興味ある方は一般的に呼ばれているウィルソンの信頼区間で調べてみてください.もちろん,JMPを使えばそんな式は存在すら忘れても良いのですね.

    通常の統計ソフトならば,nやpが小さいときは,ウィルソンの信頼区間を使うということを知っている必要があるますが,JMPだとそこを良きに計らってくれます.とはいえ,これが正規分布の母比率の信頼区間とは違うことを知らなければ,計算結果が異なることに戸惑います.そこで注が書いてあるわけです.英語版では「conputed using score intervals」と書かれているので直訳です.日本のユーザーには,スコア信頼区間ではなくウィルソンの信頼区間とした方が親切だと思うのですが,いかがでしょう.他にもWelchのt検定なんかも,統計学を前提としない間口の広さはJMPらしいですが,日本語の教科書で勉強した人には不親切ですね.その割には,「Fit Y by X」が「二変量の関係」と妙にこなれた日本語になっていたりするのが不思議なところです.聞いたところでは,芳賀先生がそのように命名されたとも聞いていいます.変な用語もありますが,まあJMPは外人だと思って付き合っていくのがいいのでしょうね.


    それではまた.

  • K値について思うこと

    一昨日に開催されたJMPの製薬業界向けWebセミナー「非臨床研究者のための統計解析入門」を聴講しました.私にとって馴染みの薄い分野なので,あまり使わない機能を紹介して頂くのは勉強になります.そもそも非臨床研究というのが聞き慣れない言葉だったのですが,話を聞くと動物などを使った医薬品の薬効や安全性の検証試験である前臨床試験のことなんですね.最近では非臨床試験と呼ばれていることは知りませんでした.非というのは良い響きではないのですが,おそらく英語のNon-clinicalを直訳したのでしょうか.工業分野の実験計画とは大きく異なって,試験の手順や分析方法,信頼性検証に至るまでがガイドラインでガチガチに固められているので,覚えることが多く大変そうです.人命が関与しているので仕方ないこととは言え,プロセスが決まった開発は自分に向かないなといつも思います.

    とはいえ,信頼性検証試験などとも関係が深い分野なので,仕事をしながらのながら聴講させていただきました.そこで「発展的なモデル>曲線のあてはめ」を紹介されていたのですが,ふと目に止まったのが,「シグモイド曲線」の中にGompertz曲線があったことです.こんなところにあったんですね.完璧に見落としてました.Gompertz曲線といえば,昨今巷ではコロナの新規感染者数の予測モデルとして有名になりました.ちなみに参考文献はこちら.あ,これアーカイブですからご注意ください. 予測モデルというよりはデータにあてはめたらフィッティングが良かったということなんですが,この発見を予測に使えると考えてしまう人もいるようです.一週間の累計を取っていたりと,このデータ特有の状況は勘案されていますが,大阪モデルで有名なK値もGompertz曲線のフィッティングで導出できるから根っこは同じです.ちなみに,K値というのは大阪大学核物理研究センターの中野先生が考案されたコロナ感染の状況を示す指標です.その定義は,先生の公開されているスライドによれば,週あたりの感染者数増加率として直近一週間の感染者数を総感染者数で割った値です.数式で書けば,a=直近一週間の新規感染者数の合計,b=一週間前の総感染者数として,K=(a-b)/a となります. JMPで描いたK値の経時変化を描いてみるとこのようになります.Days=157は7月15日です.

    確かに,感染状況の指標として使えるように見えます.冒頭の図は累積感染者数のグラフにK値を4つに順序尺度化して「重ね合わせ」してみました.感染拡大期の紫→緑から,収束期に青→赤となり,そして今また赤になっているという解釈が可能です.K値は,過去一週間の感染者数の増加率をその日の感染者数で規格化した数値ですから,K=0.01ならば,直近の一週間で総感染者数が1%増加したことを意味します.K値が感染状況の指標となるのはあたりまえのことかもしれません.このように,感染状況の指標としては何らかの意味付けができるK値ですが,一部のメディアで基本再生産数の代わりにK値が予測に使えると報道されています.これには疑義ありです.部屋の気温は温度計で分かりますが,これから寒くなるのか暑くなるのかは,時間や季節など他の情報に加えて法則性の知見が蓄えられていなければ予測不能です.

    確かに,Gompertz曲線は信頼度成長曲線とも呼ぶくらいですから,予測にも使えることは間違いありません.だからこそ,JMPのコマンドにも実装されているわけですし,システムに残っているバグの検出率を予測して,システムの品質向上と納期や人的リソースの確保の見積もりなどに,実際に使われています.ですが,この場合でも予測するのは早くても開発の中盤以降になるのが普通です.

    もちろん,タイミングが遅いほど予測精度は高くなります.開発初期に予測しても,全く使い物になりません.怖いのは予測ができてしまうことです.最悪なのは,甘い見通しを立ててそれを信じた結果,プロジェクトが立ち行かなくなってしまう事態に陥ることです.よく知られたWeibull曲線による製品の寿命予測とも状況は同じです.サンプルをほとんど破壊しないデータで寿命予測するケースがありますが,その結果は全く役に立ちません.加速試験はもちろんですが,少なくともサンプルの80%は破壊するまでデータをとらないと,やはり市場不良を引き起こす危険があります.

    そこでK値ですが,やはり状況は同じと思います.過去の感染者数に対してあてはめるには問題ありませんし,現在の状況の指標としても(他の指標と組み合わせる前提で)使えることは間違いないと思います.おそらく,予測もごく短期間(数日)程度であれば,現状の外挿として意味はあるかもしれません.ですが,報道にあったようにピークアウトのような長期的な予測には使うべきではないと考えます.経済と感染対策との間にあって,現状を楽観視できるK値による予測は,経済に重きを置く考えの人々の論拠になっています.けれども,信頼性予測のように,ここでも状況を楽観視してそれが事態を悪くしていく危険があります. 事実,この数日の感染拡大でK値のフィッティングがずれて来ています.

    驚いたのは,このズレを別の理由によって説明していることです.具体的には,複数のK値の予測曲線が重なっていて,それらを層別化すればK値の予測はうまくいっているとの説明です.これは真実かもしれません.ですが,この説明は天動説の学者の説明と同じです.観測結果が精密になるにつれて,天動説では説明しきれない惑星の動きが彼らを悩ませました.そこで持ち出したのが周転円です.更にその言い訳はプトレマイオスのエカントへと精緻化していきます.正直すごい頭脳と思いますが,それらは地動説を真実とすることであっさりと消え去りました.

    誤解のないように申し添えておきますと,私は経済とよりも感染防止に重きを置くべきと言っているわけではありません.ただ,K値にこだわり続けていると本質を見逃し,対策が遅れることは危惧しています.K値について考えることで,改めてモデリングではオーバーフィッテイングは怖いなと感じた次第です.『統計的問題解決入門』のP153(第1刷)にも書いた,宿屋の屋根に止まったカラスに紛らわされてはダメです.誤差(純誤差)を本質と分離したモデルこそ真に役立つモデルではないでしょうか.

    ではでは.


  • ポストコロナの適者生存

    5月,6月と開催されていたJMP On Airも今週で終了となりました.コロナ期間中は自宅で楽しみながら勉強させていただきました.あらためてSAS社の皆様にお礼をしたいと思います.第一回の開催ではまだオンライン視聴に慣れていなかったということもあって,少しドキドキしてスタートを待ったのは良い思い出です.1時間もヘッドセットをつけているのは辛いので,スピーカーから音を出しても邪魔にならない場所に移動して視聴していたのですが,当初は自宅の通信状態が安定せず,時々音声が途切れがちになったりして苦労しました.自宅のWifiは5GHzと2.4GHzを使い分けているのですが,調べるとその切り替えが物理的な場所によってうまくいっていないことが判明しました.幸い,その部屋には有線LANを壁に埋め込んでいたので,それを使うことができました.電気屋さんにCD管だけ埋め込んでもらって,自分で配線を引き回したのですが,当時は既に無線LANの時代でしたから,こんな手間をかける価値があるのかは悩みましたけど.実際,ほとんど使っていなかったのですが,まさかWifiの速度が有線を抜くような時代になって,有線が活躍するとは思っても見ませんでした.

    実際使ってみると,間違いなく早いし何よりも安定しているのは,今後オンライン授業をする立場としては心強く,少なくとも私の環境では有線は適者生存をしていくようです.ポストコロナ時代には変化に対応できないものは淘汰されるという趣旨のことを,(実際には言ってもいない)ダーウィンの言葉として引用した政治家がいました.変化に対応するという能動的な能力よりも,偶然を受け入れることができたという,どちらかといえば受動的な能力が,その種が生き残るには重要なんではないかと思っています.哲学的に色々難しいことを言う先生もいらっしゃいますが,要するに,運が良かったからに過ぎないのではないかと.その偶然を受け入れることで,たまたま適者となって生き延びたと言う例がたくさんあります.体色が白い蛾が,ヨーロッパにおける工業の発展に伴って発生した公害の影響で黒っぽくなっていったという現象を工業暗化といいます.19世紀後半に観察されたこの現象は,どの生物学の教科書にも載っているほど有名ですが,最近,と言っても2006年ですが,その原因が究明されました.Hof, A., Campagne, P., Rigden, D. et al. The industrial melanism mutation in British peppered moths is a transposable element. Nature 534, 102–105 (2016) によれば,ゲノム上の位置を転移するので動く遺伝子と言われるトランスポゾンの挿入が工業暗化の正体であることが突き止められました.それを変化する能力と考えることもできますが,一般的には偶発的な現象の結果として暗化した個体が発生し,それがたまたま公害で黒くなった樹木が保護色となって生き残ったのだと考えてもいいと思います.
    偶然を受け入れる能力が自分に備わっているかということを,ポストコロナ時代(未だ渦中ではありますが)の今考えています.適者生存というよりは運の良いやつが生き残る時代です.どうすれば良いのかを統計思考で考えると,試行(trial)と揺らぎが重要なのはのではないでしょうか.運はある程度制御できます.それにはとにかく試行回数を大きくし,試行を繰り返す度に行動に小さな変化を取り入れるのです.このテーマは古今東西,各所で取り上げられてます.古くは星新一のショート・ショートにもあったような記憶がありますし,手塚治虫の『ザ・クレーター』にもこのテーマの短編が収録されています.1993年のアメリカ映画『Groundhog Day』もそうでした.(今調べたら邦題は『恋はデジャ・ブ』なんてひどいのが付いてたんですね.)
    ポストコロナでは,日常でも変わったなと気づくことが多々あります.近所に深夜まで営業している割と有名な食堂があります.お客さんが10人も入らないくらい小さいお店なんですが,最近外から覗いたらその狭い店舗になんと間仕切りを設置して営業していました.その光景を見るととてもシュールです.お店に入るときに体温を計測する(される)ところも多くなったということですし,自分自身も仕事で人前でお話しするときに,非接触型の計器で体温をチェックされたりします.咳と発熱が特有な症状という感染症ですから,最も簡単な検査法ではあるのですが,信頼性はいかほどのものなんでしょうか?正常域の上限を37度あるいは37.5度としているところが多いと思います.前者は水銀体温計のデザインから,後者は法律による定義なのですけど,実はともに科学的な根拠はありません.調べれば,日本人の平均が36.89±0.34度という値があちこちに書かれています.これは田坂定孝,日新医学 44(12), 635-638, 1957が参照されているのですが,ここに示されているRangeは1σです.だから結構多くの人が上限の37.23度を超えます.このデータは古いので,現代の人の実態を反映しているのかも不明です.少なくとも老人性低体温症というのがあるくらいですから,年齢によっても基準が違うはずですし,ホルモンの日内周期の影響も受けるはずです.そもそも分布は正規分布なのかと思ってpubmedを漁ってみたら,こんな論文がありました.Myroslava Protsiv, et al, eLife 2020;9:e49555 DOI: 10.7554/eLife.49555 
    アメリカ人のデータですけど,この157年の間に体温が下がり続けているということとともに,Figure 1—figure supplement 1に分布が示されています.これによると,正規分布というよりは自然災害のリスク計算に用いられる極値分布のような裾が非対称な形状をしています.一般形は一般化極値分布(generalized extreme value distribution),属に GEVなどと呼ばれますが,体温の分布はその特殊系のガンベル型です.ガンベル分布などとも呼ばれますけど,ガンベル先生は私たちのよく知っているフィッシャー先生と一緒に極値理論を展開されたことで有名です.極値理論はリスク工学と密接に関係していて,エイヤでいうと,分布の裾の挙動に重きを置く統計理論です.乱暴ですが,平均よりも最大値を標本の代表値に据えていると言ってもいいかもしれません.リスク工学,例えば大風による被害や,河川の氾濫などの自然災害系の統計解析をされる方には馴染み深い分布かもしれません.
    JMPでは極値分布は「一変量の分布」プラットフォームに出てきます.冒頭の図は,「Iris.jmp」の「花弁の長さ」に対して「種類」の「setosa」を抜いたヒストグラムに「連続分布のあてはめ>極値」を適応したものです.因みにJMP15ではUIが少し違います.低値側の裾が少し厚くなっているのがお分かりでしょうか. この様子は箱ひげ図に赤い線で示された最頻区間と平均との関係を見てもよくわかりますね.もう一つJMPで極値分布を目にするのが「生存時間(パラメトリック)のあてはめ」プラットフォームです.因みに,分析メニューの「信頼性/生存時間分析」の下にあります.設定パネルの『分布』のところに「最小極値」と「最大極値」とがあって「対数正規」や「Weibull」の代わりに使ったりします.そういえば,信頼性分析でよく使われるWeibull分布も極値分布の一つですね.因みに,JSLでは,SEV Density(x, mu, sigma)関数が実装されています.
    私の関係する事例ではあまり極値分布が出てこないので,実際の事例に使ったことはありませんが,平均を代表値とする手法に限界を感じることが多くなってきたいので,極値理論に興味はあります.ポストコロナでは,適者生存のために今までとは違ったアプローチを取り入れていきたいので,少し勉強してみようかな.
    それでは.

  • 検査の損益に付いて考える

    先々週のブログで,PCR検査は検査というよりも,仮説検定と同じテストであるとみなすべきと書きました.もちろん,同じテストだと言っても,確率に基づいた判定ではないという点では,仮説検定とは厳密には違います.ではありますが,試行の結果(イベント)を真実と判断結果の四分割表に割り当てる手法であると緩く考えるとPCR検査の実態をより深く理解することができます.4つの領域は「あたり」と「はずれ」が半々ですが,このうちの2つの「はずれ」が意思決定手法としての評価に特に重要です.仮説検定では,それらを第一種の過誤と第二種の過誤と呼ぶわけですが,この命名では両者の違いがわかりにくいと思いませんか.英語でも,それぞれType IエラーとType IIエラー,あるいはαリスクとβリスクと言うので,分かりにくいという点では同じです.

    蛇足ですが,過誤はエラーを和訳したのだと思いますが,これは二つのあやまちの複合語です.どういう意味かと言うと,「誤ち」と「過ち」は共にあやまちですが,前者には早とちりの意味があります.高速道路の出口を誤って一つ手前で出てしまった,などと使います.一方,後者には見逃しの意味があります.あまり動詞としては使い分けないようですが,同じ例を出すと,高速道路の出口を過って降り損なった,ということになるでしょうか.第一種の過ちと書いてある参考書もありますが,厳密には間違いです.漢字を使うならば第一種の誤ちとすべきでしょうけれど,「あやまち」としておく方が無難かも知れませんね.少し脱線しました.

    日本では,アルファのaを「あわてもの」ベータのbを「ぼんやり」と紐付ける強引な覚え方もありますが,やはり,偽陽性や偽陰性という医療系で使われる用語の方が概念が明確になります.技術分野であまり使われないのは,判断の対象が陰陽と紐付けにくいものが多いからでしょうか.例えば,ある文字を「1」なのか「I」なのかを機械学習で判別することを考えると,その判断の結果は良くも悪くもないニュートラルなので,陰陽という概念と関連付けるのは困難です.

    PCR検査と仮説検定で使われる用語がそれぞれ異なるのは仕方ないとしても,両者のアナロジーをとる上で注意すべきことがあります.一つは,仮説検定では帰無仮説を棄却できなかった場合をPCR陰性と対応付ければ,仮説検定では対立仮説について何も言及しない(できない)のですが,PCR検査では陽性と判定されなければ陰性と判定されるという違いです.このことからも,PCR検査陰性の実態が理解できるのではないでしょうか.そもそも,いかなる検査でも100%陰性判定は下せないと考えるべきなのです.もう一つは,四分割表に記入される数値は,PCR検査では多数の被験者に対しての結果の割合であるのに対し,仮説検定では確率で記述することが多いということです.もちろん,両者は無限回の試行の元で一致しますが,多数回の検定を実施することはあまりないと思うので,割合と確率は違うものだということは,頭の片隅には置いておくと良いです.因みに,これが大数の法則の言わんとするところですね.

    次に,四分割表のそれぞれの領域の実態を考察していきます.
    a) 検査の結果陽性で,実際に感染していた.
    真陽性の場合で,検査が最もわかりやすい形で役に立ったケースです.インフルエンザのようにロシュ社のタミフルなどの治療薬があれば,積極的検査で早期に陽性確定する意義は大きいです.タミフルでは発症後48時間以内に服用しなければ効果がないと言われています.但し,どんな薬にも副作用は付き物なので,自然治癒に任せた方が結果としてハッピーだったというケースは全くないとは断言できません.

    b) 検査の結果は陽性だったけど,実際には感染していなかった.
    偽陽性の場合は,本人にとっては大きな不利益を被ります.無症状者であればなおのことです.二週間安静を強いられるのはまだ良いとしても,退院後に感染済みで免疫ができていると油断する人もいるだろうという懸念はあります.それと日本文化特有の要因から,プライベート面でいろいろな問題が生じることもあるようです.お店や個人などが嫌がらせを受けたり,言われのない差別をされたり,これは現在進行形で起こっていることです.いわゆる藪蛇とも違うのかもしれませんが,必要のない検査なら受けない方がハッピーだったということです.

    c) 検査陰性で実際に感染していなかった.
    真陰性の場合です.検査の結果が陰性ならば,発症していれば,少しは安心できるかもしれませんんが,症状がなければ,少なくとも今日1日は安心できるというだけのメリットとも言えます.以前のブログでも書いたように,PCR検査の信頼性を知っていれば手放しで安心はできないはずです.仮説検定と同じく陽性でないないことが示されただけで,陰性であることが示されたわけではないのです.感染していても未発症なだけかもしれないし,明日感染するかもしれません.発症者であっても,高熱が続く他の疾病もたくさんあるので高リスク者でなければ,むしろ不安に思うかもしれません.

    d) 検査陰性だけど実際には感染していた.
    この偽陰性の場合が一番怖いですね.発症者であれば行動が制限されるはずですが,無症状者であれば確実に感染を拡大させてしまいます.もっとも,高熱でバス旅行に出かけた老人などもいらっしゃいましたけど.このケースは個人にとっても家族の感染や本人の治療開始の遅れなど大きな損害を受けますが,特に社会に大きな損害をもたらします.

    上記考察から損益を数値で見積もります.因みに,非感染者が検査を受けることで感染する可能性や,感染者が検査を受けることで医療関係者に感染させてしまう可能性も考慮すべきですが,そこは今は考えていません.本来は,検査を受けるか否かの意思決定の下にPCR検査があるという二重構造を考えるべきです.法律学的に算定するのも面倒なので,それぞれのケースに平均的な個人の価値観を予想して,それをもとに利益と損を[a, b, c, d]の形式で書いてみます.それぞれのケースに値段をつけるというのも難しいところですが,いくらもらうなら(あるいは失うなら)釣り合うのかを考えてください.通常はこれらの損益を4×4行列で書き,これに混同行列をかけて期待利益を算出します.

    多くの人がa)の場合の価値を大きく見積もっていると予想しますが,バランスも考慮して仮に[10000, -10000, 100, -100000]としてみます.単位は円です.(因みに,自分の価値観ではd)の損はもっと多くしますね.)1%の有病率のもとでのPCR検査では,感染者が未発症の場合の感度が20%程度とのことなので,10000人中100人の感染者に対し陰性判定されるのが80人ですね.それから,特異度は99%としておきます.(おそらくもっと大きいとは思いますが,後述するように期待利益には大きな影響は出ないので.)とすると,残りの9900人中99人が誤まって陽性と判定されてしまうことになります.これらから.混同行列を求めます.因みに,発症当日でも偽陰性は80%なので0%にはなりません.こういう実態がようやく明らかになってきたようで,ここに良い解説があります.

    ということで,期待利益を冒頭の式から計算すると,端数を切り捨てると20-800+98-99=-781円になります.(この式では,利益と損益とが混ざってしまいましたが,修正するのも手間なのでご容赦ください.)になってしまいましたが,混同行列は有病率でも大きく変わりますし,上記では未発症者を対象としましたが,濃厚接触者で発症者であれば,更に大きく変わります.ですけれど,期待利益を計算するには些細な違いということがお分かりでしょうか.偽陰性率が決して0%にならないことを知れば,重要なのはむしろ個人の価値観です.仮に真陽性の利益を100000円とし,偽陽性の損を10000円と見積もると期待利益は119円になりますが,PCR検査の場合は個人の利益と社会の利益の双方を考えなければならないので状況は複雑です.基本は社会全体の期待利益を考えるべきです.そうなると,上述したように本来はPCR検査を受けることの損益も考慮した二重構造を考慮することになります.

    このような考察を書いてみましたが,決してPCR検査をするなということではないです.それが必要なときに直ちに結果を出せる検査体制は,クラスター対策には必須の要ですから.但し,クラスターと無関係の人にまで検査をする場合は,このようなことはどんな教科書にも書いてある考察は必要だと考えています.身近な検査に対しても同様な考察をすることで,イベントの理解が深まります.

    機械学習の分類判定の評価も勉強できますね.少し用語が違うので最初のうちは混乱するかも知れません.機械学習では感度は検出率(Recall)と呼ばれます.仮に「1」が正しい文字をどの程度「1」と分類できたかを示します.一方,紛らわしいのですが,「1」と分類した文字がどの程度本当に「1」だったのかを示すのが,精度(Precision)です.分類器の検査性能を精度で定義することもありますが,実務的には,検出率と精度の調和平均であるF値で分類器の性能とする場合が多いようです.調和平均なので二つの指標のバランスを重視していることになります.そういえばF値も色々な言葉で言われている言葉です.統一してもらえるとありがたいのですが.

    本日もお付き合い下さりありがとうございました.時間となりましたので,それではまた.


  • 検査とテスト

    今日は一日,統計リテラシーのオンライン用スライドを作成していました.オンライン向けにはアニメーションを使うのを避けたり,フォントサイズを小さめにしたりする以外にも,話の間をとるという工夫が必要です.セミナーのコンテンツは使い回す人も多いのですが,私の場合,基本的にスライドはその度に組み直します.今度はどんな相手にどんな話をしようか考えるのが楽しみの一つでもありますし,そもそも同じスライドを使っていると,自分が飽きてしまうのです.今回は,オンライン版として作り直すにあたり,せっかくなので最近の関心事からネタを拾って,コンテンツを更新しています.正直に申し上げると,コロナのことにも飽きてきた(こういうことを言うのは不謹慎ですが)ので,何か別のネタを探そうにも,未だ世間ではコロナばかりが目につきます.今週は,誰が作ったのか微陽性という変な用語まで登場してきましたね.閾値を超えたのなら陽性に違いはないのですけど,無理矢理,陰性と紐付けたいという何らかの意図があるのでしょうね.

    この変な用語にも意義はありました.この用語によって初めてあのPCR検査がアナログ的な信号処置に基づくものと知った方も多いのではないでしょうか.Real TimeのPCR検査の増幅曲線など見ていると,半導体検査などとは少し違ったものだと言うことがよくわかります.そもそも検査という日本語の訳がよくないです.英語では,半導体検査の検査はinspectionでPCR 検査の方はtestです.そう,あの仮設検定と同じtestなんです.WHOの事務局長が「検査,検査,検査」と言ったのが報道されましたが,オリジナルでは,We have a simple message for all countries: test, test, test. と言っています.

    PCR 検査と半導体検査とを比べると,擬似欠陥(第一種の過誤)や見逃し(第二種の過誤)があるということでは両者は同じですが,inspectionの方には,in(内部を)spect(見る)ion(こと)という語源からもわかるように,何かの中を(綿密に)窺うというニュアンスがあります.ですから,例えば,欠陥検査であれば,どのような種類の欠陥なのかを含めた質的,量的な結果を得ることが目的になります.因みに,誰かにinspectされるときは監査のように精査されるという意味にもなります.

    一方のtestは(厳密にはtestingとした方がいいのかもしれませんが),何らかの特性が一定の基準を超えたか否かを二値判定するものです.生徒を合格不合格に分類するのが,いわゆる学校のテストですね.半導体でもtestingは非常に重要な技術ですが,基本は動作する,しないの二値判定なのです.帰無仮説の棄却,採択による対立仮説の真偽判定としての仮設検定と同じですね.従って,ウイルス感染の陽性,陰性をある閾値で分類する手法としてのPCR検査はPCR テストと訳す方が相応しいのです.

    この語感の違いは素人でも何となくわかりますよね.PCR testがPCR検査ではなくPCR テストと訳されていれば,多くの人も一度や二度の結果で判断するのは危険だとすんなり受け入れてもらえたのではないでしょうか.検査という日本語は,テストよりも結果の信頼性が高いというイメージを誰もが持ってしまいがちです.実際は,PCR 検査では再陽性や偽陰性などの不確定な結果が珍しくないということは周知されている通りです.PCR 検査に対する認識の違いが,検査数が少ないとか,あるいは反対に不要であるとか,未だに様々な意見の衝突の火花を散らし,それをマスコミが煽って火を起こし,国民の間に分断を生じてしまったのは残念なことです.

    この感染症の最大の脅威は我々に知識の蓄積が少ないことにあります.単純にワクチンがないというだけではなく,あらゆることが未だに未知です.例えば,PCR検査の擬陽性率や擬陰性率などは,当初いくら調べても具体的な数値は見つけられませんでした.検査手法そのものの精度ではないです.PCR検査は超高感度だとは言っても,感染者からウイルスを採取できて初めて陽性になる可能性が生じるわけですから,
    超高感度故のラボ環境の汚染による偽陽性や検体採取の手技の稚拙による偽陰性まで考える必要があります.ダイヤモンド・プリンセス号ではとうに死滅したウイルスのDNAがPCR 検査で検出されましたし,インフルエンザでも,ベテランの医者と研修医では擬陰性率が大きく異なるという論文を読んだことがあります.それでも,類似の感染症の経験を踏まえて,それを元に行動できたことが,比較的マシと言える日本の状況をもたらしたのかもしれません.先人に感謝です.

    スペイン風邪の当時と違って,現代の我々には発展した医学だけでなく統計学もあるので,検査云々についてはロジカルな議論ができるはずです.ところが,当初は感度と特異度という検査の基本性能すらよくわかっていないことが障害なっていました.二つのリスクが不明ならば,その検査の結果に基づく意思決定はできません.そこで,このブログでも中国からの帰国者のデータを元に推察しました.その後,ダイヤモンド・プリンセス号をはじめとした多くの経験が蓄積してきました.ここにきて,ようやく信頼出来る値が出てきたようです.

    Variation in False-Negative Rate of Reverse Transcriptase Polymerase Chain Reaction–Based SARS-CoV-2 Tests by Time Since Exposure(RT-PCRによるSARS-CoV-2検査の偽陰性率の曝露後時間による変動)

    酵素ポリメラーゼ連鎖反応というのが,いわゆるPCRで,DNAを検出することしかできないPCRをRNA検出に適用するために開発された手法が,Reverse Transcriptase PCR 法で,RT-PCR法などととばれています.RNA を酵素による逆転写によって cDNA (cはcomplementary(相補的)の意)に変換してから,PCR 法を適用します.SARS-CoV-2も一本鎖プラス鎖RNAウイルスに分類されるRNAウイルスですから,実はPCR検査と呼ばれているのは厳密にはRT-PCR法なわけです.因みに,リアルタイムPCR(real-time PCR)というのもあって,人によってはこれをRT-PCRと呼んだりするので混乱を招いてます.リアルタイムPCR法では,DNAの複製過程における生成物の増幅蛍光をリアルタイムに測定することで定量化する手法です.このアナログ的な手法によって,野球選手の微陽性といういう判定が下されたのですが,解釈が間違っているように思います.

    計測の話はさておき,この論文のResultsをまとめると,階層ベイズモデリングを使って求めた偽陰性率は次のようになります.
    1日目 100% 95%CI(100%,100%)
    4日目 67% 95%CI(27%,94%)
    5日目(発症日) 中央値で38% 95%CI(18%,65%)
    8日目 20% 95%CI(12%,30%)
    9日目 21% (CI, 13% to 31%)
    21日目 66% (CI, 54% to 77%)

    冒頭にJMPでグラフを書いてみましたが,Figure2には毎日のデータが可視化されています.これによれば,発症前々日まではほぼ100%で,発症4日後に20%で底を打つようです.これは衝撃的なデータです.未発症者にPCR 検査しても,偽陰性100%って全く役に立たないという意味ですから.予想していたよりもPCR 検査って信頼できないですね.通常のインフルエンザでも言われていることですが,少なくとも発症前の検査は発症前日でも7割弱が見逃されてしまい,検査の精度が最も高くなるのが発症4日後とのことで,これはまだこのようなデータが,出ていなかった時からの指針と一致しています.この結果から以下のように結論するのは,おそらくほとんどの人が合理的と考えるのではないでしょうか.

    If clinical suspicion is high, infection should not be ruled out on the basis of RT-PCR alone, and the clinical and epidemiologic situation should be carefully considered.
    臨床的な疑義が大きいならば,PCR検査の結果のみで診断除外すべきではなく,臨床的および疫学的状況を注意深く検討する必要がある.

    ここの,Ruled out というのは専門用語で,識別診断の際,カルテにはR/Oなどとかかれたりしますが,その原因の可能性は低いという意味です.要するに,「発熱や咳があったら,検査陰性だったからと言って安心できないよ.もしかしたら却って手遅れになるかも.それに,動き回らないでね.できれば症状が治っても二週間くらい.」と言った感じでしょうか.この研究結果について,デザインバランスに問題があって推定は不正確であることが著者によって示唆されています.コメントにもいくつかの問題点が指摘されていて,まだ研究は続けていく必要がるようですが,何れにしても,スクリーニングポリシーの確立が急務だという認識では共通しているようです.

    統計リテラシーでは,「何のために検査をするのか?」を考える際には合理性がキーポイントになります.検査の結果と真実(感染の有無)の状況によって,4つの場合分けができます.(手元に『JMPではじめるデータサイエンス』がある方はp54に書いてあります.)PCR検査でもまさにこの状況にあります.このとき,何が起こるかを場合分けしてそれぞれの状況で何が起こるかを考えてみます.本日は,長くなってしまったので,覚えていればまた来週.

    それでは.


  • 見習うべきは台湾ではないか?

    Amzonでは,日常必需品の発送を優先しているとかで,先日,岩波新書編集部がTweetされてましたけど,書籍のような生活必需品でないものの在庫切れが目立ってきているようです.拙著も,現在両方とも在庫切れになっていて,『統計的問題解決入門』の方は,しばらく前から入荷時期未定になっていました.絶版になったわけではなく,実は,目出度く第3刷が増刷されたばかりなのです.流通が止まっている今は仕方ありませんね.連休中に勉強したいという方は,オーム社の直販を利用してくださいとのことです.第二刷では色々修正入れましたが,今回は大きな修正はありません.とはいえ,何点か変更した箇所もありますので,そのうちこの場でお知らせいたします.

    さて,今週は,先週の続きとしてAppleのMobility Trends Reportsを使って各国の比較をしてみたいと思います.


    最初に,データはAppleの特設サイトからCSVをダウンロードしてください.ここでは4月23日付けのデータを使います.先週と同じ手順でテータテーブルを作成してください.今回は,「国ごとの比較」という目的が明確なのでデータの層別化を先に実施します.具体的には「taransition_type_」の「walking」と 「geo_type_」の「country/region」のみ 取り出します.「transit」は欠けている国が多いので「walking」のみを対象としました.「driving」の方が外れ値が少ないのですが,国によって道路事情がかなり違うので,人間の活動の基本である「walking」を対象としました.何れにしても,「driving」もほぼ同じ挙動を示しているので問題はないと判断しています.更に,今回は国同士を比較したいので,「geo_type_」の「country/region」のみを抜き出します.

    そのためには「データフィルタ」で「geo_type_」と「taransition_type_」を選択して『追加』し,それぞれ「walking」「country/region」をcontrolキーを押しながら選択します.(Macではcommandキーを押します.)そして,サブセットとして別テーブルにすると63行のテーブルができるはずです.そこから「geo_type_」列と「taransition_type_」列は削除しても構いませんが,後日そのデータを見たときに,「geo_type_」と「taransition_type_」とが固定されていることがわかるように,「非表示,除外」しておく方が無難です.次に,先週と同じように「列の積み重ね」で全ての日にちの列を積み重ねます.「region」列には国の名前が入っているので,更にこれを基準として分割します.これで分析できるデーブルになりました.この次に何をやりますか?

    実は答えはありません.とにかく色々やるのが一番なのですが,自分の好きな「データ分析の周回コース」を確立しておくといいです.このデータでは,変数の正体がわかっているし,それぞれの分布はそれほど重要ではないので,いきなり多変量解析に着手してもいいと思います.因みにヒストグラムを描けば,個々の国のMobilityの度数がダブルピークとして観察できます.私の場合は,先ずは「多変量の相関」でしょうか.その上で「相関のクラスタリング」で大まかに多変量間の相関を見ます.

    どの国も活動量が時間経過とともに下がってきているという傾向は同じなので,この図のようにそれぞれはほぼ正の相関を示します.その中でも他と状況が異なっている国々が白っぽく表示されています.因みに,この図では「カラーテーマ」を「発散」の「寒色>暖色」に変更しています.完全相関の場合の色が「JMP標準』よりも濃いので,おそらく色弱の方にはこの方が見やすいはずです.

    試しに,「主成分分析」もやってみますが,やはり全体の傾向が同じなので,第一主成分のパワーが大きく(84.4%)なります.この図に示した4カ国が他と挙動が大多数の国と異なっています.日本は,対多数の中に埋もれてしまっていますが,青い矢印で示したように,第一,第二主成分の平面には乗らずやはり,少し異なった挙動と示しています.成分3が日本らしさになっているようですが,これが何を意味するのか..

    主成分分析では,何かありそうだと直観的に感じれば十分です.この違いを後続する分析で確認していきます.例えば「階層型クラスター分析」はとてもわかりやすい手法です.そのためには,今のテーブルを「転置」する必要があります.それには『転置する列』に全ての国を,『ラベル』に「日付」を割り当てて『OK』です.クラスター分析の設定は,すぐ分かると思いますので省略します.デフォルトではクラスターは7つに分割されますので,赤三角から「クラスターの保存」をしておきます.「クラスターの色分け」もしておきましょうか.

    非常に特徴的な国が見つかりましたか?因みにそれぞれのクラスターの平均の「Mobility」を折れ線グラフで描いたのが冒頭のグラフです.このグラフを参考にして,各クラスターの特徴を上げていきます.

    クラスター1(対岸の火事だったのが,飛び火して現在もロックダウン継続中)
    代表国はフランス,イタリア,この中には,ほぼ収束に向かっているニュージーランドもあります.

    クラスター2 (やはり対岸の火事と思ったか,えらい活動的だったのが数日の間に活動停止の事態に.現在もロックダウン継続中)
    代表国はスペイン.

    クラスター3 (早くから活動停止で延焼を阻止,現在も継続中)
    ここには,香港,韓国,シンガポールの3つの国しかありません.いずれの国もうまくやっていると評価されていますが,シンガポールでは外国人労働者の感染が増えて再度感染が数増えている状況です.何れにしても,うまくやってこれたのは早い段階から活動を抑えたことに加え,これらの国では日本ではできない人権侵害が功を奏したとも言えます.それを恐れた精神的なロックダウンの影響がこのクラスターのMobilityに出ているのではないでしょうか.日本の法体制のもとではこれらの国を手本にすることはできません.

    クラスター4(大丈夫なの?完全に活動停止中.)
    マカオは異彩を放っていますね.メディアがどこどこどこの国をお手本にすべきとか言ってますが,マカオこそコロナ対策では最もうまくいっている国です.おそらく,中国という国を良く知っているからこそ早い段階で対策ができたのだろうと思っています.ただ,人口も高々70万足らずの国で,産業構造も全く異なるので日本が手本とするのは難しいところです.

    クラスター5と6(対岸の火事が飛び火して活動を下げていったけど,徐々に活動開始に向かっている.)
    似ているので5と6はまとめましたが,代表国はクラスター6がロシア,オーストラリア,6がアメリカ,イギリス,ドイツといった国々です.イギリスは当初から集団免疫に言及していましたし,スウェーデンもこの戦略です.この先どうなっていきますでしょうか.

    クラスター7(何考えてんだ?)
    ここには日本と台湾しかないですが,基本的に活動量はあまり急激には下がっていません.台湾はうまくやっていると思っていたので意外でしたが,やはり中国という国を良く知っているので早い段階で独自の対策を立てられたのかなと.因みに日本と台湾を比べるとこんな感じです.

    日本は三連休中日をピークとするMobolityの増加が特徴的ですね.気が緩んだとか言われてますが,通常の週末とは変わりません.むしろ,その前後で周期が乱れていることに特徴があります.おそらく外出自粛要請を予期した,買い物に出かけた人々の影響でしょうか.日本は災害が多いこともあって,どうも私たちはアワアワしすぎのような気がします.この点,台湾は周期が乱れていないのが見事です.3月28日のピークが少しシフトしてますが,おそらく翌日の青年節(旧革命先烈紀念日)の影響でしょうか.私たち日本も台湾を見習ってもっと冷静に行動すべきではないか.そう思いましたけれど,皆さんならこのデータをどう考察しますか?機会があれば,同じデータを前にJMPユーザーでオンラインミーティングもしたいですね.

    もちろん,これは現時点での結果です.時間経過とともにまた別のクラスターを形成するはずです.今後の成り行きをデータをもとに注視していきたいと思います.

    それではまた.


  • 名古屋では外出自粛要請が無視されているのか?

    外出自粛要請の週末,いかがお過ごしでしょうか.自転車で遠乗りするには良い季節にはなってきましたが,ここ当面は控えて家に籠ることにしました.人と接触するわけではありませんが,事故でも起こせば医療リソースに負担をかけることになります.現時点では,最も避けるべき行動と判断しました.こんなときは,自宅でもできるデータ分析も趣味で良かったと思っています.最近では色々と不安な情報に振り回される毎日ですが,そんな中でもデータ分析の力を発揮して,少なくとも自分だけは正しい意思決定をしていこうと,そう考えています.皆様も,その結果をご家族や同僚など身近な人と共有して,できるだけ易しく説明してみることをお勧めします.その過程で,拙著でも強調した「視覚化」と「可視化」の違いが体感できるはずです.

    本日の記事も,基本的にJMP初心者を対象としますが,上級者の皆様とはグラフの背景を一緒に楽しんでまいりましょう.趣味のデータ分析はネタ探しから始めます.趣味ともなると電車内の広告を見てもネタになりますが,昨今ではとにかくネタには事欠きません.メディアリテラシーが身についてくると,データソースが定かでなければ,そもそも読むに値しない情報ということがわかってきますが,最近はソースが明確になっているものも多いようです.そこで,先ずは信頼できるメディアに掲載されているグラフを自分で描いてみることをお勧めします.例えば,最近こんな記事を目にしました.

    新型コロナの外出自粛要請を最も無視しているのは名古屋、アップルの移動データ解析サイトから判明

    データの出典もリンクされているので,それだけは良心的な記事と思います.

    このグラフをJMPで描いてみると,確かに名古屋のMobilityは他の都市よりも大きいようです.果たして名古屋では自粛要請が無視されているのでしょうか.このことを,JMPで確認してみます.データはAppleの特設サイトMobility Trends ReportsからCSVをダウンロードします.ここでは昨日ダウンロードしたデータを使いますが,このデータは毎日更新されているので皆さんは最新のデータをお使いください.それでは以下にグラフ作成の手順を書いていきますが,その前にどうすればいいのかを考えてみてください.

    1.ほとんどの場合,CSVは右クリックでJMPから開くのが一番早いです.特にこのような英語のデータは問題がないことが多いです.日本語のデータでは,文字化けするものもたまにありますが,その場合は,CSVを一度エクセルで開いて保存することで回避できることがあります.
    2.データがどのようにレイアウトされているかを眺めるのが次のステップです.それぞれの日付が列名になっていることがポイントですね.横軸を日付に取ってグラフを描くには,列として積み重ねたいところです.この段階で「データフィルタ」して日本の各都市のデータだけ抜き出してもいいのですが,せっかく世界各地データが記載されているので,ここでフィルタリングしてしまうのはもったいない.そこで,先にテーブルメニューから「列の積み重ね」を実施します.「2020-01-13」列以降のすべての列を『積み重ねる列』に指定して『OK』です.積み重ねない列は『すべて保持』になっていることは確認してください.
    3.「データ」列にはすべての国名や都市名ののMobilityが繋がっているので,これを分割するのですが,このデータでは注意が必要です.おそらくその地域の事情によるのだと思うのですが「taransportation_type」のうち「transit」がない地域があるのです.このデータでは,「region」の最初の「Albania」がそうです.「列の分割」は単純に列の先頭のパターンを見て,それを分割基準に採用するので,このままでは「transit」のデータが抜け落ちてしまいます.このような場合,先頭にダミーデータを付加したり,一度データを「taransportation_type」のそれぞれでフィルタリングして作成したテーブルを『結合(Join)』する手があります.「列の分割」で『基準となる列』に「region」の他に「transportation_type」を追加してもいいです.あるいは,同じ地域では「taransportation_type」間の相関が大きいので平均を採用してもいいかもしれませんし,どれか1つの「taransportation_type」を代表にとって考察してもいいでしょう.今の場合,特に日本では交通機関の利用状況に関心があるので,「transit」だけを抜き出して比較することにします.

    4.ということで「データフィルタ」で「transit」を選択して,赤三角から「サブセットの表示」して,このサブセットを「列の分割」します.設定が難しいとよく質問受けるのですが,何を分割したいのかを考えれば大丈夫です.今の場合Mobolityを国や都市毎に分割していくので,「データ」を『基準となる列』に,「region」を『分割する列』に割り当てて『OK』します.このとき「積み重ね」と違って「残りの列」は『すべて除去』がデフォルトになっていることに注意してください.『選択』にチェックを入れて,リストから「ラベル」を選択して,『選択されている列をすべて保持』をクリックしてから『OK』します.こうしてすべてのグラフ作成の元となるテーブルができたので,これを保存しておきます.
    5.それでは,日本の各都市のMobilityを比較してみます.「グラフビルダー 」でXゾーンに「ラベル」をドロップしてください.ここではデフォルトのままでいきますが,適宜変数名は変更して構いません.Yゾーンにはリストから日本の各都市を選択すればいいのですが,92列もあるのでリストの赤三角から「フィルタの表示」をするといいでしょう.この検索ボックスは文字の並びがキーとなるので「to」と打つと「Tokyo」だけでなく「Boston」なども一緒に検索されるので注意してください.「Tokyo」をYゾーンにドロップしたら,それ以外の都市を選択してYゾーンのグラフ内側の六角形の追加ゾーンにドロップしていきます.仕上げに『折れ線』ボタンをクリックすれば,冒頭に掲げたグラフになります.このグラフでは,X軸の「軸ラベル」を調整して文字が見えるように表示していますが,ラベルが間引かれていることには注意してください.

    確かに名古屋は,他と比べてMobilityが大きいようですが,そもそも「transit」が何を計測しているのか今一つ不明なので,何かを結論するのは危険です.例えば,東京と福岡では交通状況も違いますよね.東京では通勤に車を使うのはごく一部の人に限られていますし,バスでの移動はDriveになるのか,等々.いろいろな疑問は後日調べることにしましょうか.それよりも,自分で描いてみれば,このMobilityは1月13日を100%とした相対的な数値であることに気づくはずです.まだこの時期はコロナの気配は感じつつも自らの行動を制限するような状況ではなかったと思います.最初の日本での感染者が1月16日で,10人を超えたのが1月30日です.そこで,その翌日の1月31日を基準にしてこのグラフを描き直してみます.

    それには「東京」列を作成して,計算式 (:Tokyo / 128.14) * 100を入れます.以下同様にして作成したのがこの図です.

    随分印象が違いますよね.もちろん,どちらのグラフも基準が異なるだけで間違いではありません.重要なことは基準は意識的に選択できるということです.ですから,ここにデータ可視化の落とし穴があるのです.何よりも注意したいのは,私たちが意識せずに基準を選んでいる場合です.データを可視化したつもりで,正しい情報が得られないとだけでなく,こうして得られた間違った情報からの呪縛はとても強いので,この後に様々な弊害をもたらします.メディアの場合,意図的に間違ったグラフを描いて視聴者を誘導するというよりも,自分で自分に騙されて大騒ぎしているという場合もあります.更には,データの視覚化では,このように意図的に「名古屋は自粛無視している」という説得材料とすることも可能だということにも注意が必要です.

    多いこみで失敗することは誰にもあると思います.私の経験で言えば,特許を書いていて,最後の最後までロジックが破綻していることに気付かず,それまで明細書に費やした数日間がすべて無駄になったことも一回ではありません.最初の思いつきにとらわれることの怖さ,それを客観的にデータを見ることとで学べたと思っています.

    来週はこのデータを使って,別の分析をしてみようと考えています.

    それでは,また来週.

    日本や他の国はもちろんですが,ニューヨークの古き友人たちが早く安心して暮らせるようになりますように.


  • こういうときこそデータ分析しよう

    本日はこの記事に掲載されている図をJMP超入門として作成してみます.
    ここにきて「新型コロナ」の日本人感染者が爆発的に増えているワケ

    ここ暫くはあまり楽しくないネタになってしまいそうですが,どんなデータでもJMPで楽しくデータ分析してこの状況を乗り越えたいと思います.

    さて,この記事の要点は,コロナウィルスには二種類あって,そのうちの強毒株がヨーロッパで猛威を振るっているとの説を根拠に,3月24日前後にヨーロッパからの帰国者に感染が発見され,それをきっかけとして日本に感染が蔓延していったという論が展開されています.確かに図3を見る限りでは,ここに何かしらの要因があることが示唆されています.この方の言うように帰国者かもしれませんが,それを裏付けるデータは示されていません.その他にも,巷で噂されているように,ちょうどその日に決まった東京オリンピックの延期とも関係があるかもしれません.他にも,潜伏期間は1-14日で平均的には5日と言われているので,三連休で各地で賑わったことの影響もありそうです.それは脇に置いて,この図3をJMPで再現してみます.

    まずはデータから.データは,この前も使ったジャッグジャパンさんのcsvを使います.左上のcsvのリンクからダウンロードしてください.このcsvは整っているので,JMPで直接開いてOKです.さて,このデータから図3のグラフを書くにはどうしたらよいでしょうか.東京都の感染者数と日付の片対数グラフを書いて対数であてはめてみればいいのです.基本的に前回と同じ手順でいいのですが,少し考えてみてください.

    このテーブルのどの変数を使えばいいのかを最初に考えます.図3の検証が目的なので,日付けは「確定日」としておきます.「確定日」と「感染日」には数日のバイアスがあるので目的によっては注意が必要です.更に「居住都道府県」が東京である感染者数を東京都の感染者数とします.さて,このテーブルからどのようにすれば東京の累積感染者数の列を作れるでしょうか.テーブルメニューから「要約」すればいいのはすぐわかったと思います.「確定日」を『グループ化』に,そして「居住都道府県」を『サブグループ化』して『OK』です.

    「N(東京都)」列ができるので,列名を右クリックして「計算式列の新規作成>行>累積和」で「累積[N(東京都)]」列を作成してください.因みに,累積列は「統計」のところにある「Col Cumulative Sum」の計算式を使ってもいいです.興味があれば,他の都道府県も累積しておいてください.因みに「行数」列は日本全国の感染者数です.ここではこの全国と,押さえ込みが早い段階でうまくいったとされる北海道を累積しておきます.岩手県の方以外は,皆さんがお住いの都道府県を加えてみてください.

    データには外国居住者の列もできているはずです.このタイミングなので,おそらく海外からの帰国者と考えられます.多くはありませんが,合計すると昨日の時点で148名なので,これらの列を全て選択して,計算式列の新規作成の「組み合わせ>和」でまとめて,これも累計しておきます.「中華人民共和国」と「南アフリカ」を忘れないようにしてください.全部で11列あります.次に,あてはめ線の傾きに意味を持たせるために「確定日」を全国の感染者数が累計で100人を超えた2月21日を起点とした日にちを示す「Day」列を作成します.この処理は前回やったので覚えていますね.計算式は,Date Difference( 21Feb2020, :確定日, “Day” ) です.

    さて,この次にどうするかわかりますか?

    「列の積み重ね」で累計した列を積み重ねます.ここで初心者が躓くことが多いようです.グラフ不作成に迷ったら積み重ねる,と覚えておいてください.そして「二変量の関係」で「データ」を『Y』に「Day」を『X』に割り当てます.散布図のY軸のスケールを「対数」に変えます.これだけだとわかりにくいのでラベルを表示させます.そのためには,変数名を書き換えてもいいのですが,その変数の定義がわからなくなってしまうので「ラベル」列を再コード化した方が便利です.例えば,「累積[N(行数)]」を「日本」とします.「値の置換」にするのを忘れずに.それから,「ラベル」列と一番下に並ぶ最新日の行の全てを「ラベルあり」にします.散布図に現れるラベルはドラッグできますから,見やすい位置に移動させてください.

    いよいよ直線をあてはめます.片対数になっているので,赤三角から「その他のあてはめ」を選んで「自然対数」のラジオボタンをクリックしてください.このデータには0がありますからアラームが出ますが,それはスルーします.こうして得られたあてはめ線はデータ全体をあてはめていますので,「ローカルデータフィルタ」で「ラベル」を選択し,ついでに「Day」も『AND』して「Day」の上下限をそれぞれ0と9にして,あてはめ線の区間を10日間にしておきます.そうしたら,commandキー(すいませんWindowsはcontrolキーだったかな?)を押しながら「ラベル」の「東京都」をクリックします.

    はい,それでは「Day」のスライダーの真ん中を持って(手のひらのアイコンがちゃんと5本指になっているのを確認してから)左右にゆっくり動かしてみてください.あてはめ線の傾きが変わるのがよくわかるはずです.このとき,デフォルトではフィルタされた部分だけがズームされるので,今の場合は先に「Day」軸の下限を例えば-10にしておくと,スケールが固定されるというJMPの動作の特徴を覚えておくとよいでしょう.「表示」のチェックも外してください.こんな感じになればOKです.できましたか?

    因みに,冒頭に掲げたようにこのグラフはグラフビルダーでも描けますが,少なくともJMP14では対数の当てはめができないので「二変量の関係」を使いました.

    このグラフには冒頭のグラフと同じく,参考までにECDCからダウンロードしたイタリアのデータも載せておきました.ECDC (European Centre for Disease Prevention and Control:欧州疾病予防管理センター)は,スウェーデンのストックホルムに本部がある2005年にEUの伝染病の予防の専門機関として設立されました.このデータを使う場合は,「dateRep」の表示形式を『日付>m/d/y』に変えることに注意してください.必要な国をデータフィルタしてサブセット作ることをお勧めします.累積[cases]がその国の感染者数になります.後は,テーブルを「連結」すればいいのですが,このとき「一致しない行も含める」のところの「主テーブル」「結合するテーブル」の両方にチェックしてください.

    スライダーを動かしてイタリアの最新の状況を見てみると,
    Log(データ) = 10.246591 + 0.0334785*Day
    なので,爆発当初の
    Log(データ) = 3.1934324 + 0.4630768*Day
    と比べればかなり治まってきてますね.ありがたいことです.

    一方,東京のDay31から40までのあてはめ線は
    Log(データ) = 0.4262745 + 0.1494678*Day
    なので,イタリアの最新よりも傾きが大きいという状況です.これが今後イタリアの初期のようになっていくのか.それとも,北海道(最近また少し増えてきたのが心配ですが)のDay30から39まで
    Log(データ) = 4.7546329 + 0.0105227*Day
    の程度に落ち着いていくのか,見守っていきたいと思います.

    こういうと他人事のように聞こえるかもしれません.もちろん,やるべきことはやった上で,データによって客観的にリスクを見守ることが,デマに踊らされず,マスコミに煽られずに,私たちの精神を落ち着かせる一番の特効薬なのだ.そう考えています.

    それでは皆様もお大事に.本日はこれまで.