0　全体の章構成

1　一般化線型モデル(Generalized Linear Model)

リンク関数(link function)と線型予測子(linear predictor) \[g(\hat{y})=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}\]

1-1　ロジスティック関数

\[g(\hat{y})=log\left(\frac{Prob(y)}{1-Prob(y)}\right)=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}\] \[\frac{Prob(\hat{y})}{1-Prob(\hat{y})}=e^{\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}}=e^{\beta_{0}}e^{\beta_{1}x_{1}}e^{\beta_{2}x_{2}}...e^{\beta_{j}x_{j}}...e^{\beta_{p}x_{p}}\] \[Prob(\hat{y})=\frac{1}{1+e^{-(\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p})}}\]

1-2　最尤法(Maximum Liklihood Estimation)

1-3　対数オッズ比(log odds ratio)と2項ロジットモデル(binomial logit model)

　スクリプトを，以前のversionよりも少し綺麗にしました。
　読みこんだ後のデータフレイム名は短く d01 としました。

d01 <- read.csv("practice.csv")

d01$male <- c(1, 0)[d01$sex] # 男性ダミーの作成
d01$single <- c(0, 0, 1, 0, 0)[d01$q2100] # 婚姻上の地位変数から未婚ダミー作成
full <- with(d01, complete.cases(single, male, age, edu, income)) # 欠損値除外

d02 <- d01[full, c("single", "male", "age", "edu", "income")]

# glm( )関数を用いて二項分布binomial，ロジットリンクlogitを指定
result <- glm(single ~ male + age + edu + income, data=d02, family = binomial(link = "logit"))
summary(result)


Call:
glm(formula = single ~ male + age + edu + income, family = binomial(link = "logit"), 
    data = d02)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3831  -0.6923  -0.5111  -0.2816   2.3966  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.1605510  1.4740481   1.466  0.14272    
male         1.7164462  0.3807364   4.508 6.54e-06 ***
age         -0.0623619  0.0190020  -3.282  0.00103 ** 
edu         -0.0569334  0.0730136  -0.780  0.43553    
income      -0.0023085  0.0007374  -3.130  0.00175 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 344.82  on 337  degrees of freedom
Residual deviance: 306.77  on 333  degrees of freedom
AIC: 316.77

Number of Fisher Scoring iterations: 5

1-4　適合度検定(Goodness of Fit Test)と疑似決定係数(Pseudo R-squared)

　2020年6月20日現在，BaylorEdPsychパッケイジはCRANからremoveされています。
　

Package ‘BaylorEdPsych’ was removed from the CRAN repository.

Formerly available versions can be obtained from the archive.

Archived on 2020-06-08 as check problems were not corrected in time.

Please use the canonical form https://CRAN.R-project.org/package=BaylorEdPsych to link to this page.

上記のarchiveから BaylorEdPsych_0.5.tar.gz をdownloadして手動でinstallして使用する事は出来ました。
　現在検索すると，DescTools パッケイジの PseudoR2( ) が使用可能と出ます。
　ただ，私が色々と試したところ，BaylorEdPsychとDescToolsとで結果が全く異なる場合がありました。定義に従って自分で計算した結果は，BaylorEdPsychとぴったり一致します。
　どちらが正しいのか不安は有りますが，追加パッケイジ無しで自分で計算するスクリプトを公開しておきます。以下の例ではAdj.McFaddenを除くと3つの方法で一致しています。

# install.packages("ResourceSelection", repos="https://cran.ism.ac.jp")
require(ResourceSelection)
# install.packages("BaylorEdPsych", repos="http://cran.ism.ac.jp")
require(BaylorEdPsych)

疑似決定係数の定義と計算

各種の疑似決定係数の定義式は次のサイトを見ると良い。
idre UCLA, ‘FAQ: What are pseudo R-squareds?’（2019年10月22日閲覧）

LogLiklihoodをLLと表記すると，
McFadden’s pseudo R-squared: \[R^2=1-\frac{LL_{model}}{LL_{null}}\]

Cox-Snell’s pseudo R-squared: \[R^2=1-\left(\frac{L_{null}}{L_{model}}\right)^{\frac{2}{n}}\]

Nagelkerke’s pseudo R-squared: \[R^2=\frac{1-\left(\frac{L_{null}}{L_{model}}\right)^{\frac{2}{n}}}{1-{L_{null}}^{\frac{2}{n}}}\]

追加パッケイジ無しで計算する

## 自分で計算Pseudo R2
# 汎用性を高める為に， glm( )の結果オブジェクトを glm_result にcopyする
glm_result <- result # この右辺を適宜書き換える。

Lm <- glm_result$deviance/(-2)
L0 <- glm_result$null.deviance/(-2)
# Lm; L0

McFadden <- 1-Lm[1]/L0[1]     # McFadden
adjMcFadden <- 1-(Lm[1]- glm_result$rank - 1)/L0[1] # adjusted McFadden

# Cox & Snell
n_model <- glm_result$rank + glm_result$df.residual; # n_model
CoxSnell <- 1-(exp((2/n_model)*(L0[1] - Lm[1])))  # Cox & Snell

# Nagelkerke
Nagelkerke <- (1-(exp((2/n_model)*(L0[1] - Lm[1]))))/(1-exp((2/n_model)*L0[1]))

DescTools::PseudoR2(result, "all")

       McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
     0.11034136      0.08134078      0.10646341      0.16648700      0.10117853 
VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
     0.20035571      0.13461027      0.22674936      0.12519257    316.77268914 
            BIC          logLik         logLik0              G2 
   335.88791862   -153.38634457   -172.41033495     38.04798077

BaylorEdPsych::PseudoR2(result)

        McFadden     Adj.McFadden        Cox.Snell       Nagelkerke 
      0.11034136       0.07554066       0.10646341       0.16648700 
McKelvey.Zavoina           Effron            Count        Adj.Count 
      0.22674936       0.13461027       0.82840237       0.17142857 
             AIC    Corrected.AIC 
    316.77268914     316.95341203

McFadden; adjMcFadden; CoxSnell; Nagelkerke

[1] 0.1103414

[1] 0.07554066

[1] 0.1064634

[1] 0.166487

DescToolsのAdjusted McFaddenがどの様な計算なのかを探索してみる。

1-(Lm[1]- glm_result$rank - 0)/L0[1] # adjusted McFadden

[1] 0.08134078

1-(Lm[1]- glm_result$rank - 1)/L0[1] # adjusted McFadden

[1] 0.07554066

1-(Lm[1]- glm_result$rank - 2)/L0[1] # adjusted McFadden

[1] 0.06974054

大標本でのDescToolsとBaylorEdPsych

分析ケース数が数千二なるとこの二つのPseudo R2が明確にずれる様に思われたので，人為的に大きなデータセットを作成して分析してみる。

d03 <- rbind(d01, d01, d01, d01, d01, d01, d01, d01, d01, d01,
             d01, d01, d01, d01, d01, d01, d01, d01, d01, d01,
             d01, d01, d01, d01, d01, d01, d01, d01, d01, d01)
dim(d01); dim(d03)

[1] 384  43

[1] 11520    43

full <- with(d03, complete.cases(single, male, age, edu, income)) # 欠損値除外

d04 <- d03[full, c("single", "male", "age", "edu", "income")]

# glm( )関数を用いて二項分布binomial，ロジットリンクlogitを指定
result2 <- glm(single ~ male + age + edu + income, data=d04, family = binomial(link = "logit"))
summary(result2)


Call:
glm(formula = single ~ male + age + edu + income, family = binomial(link = "logit"), 
    data = d04)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3831  -0.6923  -0.5111  -0.2789   2.3966  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.1605510  0.2691231   8.028 9.90e-16 ***
male         1.7164462  0.0695126  24.693  < 2e-16 ***
age         -0.0623619  0.0034693 -17.975  < 2e-16 ***
edu         -0.0569334  0.0133304  -4.271 1.95e-05 ***
income      -0.0023085  0.0001346 -17.146  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 10344.6  on 10139  degrees of freedom
Residual deviance:  9203.2  on 10135  degrees of freedom
AIC: 9213.2

Number of Fisher Scoring iterations: 5

## 自分で計算Pseudo R2
# 汎用性を高める為に， glm( )の結果オブジェクトを glm_result にcopyする
glm_result <- result2 # この右辺を適宜書き換える。

Lm <- glm_result$deviance/(-2)
L0 <- glm_result$null.deviance/(-2)
# Lm; L0

McFadden <- 1-Lm[1]/L0[1]     # McFadden
adjMcFadden  <- 1-(Lm[1]- glm_result$rank - 1)/L0[1] # adjusted McFadden
adjMcFadden2 <- 1-(Lm[1]- glm_result$rank)/L0[1] # adjusted McFadden 2

# Cox & Snell
n_model <- glm_result$rank + glm_result$df.residual; # n_model
CoxSnell <- 1-(exp((2/n_model)*(L0[1] - Lm[1])))  # Cox & Snell

# Nagelkerke
Nagelkerke <- (1-(exp((2/n_model)*(L0[1] - Lm[1]))))/(1-exp((2/n_model)*L0[1])) 

DescTools::PseudoR2(result2, "all")

       McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
      0.1103414       0.1093747       0.1064634       0.1664870       0.1011785 
VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
      0.2003557       0.1346103       0.2267494       0.1251926    9213.1806742 
            BIC          logLik         logLik0              G2 
   9249.3018906   -4601.5903371   -5172.3100486    1141.4394230

BaylorEdPsych::PseudoR2(result2)

        McFadden     Adj.McFadden        Cox.Snell       Nagelkerke 
       0.1103414        0.1091813        0.1064634        0.1664870 
McKelvey.Zavoina           Effron            Count        Adj.Count 
       0.2267494        0.1346103        0.8284024        0.1714286 
             AIC    Corrected.AIC 
    9213.1806742     9213.1865949

McFadden; adjMcFadden; adjMcFadden2; CoxSnell; Nagelkerke

[1] 0.1103414

[1] 0.1091813

[1] 0.1093747

[1] 0.1064634

[1] 0.166487

上記の例では（結果としては良かったが？）いずれも数値が一致した。他の例でずれた原因は分からなかったが，結果としてはこの作業は計算の正確さのチェックとなった。

一般化線型モデルの特殊例としての一般線型モデル

　最後に，一般線型モデル(LM)が，一般化線型モデル(GLM)の特殊ケースである事を示しておこう。　

# LMとGLM
d01$q1700[d01$q1700 == 11] <- NA

# lm( )関数による一般線型モデル
reg1 <- lm(q1700 ~ age + edu + q1502 + income, data=d01)
summary(reg1); AIC(reg1)


Call:
lm(formula = q1700 ~ age + edu + q1502 + income, data = d01)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.9357 -1.2564  0.1654  1.2907  3.6941 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.290e+00  1.065e+00   5.904 8.83e-09 ***
age          2.246e-04  1.248e-02   0.018  0.98565    
edu          1.574e-01  5.191e-02   3.032  0.00262 ** 
q1502       -3.052e-01  5.923e-02  -5.154 4.40e-07 ***
income      -9.991e-05  3.027e-04  -0.330  0.74152    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.768 on 330 degrees of freedom
  (49 observations deleted due to missingness)
Multiple R-squared:  0.1044,    Adjusted R-squared:  0.09358 
F-statistic:  9.62 on 4 and 330 DF,  p-value: 2.275e-07

[1] 1339.341

# glm( )関数による一般化線型モデル（恒等リンク関数・正規分布）
reg2 <- glm(q1700 ~ age + edu + q1502 + income, data=d01, family = gaussian(link = "identity"))
summary(reg2); PseudoR2(reg2)


Call:
glm(formula = q1700 ~ age + edu + q1502 + income, family = gaussian(link = "identity"), 
    data = d01)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.9357  -1.2564   0.1654   1.2907   3.6941  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.290e+00  1.065e+00   5.904 8.83e-09 ***
age          2.246e-04  1.248e-02   0.018  0.98565    
edu          1.574e-01  5.191e-02   3.032  0.00262 ** 
q1502       -3.052e-01  5.923e-02  -5.154 4.40e-07 ***
income      -9.991e-05  3.027e-04  -0.330  0.74152    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 3.124805)

    Null deviance: 1151.4  on 334  degrees of freedom
Residual deviance: 1031.2  on 330  degrees of freedom
  (49 observations deleted due to missingness)
AIC: 1339.3

Number of Fisher Scoring iterations: 2

        McFadden     Adj.McFadden        Cox.Snell       Nagelkerke 
      0.10443261       0.09401082       0.30158841       0.31160894 
McKelvey.Zavoina           Effron            Count        Adj.Count 
              NA       0.10443261       0.04776119      -0.26086957 
             AIC    Corrected.AIC 
   1041.18569634    1041.36806716

偏回帰係数の推定の部分やAICは全て一致している。その他，Multiple R-squaredとMcFadden pseudo R2が一致している。

library(memisc)
memisc::mtable(reg1, reg2, summary.stats=c("sigma", "R-squared", "adj. R-squared", "F", "p", "Likelihood-ratio", "Log-likelihood", "Deviance", "AIC", "BIC", "N", "McFadden R-sq.", "Cox-Snell R-sq.", "Nagelkerke R-sq."))


Calls:
reg1: lm(formula = q1700 ~ age + edu + q1502 + income, data = d01)
reg2: glm(formula = q1700 ~ age + edu + q1502 + income, family = gaussian(link = "identity"), 
    data = d01)

==============================================
                       reg1         reg2      
----------------------------------------------
  (Intercept)          6.290***     6.290***  
                      (1.065)      (1.065)    
  age                  0.000        0.000     
                      (0.012)      (0.012)    
  edu                  0.157**      0.157**   
                      (0.052)      (0.052)    
  q1502               -0.305***    -0.305***  
                      (0.059)      (0.059)    
  income              -0.000       -0.000     
                      (0.000)      (0.000)    
----------------------------------------------
  sigma                1.768                  
  R-squared            0.104                  
  adj. R-squared       0.094                  
  F                    9.620                  
  p                    0.000        0.000     
  Log-likelihood    -663.670     -663.670     
  Deviance          1031.186     1031.186     
  AIC               1339.341     1339.341     
  BIC               1362.226     1362.226     
  N                  335          335         
  Likelihood-ratio                120.247     
  McFadden R-sq.                    0.104     
  Cox-Snell R-sq.                   0.302     
  Nagelkerke R-sq.                  0.312     
==============================================
  Significance: *** = p < 0.001;   
                ** = p < 0.01; * = p < 0.05

『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第12章　ロジスティック回帰分析

SUGINO Isamu

2021-01-23 23:40:34, Build with R4.0.2

0　全体の章構成

1　一般化線型モデル(Generalized Linear Model)

1-1　ロジスティック関数

1-2　最尤法(Maximum Liklihood Estimation)

1-3　対数オッズ比(log odds ratio)と2項ロジットモデル(binomial logit model)

1-4　適合度検定(Goodness of Fit Test)と疑似決定係数(Pseudo R-squared)

疑似決定係数の定義と計算

追加パッケイジ無しで計算する

大標本でのDescToolsとBaylorEdPsych

一般化線型モデルの特殊例としての一般線型モデル

2　ロジットモデル

2-1　多項ロジット(Multinomial Logit Model)

2-2　順序ロジット(Ordinal Logit Model)

発展1　ログリニアモデル(log-linear model; 対数線型モデル)

第12章の【練習問題】の解答

1) の答え

2) の答え

3) の答え

4) の答え

5) の答え

『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第12章 ロジスティック回帰分析

SUGINO Isamu

2021-01-23 23:40:34, Build with R4.0.2

0 全体の章構成

1 一般化線型モデル(Generalized Linear Model)

1-1 ロジスティック関数

1-2 最尤法(Maximum Liklihood Estimation)

1-3 対数オッズ比(log odds ratio)と2項ロジットモデル(binomial logit model)

1-4 適合度検定(Goodness of Fit Test)と疑似決定係数(Pseudo R-squared)

疑似決定係数の定義と計算

追加パッケイジ無しで計算する

大標本でのDescToolsとBaylorEdPsych

一般化線型モデルの特殊例としての一般線型モデル

2 ロジットモデル

2-1 多項ロジット(Multinomial Logit Model)

2-2 順序ロジット(Ordinal Logit Model)

発展1 ログリニアモデル(log-linear model; 対数線型モデル)

第12章の【練習問題】の解答

1) の答え

2) の答え

3) の答え

4) の答え

5) の答え

第12章　ロジスティック回帰分析

0　全体の章構成

1　一般化線型モデル(Generalized Linear Model)

1-1　ロジスティック関数

1-2　最尤法(Maximum Liklihood Estimation)

1-3　対数オッズ比(log odds ratio)と2項ロジットモデル(binomial logit model)

1-4　適合度検定(Goodness of Fit Test)と疑似決定係数(Pseudo R-squared)

2　ロジットモデル

2-1　多項ロジット(Multinomial Logit Model)

2-2　順序ロジット(Ordinal Logit Model)

発展1　ログリニアモデル(log-linear model; 対数線型モデル)