HADを使ってみた：8.謎の機能「曲線当てはめ」
謎の機能「曲線当てはめ」
HADはしばらく前からちょっとずつ使ってはいるんですが，この「曲線当てはめ」は，これまで一度もさわったことがありません。公開されているマニュアルにも出てこない（はずだ）。あらためて清水先生のブログを「曲線当てはめ」で検索しても何もヒットしない（2021／2／1現在）。しかし，記事をずっと読んでいくとそれらしいのが１つだけあって，
説明変数が１つの場合の一般化線形モデルにおいて，予測直線（曲線）を表示するようになった
これが，Version 14.50 へのアップデートで言及されています。記事を見てみると，ロジスティック曲線をあてはめたプロットが。なるほど，要するに1要因のロジスティック回帰ね。それをプロットに重ねてみて，結構当てはまりそう，とか，ぜんぜんダメかな，とか見たいときに使うらしい，と考えました。
とはいえこの機能をテストするのに都合の良いデータはない。さて困った，と思ったところで，1年半くらい前に買った久保拓弥先生の『データ解析のための統計モデリング入門』（通称「みどり本」）を思い出しました。たしか，ロジスティック回帰も説明されていたはずです。ということで，そのころ苦労してゲットしたサンプルデータを探しました。
現在，このデータは下記のページからダウンロード可能です。大量のデータがダウンロードされますが，今回使用したのは，第6章のサンプルデータdata4a.csv （ロジスティック回帰）です。
これは架空のデータなのですが，ある植物から8個ずつ種子をとって（N）そのうち何個が発芽するか（y）を調査した。植物の大きさ（x）や施肥条件（f:Control/Treat）で発芽率が変わるだろうか，を調べる実験，という想定のようです。もとのデータにはID変数がないので付け加えます。こんな感じになります。
従属変数の条件
施肥条件の変数もあるのですが，説明変数は１つということなので，植物の大きさ（x）で発芽数（y）を説明できるかどうか，「曲線当てはめ」を試してみますと，こんなメッセージが。
従属変数は確率である必要があるとのこと。ふむふむ。ならば，発芽数（x）を種子数（N）で割り算すればいい。数値の変換は次の手順。
計算したい変数を使用変数に設定
⇒【変数の作成】
⇒【数値変換】
⇒「変数の割り算」にチェックして，・・・？
と書きたかったのだが，割り算だけは実装されていない。残念。でも，足し算と引き算と掛け算はできるので，尺度の逆転項目の処理などはこれを使うとよい。割り算はできないようなので（できてほしい…），Excelに計算式を入力して，ratio という変数を作成した。【データ読み込み】をクリックする。
ここで，変数を追加したときには必ず【データ読み込み】ボタンをクリックする必要がある，というのが，HADを使うときのお作法。単に数値の誤りを修正したりするだけなら【データ読み込み】を押す必要はないのだが（分析にちゃんと反映される），変数を追加したときには【データ読み込み】クリックは必須。
曲線当てはめ
さて，確率を計算したので，ratio を従属変数，xを独立変数として，曲線当てはめを実行してみましょう。
使用変数に ratio, x をこの順で設定
⇒【分析】
⇒曲線当てはめ
「この順で設定」の部分はいつも通りポイントで，必ず「従属変数 独立変数」の順番。これもHADのお作法ですね。結果はこんな感じ。出力が縦に長いので，プロットの下に表示されていた数値を右側に移動させてある。なるほど，ロジスティック回帰だね。と言いたいのだが，待て待て。ちょっと変だ。
▲変なこと１：
「累積正規曲線のパラメータ」って何だ？ 数値から考えると，独立変数 x の統計量に近いが， x の統計量は M = 9.967，SD = 1.089 である。だからこれは，こんな当たり前の統計量ではない。私の知識の外だなあ。
▲変なこと２：
適合度R2乗値，-2.131。R2乗というのは文字通り「2乗」値なので，マイナスになるはずのない値である。どうしてマイナスなのだ？
いや，これは困った。まったくわからない。しかも，この分析はHAD2Rに対応していない。ヒントがない。
とりあえずロジスティック回帰を
わからないなりに何かしよう，ということで，久保「みどり本」を参照しつつ，さっきのデータを使ってロジスティック回帰分析をしてみましょう。「みどり本」では施肥条件も説明変数にしていますが，HADの分析に合わせて植物の大きさ（x）だけで分析します。ただし，Rの glm() 関数では，従属変数に確率を直接指定することはできないので，これも「みどり本」の説明の通り，y（つまり発芽"成功"） と N-y（発芽"失敗"）をマトリクスにして渡します。コードは，
summary( glm( cbind(y, N-y) ~ x, data=dat, family=binomial( ) ))
出力はこんな感じ。
でも，これではHADの出力と比較できないので，確率をY軸に設定してプロットし，出力された切片と係数を使ってロジスティック曲線を書き加えてみました。コードはこんな感じです。
(res <- summary(glm(cbind(y, N-y) ~ x, data=dat, family=binomial())))
plot(dat$x, dat$y/dat$N, pch=1, xlab="x", ylab="ratio", xlim=c(7,13))
lg <- function(z) {1 / (1+exp(-z))}
curve(lg(x*res$coefficients[2]+res$coefficients[1]), col="red", add=T)
では，両者のプロットを，軸の設定などを調整して，並べてみましょう。
おお，ぴったりです。（左がHAD，右がR）
ということで，HADの「曲線当てはめ」は，1要因のロジスティック回帰分析をして，その係数を使ってロジスティック曲線をプロットに重ねている，ということですね，きっと。
とはいえ，出力される統計量についてのヒントはまだ得られておりません。中途半端なことこの上ないのですが，ずっとここで止まっているのは嫌なので，保留にして，先に進むことにします。
そういえば
もう一つ謎だったことは，この分析メニューが，【分析】フォームの「データの要約」グループに入っていることです。ロジスティック回帰がここ？ という疑問が当然あるのですが，考えてみれば，散布図の分析のなかに直線回帰分析がすでに含まれていました。
データをプロットしたときに，これ，回帰直線引けそうだね，という感じでピッとひいてみる。同じような感覚で，別の（直線回帰なんか想定できそうにない）データをプロットしたときに，これ，ロジスティックでいけるんじゃない？ という感じでぐいっと引いてみた。それが「曲線当てはめ」なのかなあと。つまり散布図のロジスティックバージョン（？）みたいな扱い。なんていうふうに考えると，それなりに理解できそうな気がします。
ということで，今回はこれまで。