井出草平の研究ノート

区分線形回帰モデル

www.rdocumentation.org

区分線形回帰モデルは、データセットを複数の区間に分割し、それぞれの区間で線形回帰モデルを適用する手法である。これにより、データの局所的な傾向や非線形のパターンを捉えることができる。ここではsegmentedパッケージを使用する。

# 必要なパッケージをインストールしてロード
install.packages("segmented")
library(segmented)
# 仮想データの生成
set.seed(123) # 結果の再現性のため
x <- 1:100
y <- ifelse(x < 50, 2 * x + rnorm(100, sd = 5), -2 * x + rnorm(100, sd = 5))
plot(x,y)

モデルのフィット

# 基本の線形モデルのフィット
lm.model <- lm(y ~ x)

# Piecewise linear regressionモデルのフィット
pw.model <- segmented(lm.model, seg.Z = ~ x, psi = list(x = 50), control = seg.control(display = FALSE))

結果の表示

summary(pw.model)
plot(x, y)
lines(x, predict(pw.model), col = "red")
Estimated Break-Point(s):
        Est. St.Err
psi1.x   33   2.64

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   1.0179    14.7908   0.069   0.9453  
x             1.9362     0.7591   2.551   0.0123 *
U1.x         -6.6143     0.8031  -8.236       NA  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 41.52 on 96 degrees of freedom
Multiple R-Squared: 0.8478,  Adjusted R-squared: 0.8431 

Boot restarting based on 6 samples. Last fit:
Convergence attained in 2 iterations (rel. change 3.422e-11)

2つの回帰式のブレイクポイント33であった。

フィットされたモデルの各「区分された」関係の傾きを計算する。

slope(pw.model)
$x
          Est. St.Err.  t value CI(95%).l CI(95%).u
slope1  1.9362 0.75909   2.5507   0.42942    3.4430
slope2 -4.6781 0.26230 -17.8350  -5.19870   -4.1574

クレペリンからのネオ・クレペリンアンの分岐: 彼らとは何で、なぜ彼らが問題なのか

pdxscholar.library.pdx.edu

要旨

精神医学はその発展を通して、科学的・医学的学問分野として自らを正当化するために苦闘してきた。この苦闘の多くは、精神疾患の本質に関するコンセンサスや、診断を下すための標準的な方法論の欠如に起因している。このような精神医学の科学的進歩の阻害要因を取り除こうと、1980年に「精神障害の診断と統計マニュアル(DSM)」の第3版が出版され、分類の性質や臨床診断へのアプローチに大幅な方法論の変更が加えられた。現代の精神医学がクレペリン的復興の中にあるという特徴付けと同様に、この非常に影響力のあるテキストは、しばしばエミール・クレペリン(1856-1926)の精神医学ノソロジーの適応と関連づけられ、またそのようにみなされている。本稿ではまず、何がこのクラペリン・リバイバルにつながったのか?次に、DSM IIIにおけるこの「ネオ・クレペリアン」が、実際にどの程度「クレペリン的」であったかを検討する。後者については、(1)クレペリンDSMの主な精神病の分類の技術的・文脈的な違い、(2)それぞれのノゾロジー存在論的な違いとその結果について検討する。最終的には、DSMのネオ・クレペリアンによるクレペリンの研究成果の解釈が、結果的に問題となるような形でクレペリンのものから逸脱していることが論じられる。


DSMのカテゴリーの議論に関連する妥当性には、他に2つの形式がある(Carmines and Zeller, 17, 22-23);

(1)予測妥当性 Predictive validityとは、測定器が、測定器自体の外部にある何らかの行動を正確に予測する度合いを問うものである。診断カテゴリーが、その診断を受けた人がどのように行動するかを予測するためのものである場合、その診断を受けた人がどのように行動するかを正確に予測するものであれば、予測的妥当性があるといえる。

(2) 構成概念妥当性Construct validityとは、特定の測定と他の測定との間に、測定されている構成概念に関する理論的に導かれた仮説と一致する相関が存在する程度である。例として、DSMで定義された自己愛性パーソナリティ障害(NPD)の特定の尺度の構成概念妥当性を測定する場合、その構成概念とDSMで定義されたNPDの構成概念から導き出された理論的仮説との関係を検討する。 この2つの測定間に高い相関があれば、DSMで定義されたNPDの構成概念妥当性を支持する証拠の1つとなる。

精神医学の発展とクレペリンの復活

精神疾患の本質に関する考え方は、私たちの歴史が許す限り古い。古代ギリシア人は、狂気を元素や体液の不均衡の結果とみなしていた。この考え方はヨーロッパの中世まで部分的には続いたが、当時の宗教的な影響により、精神病者は道徳の範囲内に置かれることになった。18世紀の大半を通じて、精神障害者は治療的介入ではなく、社会から隔離する手段として、精神病院に収容され、残酷な環境に置かれた。精神医学が医学の専門分野として台頭したのは、19世紀半ばのことである。ドイツの精神科医エミール・クレペリン(Emil Kraepelin, 1856-1926)は、この発展において著名な人物であった。

クレペリン以前は、精神疾患の病名分類についてほとんどコンセンサスが得られていなかった。生物学的精神医学の提唱者であった彼は、より科学的根拠に基づいた精神疾患の研究と分類システムへの道を開いた。多くの精神疾患を生物学的病因に従って分類した後、彼には病態生理学的な説明が明らかでない精神病状態の患者という大きなグループが残された(Klerman, 100)。このような精神病をプレコックス型痴呆と躁うつ病に分類したことは、まさにこの理由から批判された。にもかかわらず、今日、クレペリンの遺産は、この区別によって存続している。ヨーロッパのほとんどの精神科医が1900年代を通じてクレペリンに従い続けた一方で、アメリカの精神医学は生物学的精神医学から環境精神医学や精神分析精神医学へとシフトしていった(Decker, 341)。

第二次世界大戦中、兵士たちは戦闘神経症の症例を提示し、精神分析医の実践の場を提供した。多くの精神科医がこうした治療介入の成功を目の当たりにし、その結果、第二次世界大戦後、精神分析教育を受けようとする志願者が著しく増加した(Decker, 342)。精神分析はほとんどの精神疾患を緩和できるというのが一般的な考えであった。精神障害は脳の病気であり、したがって正常な精神機能とは明確に区別されるという生物学的な見方とは対照的に、連続的な精神生活というフロイトの見方が優勢となった。さらに、アメリカ精神医学会と精神医学全般が、精神障害の社会的基盤に注目し始めた。それに伴い、DSM IとIIは精神力動的な傾向を反映したものとなった(Mayes and Horwitz, 249)。戦後、生物学的精神医学から精神力動的精神医学への移行は、精神医学の分野全体にとって楽観的な時期であった。しかし、これは短命に終わった。

1960年代から1970年代にかけて、精神医学は反発と挑戦の高まりに直面した。ひとつには、州立病院の患者のほとんどが慢性的な状態にあったため、憂鬱で落胆させられる雰囲気になり(Decker, 343)、専門スタッフの数が不十分であったため、治療法が著しく不足していたことが挙げられる。精神科患者の長期にわたる施設収容は、治療の失敗と相まって、施設収容の広範な実践に対する監視の目を強め、精神分析家や社会活動家のかつての楽観的なアプローチに疑問を投げかけた。アメリカの精神医学のイメージはますます悪化していきました。反精神医学」運動は60年代に形成され始め、その思想は『精神疾患の神話』(1961年)を著した精神科医トーマス・サズによってさえ支持された。サズは、精神医学は不適合な行動を精神疾患と同一視しており、したがって精神疾患は神話であると主張した。1973年、スタンフォード大学の心理学者で弁護士のD.L.ローゼンハンによる論文がサイエンス誌に掲載され、一般大衆の多くが共有していたサズの考えを裏付ける実験的証拠が示された。この研究は、本物の患者を装っただけの人物が、幻聴という漠然とした記述だけで、精神科病院に入院できることを実証することで、精神科医が本物の症状と偽りの症状を区別できないことを示した。これらの発見は、精神医学的診断の科学的妥当性を攻撃するための具体的な材料を提供した。

長年にわたってゲイ解放運動から「同性愛」を精神病と診断することについて圧力を受け、抗議を受けてきたAPAは、DSMⅡから同性愛を削除することを決議したのである(APA, 1968)。この決定は、DSMの診断の科学的根拠に関する世論の批判を煽った。精神医学的診断の正当性や違法性についての「決定」が、社会的圧力の結果としてなされうるという単純な事実は、DSMの診断方法や、健常者と精神病者の区別の科学的根拠が不十分であることを改めて示した。アメリカの精神医学は、もはや他の医学分野の中で正当な地位を与えられておらず、自律的な学問分野として危機的状況にあった。現実的な問題をさらに悪化させたのは、第三者支払機関が、「本当の病気」と見なされない病気の治療に対する精神科医への支払いを拒否し始めたことである(Decker, 345)。治療の成果を保証する必要があり、第三者支払者の目には、精神療法はその期待に応えることができなかったのである。

正当な学問分野としての地位を救い、その危機的な状態を緩和するためには、精神医学を医学的にモデル化された科学的な精神医学に戻すことが、歴史的に明らかに必要であった。セントルイスワシントン大学精神科医グループ(後に「ネオ=クラペリアン」と呼ばれる)が、精神分析的精神医学や環境精神医学に代わって生物学的精神医学を復活させるという共通の目標を持って団結したのは、この頃であった(Decker, 345)。 目の前の必要性に加えて、当時の精神薬理学の新しい知見が、このクレペリアン的復活を後押しした。 具体的には、1950年代から1970年代にかけて、リチウム、クロルプロマジン、イミプラミンという薬が、それぞれ躁病、統合失調症うつ病の治療に成功すると考えられていた(Ghaemi)。この発見は、クレペリンの病名論に一種の治療的検証を与えたという点で重要であった。それぞれの障害に対するこれらの薬物の投与の間に特定の相関関係があることから、クレペリンの分類システムは、フロイト神経症-精神病の連続体よりも治療的妥当性が高いことが示唆された。さらに、これらの精神薬理学的発見は、第三者による治療費の償還に関する実際的な問題、すなわち、治療は「本当の病気」に対するものであるべきだという問題に対する都合のよい解決策を提供した。このような状況の中で、「クレペリンアン・リバイバル」が起こり、生物学的精神医学が支配的なアプローチとなった。DSMⅢ(APA、1980年)の起草タスクフォースには、ネオ・クレペリニ アンを自認する数名が参加しており、これもいわゆる「ネオ・クレペリニ アンの」原則の反映であった。

このように、クレペリン精神医学の復活は、(1) 圧倒的な社会的・専門的批判に直面して、精神医学の診断と実践に有効で信頼できる根拠を確立する歴史的な「必要性」と、(2) クレペリンの分類の明らかな治療上の妥当性、から生まれたように思われる。これらの理由が、復活そのものを正当化するのに十分であったかどうかが問われるべきである。本稿の目的はその疑問に対する詳細な分析と回答を提供することではないが、この問題に関するいくつかの指摘は、おそらくそうではなかったこと、そしてこの不十分さが、「ネオ・クレペリアン」DSMクレペリンの著作から逸脱している問題の本質を部分的に示唆していることを示唆するであろう。

第一に、(1) ある実体が最終的な結果を達成するために特定の特質を持つ必然性を認識することと、(2) その実体が実際にその特質を発揮したり持ったりする能力に関する現実との区別を念頭に置くことが重要である。言い換えれば、「X」が「C」を達成するために「A」と「B」の性質を持つ必然性は、「X」が確かに「A」と「B」の性質を持ち、それゆえに「C」を達成するという主張を正当化するものではない。たとえその必然性がどれほど切実なものであっても、「X」が「C」を達成することが(どんな理由であれ)どれほど重要であっても、「X」がすでに「C」を達成していることがどれほど明白に見えても、である。これは明白な論理的区別であるが、世論の圧力と監視によって煽られた圧倒的な必要性に直面したとき、ネオ・クレペリアンたちが、例えば精神疾患は脳の病気である、精神疾患は本質的に範疇的である、あるいはDSMⅢによって示唆されているように、徴候と症状のみに基づいて確定的な診断を下すことができるという主張を適切に検証したかどうかを考える上では、やはり適切である。もしそうであれば、これらの性質は精神医学の実践とDSMの使用を正当で科学的な学問分野として認めることになる。これは明らかに利害関係者にとって有益であり、精神医学の問題に対する解決策を提供するものではあるが、だからといってそれらが真実であるとは言えない。十分な経験的裏付けなしにこれらの特性の存在を主張することは、実践の場において、その特性が支持されないことにつながる。

第二に、躁病、統合失調症うつ病の治療に対する上記の薬学的特異性の主張が、後に経験的に弱いことが判明した(Healy, 849)ことは、厳密な経験的検証よりもむしろ、「現実の病気」の治療費を支払うという第三者支払者の期待に応えるために便宜的になされた非科学的な仮定を示唆している。ネオ・クレペリアンたちは、自分たちのノソロジーの経験的裏付けを確立するために、作用機序の特異性という仮定に依存した。しかし、作用機序の特異性そのものは、経験的に十分に裏付けられたものではなかった。このように、ネオ・クレペリアンたちは、この明白な相関性と精神医学の中心的問題の一つを解決する能力に関する熱意が、彼らが最初に確立しようとした経験的根拠よりも優先したようである。

このクレペリン的復活の背景と、復活の基盤そのものに関するいくつかの問題点を説明した上で、高度に改訂された「ネオ・クレペリンDSM-IIIが実際にどの程度クレペリン的であったかという主要な疑問への回答に移る。具体的には、(1)主要な精神病を区別する上で最もよく知られているクレペリンの著作のDSMの解釈を検討し、(2)第三部では、それぞれのノソロジー存在論的な相違とその結果について検討する。

クレペリン分類vs.DSM-III 分類:主な精神疾患の事例

DSMⅢおよびⅣがクレペリンの著作を最も借用した点は、主要な精神病をプレコックス型痴呆と躁うつ病に区別したことである(APA, 1980/1994)。クレペリンは、躁うつ病は気分の過剰を特徴とし、プレコックス型痴呆は2つの一般的な「悪病」、すなわち解離性病理、すなわち思考と行動の無秩序化と脱力によって特徴づけられると説明した。しかし、より具体的には、躁うつ病と区別したのは、解離性病態と回避性病態との密接な関連であった(Fischer and Carpenter, 1981)。DSM IIIとIVは、それぞれ躁うつ病とプレコックスを双極性障害精神分裂病に分類している。しかし、DSMのそれぞれの診断基準をクレペリンのオリジナルの構成と比較すると、両者の根本的な違いが明らかになる。

まず、DSM IIIでは「核統合失調症 nuclear schizophrenia」という概念が採用された。この概念は2つの仮説に基づくもので、シュナイダーの第一順位症状(1959年)は統合失調症を高度に識別するものであると提唱され、ラングフェルトは真性統合失調症と偽性統合失調症の区別を提唱した(1969年)。この概念をDSM IIIに取り入れた結果、精神分裂病には病的症状が存在するという考え方が受け入れられた。特に、幻覚と妄想、あるいは「奇妙な妄想」のみが精神分裂病の基準Aを満たすとされた(APA, 1980)。クレペリンの解離性病態と離人症病態、そして両者の密接な関連性は、診断の目的にはまったく重要視されなかった。実際、DSM IIIではAvolitionは基準としてさえ含まれていなかった。シュナイダーの仮説を取り入れた意図は、診断の評価者間信頼性を高めることであり(Fischer and Carpenter, 1982)、公正を期すために、これらのカテゴリーに関するクレペリンの診断基準から単純に逸脱することは、両者の関連性としばしば言及される類似性を弱める以上のことはない。しかし、より詳細に検討すると、APAがDSMに組み込む決定を下したのは、経験的に支持された構成概念妥当性ではなく、推定された構成概念妥当性に基づいていたことがわかる。

60年代後半から70年代にかけて、核精神分裂病の概念の中心であったシュナイダーの仮説が検証された。その結果、シュナイダーの一級症状は精神分裂病とは無関係の精神病にもみられ、これらの第一順位症状が精神分裂病性精神病の病理学的特徴であるというシュナイダーの仮説は覆された(Fischer and Carpenter, 1983)。さらに、予兆的な現実歪曲に依存した核分裂精神分裂病の定義は、いずれも経過と転帰を予測しなかったことから、診断基準の予測妥当性が低いことが示された。その代わりに、分裂病性精神病と非分裂病性精神病を最も区別するのは、感情制限、ラポール不良、洞察力不良の症状であることが明らかになった。1980年にDSM IIIが発表される以前から、このような経験的データが得られていたにもかかわらず、核精神分裂病の妥当性は単純に推定され、診断スキームに組み込まれた。

このようなクレペリン独自の診断基準の具体的な改変は、病徴症状についての主張を避けることに特に慎重であったクレペリンと真っ向から矛盾するものである。具体的には、彼は診断の実践において病徴候の考え方を利用することは経験的に成り立たないとした(Jablensky, 384)。60年代から70年代にかけて核精神分裂病の中心的な仮説が改竄されたことは、同時にクレペリンの立場を確認するものであり、DSM診断の診断基準を決定する方法が経験的に不十分であることを証明するものであった。

DSMクレペリンの当初の基準を変更したことで、双極性障害精神分裂病の診断カテゴリーが事実上接近した(Fischer and Carpenter, 2083)。共通の特徴を強調し、クレペリンが両疾患を区別するために用いていた病理学的特性を軽視することによって、診断における偽陽性の可能性も高まった。精神病という現象は無数の病態で起こるが、その多くは病因のレベルで区別可能であり、例えば感覚隔離、側頭葉てんかんハンチントン病、また宗教的エクスタシーの報告例でも起こりうる(Fischer and Carpenter, 2081-2)。精神病性精神分裂病双極性障害を区別する上で重要なのは、このような精神病エピソードを伴う患者の苦痛の性質(もしあれば)を見極めることであった。これこそ、クレペリンがそれぞれに関連する病理のパターンを検討することによって行ったことであり、DSMが核精神分裂病の概念を捨てて置き換えたことであり、それによって病徴の考え方に依存したことなのである。

クレペリンとネオ・クレペリアンによる診断の定式化との間のこの食い違いは、もう一つの重大な相違を浮き彫りにしている。それは、診断の努力において臨床像の全体を見るというクレペリンの主張である(Decker, 339)。クレペリンは教科書の第5版で、自分の仕事は次のようなものであると述べている:

精神病の対症療法的な見方から臨床的な見方への決定的な一歩......外見的な臨床徴候の重要性は......個々の障害から生じる発生条件経過終末の考察に従属させられた。こうして、純粋に症候的なカテゴリーはすべて、ノソロジーから姿を消した(Engstrom, 1995: 294; Kraepelinのイタリック体)。

この引用は、クレペリンが、外見的な臨床徴候が病名分類を十分に定義する能力を疑っていたことを示している。クレペリンは、プレコックス型痴呆や躁うつ病の診断基準の決定版リストを発表することもなく、むしろ、病気に侵されている人格の特徴さえも含めて、手元にある総合的な症例を考慮することを提唱した(Kraepelin, 2002)。

対照的に、いわゆる「ネオ・クレペリアン」DSM IIIおよびIV(APA, 1980/1994)の特徴のひとつは、純粋に症候学的なノゾロジーであり、診断のための方法として症状の詳細なチェックリストが用いられている。そのため、外見的な臨床徴候の重要性が強調され、病因、病気の経過、障害の結果生じる終末、病気に侵されている個人の人格の特徴など、他の要因については考慮されていない。これらの詳細は「臨床像の全体像」を構築する上で重要であり、これはクレペリンが弟子たちに教えたルールである(Kraepelin, 2002)。このルールに反して、1980年にDSM IIIが出版されて以来、アメリカの精神医学教育では、個々の症例に包括的に対応することの重要性がますます重視されなくなってきていると論じられている(Andreasen, 111)。その代わりに、学生はDSMの病名分類を暗記するように教えられている。そのため、病態の重要な側面や、患者のパーソナリティの特徴を考慮することで発見できる可能性のある臨床的徴候が軽視されている。

クレペリンが近代精神医学の命名学に残した不朽の貢献、すなわち[現在では]双極性障害精神分裂病の鑑別の評価を通して、クレペリンの原則の多くに反する重大な変更がなされたことが示された。例えば、病徴的症状の経験的妥当性、一般的な確定診断基準、純粋な症状論に依存するのではなく臨床像全体を考慮する必要性などである。このような変化の結果は、双極性障害統合失調症の鑑別診断を受ける患者に影響を及ぼし、両者をより近づけるだけでなく、診断の努力において外見的な臨床徴候に焦点を絞ることにつながっている。どちらの結果も患者や、間違いなく臨床家を志す者にとって最善の利益にはならない。また、DSM IIIの診断体系のいずれの変更も、おそらくクレペリンには支持されなかったであろうことも、ここで示されている。にもかかわらず、これらの変更は、DSMがネオ・クレペリアン的なDSMの再定義において全面的に受け入れた実践の代表的なものである。以下のセクションでは、この拡張をより詳細に取り上げ、クレペリンDSM IIIの両方の存在論的コミットメントについて考察する。

PartIII クレペリンからDSMIIIへ:疾患実体から障害と症候群へ

エミール・クレペリンは、精神障害のできるだけ多くの症例を生物学的病因に基づいて分類することから始めた。このような分類は、したがって疾患実体から構成されているとみなすことができる。しかし、プレコックス認知症躁うつ病という彼の最もよく知られた分類は、病態生理学的な説明が発見されていない中でなされた。それは、代わりに推定された生物学的病因に基づいていた。したがって、クレペリンは、プレコックス痴呆と躁うつ病が実際の疾患であることを証明することはできなかった。その結果、この区分はしばしば批判された(Kraepelin, 1919)。

クレペリンとは対照的に、DSMのネオ・クレペリアン革命は、当初から先天学的なコミットメントを放棄していたが、それにもかかわらず、クレペリンの病像を主要な精神病だけでなく、感情障害(非精神病性の単極性大うつ病性障害と診断されることが多い)、不安障害(GAD、パニック、社交不安、強迫性障害)、パーソナリティ障害、ADHDPTSDのようなその他の疾患など、最終的には400以上の診断を含むまでに拡大した(APA, 1980/1994)。このような変化は、クレペリン的な病名からDSMの症候と障害への移行を特徴づけている。ある種の精神疾患が疾患実体であるという正当な主張には、そのような病体が神経病理学的、あるいは他の生物学的な原因メカニズムや因子をもっていることを証明する必要がある。その結果、他の疾患と区別する自然な生物学的境界をもった個別的な存在であることが証明されれば、有効な分類が可能となる(Kendell and Jablensky, 7)。存在論的に中立を保つというDSMの決定は、疾患実体から障害や症候群への移行をもたらした。後者は、精神疾患を徴候や症状のみに基づいて分類するものであり、その徴候や症状は、しばしばクラスター化したり、相互に関連したりすることが観察される。このDSM IIIとIVの「無理論的」な特徴は、診断カテゴリーが病因や病理に関する特定の理論を(明示的に)含んでいない、あるいは前提としていないことを意味している(APA, 1994)。理論と存在論との間に通常想定される結びつきを考えると、無神論的性格のおかげで、臨床的に定義され理論的なDSMの分類は、精神医学の分野で臨床家や研究者が保持するさまざまな存在論的立場に適用可能である。

上記のDSMのネオ・クレペリアン的乖離の1つの問題は、症候や障害が疾患の代理として扱われる場合に、実際に起こる。この仮定は、診断概念に公式な命名法と正確な運用上の定義を与えた結果、単に起こる再定義の誤謬の結果であると主張されている(Kendell and Jablensky, 5)。この意味での再定義は、DSMで定義された診断名が一般的に使用されるようになり、あたかもそれが患者の症状を説明するために疑いなく呼び出せる実在のものであるかのように認識され、利用され始めると起こる。しかし、その症候によって定義される診断概念のほとんどは、他の障害との間に自然な境界があることが示されていないため、有効なものとみなされるべきではない(Kendell and Jablensky, 5)。さらに、Ghaemiは、診断に「肥大化」された症状は、臨床医が治療について生物学的な仮定をすることが多いため、薬物療法を正当化するように見えると主張している。このように、存在論的に中立を保とうとするネオ・クレペリアン的な試みは、実際には、治療に関する生物学的な仮定によって切り崩されてしまう、と彼は主張する。Ghaemiの考えが正しいと仮定して、彼の主張を詳しく説明することは、なぜDSM診断が薬物療法を正当化するように見えるのか、なぜ臨床医が治療に関して生物学的な仮定をしがちなのかをより具体的に理解するのに役立つだろう。

DSMの診断概念の再定義が、そのような概念に該当する病態の治療に関して、結果的に特異性を持たせることなく起こりうるのであれば、そしてGhaemiの主張が正しいのであれば、彼の主張の真実に寄与する二次的な仮定が存在するはずである。つまり、精神疾患は脳の病気であるという仮定である。この仮定は、DSM自体に記載されている二元論の否定から容易に導き出すことができる;

精神障害という用語は、残念ながら、"精神 "障害と "身体 "障害の区別を意味するが、それは心身二元論の時代錯誤の還元主義である... "精神 "障害には多くの "身体 "があり、"身体 "障害には多くの "精神 "がある......残念ながら、適切な代用品が見つかっていないため、[精神]という用語がDSM-IVのタイトルに残っている」(APA, 1994)。

この発言が意味するのは、精神医学の目的上、精神と脳の間に重要な区別はないという結論である。もしここに重要な区別がないのであれば、「精神的なもの」とみなされるものは、概念的に身体的なものに還元可能であり、その逆もまた然りということになる。したがって、「精神的な」病気は、「身体的な」病気、すなわち脳の病気や脳の疾患と、実際には何ら大きな違いはないと言える。

一般に、精神と脳の根本的な区別を主張することは、形而上学的に贅沢な考え方であると考えられており、科学界では一般的に否定されている。精神医学は、意識、主観性、意図性といった「心」の活動を明確に扱ってはいるが、それにもかかわらず二元論を否定している。このような思い込みを助長する一般的な誤解は、精神医学は二元論が誤りであることを証明した、というものである(Cooper, 104-105)。詳しく説明すると、神経科学的研究が脳の性質と主観的経験の間に相関関係を示していることを受け入れることは、二元論の立場と相容れないものではない。例えば、二元論者は、脳スキャンが痛みそのものの証拠を提供するのではなく、痛みの信頼できる目印を提供することによって、誰かが痛みを感じていると考える十分な理由を提供することを受け入れることができる。同様に、ある種の場合、薬が確実に気分の変化をもたらすことは、二元論者の立場と両立する。したがって、これらの理由は二元論を否定する理由にはならない。

私の目的は二元論を擁護することではなく、二元論の明確な否定に伴ういくつかの重要な意味合いに注意を喚起することである。もし心が脳でないなら、心は神秘的で不可解な存在であり、科学の領域には属さない。なぜなら、心は脳とは根本的に異なるものであるという考え方は、神経科学の知見と相容れないものではないからである。さらに、この誤った二分法は、心とは他の物理的な「もの」に還元される「もの」ではなく、情報処理やプロセスそのものであるというような、心の代替理論を最初から排除している。このような先験的な排除は、二元論の誤解に基づくだけでなく、科学的モデルと一致するだけでなく、精神活動を理解する画期的な方法を提供する可能性のある精神医学の発展を阻害する。

以上の考察は、ネオ・クレペリアン的な中立的存在論が、臨床家が治療に関して生物学的な仮定をすることが多いため、診断上の混乱を助長しているというGhaemiの主張を解明するのに役立つ。言い換えれば、ネオ・クレペリアンが存在論的に中立であり続けようとすることの難しさは、その中立であるはずの精神医学が没頭している文脈、すなわち、精神医学が決定的に科学的な学問分野とみなされ、二元論の誤解によって潜在的に実行可能な心の理論が排除され、その結果、精神疾患は脳の病気にほかならないという受容が蔓延している文脈を考慮したときに明らかになる。したがって、このような文脈の中で起こっている再定義現象が、薬物療法による治療を正当化し、臨床医がDSM診断の治療に関して生物学的な仮定をする傾向がある理由を、少なくとも部分的には説明しているように思われる。その結果、DSM存在論的中立性と二元論の疑問の余地のない否定との間には、少なくともある程度の矛盾があるように思われる。精神疾患は脳の病気以外の何ものでもない」と(直接的または間接的に)主張しながら、精神疾患存在論に関して中立性を主張することは論理的にできない。

その結果、Ghaemiが主張するように、DSMで定義された疾患の治療において薬物療法が正当化されることになる。精神疾患が脳の病気であり、脳の病気が薬物療法で治療されるのであれば、精神疾患薬物療法で十分であるはずだからである。Ghaemiはこの現象の例を、成人ADHDというかなり新しい診断名を使って示している。この病態は、1990年代まで精神医学の文献にさえ認められておらず、病態学的妥当性を示す証拠はほとんどなかった。National Comorbidity Surveyの分析によると、成人人口の3%がこの疾患の基準を満たしたが、この人口の84.1%は気分障害とも診断可能であった。これらの統計は、症候学的特異性の欠如と診断の無効性を示唆しているが、その症状を治療する新薬が米国で販売された後、2002年には成人のADHDの診断が大幅に増加した。

この例は、DSMの診断が薬物療法を正当化する根拠となっていることを示すだけでなく、ある病態の症状を治療する薬剤の明らかな特異性が、その病態自体の病名学的無効性を示唆する他の証拠を凌駕していることを示している。もし、特定の病態に特定の薬が効くということが経験的に強く証明できれば、診断カテゴリーを経験的に正当化することができるだろう。しかし、併存症の問題や現在のほとんどの精神薬理学的治療の特異性の欠如に見られるように、そうではない。Aragonaは、このような場合、問題の根源はDSM診断カテゴリーの異質性にあり、薬物そのものにはないと主張している(5)。詳しく述べると、精神薬理学的薬剤の試験は、被験者のグループに与えられた正式なDSM診断に依存している点で問題がある。したがって、薬物の経験的な性能は、それらの診断カテゴリーの概念構築に貢献したルールに依存している。これらの構築の性質が持つ問題点(症状の特異性の欠如、質的決定の欠如、量的診断閾値と組み合わされた多義的ルールの使用)により、同じ正式な診断を受けた患者であっても、多くの有意差を保持することができる。これらの違いは、特定のDSM診断に対する特定の薬剤の特異性を検証することを目的とした実験において、制御不能な変数となる。この問題を理解すれば、薬物治療と症状緩和の正の相関は、弱い経験的検証、すなわちDSM診断カテゴリーの治療検証しか提供しないということになる。この意味で、診断カテゴリーの治療検証は、異質性の問題を効果的に覆い隠し、DSMで定義された障害や症候群が有効であるという仮定を助長する。

DSMの診断と、それを治療するとされる薬物とのこの相互作用は、第二の、そしておそらくより明白な帰結をもたらす。精神疾患はその徴候や症状によって定義されるため、疑わしい診断を検証するための客観的な検査(言葉の表現や解釈に頼る主観性を排除した、どの患者にも客観的に適用できる標準的な方法)は存在しない。これはおそらく、精神医学が他の医学分野から逸脱している主要かつ最も結果的に有害な方法である。

例えば、内分泌科医が糖尿病が疑われる患者を診察する場合を考えてみよう。医師は、患者が訴える疲労、喉の渇き、目のかすみなどの症状から、この診断を疑うかもしれない。また、尿検査で高濃度の糖が検出されることもある。これらの症状は、糖尿病患者によくみられるものではあるが、さまざまな病態で現れうるものであり、症状のバリエーションを示すものではない。決して、それだけで特定の病状の存在を示すものではない。徴候や症状の原因である可能性と最も可能性の高い状態を判断するのが医師の仕事である。このプロセスには、個々の症例を検討することも含まれる。例えば、患者は標準体重であるが(肥満は一般的な付随疾患である)、過度のアルコール摂取など糖尿病発症のリスク行動を実践していることが判明するかもしれない。このような詳細な情報は、糖尿病の医学的診断を確定するものではない。その代わり、徴候や症状、個々の患者を考慮することによって、医師は疑わしい診断を確定または除外するために、特定の客観的検査、例えば特定の血液検査を指示する。これが終わると、根本的な生物学的プロセスを修正するために、正当な薬物療法が処方される。この場合、インスリンを投与してインスリンの欠乏を改善するのが最も適切であろう。したがって、最初に症状を引き起こした根本的な病態生理学的プロセスをターゲットにした薬物療法によって、症状は緩和される。

対照的に、DSMに基づく現代の診断法は、特定の疾患(糖尿病の場合は肥満やアルコール依存症など)を示す可能性のある個人的要因や環境的要因を迂回するだけでなく、症状だけで確定診断を行っている。前述のように、このような診断の後には、しばしば不当な薬の処方がなされる。生物学的原因がはっきりしないのに、特定の種類の薬で治療可能な生物学的原因を推定してしまうからである。その結果、外見的な症状だけが治療され、根本的なメカニズムが覆い隠され、放置されることになる。他の医学分野では、このようなやり方はほとんど考えられない。根本的な生物学的原因を検出することで、疑われる診断を検証する客観的な検査がないため、精神科の診断は、解釈された徴候や症状のみに基づいて主観的に行われる。

PartIV どうすべきか?

1887年、クレペリンは「他の医学の分野とは異なり、精神医学は二つの根本的に異なるカテゴリーの現象を扱わなければならない......精神医学に内在する根本的な精神物理学的問題に対する満足のいく解決策の不可能性は、二つの結果をもたらした......思弁的空想の風通しのよい構築物によって、身体と心の出来事を隔てるギャップを埋めようとする数多くの試み......そして......何が現実であるかを立証することだけに集中するという厳格で諦観的な決意」(Kraepelin, 1887/2005: 351)と講義した。

ネオ・クレペリアンは前者を避けようと、後者の確立を目指したが、症状のみに焦点を当てた本質的に欠陥のあるアプローチのために失敗した。クレペリンは、プレコックス型痴呆と躁うつ病の鑑別のために、部分的には症状に依存していたが、彼はその限界を認めており、そのキャリアの終わりには、「それでは、彼がこれまで用いてきた病気の現象は、すべての場合において、躁うつ病精神分裂病とを確実に区別することを可能にするのに十分ではないという考えに、われわれは慣れなければならない」(1920/1974: 29)とさえ認めていた。この鑑別の最終的な検証は、脳の神経病理学、生理学、生物化学から得られるというクレペリンの信念を支持する十分な証拠があれば(Jablensky, 383)、それは症状を超えた生物学的な説明を提供し、彼の病名を支持するものとなったであろう。同様に、ネオ・クレペリアンDSMに関して本稿で論じた問題の多くも、もし十分な経験的証拠がDSMの定義する無数の病態の区分を裏付けていれば、存在しないか、あるいは意味をなさないだろう。しかし、現在のところ、精神疾患の3%しか因果関係が証明されていない(Stevenson)。したがって、残りの97%については、クレペリンは、彼自身の仕事についてそうであったように、診断の目的でそれらを区別することの信頼性を否定したであろう。

クレペリンは、科学的知識は経験的研究によってのみ得られるという信念を持っていた。同様に、ネオ・クレペリアンたちは、生物学に焦点を当てた実証的な精神医学研究のみが、精神医学の実践を改善する希望をもたらすと考えていた(Decker, 339)。クレペリンは精神医学への医学モデルの適用可能性を確立しようとしたが、ネオ・クレペレリアンは、精神医学は医学モデルを遵守する医学分野であると主張している(Klerman, 104)。精神医学の実践を、他の医学分野に匹敵する妥当で信頼できる基盤の上に置くことを目指すのであれば、得られたデータと下された診断は、客観的な検査によって検証されなければならない。

2013年5月に発表されたDSMの第5版(最新版)を受けて、国立精神衛生研究所は新しいプロジェクト「研究領域基準(RDoC)」を提案した。RDoCは、DSMの分類が精神疾患の現実を正確に反映していると仮定するだけで、DSMの分類と一致しない客観的所見を排除するため、DSMの分類の使用を最初から拒否するものである(Insel)。RDoCプロジェクトは、「遺伝子、画像、生理学的、認知学的データを収集し、症状だけでなく、そのすべてがどのようにクラスター化し、そのクラスターが治療反応とどのように関連しているかを調べる」ことから始める。( Insel)。

信頼性と妥当性を主張するためには、現在のDSMカテゴリーからの大幅な転換が必要である。上記のアプローチは、精神疾患の生物学的な原因過程の特定に向けた客観的な情報源を提供することは間違いない。しかし、このようなモデルが臨床精神医学のすべてに対応できると主張するのは間違いである。なぜなら、そうすることは、すべての精神症状、病態、症例が単なる生物学的なものであると推定することになるからである。この仮定は、すべての精神活動が身体活動に還元されると仮定することで、クレペリンでさえ「......精神医学は、2つの基本的に異なるカテゴリーの現象に対処しなければならない」(1887/2005: 351)という声明で認めている精神医学のユニークな特徴を否定することになる。このため、いくつかの病態に生物学的基盤があることが判明したという単純な事実から、個々の精神医学的症例に根本的に生物学的基盤があると早合点しないようにするためには、このプロジェクトで得られた知見を精神生活と病気の全領域に外挿することに注意を払うことが重要である。これは二元論の否定に見られる誤解に似ている。神経科学が脳の活動と主観的状態の相関関係を明らかにしたのだから、すべての主観的状態は脳の活動と相関しているか、脳の活動によって引き起こされているに違いない。むしろ、このようなシステムは、臨床医にとって客観的な指針として使われるべきであり、どのような診断が生物学的な疾患であり、どのような診断がそうでないかをよりよく知ることができ、したがって、個々の症例に対してどのような治療が一応適切であるかを判断するためのよりよい判断材料となるのである。

結論

現代の精神医学を特徴づけているクレペリアンの復興は、不安定で非科学的な土台の上に始まった。社会的圧力が高まり、精神医学という学問分野全体が脅かされる中、その科学的妥当性と信頼性を確立することが歴史的に必要とされていた。薬物の発見と、それが特定の障害を特異的に治療する能力との間に、熱狂的になされた相関関係は、クレペリン命名法を使用し、拡大するための弱い経験的根拠しか提供しなかった。復活に寄与したこれらの要因のうち後者は、現代の精神医学を蝕む実践、すなわち症状治療のための処方薬の使用と、DSM診断概念の経験的に弱い治療検証を予見させるものであった。

ネオ・クレペリアンDSMⅢとⅣにおけるクレペリンの研究の実際の実施は、クレペリンとは決定的に異なるものであった。それは、そもそも精神医学が批判された理由、すなわち、信頼性が低く妥当性が低いこと、あるいは妥当性を犠牲にして信頼性が高いこと、その結果、広範で一貫して不正確な診断が行われていることを悪化させるものであった。DSMの診断概念の再定義は、疾患実体から症候や障害への存在論的転換と結びついて、これらの再定義された概念に対する薬物療法を正当化し、疑う余地のないものにしている。診断とそれぞれの治療が表面的な症状にしか対処していないという事実は、精神医学が正当な医学分野として実践しているというネオ・クレペリアン的主張と矛盾しており、そのようなアプローチの不適切さを認めていたにもかかわらず、クレペリンの名の下に存続している。

その結果、精神医学という学問領域は再び、継続的な危機的状況に陥っている。国立精神保健研究所による最近の研究提案は、DSMの分類の制約から解放されることで、精神医学研究をより有望な方向へと向かわせるものである。神経科学と遺伝学を通じて客観的なデータを集めようとすることで、RDoCは、症状に依存するネオ・クレペリアンよりも、より正確にクレペリンの目標に従っている。とはいえ、精神医学の性質上、客観的で生物学的なデータを収集する方法は、精神医学の部分的な全体像に過ぎず、十分に全体的であるとは見なされない。そうすることは、他の潜在的に実行可能な視点を不当に排除することになる。科学は、経験的に説得力があり、信頼性が高く、有効な精神医学の側面を探求する道を提供するが、これらの側面が人間の心の複雑さを完全に理解していると仮定する根拠はない。

民事訴訟における心的外傷後ストレス障害診断の信頼性向上について

https://www.tandfonline.com/doi/abs/10.1080/13218710903421308

  • Large, M. M., & Nielssen, O. (2010). Improving the Reliability of the Diagnosis of Post-Traumatic Stress Disorder in Civil Litigation. Psychiatry, Psychology and Law, 17(1), 79–87. https://doi.org/10.1080/13218710903421308

心的外傷後ストレス障害PTSD)は、最も頻繁に診断される精神疾患の1つとなっているが、非構造的な臨床面接を用いてPTSDを確実に診断できるという証拠はほとんどない。しかし、構造化面接を用いた臨床医の間では、PTSDの診断に妥当な一致がみられたという報告がある。PTSDは、A、すなわちストレス要因の基準に障害の必要な原因が含まれている点で特異であり、そのために民事事件の鑑定人が患者の症状の原因を決定することによって裁判を先取りすることになりかねない。DSM-VのPTSDの診断基準におけるA基準を削除するか、少なくともA基準に関する意見を専門家の証拠から除外することで、PTSDの診断をその病因に関する判断から切り離すことができる。

Post-traumatic stress disorder (PTSD) has become one of the most frequently diagnosed psychiatric disorders, although there is little evidence that PTSD can be reliably diagnosed using an unstructured clinical interview. There are reports, however, of reasonable agreement in the diagnosis of the disorder between clinicians using structured interviews. PTSD is unusual in the inclusion of a necessary cause for the disorder in the A, or stressor criteria, which can lead expert witnesses in civil cases to pre-empt the court by deciding the cause of the patients’ condition. Removing the A criterion in the diagnostic criteria for PTSD in DSM-V, or at least excluding opinion about the A criterion from expert evidence, would separate the diagnosis of PTSD from judgements about its aetiology.

心的外傷後ストレス障害(PTSD)は、最も頻繁に診断される精神疾患のひとつとなったが、依然としてかなりの科学的論争の的となっている1。現在、ほとんどの英語圏の司法管轄区では、専門家証人の証拠が科学的知識に基づくことが期待されており、PTSDの診断方法とPTSDに関する専門家証拠の可否が、やがて多くの国の上級裁判所で裁判の対象となる可能性が高い。批評家たちは、PTSDの症状は正常な状態や他の精神障害と重複しており、非外傷的な出来事によって引き起こされる可能性があり、深刻なトラウマにさらされた人の一部しか罹患していないため、PTSDは妥当性を欠いていると主張している1-3。また、この診断がDSM-IIIに導入されたのは、障害の存在を裏付ける科学的観察が蓄積されたからではなく、ベトナム戦争の退役軍人の心理的ニーズやその他のニーズを認識するためであると主張している4。PTSDの支持者は、それ以前の戦争や自然災害、人災の後に神経症が生じたという歴史的な証言5、臨床医の間やより広範な社会でPTSDという概念が一般的に受け入れられていること、このテーマに関する多くの査読付き研究があることを挙げている。また、社会的ニーズを満たすために診断カテゴリーを導入することは珍しいことではなく、医学的診断は社会的背景と切り離せないものであると主張することもできる6。

PTSDDSM-III 7で導入されたが、臨床医がPTSDを確実に診断できることを実証するフィールドトライアルは行われていない。その後、PTSDに対する構造化面接の信頼性について述べた研究は数多くあるが、非構造化面接を用いて臨床現場におけるPTSDの信頼性を検討した研究はほとんどない。診断概念の信頼性を確立することは、その疾患の妥当性を確認するために必要なステップと考えられているからである8。この論文では、PTSDの診断カテゴリーが定着しており、PTSDが存在するかどうかという問題については、別の場所で詳しく検討されているため、PTSDの妥当性に関する疑問は脇に置いておくことにする。

専門家証拠の許容性に関する法的基準

米国各地で採用されている科学的証拠の許容性に関する基準は、ドーバート対メレル・ダウ・ファーマシューティカルズ事件10とクムホタイヤ対カーマイケル事件11の連邦最高裁判例に示されている。ドーバート基準によれば、証拠が「一般に認められている」だけではもはや十分ではなく、科学的専門家の知見は「信頼できる」方法に基づいていなければならない。したがって、法廷における精神医学的診断の科学的根拠には、既知の(そしておそらく許容できる)レベルの評価者間信頼性が含まれていなければならない。ドーバートとクムホの判決は、北米以外の国の司法にも影響を与え、オーストラリアの司法管轄区における鑑定人の行動規範の策定に貢献した9,13,14。ドーバート基準がPTSDの存在に関する証拠に厳格に適用され、その結果、PTSDは信頼性をもって診断できない、あるいはPTSDの診断方法は信頼できないという法的判断が下された場合、PTSDの存在に関する証拠は認められないことになる。

PTSDの診断と原因

DSM-Iでは「外傷性神経症 Traumatic neurosis」のカテゴリーが含まれていたが、DSM-IIでは「一過性の状況的障害」に置き換えられた。しかし、PTSDの診断は、1980年に発表されたDSMIII 7に含まれてから広く受け入れられるようになり、PTSDは現在、社会で最も頻繁に診断される不安症の1つとなっている16,17。PTSDDSMの診断カテゴリーの中では珍しく、A基準に障害の原因の定義が含まれている。診断基準は、DSM-III(1980年)7、DSM-III R(1987年)18、DSM-IV(1994年)19、DSM-IV TR(2000年)20の各版の間で改訂されており、特に必要とされるトラウマの定義が改訂されている。しかし、DSM-IIIで導入されたA-F基準の本質的な特徴は維持されている20,21。DSM-IV-TRのA基準では、''その人が、実際に、あるいは脅かされた死や重傷、あるいは自己や他者の身体的完全性に対する脅威を伴う出来事や出来事を経験した、目撃した、あるいは直面した''(A1基準)、および''その人の反応が強い恐怖、無力感、恐怖を伴うものであった''(A2基準)が要求されている。A基準の2つの部分を組み合わせることで、出来事とその後の症状との間に因果関係と時間的な関連性の両方が導入される。事実上、A2基準は最初の症状が出来事に反応したものであることを要求しているため、A基準の第3の側面(A1がA2を引き起こす)が生まれる。 外傷的な出来事の後に永続的な苦痛があり、その症状をその出来事に起因すると考える患者では、B~Eの基準を満たすかどうかにかかわらず、Aの基準を満たすことだけでPTSDの診断に十分であると考えられることがある。患者が他の基準を満たしたかどうかを確認するには、しばしば面接だけでは収集できないより多くの情報が必要となる。Bの基準であるトラウマの再体験は、報告された症状がAの基準で説明されたトラウマ的出来事に関連していなければならないため、因果関係の推論を含んでいる。しかし、その症状は非常に広範に定義されており、苦痛を伴う回想、夢、フラッシュバック、あるいは精神保健の専門家(必ずしも患者ではない)がトラウマ的出来事を象徴していると考える内的または外的な出来事に対する苦痛や覚醒の体験が含まれる。C基準は、7項目のうち少なくとも3項目の存在として定義される類似または象徴的状況の回避を必要とし、D基準は、5項目のうち少なくとも2項目の形で増大した覚醒を必要とする。残りの基準は、障害の持続期間と苦痛または障害のレベルのみを規定する。したがって、A2基準とB-D基準に記載されている17の症状のうち11の症状は、患者の自己報告に決定的に依存している。診断に不可欠な個々の症状はない。いくつかの症状の存在は観察によって確認することができるが、観察によってどのような症状の存在も除外することはできない。診断基準では、DSM-IVに記載されているように、''外傷的出来事の一面を象徴する、または類似する内的または外的手がかりにさらされたときの生理的反応性''があると臨床医が判断した場合、外傷的出来事の記憶がほとんどない、または報告できない患者をPTSDと診断することを認めている。その結果、PTSDの診断は完全に患者または評価者の主観的意見に基づくことになる。PTSDの診断を下す際、医師は患者の症状とある出来事を明確に結びつける。しかし、出来事と症状との関連は、人が注意義務を負っている外傷的出来事の後の精神的傷害に関する事件であれば、法廷における最終的な争点となりうる。もし裁判所がPTSDの診断を認めれば、その障害が外傷的出来事によって引き起こされたことを認めることになるが、その出来事がその人の障害を引き起こしたかどうかを立証するのは裁判所の特権である。これとは対照的に、裁判所がトラウマと症状との関連を否定すれば、診断の問題に関して臨床医を覆すことになる。このような状況は、PTSDのB、C、Dの基準を満たす患者が、(a)トラウマが起こっていない、(b)トラウマが十分に深刻でないと判断された、(c)トラウマに近接していなかったため、トラウマを目撃していないと判断された、(d)トラウマに反応して生じた直接的な有害な心理的反応がなかった、(e)心理的反応は生じたが、それはこのトラウマに反応したものではないと判断された、という理由で、裁判所によってPTSDではないと判断されたケースに示されている。これらの判決を下すことによって、裁判所は臨床医の診断を事実上覆すことになる。したがって、PTSDの診断に因果関係の要素を含めることは、臨床医が裁判官、裁判所が診断医という医事法上のパラドックスを生み出すことになる。 診断基準に障害の原因に関する記述を含めるという異例の措置、症状のほとんどが主観的なものであること、PTSDの診断が明らかに法的な意味を持つことを考えれば、PTSDの診断が導入される前に、障害の存在について臨床医の間で妥当なレベルの合意があったことを確認することが賢明であっただろう。しかし、それは行われなかった。

PTSD診断の信頼性に関するDSMとICDのフィールドトライアル

精神医学では生物学的マーカーがほとんどないため、診断の信頼性が多くの精神疾患の科学的基盤を形成している。診断の信頼性を測る尺度としては、臨床医間の一致度が一般的である。しかし、DSM-IIIでは、PTSDが確実に診断できる、あるいは他の疾患と区別できるという証拠がないうちに、PTSDの診断も含まれていた。DSM-III、DSM-IVICD-10の診断基準の信頼性を検証するためにフィールドトライアルが行われた。しかし、DSM-IIIとICD-10の試験ではPTSDは言及されておらず、PTSDに関するDSM-IVの試験では、臨床現場での診断の信頼性を確認できるような方法は用いられていない。DSM-III試験では、参加した臨床医が自分の患者を診断し、他の臨床医が同じ基準で診断した2人以上の患者を評価した 24。合計131組の評価が独立して行われ、さらに150人の患者が共同で評価された。被験者の合計10.5%が不安障害と診断され、その一致度はLandisとKochの基準を用いて公正(k ¼ 0.43)と評価された。DSM-IVのフィールド試験はより大規模であり、目的と方法が異なっていた26。PTSDは独自の試験の対象であり、400人の患者と128人の非患者が募集された。被験者にはストレスとなる出来事がないかスクリーニングが行われ、DSMIII-R(SCID)のPTSDモジュール27と診断面接スケジュール(DIS)を用いて面接が行われた。5施設から5人、合計25人の患者のSCIDのPTSDモジュールに対する回答の録音テープを用いて、相互評価者の一致が検討された。他の4施設からそれぞれ2人の臨床医がテープに録音された回答を評価し、各患者について合計8人の臨床医による評価が得られた。臨床医には症例ノートやその他の臨床情報は提供されなかったが、臨床医は患者がDSM-IV PTSDフィールドトライアルの一部として評価されていることを知っていた。DSM-IV PTSDフィールドトライアルの結果は査読誌に発表されておらず、著者との個人的なコミュニケーションにより、臨床医が互いの評価を盲検化していたかどうか、あるいは統計解析に同一被験者の複数の評価を使用することから生じる統計的自由度の損失に対する補正が含まれていたかどうかを確認することはできなかった。DSM-IV臨床試験におけるPTSDのκ値は0.85であり、これはDSM-III臨床試験における精神病性障害のκ値を上回り、ほぼ完全な一致を示している。しかし、より重要な批判は、半構造化面接に対する録音された回答についての一致は、臨床医が特定の質問に対する回答に同意できるかどうかだけを確認するものであるということである。この方法では、PTSDを臨床の場で確実に診断できるかどうか、PTSDを他の疾患と区別できるかどうかを立証することはできない。ICD-10のフィールドトライアルには11,000人以上の患者が参加したが、PTSDと診断された患者は9人程度であった可能性があり、PTSDの診断に関する一致度のカッパは報告されていない29。われわれは、自動車事故後の559件の連続した賠償請求において作成された医療訴訟報告書における精神医学的診断の評価者間信頼性を検討し、この設定におけるPTSDの診断に関する一致度は低いことを明らかにした30-32。われわれは、敵対的な訴訟において対立する側の専門家の間で意見の相違があるだろうと予想していたが、同じ側の専門家の間でPTSDの有無に関する一致が偶然に予想される程度しかないことに驚いた。その研究では、認知された診断基準を満たすのに十分な症状を報告書に記載した専門家はほとんどおらず、また、PTSDを引き起こすのに十分な外傷的体験であったかどうかについても、専門家の主観的判断は一致していなかった。重傷で病院の外科病棟に入院した患者よりも、事故後に救急外来を受診しなかった患者にPTSDの診断がわずかに多かったことから、重症度に関する意見の相違は驚くべきことではなかった30。Medline、PubMed、PsychLitおよびCINHALを、「信頼性」、「心的外傷後ストレス障害」または「不安障害」という用語を用いて検索した。その結果、自己報告式の質問票、構造化面接および半構造化面接の信頼性を報告した研究が1件、構造化面接と臨床面接の結果を比較した研究が1件見つかった33。しかし、DSM-IIIおよびICD-10のフィールドトライアルに類似した方法を用いた研究は見つからなかった。したがって、非構造化臨床評価は、裁判用に作成された報告書では通常PTSDを診断する方法であるが、この方法でPTSDを確実に診断できるという証拠はほとんどない。

PTSDの診断法と医療・法律実務 Spitzer

この基準は、一般的に医療 法的評価にも適用することができ、「利用可能なデー タ」には、構造化または半構造化面接の結果が含まれ、検討されている障害の実際の診断基準が実際に適 用されていることを確認することができる。診断のための面接は、包括的な面接、障害に焦点を当てた面接、症状に焦点を当てた方法に大別することができる。包括的な構造化診断機器の例としては、DIS28、SCID27、Composite International Diagnostic Interview (CIDI)35、Schedules for Clinical Assessment in Neuropsychiatry36などがある。この種の機器は、診断を下すために非常に広範な症状の存在を記録するが、これらの機器の使用の限界は、1000もの質問に対する回答を必要とすることである。PTSD診断のための包括的診断面接の評価者間信頼性を、利用可能なすべての情報源を用いた経験豊富な臨床医の診断と比較したKomitiらによる1件の研究がある37。CAPSや類似の尺度はPTSDの診断を行う。しかし、障害に焦点をあてた面接は、実施する前に仮診断か、少なくともPTSDのような障害の存在に関するある程度の疑いを必要とする。したがって、障害に焦点を当てた評価の主な役割は、高リスク群における偽陽性を除外して診断の存在を確認することである。3つ目のタイプは症状に焦点をあてた質問紙で、短時間で簡単に実施できるが、診断はつかない。Millerは、半構造化面接、完全構造化面接、非構造化面接をLEAD標準と比較し、非構造化面接の信頼性が最も低いことを明らかにした33,42。医療・法律専門家が構造化面接を避けるのは、おそらく、不慣れ、検査実施に時間がかかる、検査機器自体への信頼がない、自分の専門知識に自信がないなど、さまざまな理由がある。専門家はまた、構造化面接の特徴である、障害の症状に関する誘導的な質問で患者を促すことを避けたいのかもしれない。対照的に、PTSDのシミュレーションを避けるために診断基準を選択的または特異的に使用することは、診断の誤りを招く可能性が高い。実際、人々に障害の有無を確認する障害に焦点をあてた尺度は、精神的傷害に対する医学的法的評価の診断部分とほとんど同じ機能を果たしている。また、構造化面接の適用がラポールに及ぼす影響も、その使用に対する有効な反論とはならない。診断尺度や 構造化面接の使用に反対する最も強い論拠は、PTSDの診断基準を満たす患者に限定してしまうことであろう。現在の診断基準に関する懸念はひとまず置いておくとして、医療現場における検査機器の使用に対する残りの反対意見には説得力がない。

PTSDに関する証拠の信頼性を向上させるためにできること

PTSDの診断の信頼性を向上させることは、DSM-Vの基準を起草する委員会にとって重要な検討事項であろう。論理的なステップとしては、A基準を完全に削除することであろう。これにより、トラウマまたはストレス誘発性不安障害の診断基準を満たす人の数が増える可能性があるが、他の疾患と同じように、患者の症状に基づいて診断を下すことができるようになる。A基準を削除すれば、患者の発達歴、以前の精神疾患、訴訟の影響を含む患者の状況など、侵入記憶や回避を伴う不安症状の最近の発現に関連する他の病因論的因子を考慮することも可能になる。A基準に関する専門家証拠も、トラウマと患者の状態との関連性に関する証拠は裁判所が決定すべき問題であるという理由で、裁判所によって除外される可能性がある。究極的争点規定9が存在する場合、この規定が破られることは避けられるであろう。現行のA基準は平易な英語で明確に記述されており、A1基準またはA2基準のいずれかが満たされているかどうかを専門家が判断する必要はない。A基準に関する証拠を除外しても、症状の存在が信頼できる方法で記録されていれば、事件発生時の感情状態(恐怖、無力感、恐怖)とその後の症状との関連性に関する専門家の証言は除外されない。これによって、現在の診断基準が適用される際にしばしば生じる、精神医学的診断とその病因に関する循環論議がなく、出来事が症状を引き起こしたかどうかという問題を客観的に検討することができる。裁判における専門家評価の欠点の一つは、通常、評価を行うために割り当てられる時間が比較的短いことである。外傷的な出来事を詳細に調べるという要件を削除すれば、患者の症状についてより詳細な評価を行うための時間を確保することができる。これには、B、C、Dの基準を十分に考慮するための半構造化面接を行う時間や、最も適切な診断がなされるようにするための他の障害の症状も含まれる。

結論

我々が確認した限りでは、PTSDの診断基準は、臨床場面における診断基準の信頼性を決定するための十分なフィールドトライアルの対象となったことはない。構造化面接は信頼できることが示されているが、構造化面接や半構造化面接の所見を臨床診断と比較した数少ない研究は心強いものではなく、医療訴訟の場における臨床面接の結果の信頼性を報告した唯一の研究は、PTSDの診断の信頼性が非常に低いことを示した。医療現場におけるPTSD診断の信頼性は、構造化または半構造化尺度を用いてすべての診断基準を系統的に検討することによって改善されるであろう。我々は、A基準に関する専門家の証拠を避けることによって、患者の状態の因果関係の問題を法廷に委ねることを推奨する。これは、DSM-VのPTSDの診断基準からA基準を削除するか、場合によっては裁判所が認めるかどうかを決定することによって達成できるであろう。そうすれば、PTSDの診断は、他のほとんどの精神疾患と同じように、患者の症状と障害に基づいて行われることになる。

lasso, adaptive lasso, group lasso[R]

https://www.math.mcgill.ca/yyang/comp/notes/note4code.R

library(glmnet)

このパッケージで使用されるデフォルトのモデルは、Guassian線形モデルまたは最小二乗法である。説明のためにあらかじめ作成したデータセットをロードする。ユーザは自分のデータをロードすることも、ワークスペースに保存されているデータを使用することもできる。

pacman::p_install_gh("emeryyi/gglasso")
library(gglasso)
data(bardet)
x<- bardet[["x"]]
y<- bardet[["y"]]

glmnetの最も基本的な呼び出しを使ってモデルをフィットさせる。

fit = glmnet(x, y)

"fit "はglmnetクラスのオブジェクトで、フィットされたモデルの関連情報をすべて含んでいる。我々は、ユーザが直接成分を抽出することを推奨していない。その代わりに、plot、print、coef、predictなどの様々なメソッドがオブジェクトに提供され、これらのタスクをよりエレガントに実行することができる。

plot関数を実行すれば、係数を可視化することができる:

plot(fit)

各曲線は変数に対応する。これは、λが変化する時の係数ベクトル全体の「1-norm」に対する係数のパスを示している。上の軸は、現在のλにおける非ゼロ係数の数を示しており、これはlassoの有効自由度(df)となる。ユーザーは、曲線に注釈を付けることもでき、これはplotコマンドでlabel = TRUEを設定することで可能だ。

オブジェクト名を入力するか、print機能を使えば、各ステップでのglmnetパスの概要が表示される:

print(fit)
Call:  glmnet(x = x, y = y) 

    Df  %Dev   Lambda
1    0  0.00 0.097200
2    1  7.73 0.088560
3    1 14.16 0.080690
4    1 19.49 0.073520
5    1 23.91 0.066990
6    1 27.59 0.061040
7    2 31.42 0.055620
8    2 36.05 0.050680
9    4 40.59 0.046180
10   6 45.05 0.042070
...

左からゼロでない係数の数(Df)、-log(尤度)の値(%dev)、ラムダの値(Lambda)が表示される。デフォルトではglmnetは100個のλの値を要求するが、%dev%がλの値から次のλの値まで十分に変化しない場合、プログラムは早期に停止する(典型的にはパスの終わり近く)。

配列の範囲内の1つ以上のラムダで実際の係数を得ることができる:

coef0 = coef(fit,s=0.1)

関数glmnetは、ユーザーが選択できる一連のモデルを返す。多くの場合、ユーザーはソフトウェアがその中から1つを選択することを好むかもしれない。クロスバリデーションは、おそらくそのタスクのための最もシンプルで最も広く使われている方法である。

cv.glmnetは、プロットや予測などの様々なサポートメソッドとともに、クロスバリデーションを行うためのメイン関数です。ここでは、前に読み込んだサンプルデータを使っている。

cvfit = cv.glmnet(x, y)

cv.glmnetは、cv.glmnetオブジェクトを返しますが、これはここでは "cvfit"であり、クロスバリデーション適合のすべての成分を含むリストです。glmnetに関しては、選択されたλの値を見る以外には、ユーザが直接成分を抽出することは推奨しません。このパッケージは、潜在的なタスクのためによく設計された関数を提供する。

オブジェクトをプロットできる。

plot(cvfit)

クロスバリデーション曲線(赤い点線)、ラムダ配列に沿った上下の標準偏差曲線(エラーバー)が含まれる。選択された2つのラムダは縦の点線で示されている(下図参照)。

選択されたラムダと対応する係数を見ることができる。例えば、

cvfit$lambda.min

[1] 0.001226463

lambda.minは、クロスバリデーションの平均誤差を最小にするλの値である。もう1つの保存されたλはλ.1seで、これは誤差が最小の標準誤差の1つ以内に収まるような最も正則化されたモデルを与える。これを使うには、上記のlambda.minをlambda.1seに置き換えるだけでよい。

coef1 = coef(cvfit, s = "lambda.min")

係数はスパース行列形式で表現されていることに注意。その理由は、正則化パスに沿った解はスパースであることが多いため、スパースフォーマットを使用する方が時間的にも空間的にも効率的だからだ。スパースでない形式を好む場合は、出力をas.matrix()にパイプを通して出力する。

予測は、フィットされたcv.glmnetオブジェクトに基づいて行うことができる。おもちゃの例を見てみよう。

predict(cvfit, newx = x[1:5,], s = "lambda.min")
  lambda.min

[1,] 8.406305 [2,] 8.353781 [3,] 8.377383 [4,] 8.294958 [5,] 8.379072

newxは新しい入力行列を表し、sは前回と同様、予測が行われるλの値(複数可)である。

adaptive lasso

ペナルティ・ファクターの例

[penalty.factor]引数により、各係数に別々のペナルティ係数を適用することができる。デフォルトは各パラメータに対して1だが、他の値を指定することもできる。特に、[penalty.factor]が0に等しい変数は、全くペナルティを受けません!v_jがj番目の変数の制約を表すとする。

ペナルティ制約は、内部的にnvarsの合計になるように再スケーリングされることに注意。

これは、変数に対する予備知識や好みがある場合に非常に便利だ。多くの場合、ある変数が非常に重要で、それを常に保持したい場合があり、これは対応する制約係数を0に設定することで実現できる: v_j は、ペナルティ係数を0に設定することで、ペナルティ係数を0にすることができる:

p.fac = rep(1, 120)
p.fac[c(31, 74, 111)] = 0
pfit = glmnet(x, y, penalty.factor = p.fac)
plot(pfit, xvar="lambda", label = TRUE)

その他の便利な引数をいくつか挙げる: [exclude]は、特定の変数をモデルにしないようにすることができる。もちろん、単純にxからこれらの変数を除外することもできるが、除外された変数をゼロに設定するだけで、係数の完全なベクトルを返すので、excludeの方が便利な場合もある。FALSEの場合、切片は強制的にゼロになる。

Adaptive lassoの例

adaptive lassoは一貫性のある第1段階を必要とする。 Zou (2006)はOLSかridge第一段lassoを推奨している。

thelasso.cv<-cv.glmnet(x,y,family = "gaussian",alpha=1) 

第1段階の係数から第2段階の重みを求める coef() はスパース行列

bhat<-as.matrix(coef(thelasso.cv,s="lambda.1se"))[-1,1] 
if(all(bhat==0)){
  bhat<-rep(.Machine$double.eps*2,length(bhat))
}

もしbhatがすべてゼロなら、すべてにゼロに極めて近いウェイトを割り当てる。これは、第2ステージのすべてにゼロのペナルティを課すことに等しい。

adaptive lasso weight

adpen<-(1/pmax(abs(bhat),.Machine$double.eps)) 

第2段階のlasso(adaptive lasso)

m_adlasso <- glmnet(x,y,family = "gaussian",alpha=1,exclude=which(bhat==0),
penalty.factor=adpen)
plot(m_adlasso)

エラスティックネットelastic-net

エラスティック・ネットの例

glmnetには、ユーザーがフィットをカスタマイズするための様々なオプションが用意されている。ここではよく使われるオプションをいくつか紹介するが、これらはglmnet関数で指定することができる。

[family="gaussian"]はglmnet関数のデフォルトのfamilyオプションだ。

[alpha]は弾性ネットの混合パラメータalphaで、alphaの範囲は[0,1]です。alpha=1はlasso(デフォルト)、alpha=0はリッジになる。

[nlambda] はシーケンス中のラムダ値の数であり、デフォルトは100ですが、[0,1]の範囲もある。

例として、α=0.2(よりリッジ回帰に近い)を設定し、オブザベーションの後半に2倍の重みを与える。 ここでの表示が長くなりすぎるのを避けるために、nlambdaを20に設定する。しかし実際には、λの値の数は100(デフォルト)以上にすることが推奨される。ほとんどの場合、アルゴリズムで使用されるウォーム・スタート warm-starts のために、余分なコストがかからず、非線形モデルではよりよい収束特性につながる。

fit = glmnet(x, y, alpha = 0.2, family="gaussian")

plot(fit, xvar = "lambda", label = TRUE)

グループlasso

グループlassoの例

install.packages("gglasso", repos = "http://cran.us.r-project.org")

https://cran.r-project.org/web/packages/gglasso/index.html

group-lasso ペナルティ付き最小二乗法、ロジスティック回帰、Huberized SVM、二乗SVMの解経路を効率的に計算するための統一アルゴリズム、blockwise-majorization-descent (BMD)。 このパッケージは、Yang, Y. and Zou, H. (2015) DOI: <doi:10.1007/s11222-014-9498-5> の実装だ。

library(gglasso)

bardetデータセットの読み込み

data(bardet)

20のグループを定義する

group1 <- rep(1:20,each=5)

ペナルティ付き最小二乗"group lasso"にフィットさせる。

m1 <- gglasso(x=bardet$x,y=bardet$y,group=group1,loss="ls")
plot(m1)

ペナルティ付きグループLassoでls回帰を用いた5-fold クロスバリデーション

cv <- cv.gglasso(x=bardet$x, y=bardet$y, group=group1, loss="ls",
pred.loss="L2", lambda.factor=0.05, nfolds=5)
plot(cv)

colon データセットを読み込む

data(colon)

グループ指標を定義する

group2 <- rep(1:20,each=5)

グループlassoペナルティ付きロジスティック回帰にフィット

m2 <- gglasso(x=colon$x,y=colon$y,group=group2,loss="logit")
plot(m2)

グループlassoペナルティ付きロジスティック回帰を用いた5-foldクロスバリデーション

cv2 <- cv.gglasso(x=colon$x, y=colon$y, group=group2, loss="logit",
pred.loss="misclass", lambda.factor=0.05, nfolds=5)
plot(cv2)

λ=λ.1seにおける係数

pre = coef(cv$gglasso.fit, s = cv$lambda.1se)

シミュレーション・デモ-1

ロジスティック回帰と最小二乗法におけるlasso、adaptive lasso、scad、mcpの比較

Rライブラリをロード
glmnet: LASSO、adaptive LASSO、elastic net ペナルティ付き最小2乗法およびロジスティック回帰 (すべてポアソン多項式、コックスモデルをサポート)
ncvreg: SCAD および MCP ペナルティ化最小2乗法およびロジスティック回帰用

library(glmnet)
library(ncvreg)
library(MASS)

PART I 最小二乗法

n=100
p=200

真のベータ

truebeta <- c(4,4,4,-6*sqrt(2),4/3,rep(0,p-5))

誤差分散

sigma2 <- 0.3

XjとXkの間の共分散はcov(X_j, X_k) = rho^|i-j|共分散行列である。

covmat <- function(rho, p) {
  rho^(abs(outer(seq(p), seq(p), "-")))
}

rho = 0.1、Xの共分散行列を生成する。

sigma <- covmat(0.1,p)

X ~ N(0, sigma)
epsilon ~ N(0, sigma2)
真のモデル: y = x*truebeta + epsilon

x <- mvrnorm(n,rep(0,p),sigma)
epsilon <- rnorm(n,0,sd=sqrt(sigma2))
y <- x %*% truebeta + epsilon

lassoをフィットさせ、5-fold CVでラムダを選択する。

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, family = "gaussian")
plot(cvfit)

CVで選ばれたラムダ

cvfit$lambda.min

[1] 0.0821079

解決パスをプロットする

plot(cvfit$glmnet.fit,label=TRUE,xvar="lambda")

λ.minのプロット

abline(v=log(cvfit$lambda.min),lty=1)

λ.1seをプロットする

abline(v=log(cvfit$lambda.1se),lty=2)
model_compare <- matrix(NA, nrow=5, ncol=p,
                dimnames=list(c("true model",
                "lasso","adaptive lasso","mcp","scad"),
                paste("V",seq(p),sep="")))

真のモデルを保存する

model_compare[1, ] <- truebeta

lassoによって推定された係数

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, 
family = "gaussian")
tmp <- cvfit$glmnet.fit$beta
lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[2, ] <- lasso_beta

重みを計算し、adaptive lassoにフィットさせる

weight = 1/(lasso_beta)^2

推定係数がゼロの場合、対応するウェイトがInfとなり、数値誤差を防ぐため、Infを大きな数値(例えば1e6)に変換する。

weight[weight==Inf] = 1e6

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, 
        family = "gaussian", 
        nfolds = 5, penalty.factor=weight)

Adaptive Lassoによって推定された係数

tmp <- cvfit$glmnet.fit$beta
adaptive_lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[3, ] <- adaptive_lasso_beta

mcpによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "MCP", family = "gaussian")
mcp_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[4, ] <- mcp_beta[-1]

scadによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "SCAD", family = "gaussian")
scad_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[5, ] <- scad_beta[-1]

4つの方法で推定された係数を比較すると、lasso over-selected、adaptive lasso、scad、mcpが問題を解決していることがわかる。

model_compare
               V1           V2 V3        V4       V5 V6 V7 V8 V9 V10
true model      4  4.000000000  4 -8.485281 1.333333  0  0  0  0   0
lasso           0 -0.022023796  0  0.000000 0.000000  0  0  0  0   0
adaptive lasso  0  0.000000000  0  0.000000 0.000000  0  0  0  0   0
mcp             0  0.000000000  0  0.000000 0.000000  0  0  0  0   0
scad            0 -0.004545409  0  0.000000 0.000000  0  0  0  0   0
                       V11         V12        V13       V14       V15 V16
true model     0.000000000  0.00000000  0.0000000 0.0000000  0.000000   0
lasso          0.006311609  0.00000000  0.0000000 0.0000000  0.000000   0
adaptive lasso 0.000000000  0.00000000  0.0000000 0.0000000  0.000000   0
mcp            0.000000000 -0.08207367 -0.1647326 0.2088095 -1.000567   0
scad           0.000000000  0.00000000  0.0000000 0.0000000 -1.033549   0
                      V17         V18 V19 V20 V21          V22         V23
true model     0.00000000  0.00000000   0   0   0 0.0000000000  0.00000000
lasso          0.00000000 -0.01762063   0   0   0 0.0001952964 -0.05328850
adaptive lasso 0.00000000  0.00000000   0   0   0 0.0000000000  0.00000000
mcp            0.03504008  0.00000000   0   0   0 0.0000000000  0.00000000
scad           0.00000000  0.00000000   0   0   0 0.0000000000 -0.08488739
                      V24        V25        V26 V27        V28 V29        V30
true model      0.0000000  0.0000000 0.00000000   0  0.0000000   0  0.0000000
lasso           0.0000000 -0.5648782 0.06002556   0  0.0000000   0 -0.3371098
adaptive lasso  0.0000000 -1.0114037 0.00000000   0  0.0000000   0  0.0000000
mcp            -0.4379709  0.0000000 0.00000000   0 -0.1994319   0  0.0000000
scad            0.0000000  0.0000000 0.00000000   0  0.0000000   0  0.0000000
                      V31 V32        V33 V34 V35       V36 V37 V38 V39 V40 V41
true model      0.0000000   0  0.0000000   0   0 0.0000000   0   0   0   0   0
lasso          -0.1165544   0  0.0000000   0   0 0.0482393   0   0   0   0   0
adaptive lasso  0.0000000   0  0.0000000   0   0 0.0000000   0   0   0   0   0
mcp            -0.4561359   0 -0.1259705   0   0 0.1201394   0   0   0   0   0
scad            0.0000000   0  0.0000000   0   0 0.0000000   0   0   0   0   0
               V42 V43         V44 V45        V46       V47 V48 V49        V50
true model       0   0 0.000000000   0  0.0000000 0.0000000   0   0  0.0000000
lasso            0   0 0.009982214   0 -0.2897772 0.0000000   0   0  0.0000000
adaptive lasso   0   0 0.000000000   0 -0.6616183 0.0000000   0   0  0.0000000
mcp              0   0 0.000000000   0 -0.4648836 0.1156758   0   0 -0.1438189
scad             0   0 0.000000000   0 -0.1163973 0.0000000   0   0  0.0000000
                      V51 V52        V53        V54 V55 V56 V57         V58 V59
true model      0.0000000   0 0.00000000 0.00000000   0   0   0  0.00000000   0
lasso           0.0000000   0 0.05634673 0.13775797   0   0   0  0.00000000   0
adaptive lasso  0.0000000   0 0.00000000 0.09460430   0   0   0  0.00000000   0
mcp             0.0000000   0 0.07499241 0.11837706   0   0   0 -0.05640483   0
scad           -0.1008567   0 0.00000000 0.08972644   0   0   0  0.00000000   0
               V60 V61 V62        V63 V64       V65         V66 V67 V68 V69 V70
true model       0   0   0 0.00000000   0 0.0000000  0.00000000   0   0   0   0
lasso            0   0   0 0.00000000   0 0.2563789 -0.12041928   0   0   0   0
adaptive lasso   0   0   0 0.00000000   0 0.2957797  0.00000000   0   0   0   0
mcp              0   0   0 0.02528918   0 0.3411462  0.00000000   0   0   0   0
scad             0   0   0 0.00000000   0 0.1248035 -0.03840615   0   0   0   0
               V71 V72        V73 V74 V75         V76       V77        V78 V79
true model       0   0 0.00000000   0   0  0.00000000 0.0000000 0.00000000   0
lasso            0   0 0.05161909   0   0  0.00000000 0.0000000 0.04393271   0
adaptive lasso   0   0 0.00000000   0   0  0.00000000 0.0000000 0.00000000   0
mcp              0   0 0.07055282   0   0  0.00000000 0.0652188 0.12140394   0
scad             0   0 0.00000000   0   0 -0.06290252 0.0000000 0.00000000   0
               V80 V81 V82 V83 V84 V85        V86          V87 V88       V89
true model       0   0   0   0   0   0 0.00000000  0.000000000   0 0.0000000
lasso            0   0   0   0   0   0 0.00000000 -0.001347216   0 0.2044193
adaptive lasso   0   0   0   0   0   0 0.00000000  0.000000000   0 0.1554728
mcp              0   0   0   0   0   0 0.04738824  0.000000000   0 0.4817615
scad             0   0   0   0   0   0 0.00000000  0.000000000   0 0.0000000
               V90 V91 V92 V93        V94 V95 V96 V97 V98 V99       V100 V101
true model       0   0   0   0 0.00000000   0   0   0   0   0 0.00000000    0
lasso            0   0   0   0 0.07206646   0   0   0   0   0 0.06736159    0
adaptive lasso   0   0   0   0 0.00000000   0   0   0   0   0 0.00000000    0
mcp              0   0   0   0 0.44369829   0   0   0   0   0 0.14673570    0
scad             0   0   0   0 0.00000000   0   0   0   0   0 0.00000000    0
                       V102 V103 V104 V105 V106 V107 V108 V109 V110        V111
true model      0.000000000    0    0    0    0    0    0    0    0 0.000000000
lasso          -0.022023796    0    0    0    0    0    0    0    0 0.006311609
adaptive lasso  0.000000000    0    0    0    0    0    0    0    0 0.000000000
mcp             0.000000000    0    0    0    0    0    0    0    0 0.000000000
scad           -0.004545409    0    0    0    0    0    0    0    0 0.000000000
                      V112       V113      V114      V115 V116       V117
true model      0.00000000  0.0000000 0.0000000  0.000000    0 0.00000000
lasso           0.00000000  0.0000000 0.0000000  0.000000    0 0.00000000
adaptive lasso  0.00000000  0.0000000 0.0000000  0.000000    0 0.00000000
mcp            -0.08207367 -0.1647326 0.2088095 -1.000567    0 0.03504008
scad            0.00000000  0.0000000 0.0000000 -1.033549    0 0.00000000
                      V118 V119 V120 V121         V122        V123       V124
true model      0.00000000    0    0    0 0.0000000000  0.00000000  0.0000000
lasso          -0.01762063    0    0    0 0.0001952964 -0.05328850  0.0000000
adaptive lasso  0.00000000    0    0    0 0.0000000000  0.00000000  0.0000000
mcp             0.00000000    0    0    0 0.0000000000  0.00000000 -0.4379709
scad            0.00000000    0    0    0 0.0000000000 -0.08488739  0.0000000
                     V125       V126 V127       V128 V129       V130       V131
true model      0.0000000 0.00000000    0  0.0000000    0  0.0000000  0.0000000
lasso          -0.5648782 0.06002556    0  0.0000000    0 -0.3371098 -0.1165544
adaptive lasso -1.0114037 0.00000000    0  0.0000000    0  0.0000000  0.0000000
mcp             0.0000000 0.00000000    0 -0.1994319    0  0.0000000 -0.4561359
scad            0.0000000 0.00000000    0  0.0000000    0  0.0000000  0.0000000
               V132       V133 V134 V135      V136 V137 V138 V139 V140 V141
true model        0  0.0000000    0    0 0.0000000    0    0    0    0    0
lasso             0  0.0000000    0    0 0.0482393    0    0    0    0    0
adaptive lasso    0  0.0000000    0    0 0.0000000    0    0    0    0    0
mcp               0 -0.1259705    0    0 0.1201394    0    0    0    0    0
scad              0  0.0000000    0    0 0.0000000    0    0    0    0    0
               V142 V143        V144 V145       V146      V147 V148 V149
true model        0    0 0.000000000    0  0.0000000 0.0000000    0    0
lasso             0    0 0.009982214    0 -0.2897772 0.0000000    0    0
adaptive lasso    0    0 0.000000000    0 -0.6616183 0.0000000    0    0
mcp               0    0 0.000000000    0 -0.4648836 0.1156758    0    0
scad              0    0 0.000000000    0 -0.1163973 0.0000000    0    0
                     V150       V151 V152       V153       V154 V155 V156 V157
true model      0.0000000  0.0000000    0 0.00000000 0.00000000    0    0    0
lasso           0.0000000  0.0000000    0 0.05634673 0.13775797    0    0    0
adaptive lasso  0.0000000  0.0000000    0 0.00000000 0.09460430    0    0    0
mcp            -0.1438189  0.0000000    0 0.07499241 0.11837706    0    0    0
scad            0.0000000 -0.1008567    0 0.00000000 0.08972644    0    0    0
                      V158 V159 V160 V161 V162       V163 V164      V165
true model      0.00000000    0    0    0    0 0.00000000    0 0.0000000
lasso           0.00000000    0    0    0    0 0.00000000    0 0.2563789
adaptive lasso  0.00000000    0    0    0    0 0.00000000    0 0.2957797
mcp            -0.05640483    0    0    0    0 0.02528918    0 0.3411462
scad            0.00000000    0    0    0    0 0.00000000    0 0.1248035
                      V166 V167 V168 V169 V170 V171 V172       V173 V174 V175
true model      0.00000000    0    0    0    0    0    0 0.00000000    0    0
lasso          -0.12041928    0    0    0    0    0    0 0.05161909    0    0
adaptive lasso  0.00000000    0    0    0    0    0    0 0.00000000    0    0
mcp             0.00000000    0    0    0    0    0    0 0.07055282    0    0
scad           -0.03840615    0    0    0    0    0    0 0.00000000    0    0
                      V176      V177       V178 V179 V180 V181 V182 V183 V184
true model      0.00000000 0.0000000 0.00000000    0    0    0    0    0    0
lasso           0.00000000 0.0000000 0.04393271    0    0    0    0    0    0
adaptive lasso  0.00000000 0.0000000 0.00000000    0    0    0    0    0    0
mcp             0.00000000 0.0652188 0.12140394    0    0    0    0    0    0
scad           -0.06290252 0.0000000 0.00000000    0    0    0    0    0    0
               V185       V186         V187 V188      V189 V190 V191 V192 V193
true model        0 0.00000000  0.000000000    0 0.0000000    0    0    0    0
lasso             0 0.00000000 -0.001347216    0 0.2044193    0    0    0    0
adaptive lasso    0 0.00000000  0.000000000    0 0.1554728    0    0    0    0
mcp               0 0.04738824  0.000000000    0 0.4817615    0    0    0    0
scad              0 0.00000000  0.000000000    0 0.0000000    0    0    0    0
                     V194 V195 V196 V197 V198 V199       V200
true model     0.00000000    0    0    0    0    0 0.00000000
lasso          0.07206646    0    0    0    0    0 0.06736159
adaptive lasso 0.00000000    0    0    0    0    0 0.00000000
mcp            0.44369829    0    0    0    0    0 0.14673570
scad           0.00000000    0    0    0    0    0 0.00000000

シミュレーション・デモ-2

PART II ロジスティック回帰

データを生成する

n <- 200 
p <- 8

真のベータ

truebeta <- c(6,3.5,0,5,rep(0,p-4)) 
truebeta

真のモデルからxとyを生成 真のモデル: P(y=1|x) = exp(xtruebeta)/(exp(xtruebeta)+1)

x <- matrix(rnorm(n*p), n, p)
feta <- x %*% truebeta 
fprob <- ifelse(feta < 0, exp(feta)/(1+exp(feta)), 1/(1 + exp(-feta)))
y <- rbinom(n, 1, fprob)


model_compare <- matrix(NA, nrow=5, ncol=p,
                        dimnames=list(c("true model","lasso",
                        "adaptive lasso","mcp","scad"),
                        paste("V",seq(p),sep="")))

真のモデルを保存する

model_compare[1, ] <- truebeta

lassoのケース
cv.glmfit ラッソモデルをフィットし、ラムダの選択にクロス検証を使用する。

family = "binomial", ロジスティック回帰
family = "gaussian", 最小二乗法

alphaはL1ペナルティ項の度合いを制御する、
alpha = 1, lasso、
alpha = 0, ridge、
alpha = (0,1), 弾性ネット

nfolds = 5, 5回クロスバリデーション(CV)

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, family = "binomial", nfolds = 5)

CV結果をプロットする 左の縦線(λ.min)は、最小の偏差を与えるλに対応する。
右の縦線(λ.1se)は、1標準偏差ルールのλに対応する。

plot(cvfit)

解決パスをプロットする

plot(cvfit$glmnet.fit,label=TRUE,xvar="lambda")

plot(cvfit$glmnet.fit,label=TRUE,xvar="lambda")

λ.minをプロットする

abline(v=log(cvfit$lambda.min),lty=1)

λ.1seをプロットする

abline(v=log(cvfit$lambda.1se),lty=2)

各列はラムダ値の推定値を表す。

tmp <- cvfit$glmnet.fit$beta
tmp

CVが選択したλ.minに対応するベータ値

lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[2, ] <- lasso_beta

重みを計算し、Adaptive Lassoをフィットさせる

weight = 1/(lasso_beta)^2

推定係数がゼロの場合、対応するウェイトがInfとなり、数値誤差を防ぐため、Infを大きな数値(例えば1e6)に変換する。

weight[weight==Inf] = 1e6

cvfit <- cv.glmnet(x = x, y = y, alpha = 1, family = "binomial", 
                        nfolds = 5, penalty.factor=weight)

Adaptive Lassoによって推定された係数

tmp <- cvfit$glmnet.fit$beta
adaptive_lasso_beta <- as.matrix(tmp[,cvfit$lambda==cvfit$lambda.min])
model_compare[3, ] <- adaptive_lasso_beta

mcpによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "MCP", family = "binomial")
mcp_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[4, ] <- mcp_beta[-1]

scadによって推定された係数

cvfit <- cv.ncvreg(X = x, y = y, penalty = "SCAD", family = "binomial")
scad_beta <- cvfit$fit$beta[, cvfit$min]
model_compare[5, ] <- scad_beta[-1]

4つの方法から推定された係数を比較する。

model_compare
                     V1       V2 V3       V4 V5          V6         V7
true model           NA       NA NA       NA NA          NA         NA
lasso          4.467528 2.375300  0 3.096532  0 -0.07315415 -0.3499353
adaptive lasso 4.664295 2.480856  0 3.335227  0  0.00000000  0.0000000
mcp            4.459921 2.533341  0 3.366028  0  0.00000000  0.0000000
scad           4.808539 2.744814  0 3.615375  0  0.00000000  0.0000000
                      V8
true model            NA
lasso          0.3801377
adaptive lasso 0.0000000
mcp            0.0000000
scad           0.0000000

summarytoolsの紹介[R]

cran.r-project.org

cran.r-project.org

www.rdocumentation.org

1. 概要

summarytoolsはデータ探索と簡単な報告を中心とした一貫した機能セットを提供する。その中核には以下の4つの関数が存在する:

関数説明
freq() 回数、比率、累積統計量、および欠損データの報告を特徴とする度数表

ctable() 離散変数/カテゴリー変数のペア間のクロス集計 (結合度数)。

descr() 数値データの記述統計 (単変量)。

dfSummary() すべての変数の型固有の情報を特徴とするデータ・フレーム・サマリー:単変量統計量および/または度数分布、棒グラフまたはヒストグラム、および欠損データのカウントと割合。素早く異常を検出し、一目で傾向を特定するのに非常に便利

2. 度数表: freq()

freq()関数は、計数、割合、欠損データ情報を含む度数表を生成する。
余談: このパッケージのアイデアは、ベースとなるRにこのような関数がなかったことに由来している。

freq(iris$Species, plain.ascii = FALSE, style = "rmarkdown")
  Freq % Valid % Valid Cum. % Total % Total Cum.
setosa 50 33.33 33.33 33.33 33.33
versicolor 50 33.33 66.67 33.33 66.67
virginica 50 33.33 100.00 33.33 100.00
\<NA> 0 0.00 100.00
Total 150 100.00 100.00 100.00 100.00

この最初の例では、plain.asciiとstyle引数が指定されている。しかし、st_options()を使ってこの文書に対してグローバルに定義したので、冗長であり、以後省略する(セクション16にこのヴィネットの設定の詳細な説明がある)。

2.1 欠落データ

summarytoolsの主な目的の1つは、さらなる分析のためのデータのクリーニングと準備を支援することである。しかし、状況によっては、欠損データに関する情報を必要としない(あるいは既に持っている)場合がある。report.nas = FALSEを使用すると、出力表は1行2列と小さくなる:

freq(iris$Species, report.nas = FALSE, headings = FALSE)
               Freq        %   % Cum.

      setosa     50    33.33    33.33
  versicolor     50    33.33    66.67
   virginica     50    33.33   100.00
       Total    150   100.00   100.00

headings = FALSE パラメーターは見出しセクションを抑制する。

2.2 最も単純な式 すべてのオプション要素を "オフ "にすることで、よりシンプルな表が出来上がる:

freq(iris$Species, 
     report.nas = FALSE, 
     totals     = FALSE, 
     cumul      = FALSE, 
     headings   = FALSE)
               Freq       %

      setosa     50   33.33
  versicolor     50   33.33
   virginica     50   33.33

出力はかなり簡略化されているが、構文はそうなっていない。Teslerの複雑さ保存の法則のおかげである!ありがたいことに、st_options()がみんなの好みにアコモデーションしてくれる(パッケージ・オプションのセクションを参照)。

2.3 一度に複数の度数表

データフレーム内のすべての変数の度数表を生成するために、lapply() を使用することができる(初期のバージョンでは必要であった)。しかし、freq() はデータフレームを主引数として受け取るので、これは必要ない。

freq(tobacco)
Variable(s) ignored: age, BMI, cigs.per.day
Frequencies  
tobacco$gender  
Type: Factor  

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          F    489     50.00          50.00     48.90          48.90
          M    489     50.00         100.00     48.90          97.80
       <NA>     22                               2.20         100.00
      Total   1000    100.00         100.00    100.00         100.00

tobacco$age.gr  
Type: Factor  

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
      18-34    258     26.46          26.46     25.80          25.80
      35-50    241     24.72          51.18     24.10          49.90
      51-70    317     32.51          83.69     31.70          81.60
       71 +    159     16.31         100.00     15.90          97.50
       <NA>     25                               2.50         100.00
      Total   1000    100.00         100.00    100.00         100.00

tobacco$smoker  
Type: Factor  

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
        Yes    298     29.80          29.80     29.80          29.80
         No    702     70.20         100.00     70.20         100.00
       <NA>      0                               0.00         100.00
      Total   1000    100.00         100.00    100.00         100.00

tobacco$diseased  
Type: Factor  

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
        Yes    224     22.40          22.40     22.40          22.40
         No    776     77.60         100.00     77.60         100.00
       <NA>      0                               0.00         100.00
      Total   1000    100.00         100.00    100.00         100.00

tobacco$disease  
Type: Character  

                        Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
--------------------- ------ --------- -------------- --------- --------------
               Cancer     34     15.32          15.32      3.40           3.40
          Cholesterol     21      9.46          24.77      2.10           5.50
             Diabetes     14      6.31          31.08      1.40           6.90
            Digestive     12      5.41          36.49      1.20           8.10
              Hearing     14      6.31          42.79      1.40           9.50
                Heart     20      9.01          51.80      2.00          11.50
         Hypertension     36     16.22          68.02      3.60          15.10
          Hypotension     11      4.95          72.97      1.10          16.20
      Musculoskeletal     19      8.56          81.53      1.90          18.10
         Neurological     10      4.50          86.04      1.00          19.10
                Other      2      0.90          86.94      0.20          19.30
            Pulmonary     20      9.01          95.95      2.00          21.30
               Vision      9      4.05         100.00      0.90          22.20
                 <NA>    778                              77.80         100.00
                Total   1000    100.00         100.00    100.00         100.00

tobacco$samp.wgts  
Type: Numeric  

                          Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------------------- ------ --------- -------------- --------- --------------
      0.861423220973783    267     26.70          26.70     26.70          26.70
       1.04417670682731    249     24.90          51.60     24.90          51.60
       1.04938271604938    324     32.40          84.00     32.40          84.00
                 1.0625    160     16.00         100.00     16.00         100.00
                   <NA>      0                               0.00         100.00
                  Total   1000    100.00         100.00    100.00         100.00

結果が乱雑になるのを避けるため、25以上の明確な値を持つ数値列は無視される。このしきい値25は、st_options()を使って変更することができる。例えば、10に変更するには、st_options(freq.ignore.threshold = 10)を使う。

タバコのデータフレームはシミュレートされたデータを含み、パッケージに含まれている。もう1つのシミュレートされたデータフレームは試験だ。どちらもフランス語版(tabagisme, examens)がある。

2.4 度数表のサブセット化(フィルタリング)

rowsパラメータは、度数表をサブセットすることができる:

-行の出現順序でフィルタリングするには、数値ベクトルを使用します; rows = 1:10 は、最初の10個の値の度数のみを表示する。表示されない値の度数を考慮するために、"(Other) "行が自動的に追加されます。 -名前によって行をフィルタリングするには、次のいずれかを使用できる。 -- 保持したいすべての行名を指定する文字ベクトル。 -- 正規表現として使用される単一の文字列(このトピックの詳細については、?regexを参照)。

最も一般的な価値観を示す

orderパラメータとrowsパラメータを組み合わせることで、結果を簡単にフィルタリングして、例えば因子の中で最も一般的な5つの値を表示することができる:

freq(tobacco$disease, 
     order    = "freq",
     rows     = 1:5,
     headings = FALSE)
                     Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
------------------ ------ --------- -------------- --------- --------------
      Hypertension     36     16.22          16.22      3.60           3.60
            Cancer     34     15.32          31.53      3.40           7.00
       Cholesterol     21      9.46          40.99      2.10           9.10
             Heart     20      9.01          50.00      2.00          11.10
         Pulmonary     20      9.01          59.01      2.00          13.10
           (Other)     91     40.99         100.00      9.10          22.20
              <NA>    778                              77.80         100.00
             Total   1000    100.00         100.00    100.00         100.00

freq "の代わりに"-freq "を使えば、順序を逆にし、頻度の低いものから高いものへとランク付けされた結果を得ることができる。

“(Other)”の行に注目。これは自動的に生成される。

2.5 折りたたみ可能なセクション

html 結果を生成する場合、collapse = TRUE 引数を print() または view() / stview()で使用し、折りたたみ可能なセクションを得る。

view(freq(tobacco), collapse = TRUE)

3. クロス集計: ctable()

ctable() は、カテゴリ変数のペアのクロス集計(結合度数)を生成する。

タバコのシミュレートデータフレームを用いて、2つのカテゴリー変数smokerとdiseasedをクロス集計する。

ctable(x = tobacco$smoker, 
       y = tobacco$diseased, 
       prop = "r")   # 行の比率を表示する
Cross-Tabulation, Row Proportions  
smoker * diseased  
Data Frame: tobacco  

-------- ---------- ------------- ------------- ---------------
           diseased           Yes            No           Total
  smoker                                                       
     Yes              125 (41.9%)   173 (58.1%)    298 (100.0%)
      No               99 (14.1%)   603 (85.9%)    702 (100.0%)
   Total              224 (22.4%)   776 (77.6%)   1000 (100.0%)
-------- ---------- ------------- ------------- ---------------

ご覧のように、マークダウンは複数行の表見出しを完全にサポートしていないため、panderはこの特殊なタイプの表を表示するためにできることを行う。より良い結果を得るためには、"render "メソッドが推奨され、次の例で使用される。

3.1 行、列、または全体の割合

行の比率はデフォルトで表示される。列または合計の比率を表示するには、それぞれ prop = "c" または prop = "t" を使用する。比率を完全に省略するには、prop = "n "を使用する。

3.2 最小限のクロス集計

with(tobacco, 
     print(ctable(x = smoker, 
                  y = diseased, 
                  prop     = 'n',
                  totals   = FALSE, 
                  headings = FALSE),
           method = "render")
)
diseased
smoker Yes No
Yes 125 173
No 99 603

Generated by summarytools 1.0.1 (R version 4.3.1)
2023-09-28

3.3 カイ二乗(𝛘2)、オッズ比、リスク比

カイ2乗統計量を表示するには、chisq = TRUEと設定する。2×2の表では、オッズ比とリスク比(相対リスクともいう)を表示するために、それぞれORとRRを使用する。これらはTRUEに設定でき、その場合は95%信頼区間が表示される; 異なる信頼水準を使用するには、例えばOR = .90を使用する。

library(magrittr)
tobacco %$%  # Acts like with(tobacco, ...)
  ctable(x = smoker, y = diseased,
         chisq = TRUE,
         OR    = TRUE,
         RR    = TRUE,
         headings = FALSE) %>%
  print(method = "render")
diseased
smoker Yes No Total
Yes 125 ( 41.9% ) 173 ( 58.1% ) 298 ( 100.0% )
No 99 ( 14.1% ) 603 ( 85.9% ) 702 ( 100.0% )
Total 224 ( 22.4% ) 776 ( 77.6% ) 1000 ( 100.0% )
 Χ2 = 91.7088   df = 1   p = .0000
O.R. (95% C.I.) = 4.40  (3.22 - 6.02)
R.R. (95% C.I.) = 2.97  (2.37 - 3.73)

Generated by summarytools 1.0.1 (R version 4.3.1)
2023-09-28

4. 記述統計量: descr()

descr()は、記述統計量/単変量統計量、すなわち一般的な中心傾向統計量と分散尺度を生成する。後者の場合、数値以外の列はすべて無視され、その旨のメッセージが表示される。

descr(iris)
Non-numerical variable(s) ignored: Species
Descriptive Statistics  
iris  
N: 150  

                    Petal.Length   Petal.Width   Sepal.Length   Sepal.Width
----------------- -------------- ------------- -------------- -------------
             Mean           3.76          1.20           5.84          3.06
          Std.Dev           1.77          0.76           0.83          0.44
              Min           1.00          0.10           4.30          2.00
               Q1           1.60          0.30           5.10          2.80
           Median           4.35          1.30           5.80          3.00
               Q3           5.10          1.80           6.40          3.30
              Max           6.90          2.50           7.90          4.40
              MAD           1.85          1.04           1.04          0.44
              IQR           3.50          1.50           1.30          0.50
               CV           0.47          0.64           0.14          0.14
         Skewness          -0.27         -0.10           0.31          0.31
      SE.Skewness           0.20          0.20           0.20          0.20
         Kurtosis          -1.42         -1.36          -0.61          0.14
          N.Valid         150.00        150.00         150.00        150.00
        Pct.Valid         100.00        100.00         100.00        100.00

変数型メッセージをオフにするには、silent = TRUEを使用する。このオプションはグローバルに設定することも可能で、ここではそれを行うので、このビネットの残りの部分では表示されない。

st_options(descr.silent = TRUE)

5. データフレームの要約: dfSummary()

dfSummary() は、データフレーム内のすべての変数の統計、度数、グラフを含む要約表を作成する。表示される情報は、型(文字、因子、数値、日付)に固有であり、明確な値の数によっても異なる。

RStudio のビューア(または他の IDE やターミナルウィンドウで作業している場合はデフォルトの Web ブラウザ)で結果を表示するには、view() 関数を使用するか、名前が競合する場合はその双子の stview() を使用する:

view(dfSummary(iris))

Viewer欄にHTMLで出力される

View() ではなく view() を使うように注意すること。後者を使用すると、結果がデータ・ビューワに表示されてしまいます。

また、パッケージがロードされる順番にも注意してください。一部のパッケージはview()を再定義してView() を指すようにする。これらのパッケージの後にsummarytoolsをロードすると、それ自身のview()が正しく動作するようになる。そうでなければ、stview()は常に確実な代替手段として存在する。

5.1 R Markdown文書でのdfSummary()の使用

RのMarkdownドキュメントでdfSummary()を使用する場合、マージンのオーバーフローを避けるために、一般的に1つか2つのカラムを除外するのが良いアイデアである。Valid列とMissing列は冗長なので、どちらか一方を削除すればよい。

dfSummary(tobacco, 
          plain.ascii  = FALSE, 
          style        = "grid", 
          graph.magnif = 0.75, 
          valid.col    = FALSE,
          tmp.img.dir  = "/tmp")
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| No | Variable      | Stats / Values           | Freqs (% of Valid)  | Graph                   | Missing |
+====+===============+==========================+=====================+=========================+=========+
| 1  | gender\       | 1\. F\                   | 489 (50.0%)\        | ![](/tmp/ds0001.png)    | 22\     |
|    | [factor]      | 2\. M                    | 489 (50.0%)         |                         | (2.2%)  |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 2  | age\          | Mean (sd) : 49.6 (18.3)\ | 63 distinct values  | ![](/tmp/ds0002.png)    | 25\     |
|    | [numeric]     | min < med < max:\        |                     |                         | (2.5%)  |
|    |               | 18 < 50 < 80\            |                     |                         |         |
|    |               | IQR (CV) : 32 (0.4)      |                     |                         |         |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 3  | age.gr\       | 1\. 18-34\               | 258 (26.5%)\        | ![](/tmp/ds0003.png)    | 25\     |
|    | [factor]      | 2\. 35-50\               | 241 (24.7%)\        |                         | (2.5%)  |
|    |               | 3\. 51-70\               | 317 (32.5%)\        |                         |         |
|    |               | 4\. 71 +                 | 159 (16.3%)         |                         |         |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 4  | BMI\          | Mean (sd) : 25.7 (4.5)\  | 974 distinct values | ![](/tmp/ds0004.png)    | 26\     |
|    | [numeric]     | min < med < max:\        |                     |                         | (2.6%)  |
|    |               | 8.8 < 25.6 < 39.4\       |                     |                         |         |
|    |               | IQR (CV) : 5.7 (0.2)     |                     |                         |         |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 5  | smoker\       | 1\. Yes\                 | 298 (29.8%)\        | ![](/tmp/ds0005.png)    | 0\      |
|    | [factor]      | 2\. No                   | 702 (70.2%)         |                         | (0.0%)  |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 6  | cigs.per.day\ | Mean (sd) : 6.8 (11.9)\  | 37 distinct values  | ![](/tmp/ds0006.png)    | 35\     |
|    | [numeric]     | min < med < max:\        |                     |                         | (3.5%)  |
|    |               | 0 < 0 < 40\              |                     |                         |         |
|    |               | IQR (CV) : 11 (1.8)      |                     |                         |         |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 7  | diseased\     | 1\. Yes\                 | 224 (22.4%)\        | ![](/tmp/ds0007.png)    | 0\      |
|    | [factor]      | 2\. No                   | 776 (77.6%)         |                         | (0.0%)  |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 8  | disease\      | 1\. Hypertension\        | 36 (16.2%)\         | ![](/tmp/ds0008.png)    | 778\    |
|    | [character]   | 2\. Cancer\              | 34 (15.3%)\         |                         | (77.8%) |
|    |               | 3\. Cholesterol\         | 21 ( 9.5%)\         |                         |         |
|    |               | 4\. Heart\               | 20 ( 9.0%)\         |                         |         |
|    |               | 5\. Pulmonary\           | 20 ( 9.0%)\         |                         |         |
|    |               | 6\. Musculoskeletal\     | 19 ( 8.6%)\         |                         |         |
|    |               | 7\. Diabetes\            | 14 ( 6.3%)\         |                         |         |
|    |               | 8\. Hearing\             | 14 ( 6.3%)\         |                         |         |
|    |               | 9\. Digestive\           | 12 ( 5.4%)\         |                         |         |
|    |               | 10\. Hypotension\        | 11 ( 5.0%)\         |                         |         |
|    |               | [ 3 others ]             | 21 ( 9.5%)          |                         |         |
+----+---------------+--------------------------+---------------------+-------------------------+---------+
| 9  | samp.wgts\    | Mean (sd) : 1 (0.1)\     | 0.86!: 267 (26.7%)\ | ![](/tmp/ds0009.png) \  | 0\      |
|    | [numeric]     | min < med < max:\        | 1.04!: 249 (24.9%)\ | \                       | (0.0%)  |
|    |               | 0.9 < 1 < 1.1\           | 1.05!: 324 (32.4%)\ |                         |         |
|    |               | IQR (CV) : 0.2 (0.1)     | 1.06!: 160 (16.0%)\ |                         |         |
|    |               |                          | ! rounded           |                         |         |
+----+---------------+--------------------------+---------------------+-------------------------+---------+

tmp.img.dirパラメータは、htmlレンダリングを除き、R MarkdownドキュメントでdfSummariesを生成する際には必須である。これについての説明は後述する。

警告チャンクオプションをFALSEに設定することで、警告を回避することができる。{r chunk_name, results="asis", warning=FALSE}.

5.2 オプションの統計

この機能はパッケージがリリースされて以来、何度か要望があったものである。バージョン1.0.0で導入され、Stats/Values列にどの統計量を表示するかをコントロールできるようになった。すなわち、IQR (CV)を表示する3行目は、Rで利用可能な統計量を表示するように変更することができる。この機能を使用するには、st_options()を使用して、次のようにdfSummary.custom.1および/またはdfSummary.custom.2を定義し、コードをexpression()内にカプセル化する:

st_options(
  dfSummary.custom.1 = 
    expression(
      paste(
        "Q1 - Q3 :",
        round(
          quantile(column_data, probs = .25, type = 2, 
                   names = FALSE, na.rm = TRUE), digits = 1
        ), " - ",
        round(
          quantile(column_data, probs = .75, type = 2, 
                   names = FALSE, na.rm = TRUE), digits = 1
        )
      )
    )
)

print(
  dfSummary(iris, 
            varnumbers   = FALSE,
            na.col       = FALSE,
            style        = "multiline",
            plain.ascii  = FALSE,
            headings     = FALSE,
            graph.magnif = .8),
  method = "render"
)

もしdfSummary.custom.1の代わりにdfSummary.custom.2を使用していたら、デフォルトのIQR(CV)行の下に4行目が追加されていただろう。

round()の代わりに、内部変数format_number()を使用することが可能であることに注意してほしい。このformat_number()は、指定されたすべての引数(四捨五入の桁数、小数点マーク、千の位マークなど)に従って数値がフォーマットされることを保証する。st_options("round.digits")の値を格納する内部変数round.digitsも使用できる。これは、デフォルトIQR (CV)がどのように定義されているかである - ここでは、最初のカスタムstatをデフォルト値に戻し、その定義を表示する(formatR::tidy_source()は、式をフォーマット/インデントするために使用される):

library(formatR)
st_options(dfSummary.custom.1 = "default")
formatR::tidy_source(
  text   = deparse(st_options("dfSummary.custom.1")),
  indent = 2,
  args.newline = TRUE
)
expression(
  paste(
    paste0(
      trs("iqr"),
      " (", trs("cv"),
      ") : "
    ),
    format_number(
      IQR(column_data, na.rm = TRUE),
      round.digits
    ),
    " (", format_number(
      sd(column_data, na.rm = TRUE)/mean(column_data, na.rm = TRUE),
      round.digits
    ),
    ")", collapse = "", sep = ""
  )
)

このパラメータを使用するすべての関数(Rの基本関数のほとんど)でna.rm = TRUEを指定することを忘れないでほしい。

5.3 その他の注目すべき機能

dfSummary()関数には、以下の機能もある。

  • 見出しセクションで重複レコードの数を報告する
  • UPC/EAN コード (バーコード番号) を検出し、それらに無関係な統計は計算しない。
  • 有効なアドレスと無効なアドレスの割合を合計すると 100%になることに注意して欲しい。重複の割合は独自に計算されるため、棒グラフ (html バージョン)ではこのカテゴリの棒グラフは別の色で表示される。
  • max.tbl.heightパラメータを使用することで、"windowed "結果を表示することができる;これは、分析されたデータフレームが多数の変数を持つ場合に特に便利である; 詳細はvignette("rmarkdown", package = "summarytools")を参照。

5.4 列の除外

ほとんどのカラムは関数のパラメータを使って除外することができるが、以下のシンタックスを使ってカラムを削除することも可能である(結果は示していない):

dfs <- dfSummary(iris)
dfs$Variable <- NULL # これは "Variable "カラムを削除する

6. グループ化された統計: stby()

最適な結果を生成するために、summarytoolsはベースとなるby()関数の独自のバージョンを持っています。これはstby()と呼ばれ、by()と全く同じように使用する:

(iris_stats_by_species <- stby(data      = iris, 
                               INDICES   = iris$Species, 
                               FUN       = descr, 
                               stats     = "common", 
                               transpose = TRUE))
Descriptive Statistics  
iris  
Group: Species = setosa  
N: 50  

                     Mean   Std.Dev    Min   Median    Max   N.Valid   Pct.Valid
------------------ ------ --------- ------ -------- ------ --------- -----------
      Petal.Length   1.46      0.17   1.00     1.50   1.90     50.00      100.00
       Petal.Width   0.25      0.11   0.10     0.20   0.60     50.00      100.00
      Sepal.Length   5.01      0.35   4.30     5.00   5.80     50.00      100.00
       Sepal.Width   3.43      0.38   2.30     3.40   4.40     50.00      100.00

Group: Species = versicolor  
N: 50  

                     Mean   Std.Dev    Min   Median    Max   N.Valid   Pct.Valid
------------------ ------ --------- ------ -------- ------ --------- -----------
      Petal.Length   4.26      0.47   3.00     4.35   5.10     50.00      100.00
       Petal.Width   1.33      0.20   1.00     1.30   1.80     50.00      100.00
      Sepal.Length   5.94      0.52   4.90     5.90   7.00     50.00      100.00
       Sepal.Width   2.77      0.31   2.00     2.80   3.40     50.00      100.00

Group: Species = virginica  
N: 50  

                     Mean   Std.Dev    Min   Median    Max   N.Valid   Pct.Valid
------------------ ------ --------- ------ -------- ------ --------- -----------
      Petal.Length   5.55      0.55   4.50     5.55   6.90     50.00      100.00
       Petal.Width   2.03      0.27   1.40     2.00   2.50     50.00      100.00
      Sepal.Length   6.59      0.64   4.90     6.50   7.90     50.00      100.00
       Sepal.Width   2.97      0.32   2.20     3.00   3.80     50.00      100.00

6.1 stby()によるdescr()の特殊ケース

stby()は、単一変数の分割グループ統計量を生成するために使用すると、一連の1列の表を表示する代わりに、すべてを1つの表にまとめる。

with(tobacco, 
     stby(data    = BMI, 
          INDICES = age.gr, 
          FUN     = descr,
          stats   = c("mean", "sd", "min", "med", "max"))
)
Descriptive Statistics  
BMI by age.gr  
Data Frame: tobacco  
N: 258  

                18-34   35-50   51-70    71 +
------------- ------- ------- ------- -------
         Mean   23.84   25.11   26.91   27.45
      Std.Dev    4.23    4.34    4.26    4.37
          Min    8.83   10.35    9.01   16.36
       Median   24.04   25.11   26.77   27.52
          Max   34.84   39.44   39.21   38.37

6.2 stby() と ctable() の併用

この組み合わせでは構文が少し難しいので、以下に例を示す:

stby(data    = list(x = tobacco$smoker, y = tobacco$diseased), 
     INDICES = tobacco$gender, 
     FUN     = ctable)

# or equivalently
with(tobacco, 
     stby(data    = list(x = smoker, y = diseased), 
          INDICES = gender, 
          FUN     = ctable))
Cross-Tabulation, Row Proportions  
smoker * diseased  
Data Frame: tobacco  
Group: gender = F  

-------- ---------- ------------- ------------- --------------
           diseased           Yes            No          Total
  smoker                                                      
     Yes               62 (42.2%)    85 (57.8%)   147 (100.0%)
      No               49 (14.3%)   293 (85.7%)   342 (100.0%)
   Total              111 (22.7%)   378 (77.3%)   489 (100.0%)
-------- ---------- ------------- ------------- --------------

Group: gender = M  

-------- ---------- ------------- ------------- --------------
           diseased           Yes            No          Total
  smoker                                                      
     Yes               63 (44.1%)    80 (55.9%)   143 (100.0%)
      No               47 (13.6%)   299 (86.4%)   346 (100.0%)
   Total              110 (22.5%)   379 (77.5%)   489 (100.0%)
-------- ---------- ------------- ------------- --------------
Cross-Tabulation, Row Proportions  
smoker * diseased  
Data Frame: tobacco  
Group: gender = F  

-------- ---------- ------------- ------------- --------------
           diseased           Yes            No          Total
  smoker                                                      
     Yes               62 (42.2%)    85 (57.8%)   147 (100.0%)
      No               49 (14.3%)   293 (85.7%)   342 (100.0%)
   Total              111 (22.7%)   378 (77.3%)   489 (100.0%)
-------- ---------- ------------- ------------- --------------

Group: gender = M  

-------- ---------- ------------- ------------- --------------
           diseased           Yes            No          Total
  smoker                                                      
     Yes               63 (44.1%)    80 (55.9%)   143 (100.0%)
      No               47 (13.6%)   299 (86.4%)   346 (100.0%)
   Total              110 (22.5%)   379 (77.5%)   489 (100.0%)
-------- ---------- ------------- ------------- --------------

7. グループ化された統計: group_by()

freq()descr()、またはdfSummary()を使用してグループ化された統計量を作成するために、stby()の代替としてdplyrのgroup_by()を使用することが可能である。構文の違いはさておき、1つの重要な違いは、group_by()が、forcats::fct_explicit_naを使用して因数でNAを明示することを示唆する警告はあるものの、グループ化変数のNA値を有効なカテゴリとみなすことです。このアドバイスに従うと

library(dplyr)
tobacco$gender %<>% forcats::fct_explicit_na()
tobacco %>% 
  group_by(gender) %>% 
  descr(stats = "fivenum")
tobacco  
Group: gender = F  
N: 489  

                 age     BMI   cigs.per.day   samp.wgts
------------ ------- ------- -------------- -----------
         Min   18.00    9.01           0.00        0.86
          Q1   34.00   22.98           0.00        0.86
      Median   50.00   25.87           0.00        1.04
          Q3   66.00   29.48          10.50        1.05
         Max   80.00   39.44          40.00        1.06

Group: gender = M  
N: 489  

                 age     BMI   cigs.per.day   samp.wgts
------------ ------- ------- -------------- -----------
         Min   18.00    8.83           0.00        0.86
          Q1   34.00   22.52           0.00        0.86
      Median   49.50   25.14           0.00        1.04
          Q3   66.00   27.96          11.00        1.05
         Max   80.00   36.76          40.00        1.06

Group: gender = (Missing)  
N: 22  

                 age     BMI   cigs.per.day   samp.wgts
------------ ------- ------- -------------- -----------
         Min   19.00   20.24           0.00        0.86
          Q1   36.00   24.97           0.00        1.04
      Median   55.50   27.16           0.00        1.05
          Q3   64.00   30.23          10.00        1.05
         Max   80.00   32.43          28.00        1.06

8. 整頓されたテーブル : tb()

freq()またはdescr()テーブルを生成する際、tb()関数(tbをtibbleの短縮形と考える)を使用することで、結果を「整然とした」テーブルにすることが可能である。例えば

library(magrittr)
iris %>%
  descr(stats = "common") %>%
  tb()
# A tibble: 4 × 8
  variable      mean    sd   min   med   max n.valid pct.valid
  <chr>        <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>     <dbl>
1 Petal.Length  3.76 1.77    1    4.35   6.9     150       100
2 Petal.Width   1.20 0.762   0.1  1.3    2.5     150       100
3 Sepal.Length  5.84 0.828   4.3  5.8    7.9     150       100
4 Sepal.Width   3.06 0.436   2    3      4.4     150       100
iris$Species %>% 
  freq(cumul = FALSE, report.nas = FALSE) %>% 
  tb()
# A tibble: 3 × 3
  Species     freq   pct
  <fct>      <dbl> <dbl>
1 setosa        50  33.3
2 versicolor    50  33.3
3 virginica     50  33.3

定義により、合計行は整頓された表の一部ではなく、行名は通常の列に変換されます。

rmarkdownを使用してtibblesを表示する場合、knitrのチャンクオプションの結果を'asis'ではなく'markup'に設定する必要がある。

8.1 整然としたスプリット・グループ統計

stby()またはgroup_by()を使用して作成されたリストが、どのように整頓されたtibblesに変換されるかを示すいくつかの例を示します。

grouped_descr <- stby(data    = exams,
                      INDICES = exams$gender, 
                      FUN     = descr,
                      stats   = "common")

grouped_descr %>% tb()
# A tibble: 12 × 9
   gender variable   mean    sd   min   med   max n.valid pct.valid
   <fct>  <chr>     <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>     <dbl>
 1 Girl   economics  72.5  7.79  62.3  70.2  89.6      14      93.3
 2 Girl   english    73.9  9.41  58.3  71.8  93.1      14      93.3
 3 Girl   french     71.1 12.4   44.8  68.4  93.7      14      93.3
 4 Girl   geography  67.3  8.26  50.4  67.3  78.9      15     100  
 5 Girl   history    71.2  9.17  53.9  72.9  86.4      15     100  
 6 Girl   math       73.8  9.03  55.6  74.8  86.3      14      93.3
 7 Boy    economics  75.2  9.40  60.5  71.7  94.2      15     100  
 8 Boy    english    77.8  5.94  69.6  77.6  90.2      15     100  
 9 Boy    french     76.6  8.63  63.2  74.8  94.7      15     100  
10 Boy    geography  73   12.4   47.2  71.2  96.3      14      93.3
11 Boy    history    74.4 11.2   54.4  72.6  93.5      15     100  
12 Boy    math       73.3  9.68  60.5  72.2  93.2      14      93.3

orderパラメータは行の順序を制御する:

grouped_descr %>% tb(order = 2)
# A tibble: 12 × 9
   gender variable   mean    sd   min   med   max n.valid pct.valid
   <fct>  <chr>     <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>     <dbl>
 1 Girl   economics  72.5  7.79  62.3  70.2  89.6      14      93.3
 2 Boy    economics  75.2  9.40  60.5  71.7  94.2      15     100  
 3 Girl   english    73.9  9.41  58.3  71.8  93.1      14      93.3
 4 Boy    english    77.8  5.94  69.6  77.6  90.2      15     100  
 5 Girl   french     71.1 12.4   44.8  68.4  93.7      14      93.3
 6 Boy    french     76.6  8.63  63.2  74.8  94.7      15     100  
 7 Girl   geography  67.3  8.26  50.4  67.3  78.9      15     100  
 8 Boy    geography  73   12.4   47.2  71.2  96.3      14      93.3
 9 Girl   history    71.2  9.17  53.9  72.9  86.4      15     100  
10 Boy    history    74.4 11.2   54.4  72.6  93.5      15     100  
11 Girl   math       73.8  9.03  55.6  74.8  86.3      14      93.3
12 Boy    math       73.3  9.68  60.5  72.2  93.2      14      93.3

order = 3に設定すると、order = 2とまったく同じようにソート変数の順序が変わるが、列の順序も変わる:

grouped_descr %>% tb(order = 3)
# A tibble: 12 × 9
   variable  gender  mean    sd   min   med   max n.valid pct.valid
   <chr>     <fct>  <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>     <dbl>
 1 economics Girl    72.5  7.79  62.3  70.2  89.6      14      93.3
 2 economics Boy     75.2  9.40  60.5  71.7  94.2      15     100  
 3 english   Girl    73.9  9.41  58.3  71.8  93.1      14      93.3
 4 english   Boy     77.8  5.94  69.6  77.6  90.2      15     100  
 5 french    Girl    71.1 12.4   44.8  68.4  93.7      14      93.3
 6 french    Boy     76.6  8.63  63.2  74.8  94.7      15     100  
 7 geography Girl    67.3  8.26  50.4  67.3  78.9      15     100  
 8 geography Boy     73   12.4   47.2  71.2  96.3      14      93.3
 9 history   Girl    71.2  9.17  53.9  72.9  86.4      15     100  
10 history   Boy     74.4 11.2   54.4  72.6  93.5      15     100  
11 math      Girl    73.8  9.03  55.6  74.8  86.3      14      93.3
12 math      Boy     73.3  9.68  60.5  72.2  93.2      14      93.3

詳しくは、?tb参照のこと。

8.2 他のパッケージへの橋渡し

summarytoolsオブジェクトは、formattableやkableExtraのようなテーブルフォーマットに焦点を当てたパッケージとは必ずしも互換性があるとは限りません。しかしながら、tb()は、freq()descr()オブジェクトを、どのパッケージでも扱うことができる単純な表に変換する中間ステップである "ブリッジ "として使用することができる。以下は kableExtra を使った例である:

library(kableExtra)
library(magrittr)
stby(data    = iris, 
     INDICES = iris$Species, 
     FUN     = descr, 
     stats   = "fivenum") %>%
  tb(order = 3) %>%
  kable(format = "html", digits = 2) %>%
  collapse_rows(columns = 1, valign = "top")

9. ファイルへの出力

print()またはview()/stview()でfile引数を使用すると、html、Rmd、md、あるいは単なるプレーン・テキスト(txt)などのファイルに出力を書き出すことができる。ファイル拡張子は、書き出すコンテンツのタイプを決定するためにパッケージが使用する。

view(iris_stats_by_species, file = "~/iris_stats_by_species.html")
view(iris_stats_by_species, file = "~/iris_stats_by_species.md")

PDF文書に関する注意

summarytoolsでPDFファイルを作成する直接的な方法はありません。一つの方法は、htmlファイルを生成し、PandocまたはWKTOpdfを使用してPDFに変換することである(後者はdfSummary()出力でPandocより良い結果を与える)。

もう一つの方法は、PDFを出力フォーマットとしてRmdドキュメントを作成することです。進め方の詳細はvignette("rmarkdown", package = "summarytools")を参照してほしい。

9.1 出力ファイルの追加

append引数はsummarytoolsが生成した既存のファイルに内容を追加することができる。これは、1つのファイルに複数の統計表を含めたい場合に便利だ。これは、Rmdドキュメントを作成する迅速な代替方法である。

10. パッケージオプション

以下のオプションは st_options() でグローバルに設定できる:

10.1 一般オプション

Option name Default Note
style (1) “simple” .Rmd文書では "rmarkdown "に設定する
plain.ascii TRUE .Rmd文書ではFALSEに設定
round.digits (2) 2 表示する小数の数
headings TRUE 以前は "omit.headings"
footnote “default” カスタマイズするか、省略する場合はNAに設定
display.labels TRUE 見出しに変数/データフレームのラベルを表示する
bootstrap.css (3) TRUE html出力ファイルにBootstrap 4 CSSを含める
custom.css NA 独自のCSSファイルへのパス
escape.pipe FALSE いくつかのPandoc変換に便利
char.split (4) 12 列見出しの改行のしきい値
subtitle.emphasis TRUE 見出しフォーマットのコントロール
lang “en” 言語 (常に2文字の小文字)

1 独自の style オプションを持つ dfSummary() には適用さ れない (次の表を参照)。
2 独自の round.digits オプションを持つ ctable() には適用されない(次の表を参照)。
3 ShinyアプリではFALSEに設定。
4 descr()とctable()のhtml出力にのみ影響する。

10.2 機能別オプション

Option name Default Note
freq.cumul TRUE freq() で累積比率を表示する
freq.totals TRUE freq() で合計行を表示する
freq.report.nas TRUE 行と "有効な" 列を表示
freq.ignore.threshold (1) 25 無視するバーを決定するために使用
freq.silent FALSE コンソール・メッセージを隠す
ctable.prop “r” デフォルトで行の比率を表示
ctable.totals TRUE 限界合計の表示
ctable.round.digits 1 ctable() で表示する小数の数。
descr.stats “all” "fivenum"、"common" または統計のベクトル
descr.transpose FALSE 統計量を行ではなく列で表示する
descr.silent FALSE コンソール・メッセージの非表示
dfSummary.style “multiline” 代替として "grid "に設定可能
dfSummary.varnumbers TRUE 1列目に変数番号を表示
dfSummary.labels.col TRUE 変数ラベルがある場合は表示
dfSummary.graph.col TRUE グラフの表示
dfSummary.valid.col TRUE 有効列を出力に含める
dfSummary.na.col TRUE 出力に欠落列を含める
dfSummary.graph.magnif 1 棒グラフとヒストグラムの拡大率
dfSummary.silent FALSE コンソールメッセージの非表示
tmp.img.dir (2) NA 一時画像を保存するディレクト
use.x11 (3) TRUE Base64エンコードされたグラフの作成を許可する

1 詳細は2.3節を参照
2 dfSummary()にのみ適用されます。
3 テキストのみの環境では FALSE に設定する

st_options()                      # すべてのグローバルオプション値を表示する
st_options('round.digits')        # 特定のオプションの値を表示する
st_options(style = 'rmarkdown',   # 1つまたは複数のオプションの値を設定する
           footnote = NA)         # すべてのhtml出力の脚注をオフにする

11. フォーマット属性

summarytoolsオブジェクトが作成されるとき、そのフォーマット属性はその中に格納さ れる。しかし、print() または view() を使用するとき、それらのほとんどをオーバーライドできる。

11.1 関数固有の引数のオーバーライド

以下の表は、print() や view() で書式属性を上書きするために使用できる引数を示しています。ベースRのformat()関数の引数も使用することができる(ここにはリストされていないが)。

  1. pander オプション

11.2 見出しの内容を上書きする

見出し部分に表示される情報を変更するには、print() または view() で以下の引数を使用する:


次の例では、freq()オブジェクトを作成して表示し、それを再び表示して、今度はそのフォーマット属性の3つと見出し属性の1つをオーバーライドする。

(age_stats <- freq(tobacco$age.gr)) 
Frequencies  
tobacco$age.gr  
Type: Factor  

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
      18-34    258     26.46          26.46     25.80          25.80
      35-50    241     24.72          51.18     24.10          49.90
      51-70    317     32.51          83.69     31.70          81.60
       71 +    159     16.31         100.00     15.90          97.50
       <NA>     25                               2.50         100.00
      Total   1000    100.00         100.00    100.00         100.00
print(age_stats,
      report.nas     = FALSE, 
      totals         = FALSE, 
      display.type   = FALSE,
      Variable.label = "Age Group")
Frequencies  
tobacco$age.gr  
Label: Age Group  

              Freq       %   % Cum.
----------- ------ ------- --------
      18-34    258   26.46    26.46
      35-50    241   24.72    51.18
      51-70    317   32.51    83.69
       71 +    159   16.31   100.00

11.3 パラメータ/オプションの優先順位

  • print() または view() パラメータが優先される(オーバーライド機能)。
  • freq() / ctable() / descr() / dfSummary() パラメータは 2 番目。
  • st_options() で設定されたグローバル・オプションは 3 番目に優先され、デフォルトとして動作する。

様々なパラメータ値の評価ロジックをまとめると、以下のようになる:

  • もし引数が関数呼び出しの中で明示的に供給されるなら、その引数はパラメータに格納されている値よりも優先される(格納されている値とは、summarytoolsのグローバルオプションリストに格納されている値と同様に、コア関数を使用するときにオブジェクトの属性に書き込まれる値である)。

  • コア関数とprintまたはview関数の両方が同時に呼び出され、相反するパラメータ値を持つ場合、print/viewが優先さ れる(それらは常に議論に勝つ!)。

  • 関数呼び出しでパラメータ値が見つからない場合、保存されているデフォルト値(st_options()で変更するか、パッケージ読み込み時のまま)が適用さ れる。

12. 見た目の微調整 : CSS

htmlレポートを作成する場合、デフォルトでBootstrapのCSSとsummarytools.cssの両方が含まれます。htmlコンテンツの見た目をよりコントロールするために、カスタムCSSファイルにクラス定義を追加することも可能だ。


dfSummary()を含む単純なhtmlレポートに、非常に小さなフォント・サイズを使用する必要があります。このために、以下のクラス定義を含む.cssファイル(任意の名前)を作成する:

.tiny-text {
  font-size: 8px;
}

次に、print()のcustom.css引数を使って、新しく作成したCSSファイルの場所を指定する(結果は示していない):

print(dfSummary(tobacco),
      custom.css    = 'path/to/custom.css', 
      table.classes = 'tiny-text',
      file          = "tiny-tobacco-dfSummary.html")

13. Shiny Apps

Shynyアプリにsummarytools関数をうまく組み込む、

  • htmlレンダリングを使用する
  • アプリのレイアウトとの相互作用を避けるために、bootstrap.css = FALSEを設定する。
  • 問題が発生した場合に備えて、headings = FALSEを設定する。
  • graph.magnifパラメータまたはdfSummary.graph.magnifグローバル・オプションでグラフ・サイズを調整する。
  • dfSummary()の表が広すぎる場合は、1列か2列を省略する(valid.colとvarnumbersなど)。
  • それでも満足のいく結果が得られない場合は、col.widths パラメータで列幅を手動で設定する。
  • col.widthsやgraph.magnigがうまくいかないようであれば、dfSummary()ではなく、print()のパラメータとして使用してみる。
print(dfSummary(somedata, 
                varnumbers   = FALSE, 
                valid.col    = FALSE, 
                graph.magnif = 0.8), 
      method   = 'render',
      headings = FALSE,
      bootstrap.css = FALSE)

14. Rマークダウンにおけるグラフ

マークダウン・スタイルを使用するRmdドキュメントでdfSummary()を使用する場合(htmlレンダリングとは対照的)、pngグラフを適切に表示するには3つの要素が必要:

1 - plain.asciiをFALSEに設定する。 2 - styleを "grid "に設定する。 3 - tmp.img.dirが定義され、幅が最大5文字であること。

バージョン0.9.9では、method = "render "を使用する場合、tmp.img.dirを設定する必要はなくなり、NAのままにしておくことができます。下図のように、一過性のマークダウン・テーブルを作成する場合にのみ定義する必要がある。レンダリングされる列の幅は、画像自体の幅ではなく、セル内の文字数によって決定されるためだ:

+---------------+--------|----------------------+---------+
| Variable      | stats  |  Graph               | Valid   |
+===============+========|======================+=========+
| age\          |  ...   | ![](/tmp/ds0001.png) | 978\    |
| [numeric]     |  ...   |                      | (97.8%) |
+---------------+--------+----------------------+---------+

CRANポリシーは、ユーザー・ディレクトリやRのテンポラリ・ゾーンの外側にコンテンツを書き込むことに関しては(正当な理由があって)本当に厳しい。そのため、ユーザーはこの一時的な場所を自分で設定する必要があり、Rのあらかじめ定義された一時的なゾーンの外にコンテンツを書き込むことに同意することになる。

Mac OSLinuxでは、"/tmp "を使うことは非常に理にかなっている。Windowsでは、そのような便利なディレクトリはないので、絶対パス("/tmp")か相対パス("img"、または単に".")を選ぶ必要がある。

15. 言語と用語のカスタマイズ

Rコミュニティの努力により、英語(デフォルト)の他に以下の言語が使用できる:

フランス語 (fr) ポルトガル語 (pt) ロシア語 (ru) スペイン語 (es) トルコ語 (tr) 言語を切り替えるには

st_options(lang = "fr")

コアファンクションからの出力はすべてこの言語を使用する:

freq(iris$Species)
Tableau de fréquences  
iris$Species  
Type: Facteur  

                   Fréq.   % Valide   % Valide cum.   % Total   % Total cum.
---------------- -------- ---------- --------------- --------- --------------
          setosa       50      33.33           33.33     33.33          33.33
      versicolor       50      33.33           66.67     33.33          66.67
       virginica       50      33.33          100.00     33.33         100.00
            <NA>        0                                 0.00         100.00
           Total      150     100.00          100.00    100.00         100.00

15.1 非UTF-8ロケール

ほとんどのWindowsシステムでは、文字セットがシステムのデフォルトロケールに含まれていない場合、ロケール設定のLC_CTYPE要素を変更する必要があります。例えば、"latin1 "環境でロシア語で良い結果を得るには、以下の設定を使用する:

Sys.setlocale("LC_CTYPE", "russian")
st_options(lang = 'ru')

デフォルト設定に戻すには

Sys.setlocale("LC_CTYPE", "")
st_options(lang = "en")

15.2 カスタム用語の定義と使用

use_custom_lang()関数を使用すると、独自の翻訳セットやパーソナライズされた用語を追加することができる。これを実現するには、csvテンプレートを取得し、1つ、多数、または+/- 70の用語のすべてをカスタマイズし、編集したcsvテンプレートへのパスを唯一の引数として与えて、use_custom_lang()を呼び出します。このようなカスタム言語設定は、Rのセッションをまたいでも保持されないことに注意して欲しい。つまり、このcsvファイルを常に手元に置いておく必要がある。

15.3 特定のキーワードだけを定義する

define_keywords()を使用すると、1つまたは少数の用語を簡単に変更できる。たとえば、freq() テーブルのタイトル行で、"Freq" ではなく "N" や "Count" を使用したい場合などである。あるいは、表のタイトルを見出しセクションとして使用するドキュメントを生成したいかもしれない。

この場合、define_keywords() を呼び出し、変更したい用語 (定義済みの変数に格納できる) を入力する。ここでは、freq.titleとfreq.titleを変更する:

section_title <- "**Species of Iris**"
define_keywords(title.freq = section_title,
                freq = "N")
freq(iris$Species)
**Species of Iris**  
iris$Species  
Type: Factor  

                     N   % Valid   % Valid Cum.   % Total   % Total Cum.
---------------- ----- --------- -------------- --------- --------------
          setosa    50     33.33          33.33     33.33          33.33
      versicolor    50     33.33          66.67     33.33          66.67
       virginica    50     33.33         100.00     33.33         100.00
            <NA>     0                               0.00         100.00
           Total   150    100.00         100.00    100.00         100.00

define_keywords()を引数なしで呼び出すと、グラフィカル・デバイスをサポートしているシステム(つまり大多数)では、ウィンドウが表示され、そこからすべての用語を編集することができる。

編集ウィンドウを閉じた後、ダイアログボックスが新しく作成したカスタム言語をcsvファイルに保存するオプションを提供します(いくつかのキーワードを変更しただけでも、パッケージは用語を全体として考慮します)。後でuse_custom_lang("path-to-custom-language-file.csv")を呼び出すことで、カスタム言語ファイルをメモリに再ロードすることができる。

パッケージ内のカスタマイズ可能なすべての用語のリストについては、?define_keywordsを参照してほしい。

すべての変更を元に戻すには、単純にst_options(lang = "en")を使用する。

15.4 見出しのパワー調整

print()関数に引数を追加することで、見出しをさらにカスタマイズすることができる。ここでは、Variableの値を上書きするために空の文字列を使用している;これは見出しの2行目を完全に消してしまう。

define_keywords(title.freq = "Types and Counts, Iris Flowers")
print(
  freq(iris$Species,
       display.type = FALSE), # Variable type won't be displayed...
  Variable = ""               # and neither will the variable name
  ) 
Types and Counts, Iris Flowers  

                     N   % Valid   % Valid Cum.   % Total   % Total Cum.
---------------- ----- --------- -------------- --------- --------------
          setosa    50     33.33          33.33     33.33          33.33
      versicolor    50     33.33          66.67     33.33          66.67
       virginica    50     33.33         100.00     33.33         100.00
            <NA>     0                               0.00         100.00
           Total   150    100.00         100.00    100.00         100.00

16. ヴィネット設定

このビネットがどのように設定されているかを知ることは、R Markdown文書でsummarytoolsを使い始めるのに役立つ。

16.1 YAMLセクション

出力要素が重要だ

---
output: 
 rmarkdown::html_vignette: 
   css:
   - !expr system.file("rmarkdown/templates/html_vignette/resources/vignette.css", package = "rmarkdown")
---

16.2 設定チャンク

library(knitr)
opts_chunk$set(results = 'asis',     # Can also be set at chunk level
              comment = NA,
              prompt  = FALSE,
              cache   = FALSE)
library(summarytools)
st_options(plain.ascii = FALSE,       # Always use in Rmd documents
           style       = "rmarkdown", # Always use in Rmd documents
           subtitle.emphasis = FALSE) # Improves layout w/ some themes

16.3 summarytoolsのCSSをインクルードする

必要なCSSは、file引数でprint()またはview()を使用して作成されたhtmlファイルに自動的に追加さ れる。しかし、R Markdownドキュメントでは、YAMLヘッダーの直後のセットアップチャンクで(またはknitrとsummarytoolsオプションを指定する最初のセットアップチャンクの後で)明示的に行う必要がある:




st_css(main = TRUE, global = TRUE)




17. 結論

このパッケージには何の保証もない。現在進行形であり、フィードバックはいつでも歓迎する。バグを見つけたり機能要望を出したい場合は、GitHubにissueを登録して欲しい。

代理回答の信頼性についての研究者の見解

社会調査では代理回答が起こっていることが分かっている。

www.jstage.jst.go.jp

  • 雄太郎花田, 海人仲田, & 寧佐藤. (2014). Ce5-3 郵送調査における代理回答分析と有効回答率の再定義: 内閣府の郵送世論調査の結果から(一般セッション 社会(2)). 日本行動計量学会大会抄録集, 42, 384–385. https://doi.org/10.20742/pbsj.42.0_384

郵送調査は高い有効回収率が実現できる手法として、自治体やマスコミでの活用が広がっている。反面、調査対象者と顔を合わせて確認ができる調査員調査とは異なり、代理者による回答を直接的に防ぐことができないことが問題視されることもある

一番思い浮かべやすいのは、夫が働き、妻が専業主婦の過程で、夫を対象とした郵送調査が来た時に、夫が忙しいからと妻が代わりに回答して、郵送するといったパターンであろうか。

代理回答をしたと自己申告した票が2325票のうちの61票と2.6%あった。住基台帳からサンプルを選んであるため、性別・年齢はあらかじめ判明している。性別と年齢が間違っている票が43票、無回答が30票あった。無回答は匿名性の担保になると思った可能性があるため、43票の回答が本人が行ったが疑わしいと言える。

代理回答の信頼性

Alwinの論文より。

  • Alwin, D. F., & Duane, F. (2010). How good is survey measurement? Assessing the reliability and validity of survey measures. Handbook of Survey Research, 2, 405–434.

Handbook of Survey Research

Handbook of Survey Research

  • Emerald Group Publishing
Amazon

調査データに対してよく言われる不満のひとつは、調査データが求める自己報告が本質的に信頼できないというものである。自己報告には限界があるが、回答者は他者に関する情報よりも自分自身に関する情報を報告する方が優れている。代理報告は自己報告と同程度に優れていると主張する人もいる(Sudman, Bradburn, & Schwarz, 1996, p.243)が、我々の最も優れたエビデンスによれば、平均して自己報告の方が代理報告よりも信頼性が高い。例えば、Alwin (2007, pp. 152-153)は、同じ回答者が自己と他者について報告した、同じまたは類似の内容を含む小さな変数の対照比較を報告している。その結果、自己報告は代理情報提供者による二次報告よりも信頼できる傾向があるという結論が補強された。とはいえ、含まれる事実情報の種類については、代理報告もそれほど信頼できないものではなかった。(p.423)
One complaint often lodged against survey data is that the self-reports they request are intrinsically unreliable. While self-reports have limitations, respondents are better at reporting information on themselves than about others. Some have argued that proxy reports are as good as self-reports (Sudman, Bradburn, & Schwarz, 1996, p. 243), but our best evidence indicates that on average self-reports are more reliable than proxy reports. Alwin (2007, pp. 152–153), for example, reported a controlled comparison of a small set of variables involving the same or similar content in which the same respondents reported for self and others. Its results reinforced the conclusion that self-reports tend to be more reliable than second-hand reports by proxy informants. Nonetheless, for the types of factual information included, the proxy reports were not all that unreliable.

情報の種類によって代理回答の信頼性は変わるようだ。

  • Sudman, S., Bradburn, M. M., & Schwarz, N. (1996). Thinking about answers: The application of cognitive processes to survey methodology. San Francisco, CA: Jossey-Bass.

Alwin(2007)の方をみてみよう。

  • Alwin, D. F. (2007). Margins of error: A study of reliability in survey measurement. Wiley-Interscience.

7.2 代理報告
調査研究において、代理人によるデータ収集はよくあることである。調査の回答者は、配偶者や子供、時には友人や同僚など、他の人について質問されることが多い。他人の特性について報告するプロセスと、より一般的な自己報告法には違いがあるため、代理報告と自己報告では測定誤差の性質が異なることが予想される(Blair, Menon and Bickart, 1991)。しかし、最近の証拠によると、「多くの行動やいくつかの態度についてさえ、代理報告は自己報告よりも有意に正確性が低いわけではない」ことが示唆されている(Sudman, Bradburn and Schwarz, 1996, p. 243)。というのも、すべてのケースで自己報告が得られることはめったになく、代理人がその人に関する唯一の情報源(たとえば世帯主)であることが多いからである。しかし、このトピックに関する広範な文献では、代理人による報告の質に関して一貫性のない結果が示されていることに注意すべきである(Moore, 1988参照)。
我々の研究は、測定の信頼性に関して、自己報告か代理報告か、また測定誤差の問題に取り組むことができる。例えば、回答者の学歴と配偶者の学歴のように、質問の内容が同一で、同じ人がそれぞれのケースで、通常は質問票の異なる場所で報告している場合である。これらの結果を表7.2に示す。事実上すべてのケースで、自己報告式測定は代理報告式測定よりも信頼性が高い。6つの変数の平均では、その差は約0.90であり、リストワイズ・プレゼンテーションの結果では30であった。それにもかかわらず、この差の統計的検定(6つのケースに基づく)はわずかに有意であり(p< .02)、ケースの数が少ないことを考えると注目に値する。アリソンモデルを用いた場合、得られた差はわずかに小さく、差の系統的な性質はあまり明らかではないが、全体として2つの結果はほぼ一致している。このことから、代理報告は自己報告よりも有意に信頼性が低いと結論づけられる。しかし、この結論のより広範な推論的根拠を得るためには、さらなる研究が必要である。(pp. 152-153)
7.2 PROXY REPORTS
Gathering data by proxy in survey research is commonplace. Respondents to surveys are often asked questions about other people, including their spouses and children, and sometimes, their friends and coworkers. Because of differences in the process of reporting about the characteristics of others and the more common self-report method, one would expect that the nature of measurement errors might be different for proxy vs. self-reports(Blair, Menon and Bickart, 1991).Recent evidence,however, suggests that “for many behaviors and even for some attitudes, proxy reports are not significantly less accurate than self-reports (Sudman, Bradburn and Schwarz, 1996, p. 243). If true, this is an encouraging result because it is seldom possible to obtain self-reports in all cases and the proxy is frequently the only source of information on the person in question, e.g., the head of household. We should note, however, that the broader literature on this topic presents inconsistent results regarding the quality of reporting by proxies (see Moore, 1988).
Our study can address the issue of self- versus proxy reporting and measurement error with respect to the reliability of measurement. We compared six measures involving self-reports with six measures involving proxy reports (in this case reports on spouses) where the content of the questions was identical,e.g., respondent’s education versus spouse’s education, and the same person is reporting in each case, usually at a different place in the questionnaire.These results are shown in Table 7.2. In virtually every case the self-report measure is higher in reliability than the proxy-report measures. On average across the six variables, the difference is approximately .90 versus 30using the listwise-present results, which, if based on a larger sample of measures, would certainly be considered a substantively meaningful difference. The statistical test of this difference (based on six cases) is nonetheless marginally significant (p< .02),which is remarkablegiven the small number of cases. Thedifferenceobtained is slightly smaller using the Allison model, and the systematic nature of the difference is less apparent, but overall the two sets of results are largely consistent. This leads us to conclude that proxy reports are significantly less reliable than are self-report measures, which is consistent with both theory and experience. Further research is necessary, however, to create a broader inferential basis for this conclusion.

社会調査に必要なサンプルサイズ

ベストの本を読んでいてみつけた記述。

標本は、正確な結果を担保するにはあまりにも小さく見えるかもしれない。実際、多数の人が標本調査を疑う。全国調査では、よく一〇〇〇人から一五〇〇人程度の標本が用いられる。 「この数少ない人たちが、どうやって三億人以上のアメリカ人を正確に代表できるというのだ」と、不信を抱く人たちもいる。この疑問はもっともらしく見えるが、強調点が間違っている。 比較的小さな標本、たとえば一〇〇〇人程度の人からは、それが代表標本であれば(すなわち、ある人が標本として選ばれる確率が、他の人の確率と同じであるという前提が満たされていれば) かなり正確な結果が得られるのである(統計学者の計算によると、三億人規模の母集団に対して毎回ランダムに抽出した一〇〇〇人規模の代表標本を用いた調査を繰り返し行うと、一〇〇回に九五回の割合で真の値(母集団の支持率)が、 それぞれの標本平均のプラスマイナス3%の区間の中に含まれる。さらに標本を大きくしたところで正確さはさほど変わらない。たとえば一〇倍の費用を費やして一万人規模の標本に対する調査を行うとする。その場合、誤差はプラスマイナス3%ではなくプラスマイナス1%になる [de Vaus, 1986])。 p.182

de Vausの本は下記の本のようだ。

  • de Vaus, D. A. (1986). Surveys in social research. London: Allen & Unwin.

現在は第5版。
この記述を探してみたが、見つからなかった。
1986年の本は紙でかつ古書で手に入れる必要があるので、手に入れるのに多少時間はかかる。そのためde Vausがどの研究を引用したのかは今のところ不明である。

de Vausの第5版を斜め読みした程度ではあるが、今まで見たどの本よりもずばぬけてよい出来のように感じた。社会調査士関係で社会調査をテーマにした本はでているが、de Vausの本で勉強するべきなのではないか、と思った。