Misc

業務中に発生した、研究デザインや解析に関する雑多なやり取りの記録。

■　抄読会でのコメント　初回発症しか取らないのは？　に付いて補足

　Coxモデルなどの追跡データを扱う解析では、発症までの時間（生存時間）や、生存時間から導かれる関数をアウトカムとしてモデル化することが多いです。そのため臨床研究ではイベントの初回発症までの時間を記録し、その時点で打ち切りにするという処理が一般的です。よく見るカプランマイヤー曲線もそういうものですね。

　別コメントでも有りましたが、一回イベントを起こすケースと、イベントを繰り返すケースではそもそも患者特性が違う可能性もあるので、その面でも初回のみに揃えることは利点でもあります。

　繰り返しイベント・再発を含めて扱いたい場合には、さきほど話題に出た frailty model や、multistate model, Andersen-Gill model, PWP model といった（マニアックな）方向に足を踏み入れる必要があります。心血管病、がん領域、発作や増悪をアウトカムとする場合には再発を扱うことがあります。https://pubmed.ncbi.nlm.nih.gov/29314447/　個人的な印象では心不全の研究が多いですね。ただ、手法のとっつきにくさに加えて、解釈が難しかったりで実りが大きくないので、今後もあまり広まらないかも。 

■　一般化傾向スコア(Generalized propensity score, by Hirano & Imbens)

https://www.math.mcgill.ca/dstephens/PSMMA/Articles/HIrano-Imbens-2004.pdf

https://www.nber.org/system/files/working_papers/t0237/t0237.pdf

解析自体は単純で、多項ロジスティックで各群(割付)に入る確率を推定するだけです。どちらかといえば、通常診療のデータベース研究で、処方薬の処方量（少用量、中用量、大用量）ごとの背景分布を揃えたい、といったモチベーションで使う解析ですね。 

■　因果推論から二重頑健推定まで

■　因果推論の目的

　臨床研究の目的は、曝露とアウトカムの関連を明らかにすることです。このとき概念的な"真の関連"を有限のサンプルから予想するのが推定です。ランダム化比較試験(RCT)のエビデンスが重視されるのは、介入群と統制群で見られたアウトカムの差が、平均因果効果（"真の関連"の集団平均）と一致する＝バイアスのない推定値（不偏推定値）が得られるからです。

　一方観察研究では、曝露とアウトカムの関連に影響を与える因子、すなわち交絡因子を調整（モデリング）などの手法でコントロールしますが、「コントロールした交絡因子が本当に正しい交絡因子のセットであるか」を確かめられないので通常のモデリングから得られた結果は"真の関連"から幾分かズレた（＝バイアスのかかった）推定値となります。

　最近流行りの因果推論は、観察研究のデータを用いながら、可能な限りバイアスを減らして推定値を真値に近づけようとする一連の試みを指します。

■　傾向スコア

　傾向スコアは、因果推論の一手法としてRosenbaumとRubinによって開発・実装されたものです。

　RCTと観察研究の大きな違いは、前者は曝露因子（割り当て）が共変量の影響を受けない一方で、後者は曝露因子と共変量の間に相関があります。曝露因子と共変量の間に相関があると、推定値（回帰モデルであれば回帰係数）と平均因果効果が一致しなくなります。そこで、共変量と曝露因子が独立になるような関数b(x)を考えます。傾向スコアは、この性質をもつ関数のひとつです。傾向スコアは「共変量を条件づけたときに、曝露X=1となる確率変数e」と定義されます。ざっくり言い換えると共変量を右辺、曝露因子を左辺において、共変量から曝露因子を予測する回帰式がこれにあたります。独立変数に対して傾向スコアで条件付けることで、独立変数と共変量の相関を打ち消すことができます。すると、推定値が平均因果効果と一致する（＝不偏推定値）ようになります。

　傾向スコアはもともとマッチングや層別化の手法としてよく利用されましたが、サンプルサイズや誤差を犠牲にしないために、傾向スコアの逆数を用いた重み付け（IPW推定量）が回帰モデルのなかで頻用されるようになっています。

■　二重頑健推定(Doubly robust estimation)

　傾向スコアの問題点として、傾向スコアを作るときに回帰式にふくめた共変量のセットが正しくないと、不偏推定量が得られない可能性があります。そこで傾向スコア（というかIPW推定量）に加えて、共変量とアウトカムの回帰式を組み合わせた二重頑健推定量が考案されました。

　二重頑健推定量の性質として、傾向スコア（つまり"共変量→曝露因子"の回帰モデル）か、アウトカムモデル（"共変量→アウトカム"の回帰モデル）のいずれかが正しく指定されていれば、不偏推定量が得られるという特徴があります。傾向スコアがうまく作られていない場合の保険がかけられている感じです。

　（※ ただし傾向スコアもアウトカムモデルも、共変量のセットが正しいかどうかを精査するのは困難なので、本当に保険として役立つのかは個人的には疑問ですが・・・）

Further reading:

　https://qiita.com/saltcooky/items/85124a1c5a360f822ff6

　https://www.slideshare.net/shuyo/observational-study-causality-52737410

　星野崇宏：調査観察データの統計科学 (岩波書店)

■　相加モデル(Additive model)と相乗モデル(Multiplicative model)

曝露因子がアウトカムに与える影響の大きさを考えるときに、アウトカムの差で評価するAdditive model （相加的モデル）とアウトカムの比で評価する Multiplicative interaction（相乗的モデル）があります。

　たとえば　10年間の外傷発生割合を 多量飲酒なし＝2%　多量飲酒あり＝ 6%とすると、多量飲酒によって外傷発生が2%→6%と4ポイント増えます（リスク差 4%）。比は3倍です。

　コホート研究ではロジスティックモデル・ポワソンモデル・ハザードモデルなど、アウトカムの比（オッズ比やリスク比）を考えることが多いので、相乗的効果が議論されることが多いです。ただし、公衆衛生の観点では集団全体に対しての寄与の大きさを考えるためには、相加的効果も重視されます。コホート研究でも、アウトカムの発生率を連続量とみなして重回帰モデルを適用することで、共変量調整をした発生率の差（リスク差）を求めることができます。

Further reading:

https://www.krsk-phs.com/entry/interaction 

■　変化量のモデリング

変化量の解析の場合、

(1) Y_follow = X + 共変量 + Y_baseline 

(2) Y_follow - Y_baseline = X + 共変量

(3) Y_follow  -Y_baseline= X + 共変量 + Y_baseline 

モデル(1)は、特にRCTのような場面ではYを正しく推定(不偏推定量)できるようです。一方で観察研究の場合、ベースライン値とフォローアップ値に因果関係が仮定できるなら(1)がいいとか、時間に依存しない未測定の交絡が存在するときは差をとる(2,3)のがが良いとか、いろいろディスカッションがあります。

(https://stats.stackexchange.com/questions/15713/is-it-valid-to-include-a-baseline-measure-as-control-variable-when-testing-the-e

　http://www.statisticalhorizons.com/wp-content/uploads/Allison.SM90.pdf) 

■　乳がんの予後予測モデル

https://www.bmj.com/content/381/bmj-2022-073800/

Development and internal-external validation of statistical and machine learning models for breast cancer prognostication: cohort study

ラストオーサーのJ Hippisley-Coxや G Collinsは、予測モデル開発やその方法論の界隈では有名な人たちです。

対象者：UKのGeneral practitionerデータベース(QResearch)

アウトカム：乳がん関連死亡

解析：　Cox モデル、競合リスクモデル、XGBoost, ニューラルネットワーク(NN)の4種類比較

■　Coxモデルと競合リスクモデルでは、emergency presentation の推定値が変わったり、CKDや降圧薬が選ばれなくなったりしていて、競合リスクモデルでは、たしかにがん以外の死亡の影響を除去しているようです。ただ、単純なC統計量やキャリブレーションでは Coxモデルのほうが良く(Table 2)、競合リスクを加味したdecision curve analysisだけで競合リスクモデルのほうが良い（Fig 8）など、競合リスクを考慮するメリットはあまり大きくは無いようです。

■　XGBoost は決定木(regression & classification tree)を応用した手法で、分類・予測を行う教師あり学習で代表的な手法です。特徴として、決定木を[直列に]つなげるようなステップを踏みます。集団から小サンプルを抽出して決定木を当てはめる→次の小サンプルを抽出して、前のステップの決定木(のパラメータ)を修正する…という方法を繰り返していくことでパラメータを収束させていきます。似たような手法としてランダムフォレストがありますが、ランダムフォレストは、決定木を[並列に]たくさん作って、それをメタ解析みたく統合する方法です。→　https://smart-hint.com/ml/gbdt-image/　　　

　XGBoostでは、特に中～高リスク者において予測確率を過大に推定する傾向が見られました(Fig6の上のキャリブレーションプロットで、お辞儀しているのが 実測確率＜予測確率 になっている部分)。XGBoostで選択された変数の重要度（Supplementary table 9）を見てみると、 がんのStage 4、がんのStage1 、Emergency presentation・・・と続いています。このことから、ステージが進行していたら、有無を言わせず高い予測確率を返すモデルになっていると予想されます。特にXGBoostは決定木をベースとしており、Coxモデルのように変数を線形結合していないので、決定木の最初に判定につかう変数の寄与が大きくなってしまい、結果的に 進行がん患者での死亡率の過大推定が起きていると思われます。

■　NNは、入力を線形変換して出力する計算機＝"ノード"が、たくさんネットワーク状につながったモデルです。神経を模しているとか何とかという話はどこかで聞いたことがあると思います。https://jp.mathworks.com/discovery/neural-network.html　　本研究でのNNの中身についてはパラメータが見えないのでよくわかりませんが、XGBoost と同じように、おそらくがんのステージの寄与が大きいと思います。本研究では、比較的発症リスクが低めの人では過大評価、高めの人では過小評価される傾向がみられていますね(Fig6)。Fig 6の下段のPredicted Risk を見るとPredicted risk が 0.05, 0.25のあたりで飛び出たバーが見えます。あくまでも予想ですが、ステージ1, 4 がそれぞれに対応しているような気がします。 Stage 4の多くの人が 0.25あたりの特定の予測値に一律で決められてしまうような、極端な挙動をしているのかもしれません。実際、ステージ別に評価した場合 (Supplementary Fig 6)に、Stage 1とStage 4ではどちらもNNの精度が飛び抜けて悪いので、Stage 1、Stage4の人では正しく判別ができていないようです。

■　なお、集団全体での予測値と観測値の差を表すCalibration-in-the-largeや、キャリブレーションプロットの傾きを表す calibration slope は、線形モデルよりも機械学習モデルのほうがより1に近く、キャリブレーションが取れているという結果でした。言い換えると、集団全体としての絶対リスクは、機械学習モデルでもずれていなかったといえます。

以上を踏まえると、がんのステージが死亡率に与える寄与が大きいため、機械学習モデルではステージの変数の影響が強いがために、予測値の上振れ・下振れが起きてしまって予測精度が悪かったものと思われました。

ただ、著者らが It should be noted ...と書いているように、あくまでも「QResearch」というGPの日常診療データベース（所謂RWD）に当てはめた場合の結果なので、異なるデータベースでは結果は違うかもしれません。たとえば院内がん登録のように、がん治療の情報をもう少し増やすことで、相対的にステージの変数の寄与が小さくなって、機械学習モデルが優位になるかもしれません。あるいは単施設で腫瘍の画像データを使いたければ、単純な線形モデルは使えません。目的に応じて、適した手法を検証していくのが必要のようです。 

■　βエラー

検定や推定そのものはβエラーを考慮しているものではありませんが、信頼区間が広いような場合にはだいたいイベント発生数（アウトカム）が少ないので確かに一定の参考にはなります。この場合「ギリギリで１をまたぐかどうか」はβエラーとはあまり関係なく、どちらかというと「信頼区間が広いかどうか」が重要ですね。信頼区間の広さは結局のところ発症数に影響されますので、区間推定を介さずとも発症数が多いかどうかで信頼できるか（βエラーを起こしていないか）を考えることが多いと思います。

また、α・βエラーを考える時にはアウトカムの数だけでなく曝露の例数も考えます。たとえばタバコと死亡率の関連をみたいと考えたとき、日本では喫煙者が比較的多いので関連を検出しやすく、また棄却もしやすくなります。一方で、大麻と死亡率の関連を見たい場合、日本では大麻使用者が極めて少ない（し申告もしない）ので無作為に10万人集めても使用有りと答えるひとが一桁もいないかもしれません。たとえば3人いたとして、追跡期間中の死亡者が1人だった場合と2人だった場合で、推定されるリスク比は大きく動きますが、偶然や他の疾患の影響も大きいのでリスク比の数字自体あまり信頼できませんね。なので、こういう研究は曝露人口がそれなりに多い国でないとそもそも研究が難しいといえます。 

Q. 論文でtable1に平均値と標準偏差を提示することが多いと思うのですが、このときの標準偏差にも標本分散ではなく不偏分散を使うのはなぜでしょうか？ 

>

あまり意識して使い分けてなかったですが、SASやSPSS含めてだいたいの統計ソフトは記述統計の計算で不偏分散をデフォルトにしていますね。

　質問の回答として、明確なものではない私の推測ですが、おそらく（１）実用上は標本分散と不偏分散でほとんど値は変わらないのと、（２）標本分散を使うことの意義があまり大きくないというのがあるかもしれません。

　まず、標準偏差を標本分散と不偏分散それぞれから計算した場合、その違いはn=10で5%程度、n=100あれば0.5%とほとんど無視できるレベルです。

　そのうえで、標本分散自体は、国勢調査のような全数調査（調査対象集団＝母集団）では数値の意味づけができますが、そうでない場合、つまり標本抽出をした場合には「母分散の偏った(biased)推定値」なので使い所が難しいと思います。なので、統計学者は不偏分散を指して"分散”と呼ぶのがほとんどのようです。統計の教科書的には標本分散と不偏分散を分けていますが、実用上出てくるのは後者がほとんどだと思います。

　論文のtable 1は、確かにその標本の分布を表すために示しますが、t検定をつかって群間比較をしたりもしますよね。その時に比較しているのは、あくまでも標本から推測される母集団同士の平均値です。（なので、t検定は２標本の不偏分散を使いますよね）　なので、標本の分布を示すという行為にも背景の母集団を常に念頭に置いているとも考えられます。このような検定との一貫性も含めて、不偏分散やそのルートをとった標準偏差が使われているのではないかなと思います。

　"標本分散"という日本語がややこしさの一因でもありますね。標本分散は"標本で使う分散"ではありません。あくまでも分散という指標　1/n * Σ(xi-m)^2　があって、それを母集団で調べられたら母分散と呼びますが、"標本に同じ式を無理くり当てはめたもの"が標本分散です。 むしろ不偏分散のほうが"標本(で使う)分散"と言い換えたほうが良いかも。

http://www.medicine.mcgill.ca/epidemiology/hanley/IntMedResidents/SD-SE.pdf

Fig 1で整理していますが、σ: SD of Population (母標準偏差)と SD: Estimate of Population SD from sample (標本から推測した母標準偏差)と記載されてます。前者が、教科書でいう標本標準偏差と同じ、後者が不偏標準偏差と同じ意味です。このような表記のほうが混乱が少ないような気がします。 

■　欠測値

MCAR出ない限りは基本的にバイアスを生じますので、本来は少なくても補完ができたほうが望ましいです。研究室の場合は、多重代入などが一般的になる前から、多変量解析ではリストワイズ除去をしているので古い論文と方法を合わせるという意味合いと、欠測例が0.5%も無いのでほとんど無視できるだろう、という判断で除いてるわけです。

欠測値については「どのように処理したのか」を明記するようSTROBE statementでも言及されているので、自分が査読するときには、まずこの記述があるかを確認しています。変数にもよりますが、一般的に使われる検査値などの共変量については、多くのコホート研究では5%程度、リアルワールドデータのように診療録など事前に標準化されていないデータを使った場合は10%位の欠測値を報告しているものが多い印象です。あくまで個人的な感覚ですが、重要そうな変数が欠測頻度10%程度あれば多重補完してみたら？と提案することがあります。

また、これも経験的に言われていることですが、一つの変数で欠測がおおむね20%を超えるような場合は、処理していても慎重に読んだほうがよいですね。要は、多重代入法を使ったとしても誤差(95%CI)が結構広がってしまうので、有意か有意でないかの判断に影響する可能性があがります。どちらにしてもリストワイズより多重代入のほうがベターですが、そもそも極端に欠測が多いデータは、本来はデザインの時点で極力減らすための努力が必要ですね。

（https://www.slideshare.net/koichirogibo/ss-32418161　では25%以上は怪しいと言われてます） 

Q. 

①マッチマージとは？

②proc regとproc glmの使い分けは？

③残差分析の見方は？

>

1)について、マージは実際 IDなどで使う場面が多いですね。

[DATA SBP_SAMPLE]

ID | SBP

1 | 121

2 | 140

3 | 155

[DATA BPM_SAMPLE]

ID | BPM

1 | 85

3 | 61

という２つのデータがあるとき

DATA ALL; MERGE BPM_SAMPLE;

By ID;

とマージすると下のデータセットができます。

[DATA ALL]

ID | SBP | BPM

1 | 121 | 85

2 | 140 | .

3 | 155 | 61

一方、今回テキストで使われたマッチマージは、全員に同じ数字をマージさせるために使うテクニックです。

たとえば、SBP 120以上、140以上、160以上というカットオフ値をデータに付与するために、データセットを準備してみます。

DATA CUTOFF;

    input mv1 CO1 CO2 CO3;

    cards;

    1 120 140 160

    ;

で下のデータセットができます。

[DATA CUTOFF]

mv| CO1 | CO2 | CO3

1 | 120 | 140 | 160

これとさきほどの DATA ALL をくっつける場合、

DATA ALL; SET ALL;

    mv=1;

DATA ALL_MERGED; MERGE ALL CUTOFF;

    by mv;

とすると、データセットの横にカットオフとなる基準値を貼り付けることができます。

[DATA ALL_MERGED]

ID | SBP | BPM | mv | CO1 | CO2 | CO3

1 | 121 | 85 | 1 | 120 | 140 | 160

2 | 140 | . | 1 | 120 | 140 | 160

3 | 155 | 61 | 1 | 120 | 140 | 160

もし両方のデータセットにmvが存在しない状態で（＝byステートメントを使わないで）マージさせると・・・

DATA ALL_MERGED; MERGE ALL CUTOFF;

↓のように、一行目にしかマージされません。

[DATA ALL_MERGED]

ID | SBP | BPM | CO1 | CO2 | CO3

1 | 121 | 85 | 120 | 140 | 160

2 | 140 | . | . | . | . 

3 | 155 | 61 | . | . | . 

なので、mvは、全ての行に繰り返し貼りつけるために使われます。

ただし、あるカットオフ値を超えるかどうかは if thenで   if SBP >=140 then HT=1;のように書くことが多いですし、可読性が高いです。また教科書のように３分位や４分位で分ける場合も、proc rankを使うことが多いので、現在は↑のように定数を使ったマージをする機会はあまりないかもです。

2) についてですが、 GLM は一般線形モデル(general linear model)のためのプロシジャなので、アウトカムの確率分布に正規分布を仮定した検定・線形モデルの汎用プロシジャです。つまり、t検定・分散分析 & 共分散分析・単回帰分析&重回帰分析をカバーしています。

　t検定にはPROC TTEST、分散分析には PROC ANOVA, 回帰分析には PROC REGというそれぞれ専用のプロシジャがありますが、どれもPROC GLMで再現可能です。たとえば、PROC GLMで、右辺に連続量を入れてSOLUTIONをオプションをつけると推定値が出てきますが、これはPROC REGで出てくる回帰係数と同じものです。

　使い分けとしては、使いたいオプションによって書き分けることができます。PROC REGでは、残差分析のオプションが豊富ですし ODS Graphics を指定するとデフォルトで残差分析のグラフをたくさん出してくれます。また、多重共線性を調べるためのVIF (分散拡大係数)を計算するオプションが使えるといった利点があります。私は VIFチェックのために良く使っています。

　一方で、PROC REGは（というか一般的に重回帰モデルは）両辺が連続量のモデルを前提につくられていたので、classステートメントを使うことができません。なので、カテゴリ変数を入れる場合には事前にダミー変数を入れる必要があります。

たとえば diab_c (nonDM=0, prediabetes=1, DM=2)という変数があったときにPROC GLMなら↓のようにすれば非糖尿病者に対する前糖尿病・糖尿病者の回帰係数を出すことができます。

PROC GLM;

    CLASS diab_c (ref=first);

    MODEL Y = diab_c/solution ss3;

一方、PROC REGで同じ解析を再現するなら下のようになります。

DATA SAMPLE; SET SAMPLE;

if diab_c=0 then do; diab_d1=0; diab_d2=0; end;

if diab_c=1 then do; diab_d1=1; diab_d2=0; end;

if diab_c=2 then do; diab_d1=0; diab_d2=1; end;

PROC REG;

    MODEL Y = diab_d1 diab_d2;

ちなみに、アウトカムの(誤差)分布を正規分布以外に拡張したモデル(= 一般化線形モデル General"ized" linear model)は、SASではPROC GENMODが対応しています。t検定・分散分析 & 共分散分析・単回帰分析&重回帰分析　に加えて、ロジスティック回帰・ポワソン回帰 などを一つのプロシジャで書くことができます。

3-4) ですが、最小二乗法の性質上、予測値と残差の間に相関が無いように係数が決まります。

なので　proc reg; model PRED=RESID; は必ず無関連 (β=0, p=1.00)になります。

むしろ予測値を横軸・残差を縦軸にとったプロットを視覚的にみて、↓のようになんらかの傾向がないことを確認します。

https://corvus-window.com/all_residual-plot/

データ数が極端に多いと視覚的な判断も難しくなりますが、データが十分に多ければ外れ値の影響も受けにくくなり、正規性から逸脱していても推定値はロバストになるので、残差分析じたいそれほど重視しなくても良いと思います。

最後のβエラーについて補足すると、検定や推定そのものはβエラーを考慮しているものではありませんが、信頼区間が広いような場合にはだいたいイベント発生数（アウトカム）が少ないので確かに一定の参考にはなります。この場合「ギリギリで１をまたぐかどうか」はβエラーとはあまり関係なく、どちらかというと「信頼区間が広いかどうか」が重要ですね。信頼区間の広さは結局のところ発症数に影響されますので、区間推定を介さずとも発症数が多いかどうかで信頼できるか（βエラーを起こしていないか）を考えることが多いと思います。

また、α・βエラーを考える時にはアウトカムの数だけでなく曝露の例数も考えます。たとえばタバコと死亡率の関連をみたいと考えたとき、日本では喫煙者が比較的多いので関連を検出しやすく、また棄却もしやすくなります。一方で、大麻と死亡率の関連を見たい場合、日本では大麻使用者が極めて少ない（し申告もしない）ので無作為に10万人集めても使用有りと答えるひとが一桁もいないかもしれません。たとえば3人いたとして、追跡期間中の死亡者が1人だった場合と2人だった場合で、推定されるリスク比は大きく動きますが、偶然や他の疾患の影響も大きいのでリスク比の数字自体あまり信頼できませんね。なので、こういう研究は曝露人口がそれなりに多い国でないとそもそも研究が難しいといえます。 

Q.

①Cochran-Armitage検定は一度回帰直線を描いてから，その理論値と観測値を使ってχ二乗分布に基づいた分散分析を行い，直線性を検定するという理解でよろしいでしょうか．

②また，上記のように最後のステップがχ二乗分布による分散分析なので，P値でしか表現しえないと考えてもいいでしょうか？ 

>

Cochran Armitage testの解説も兼ねて回答です。

　P値のみ表現する検定としてよく使われるのが傾向性検定で、なかでも割合（２値変数）の傾向性を検定する Cochran-Armitage検定、連続量の傾向性を検定するJonckheer-Terpstra検定は比較的使用頻度が高いです。特にゲノムの症例対照研究でアレル保有数(0,1,2)ごとの疾病割合を示したりするときに慣習的に使われます。

Cochran-Armitage 検定は、割合 (p) にたいして回帰直線を引いて、その傾きが有意かどうかを検定します。

　p: 疾病頻度

　s: 説明変数のクラス (たとえば 4分位で分けた群に0, 1, 2, 3を割り当てたもの)　

のとき、 p = α + βs 　となる直線を最小二乗法で引き、

 　帰無仮説 H0 : p(s=0) = p(s=1)  = p(s=2)  = p(s=3) 　→　β＝0

　対立仮説H1:  p(s=0) ≦ p(s=1)  ≦ p(s=2)  ≦ p(s=3) 　（もしくはp(s=0) ≧ p(s=1)  ≧ p(s=2)  ≧ p(s=3)）　→ β≠0

を検定しています。単調増加ないしは単調減少傾向の検定といわれるものですね。

検定の仕方としては、

（１）モデルのばらつきを平方和に分解していき、検定統計量z^2についてχ二乗検定するもの ：z^2 ~ Chisquare(df=1)  ↓と

　https://cochineal19.hatenablog.com/entry/2020/12/13/003731

　　（＝先生が書かれている分散分析的なアプローチ）

（２）代数的にスコア統計量(z)を導いて、正規分布を仮定した検定を行うもの： z ~ Normal

があります。どちらも式変形すると等価になりますが、SASの計算上は（２）が使われているので、PROC FREQの trendオプションで検定したときにはχ二乗値ではなくZ値が表示されます。

さらに、この検定結果はロジスティック回帰での「包括帰無仮説 : BETA=0 の検定」で出てくるスコア検定のものとも一致します。

(各回帰係数の推定値の検定にはスコア検定ではなくWald検定が使われているので、この値とは少しだけズレます)

基本的に回帰直線の傾きβを検定しているので、βの大きさを解釈することもできないことはないです。

ただし、この場合には「X(群)の番号が1つ上がるときにyの割合がどれほど動くか」を見ていることになりますが、正直なところこの値は解釈する意味があまりありません。

たとえば血圧値を４分位したときにQ1群→Q2群、Q2群→Q3群、Q3群→Q4群の間で、疾病頻度が平均的にどれほど増加（減少）するかが推定値βの解釈ですが、Q1やQ2という区分けが集団依存なので重要な情報ではないのです。

ということで、先生の質問ですが

①　あってます。計算上はさらに式変形した値が使われています。

②　方法論的な問題よりも、むしろ傾向性検定の目的がそもそも推定値を求めることに意味を置いていない（単調増加/単調減少の有無にのみ興味がある）から、と考えています。

なお、↑で書いたようにCA検定とロジスティック回帰の傾きβのWald検定は、ほんのすこしだけ値がズレるのですが、見ているものはほぼほぼ同じなので久山町研究や他の疫学コホートでは傾向性検定としてロジスティック回帰を使っていることも多いです。この場合も、群の番号１上昇ごとのβやオッズ比を出すことは出来ますが、解釈する意味が薄いのでp値のみを報告するのが一般的です。 

勉強会で一般化線形モデルまで話題が出るようになって、いいですね。

さて、質問についてですが、大きく次の２つですね。

１）proc logistic で調整済みの割合を出せるか。

２）一般化線形モデルで二項分布＋ロジットリンク関数を指定した場合と、ロジスティック回帰モデルの違いはなにか。

１）proc logisticからロジスティックモデルのパラメータ推定値が出ているので、そこからエクセルで計算すれば割合の予測値を出せます。

ln(p/{1-p})=Σβ  から変形して

 p = 1/(1+e^-{Σβ})　に当てはめればいいですね。

添付ファイルでは、丹後本の低体重データを使って、

 Y=低体重(lwt) X=喫煙(smoke) の単変量ロジスティック回帰を回してみました。

proc logistic descending;

    model low = smoke;

切片=-1.087, β(smoke)=0.704 なので

smoke =0 の場合 p=1/(1+exp^-(-1.087+0*0.704)) =0.2521,

smoke =1 の場合 p=1/(1+exp^-(-1.087+1*0.704)) =0.4054となります。

調整していないので、

proc freq; 

     tables low*smoke;

で出したクロス表での割合とほぼ一致しますね。調整する場合も、共変量の平均値を proc meansで出しておいて、共変量はβ×平均値を代入すれば算出ができます。

※ちなみに、proc logistic では予測平均値を直接出せないと思っていたのですが、↑で先生が調べられたようにlsmeansが実装されていたので、予測平均値を出すこともできました。ただ、このときにクラス変数のパラメータの置き方が REF型(param=ref)ではなくGLM型(param=glm)になったりして、少し書きなれる必要があります。

２）両者は、等しい(equivalent)モデルです。

添付には↑と同じデータをgenmodでも回した結果が出ています。

proc genmod descending;

class smoke;

model low  = smoke/dist=bin link=logit;

lsmeans smoke/ilink;

定式化が違うため、切片などが変わっていますが、結果的に得られる推定値・予測値は全く同じ値になります。

↑のように、logistic を使っても出せることは出せますが、parameterizationがGLM型に変わるので、同じモデルでもオッズ比を出すときと予測割合を出すときで推定値が変わります。

一方、genmod を使うとロジスティック回帰も線形回帰も一緒に扱えるので、プログラムを後から書き直すときに若干楽だったりします。（→proc glmとproc logisticの文法の違いとかを考えないで済む）

添付には 修正ポワソンを同じデータに対してかけた結果も載せています。割合の予測値は 0.2522, 0.4054とほぼ同じ結果が得られていますが、仮定している分布が二項分布ではなくポワソン分布なので、途中の推定値は全然違いますね。

◆ Retrospective studyの意味

　"retrospective study"の語義・定義は、大事な部分ですね。基本的な考え方は手島先生が書いてくださったとおりです。

この用語の使い方は、疫学の界隈でも一貫したものではなく、混乱を招いているものです。

　症例対照研究は、アウトカムを生じた人を集めて、それに対応する対照群を集めて、曝露の情報を調べることになります。対象者をサンプリングするための基準がアウトカム発生有無を起点とするので、アウトカムから曝露に遡る＝時間的に後ろから前に戻っていく、ということでretrospectiveと言われます。これはイメージがしやすいですね。

　回顧的コホートは、あくまでも曝露の情報に基づいて対象者を集めるので、時間的な順序は曝露→アウトカム、と前向きコホートと同様です。ただし、ここで問題になるのは、「曝露とアウトカムの情報が必ずしも独立していない」という点です。実臨床では、たとえばある疾病が疑われるときには、診断を確定させるために画像検査をするけど、そうでなければ画像検査はしない、といった取捨選択がされています。そうすると、疾病（アウトカム）の確率が高い人では検査（曝露）の情報がとれていて、反対に確率が低い人では曝露の情報に欠ける、といった情報の偏りが生じます。アウトカムによって曝露が左右されてしまうので、そのまま解析すると歪んだ(biased)結果が得られてしまいます。あるいは、頻回に詳細な検査をできた人では確定診断が付けやすく、検査をできなかった人では診断がつけにくい、といったように曝露がアウトカムの測定精度に影響する場合も往々にしてあります。

　回顧的コホート研究を”後ろ向き=retrospective”と称しているのは、既にアウトカムが記録済みの過去のデータベースを用いるので、それを「解析するひと」にとって過去の情報という意味ですが、正直なところ、これをretrospective と呼ぶのはナンセンスな面があります。Porta et al『Dictionary of Epidemiology』といった世界的に広く使われている書籍がこの定義を紹介しているのが問題ですが。。。たとえば、久山町研究の1960年代のデータは、解析する我々からすると過去のデータベースで、曝露やアウトカムの測定に関わった人は誰もいませんが、まぎれもなく前向きコホートですね。Rothman の『Modern Epidemiology』第４版では、Prospective vs Retrospective と題して、この定義のブレを解説しています。（添付）

　回顧的コホートと比べると、前向きコホート研究の場合には、対象者を設定して、曝露を測定し、アウトカムを測定するという順序で調べていきますし、アウトカムも事前に定義した方法で測定することが求められるので、曝露因子の精度にアウトカムが影響したり、アウトカムの精度に曝露因子が影響したりすることは少ないと言えます。ただし前向きコホート研究であっても、アウトカムの診断をつけるときに曝露の情報を加味してしまったりすると、結果にバイアスがかかります。たとえば、脳卒中での死亡例の病型を決めたいとき「この人は血圧は低かったけどヘビースモーカーだったから脳梗塞だ」、「昔から血圧が高かったから脳出血としよう」としてしまうと、喫煙と脳梗塞・高血圧と脳出血の関連を解析したときに、見かけ上の関連が強くなります。一方で、回顧的コホートでも死亡のようなハードアウトカムであれば、曝露因子によって変わることは稀です。このように、単にデータの取得方法が前向きか、後ろ向きかという区分は、バイアスの懸念が大きいかどうかと必ずしも対応できていません。それよりも「曝露とアウトカムが、互いの測定（有無や精度）に影響をしうるか」を考えるほうが、デザインの質を適切に評価できるとModern Epidemiologyでは述べられています。

　実際、論文でも retrospective studyという言葉は様々に用いられているのが現状ですので、この言葉にあまり惑わされず、曝露・アウトカムの測定方法を読んでバイアスの程度を想像するのが重要です。

◆コホート内症例対照研究の使いどころ

　コホート内で稀な疾病を研究する場合にも、可能であればコホート全体で解析した方が望ましいです。

たとえば地域住民1万人のコホートで20人の症例が発症したとしましょう。症例20人・対照20人に絞るよりは、症例20人・対照9980人で解析した方が、対照群の人数が多いのでオッズ比などを計算したときに誤差（95%CI）がより小さく推定できるはずです。また、症例20人・対照20人に絞った場合、症例群にマッチングをさせるので、解析対象者40人の背景特性は症例20人に似通ってしまい、地域住民からは偏った集団になり、一般化可能性が下がります。一方、対照群9980人を残しておけば、解析対象集団は元の地域住民1万人と一致するので、一般化可能性を保てます。

　では、コホート内症例対照研究が必要なタイミングは？というと曝露の測定リソース（主にお金）が足りない場合です。「9980人分の保存試料全部を測定するには予算が足りない」とか「全例測るのはもったいない」とか「全員にアンケートを取りに行く人件費が出せない」といった場合にコホートから部分集団を抽出して測定をするわけです（たとえばhttps://pubmed.ncbi.nlm.nih.gov/38497939/）。あとは、大規模データベース研究の場合に、十分に症例・対照の人数がいるので（誤差が十分小さいので）、解析時のコンピュータの計算負荷を抑えるためにあえて人数を絞るケースもあります（https://www.bmj.com/content/381/bmj-2022-072770.long）。

　なお、コホート内で症例・対照をマッチさせると、上記のように元のコホート集団から偏った集団（発症集団の背景特性）になってしまうのが問題でした。それを乗り越えるため、ケースコホート研究といった手法が提案されて、広まりつつあります。ケースコホートでは、対照群を抽出するときに、症例とマッチングをするのではなく、元のコホート集団からランダムに抽出をします。そうすることで、対照群が元のコホート集団の特性を反映するので、オッズ比等を計算するときにバイアスがかかりにくく、一般化可能性を高めることができます。