ロジスティック回帰２

前回は、ロジスティック回帰モデルの概念を説明した上で、カブトガニの甲羅の幅を説明変数とする、単項ロジスティック回帰モデルを当てはめました。で、当てはめた上で、予測値を計算して、グラフにプロットしました。

とりあえず、前回で、データにモデルを当てはめることは出来ました。
でも、データを使って計算できる事と、計算に用いたモデルが、データに当てはまってる事とは、全く別問題です。
言い換えれば、数字とデータ形式さえ合えば、計算することなんて簡単です。なんせ、計算するのは、人でなくてコンピューターですもん。
と言うわけで、当てはめた統計モデルが、データと整合性があるかを、今回はチェックしてみます

長々と書きましたが、簡単に言えば、回帰分析でいう回帰診断とか、残差分析と同じことです。

まず、前回の計算結果は、こんな感じでした
ただし、今回は、corパラメータを指定して、coefficient内の相関係数も出力しています

===============================================================
>summary(crab.result, cor=T)


Call:
glm(formula = crab.data$satell.flag ~ crab.data$width, family = binomial)

Deviance Residuals:
Min 1Q Median 3Q Max
-2.0281 -1.0458 0.5480 0.9066 1.6941

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.3508 2.6287 -4.698 2.62e-06 ***
crab.data$width 0.4972 0.1017 4.887 1.02e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 225.76 on 172 degrees of freedom
Residual deviance: 194.45 on 171 degrees of freedom
AIC: 198.45

Number of Fisher Scoring iterations: 4

Correlation of Coefficients:
(Intercept)
crab.data$width -1.00

===============================================================

Deviance Residuals:
Null deviance

が逸脱度を表します
また、

crab.result%reiduals

ここで、crab.resultは、関数glm()の結果が入っているオブジェクト（変数です）
コマンドを入力すれば、個別データの残差についてデータを得ることが出来ます。が、今回は、まだ不勉強のため、パスです

と言うわけで、今回計算するのは、予測値の９５％信頼区間です
つまり、前回計算した予測値が、９５％の確率で、どの範囲でばらつくかを、グラフ化しようと言うわけです。予測値の区間推定を考えるのは、Coefficientsたちが、推計値だったことを考えれば、当然ですよね。

まず、大標本下での線型予測子の標準誤差（ASE)は、下式になります

Var(hat(α)+hat(β))=Var(hat(α)+x^2Var(β)+2x*Cov(hat(α),hat(β))

*1 hat()は推計値　
*2 x 独立変数 を表します
*3 Cov(α,β)= r * sqrt(Var(α))*sqrt(Var(β))
ただし、rは相関係数

早速、線型予測子、つまり、オッズの95％信頼区間を計算しましょう。
まず、coefficientの分散を変数に叩き込んでおきます

> var.a <- 2.6287*2.6287
> var.b <- 0.1017*0.1017

続いて、Coefficient間の相関係数も変数に入れておきます

> r.coef <- -1

準備が整いましたので、先ほど示した線型予測子の標準誤差(ASE)を計算しましょう

> ASE.linear.predictor <- var.a+(crab.data$width^2)*var.b+(c*crab.data$width*sqrt(var.a*var.b))

ごちゃごちゃしてますが、掛け算と足し算しかありません
小学生でも、計算できます
さて、つづきまして、線型予測子の標準誤差を用いて、線型予測子の95％信頼区間を計算しましょう

ASE(Asymptotic Standard Error)は、大標本下で正規分布に従うので、95％信頼区間の信頼係数は、1.96です。そんなわけで、上限と下限の信頼区間は

> upper95.logit.ci <- crab.result$linear.predictor+(1.96*sqrt(ASE.linear.predictor))

> lower95.logit.ci <- crab.result$linear.predictor-(1.96*sqrt(ASE.linear.predictor))

で計算できます
ちなみに、crab.result$linear.predictorと指定するれば、ワザワザ計算せずとも、Ｒが、独立変数に対する線型予測子を返してくれます。

さて、線型予測子の信頼区間が求まりました
でも、線型予測子は、オッズの信頼区間です
オッズに対して、線型モデルを当てはめたのが、ロジスティック回帰でしたから。
でも、オッズのままでは、モデルの解釈が難しいですよね
だから、ここから、前使った式で、オッズから、確率に計算しなおします。計算式は、確率を独立変数の関数として考えた、こんな感じの式でしたね

π(x)=exp(α+βx)/(1+exp(α+βx))

さて、α+βxは、いままでに計算済みです
そこで、上の式に、さっき計算した２つの変数

upper95.logit.ci
lower95.logit.ci

を叩き込んで、確率の信頼区間を計算しましょう
それが、こんな感じです

> upper95.ci <- exp(upper95.logit.ci)/(1+exp(upper95.logit.ci))
> lower95.ci <- exp(lower95.logit.ci)/(1+exp(lower95.logit.ci))

言うまでもなく、この二つが、上限と下限の信頼区間のベクトルです


さて、つまんねえ計算が終わったので、上の結果を、予測確率と共に、グラフにしてみましょう。もちろん、x軸は、カブトガニの甲羅の幅ですよ

最初に、下限の信頼区間

> plot(crab.data$width[width.order],lower95.ci[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",col="blue",type="l")

上書きを許可
> par(new =T)

次に、上限の信頼区間

> plot(crab.data$width[width.order],upper95.ci[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",col="red",type="l")

上書きを許可

> par(new =T)

つづいて、予測確率をプロット

> plot(crab.data$width[width.order],crab.result$fitted.values
[width.order],xlim=c(21,34),ylim=c(0,1),xlab="Crab.Width",ylab="Probability",type="l",
+ main ="Ex. of Logistic Regression/ 95% of C.I.")

上書きを許可しておいて
> par(new =T)

最後に、カテゴリー化したデータから計算したサンプル確率をプロット
>plot(crab.category$width.cate,crab.category$SampleRatio,xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",col="red")

で、出てきたのが、上のグラフです

赤：上限の９５％信頼区間
青：下限の９５％信頼区間
黒：予測確率
赤の点：カテゴリー化したデータから計算したサンプル確率

なんだか、えらい手間をかけたワリには、簡単なグラフですね
すでに、計算してあればいいのに、なんて思ってしまいますね

次回は、Devianceとかresidualsあたりを説明する予定です


=======================================================
なんだか、えらく長くなった

（参考文献)
カテゴリカルデータ解析入門
Alan Agresti
サイエンティスト社

のp.151あたりを計算

(感想)
データハンドリングが大変
アクセスとかエクセルなら、簡単なのになんて思ったりします
あと、Coefficient間の相関係数が、どうすれば出せるのかが
中々分からなかった。
やったことは大したことないのに、コマンドを調べたりするのは
みょーーーーに、大変だったのでした

なんか、先に、重回帰分析とかを、試してみればよかったかも
そんな気がしている

いずれにしても、まだまだ勉強中の身ですね
=======================================================
前回は、ロジスティック回帰モデルの概念を説明した上で、カブトガニの甲羅の幅を説明変数とする、単項ロジスティック回帰モデルを当てはめました。で、当てはめた上で、予測値を計算して、グラフにプロットしました。

とりあえず、前回で、データにモデルを当てはめることは出来ました。
でも、データを使って計算できる事と、計算に用いたモデルが、データに当てはまってる事とは、全く別問題です。
言い換えれば、数字とデータ形式さえ合えば、計算することなんて簡単です。なんせ、計算するのは、人でなくてコンピューターですもん。
と言うわけで、当てはめた統計モデルが、データと整合性があるかを、今回はチェックしてみます

長々と書きましたが、簡単に言えば、回帰分析でいう回帰診断とか、残差分析と同じことです。

まず、前回の計算結果は、こんな感じでした
ただし、今回は、corパラメータを指定して、coefficient内の相関係数も出力しています

===============================================================
>summary(crab.result, cor=T)


Call:
glm(formula = crab.data$satell.flag ~ crab.data$width, family = binomial)

Deviance Residuals:
Min 1Q Median 3Q Max
-2.0281 -1.0458 0.5480 0.9066 1.6941

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.3508 2.6287 -4.698 2.62e-06 ***
crab.data$width 0.4972 0.1017 4.887 1.02e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 225.76 on 172 degrees of freedom
Residual deviance: 194.45 on 171 degrees of freedom
AIC: 198.45

Number of Fisher Scoring iterations: 4

Correlation of Coefficients:
(Intercept)
crab.data$width -1.00

===============================================================

Deviance Residuals:
Null deviance

が逸脱度を表します
また、

crab.result%reiduals

ここで、crab.resultは、関数glm()の結果が入っているオブジェクト（変数です）
コマンドを入力すれば、個別データの残差についてデータを得ることが出来ます。が、今回は、まだ不勉強のため、パスです

と言うわけで、今回計算するのは、予測値の９５％信頼区間です
つまり、前回計算した予測値が、９５％の確率で、どの範囲でばらつくかを、グラフ化しようと言うわけです。予測値の区間推定を考えるのは、Coefficientsたちが、推計値だったことを考えれば、当然ですよね。

まず、大標本下での線型予測子の標準誤差（ASE)は、下式になります

Var(hat(α)+hat(β))=Var(hat(α)+x^2Var(β)+2x*Cov(hat(α),hat(β))

*1 hat()は推計値　
*2 x 独立変数 を表します
*3 Cov(α,β)= r * sqrt(Var(α))*sqrt(Var(β))
ただし、rは相関係数

早速、線型予測子、つまり、オッズの95％信頼区間を計算しましょう。
まず、coefficientの分散を変数に叩き込んでおきます

> var.a <- 2.6287*2.6287
> var.b <- 0.1017*0.1017

続いて、Coefficient間の相関係数も変数に入れておきます

> r.coef <- -1

準備が整いましたので、先ほど示した線型予測子の標準誤差(ASE)を計算しましょう

> ASE.linear.predictor <- var.a+(crab.data$width^2)*var.b+(c*crab.data$width*sqrt(var.a*var.b))

ごちゃごちゃしてますが、掛け算と足し算しかありません
小学生でも、計算できます
さて、つづきまして、線型予測子の標準誤差を用いて、線型予測子の95％信頼区間を計算しましょう

ASE(Asymptotic Standard Error)は、大標本下で正規分布に従うので、95％信頼区間の信頼係数は、1.96です。そんなわけで、上限と下限の信頼区間は

> upper95.logit.ci <- crab.result$linear.predictor+(1.96*sqrt(ASE.linear.predictor))

> lower95.logit.ci <- crab.result$linear.predictor-(1.96*sqrt(ASE.linear.predictor))

で計算できます
ちなみに、crab.result$linear.predictorと指定するれば、ワザワザ計算せずとも、Ｒが、独立変数に対する線型予測子を返してくれます。

さて、線型予測子の信頼区間が求まりました
でも、線型予測子は、オッズの信頼区間です
オッズに対して、線型モデルを当てはめたのが、ロジスティック回帰でしたから。
でも、オッズのままでは、モデルの解釈が難しいですよね
だから、ここから、前使った式で、オッズから、確率に計算しなおします。計算式は、確率を独立変数の関数として考えた、こんな感じの式でしたね

π(x)=exp(α+βx)/(1+exp(α+βx))

さて、α+βxは、いままでに計算済みです
そこで、上の式に、さっき計算した２つの変数

upper95.logit.ci
lower95.logit.ci

を叩き込んで、確率の信頼区間を計算しましょう
それが、こんな感じです

> upper95.ci <- exp(upper95.logit.ci)/(1+exp(upper95.logit.ci))
> lower95.ci <- exp(lower95.logit.ci)/(1+exp(lower95.logit.ci))

言うまでもなく、この二つが、上限と下限の信頼区間のベクトルです


さて、つまんねえ計算が終わったので、上の結果を、予測確率と共に、グラフにしてみましょう。もちろん、x軸は、カブトガニの甲羅の幅ですよ

最初に、下限の信頼区間

> plot(crab.data$width[width.order],lower95.ci[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",col="blue",type="l")

上書きを許可
> par(new =T)

次に、上限の信頼区間

> plot(crab.data$width[width.order],upper95.ci[width.order],xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",col="red",type="l")

上書きを許可

> par(new =T)

つづいて、予測確率をプロット

> plot(crab.data$width[width.order],crab.result$fitted.values
[width.order],xlim=c(21,34),ylim=c(0,1),xlab="Crab.Width",ylab="Probability",type="l",
+ main ="Ex. of Logistic Regression/ 95% of C.I.")

上書きを許可しておいて
> par(new =T)

最後に、カテゴリー化したデータから計算したサンプル確率をプロット
>plot(crab.category$width.cate,crab.category$SampleRatio,xlim=c(21,34),ylim=c(0,1),xlab="",ylab="",col="red")

で、出てきたのが、上のグラフです

赤：上限の９５％信頼区間
青：下限の９５％信頼区間
黒：予測確率
赤の点：カテゴリー化したデータから計算したサンプル確率

なんだか、えらい手間をかけたワリには、簡単なグラフですね
すでに、計算してあればいいのに、なんて思ってしまいますね

次回は、Devianceとかresidualsあたりを説明する予定です


=======================================================
なんだか、えらく長くなった

（参考文献)
カテゴリカルデータ解析入門
Alan Agresti
サイエンティスト社

のp.151あたりを計算

(感想)
データハンドリングが大変
アクセスとかエクセルなら、簡単なのになんて思ったりします
あと、Coefficient間の相関係数が、どうすれば出せるのかが
中々分からなかった。
やったことは大したことないのに、コマンドを調べたりするのは
みょーーーーに、大変だったのでした

なんか、先に、重回帰分析とかを、試してみればよかったかも
そんな気がしている

いずれにしても、まだまだ勉強中の身ですね
=======================================================