今回はロジスティック回帰編です。
今回もRは組み込みのタイタニックのデータを使ってみました。

data.Titanic

titanic<- data.frame(Titanic)
head(titanic)

> head(titanic)
  Class    Sex   Age Survived Freq
1   1st   Male Child       No    0
2   2nd   Male Child       No    0
3   3rd   Male Child       No   35
4  Crew   Male Child       No    0
5   1st Female Child       No    0
6   2nd Female Child       No    0

生存したかどうかの2値データの解析なのでロジスティック回帰をやっていきます。
分析によさげなカテゴリカルデータがあるものの頻度データになっているのでどうしようか思っていたらいい記事がありました。

やっぱこういう記事あるんですね。
glmにweightsという引数があるのを初めて知りました。
集計データでは、weightsに頻度を入れればそのまま使えるんですね。

Classが4つのカテゴリに分けられた変数であるため、これについて多重比較します。
まずはモデルをチェック。

glm

fit<-glm(Survived~Class+Sex+Age,weights=Freq,family=binomial,data=titanic)
summary(fit)

> summary(fit)

Call:
glm(formula = Survived ~ Class + Sex + Age, family = binomial, 
    data = titanic, weights = Freq)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.6853     0.2730   2.510   0.0121 *  
Class2nd     -1.0181     0.1960  -5.194 2.05e-07 ***
Class3rd     -1.7778     0.1716 -10.362  < 2e-16 ***
ClassCrew    -0.8577     0.1573  -5.451 5.00e-08 ***
SexFemale     2.4201     0.1404  17.236  < 2e-16 ***
AgeAdult     -1.0615     0.2440  -4.350 1.36e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2769.5  on 23  degrees of freedom
Residual deviance: 2210.1  on 18  degrees of freedom
AIC: 2222.1

Number of Fisher Scoring iterations: 5

ただこれだとサンプル数が少ない分Dispersion parameterがちゃんと計算できないですね。
までも一応多重比較をやってみます。

多重比較

library(multcomp)
multicomp<-glht(fit,mcp(Class="Tukey"))
summary(multicomp)

> summary(multicomp)

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = Survived ~ Class + Sex + Age, family = binomial, 
    data = titanic, weights = Freq)

Linear Hypotheses:
                Estimate Std. Error z value Pr(>|z|)    
2nd - 1st == 0   -1.0181     0.1960  -5.194  < 1e-04 ***
3rd - 1st == 0   -1.7778     0.1716 -10.362  < 1e-04 ***
Crew - 1st == 0  -0.8577     0.1573  -5.451  < 1e-04 ***
3rd - 2nd == 0   -0.7597     0.1764  -4.308 0.000115 ***
Crew - 2nd == 0   0.1604     0.1738   0.923 0.790536    
Crew - 3rd == 0   0.9201     0.1486   6.192  < 1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

とりあえず実行できました。
調べてみると頻度に応じてデータを増やすやり方もわかりました。

同じようにやってみます。

頻度に応じてデータを複製

titanic2 <- data.frame(lapply(titanic,function(i){rep(i,titanic[,5])}))
#Freqの列を削除
titanic2 <- titanic2[,-5]
head(titanic2)

> head(titanic2)
  Class  Sex   Age Survived
1   3rd Male Child       No
2   3rd Male Child       No
3   3rd Male Child       No
4   3rd Male Child       No
5   3rd Male Child       No
6   3rd Male Child       No

いい感じです。
もう一度glmをやってみます。

fit2<-glm(Survived~Class+Sex+Age,family=binomial,data=titanic2)
summary(fit2)

> summary(fit2)

Call:
glm(formula = Survived ~ Class + Sex + Age, family = binomial, 
    data = titanic2)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.6853     0.2730   2.510   0.0121 *  
Class2nd     -1.0181     0.1960  -5.194 2.05e-07 ***
Class3rd     -1.7778     0.1716 -10.362  < 2e-16 ***
ClassCrew    -0.8577     0.1573  -5.451 5.00e-08 ***
SexFemale     2.4201     0.1404  17.236  < 2e-16 ***
AgeAdult     -1.0615     0.2440  -4.350 1.36e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2769.5  on 2200  degrees of freedom
Residual deviance: 2210.1  on 2195  degrees of freedom
AIC: 2222.1

Number of Fisher Scoring iterations: 4

Dispersion parameterもほぼ1になりそうですね。
一応計算します。

> sum(residuals(fit2,type ="pearson")^2)/fit2$df.residual
[1] 1.023531

performanceパッケージで検定してみます。

library(performance)
check_overdispersion(fit2)

> check_overdispersion(fit2)
 エラー: Overdispersion checks cannot be used for Bernoulli models.

ゼロイチデータでは使えないっぽいです。
前回と同じくどの分布が適切なのかランダムフォレストのアルゴリズムで調べてくれる関数を利用してみます。

> check_distribution(fit2)
# Distribution of Model Family

Predicted Distribution of Residuals

         Distribution Probability
 poisson (zero-infl.)         22%
               cauchy         19%
               normal         19%

Predicted Distribution of Response

  Distribution Probability
     bernoulli         88%
 beta-binomial          9%
      binomial          3%

問題なさそうですね。
最初の重み付けしたglmと違いがないか、念のためこちらでも多重比較を実行します。

multicomp2<-glht(fit2,mcp(Class="Tukey"))
summary(multicomp2)

> summary(multicomp2)

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = Survived ~ Class + Sex + Age, family = binomial, 
    data = titanic2)

Linear Hypotheses:
                Estimate Std. Error z value Pr(>|z|)    
2nd - 1st == 0   -1.0181     0.1960  -5.194   <1e-04 ***
3rd - 1st == 0   -1.7778     0.1716 -10.362   <1e-04 ***
Crew - 1st == 0  -0.8577     0.1573  -5.451   <1e-04 ***
3rd - 2nd == 0   -0.7597     0.1764  -4.308   <1e-04 ***
Crew - 2nd == 0   0.1604     0.1738   0.923    0.791    
Crew - 3rd == 0   0.9201     0.1486   6.192   <1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

特に変わらなくて安心しました。

前回と同じようにグラフに結果を描画してみます。
生存しているかどうかがYesとNoの文字列なので1と0に変換します。
ゼロイチデータなので今回はバブルチャートにしてみました。

グラフに描画

alphabet = cld(multicomp2,decreasing=F)
alphabet

library(tidyverse)
titanic2 <- dplyr::mutate(titanic2, Survived = str_replace_all(Survived, pattern = c("No" = "0", "Yes" = "1")))
titanic2$Survived<-as.numeric(titanic2$Survived)
abc = alphabet[["mcletters"]][["Letters"]]
plot = 
  ggplot(data = titanic2, aes(x = Class, y = Survived))+
  geom_count()+
  stat_summary(geom = 'text', label =abc, fun = max, vjust = -1)+
  theme_classic()
print(plot)