Tsuruta

このウェブページは、以下の書籍（以下、テキストという）を R で実行してみた記録です。テキストは STATA や SPSS に触れていて、web 付録にはエクセルファイルがありますが、R については触れていません。web 付録をダウンロード・展開した上で、試してみてください。

すべての医療系学生・研究者に贈る独習統計学応用編24講 ―分割表・回帰分析・ロジスティック回帰―

https://www.asakura.co.jp/detail.php?book_code=12217

鶴田陽和(著)

定価 3,850 円（本体 3,500 円＋税）

A5判／248ページ
刊行日：2016年10月25日
ISBN：978-4-254-12217-6　C3041

Amazon

第4講

定義・説明なしに自由度が登場している。例えば、10人を観察して身長を測る場合、その値は何にも影響を受けない (A さんの身長は B に影響を受けない) ため、値は自由である。この場合、自由度はサンプルサイズと同じく 10 である。

一方で、平均値を決定した場合、9人目までは自由に値を取れるが、10人目は平均に合わせて値を決めなければならない。この場合は自由度が 9 になる。

p.26 図4.4 は、web 付録の付録_基本編.xlsx の「§4_t分布のグラフ」シートによると、\(\infty\) は現実には 100 である。

同様の図を R で描くには、データは特に必要ではない。

curve(dt(x, df=100), from=-4, to=4)
curve(dt(x, df=3), from=-4, to=4, add=TRUE)
curve(dt(x, df=1), from=-4, to=4, add=TRUE)

p. 27 の Note 1 について、テキストでは結果が書かれていないが、Excel で計算すると 2.776445105 が得られるはずである。

これを同じことを R で行うには、以下のとおり。

qt(0.975, 4)

## [1] 2.776445

第6講

p.37 の図6.2 は、t 分布と同様です。

curve(dchisq(x, df = 1), from = 0, to = 20)
curve(dchisq(x, df = 2), from = 0, to = 20, add=TRUE)
curve(dchisq(x, df = 3), from = 0, to = 20, add=TRUE)
curve(dchisq(x, df = 4), from = 0, to = 20, add=TRUE)
curve(dchisq(x, df = 10), from = 0, to = 20, add=TRUE)

第13講

library(readxl)
df13 <- read_excel("Appendix/計算例/付録_回帰分析編.xlsx", 
    sheet = "図13.1左",
    range = "C7:D27")

library(ggplot2)
gg1 <- ggplot(df13, aes(X, Y1)) +
  geom_point() +
  geom_abline()
gg1

cor(df13$X, df13$Y1)

## [1] 0.8926096

library(ggplot2)
gg1 <- ggplot(df13, aes(X, Y1)) +
  geom_point() +
  annotate(geom = "text", label = "r = 0.893", x =5, y = -5)
gg1

library(ggplot2)
gg1 <- ggplot(df13, aes(X, Y1)) +
  geom_point() +
  annotate(geom = "text", label = paste("r =", cor(df13$X, df13$Y1)), x =5, y = -5)
gg1

library(ggplot2)
gg1 <- ggplot(df13, aes(X, Y1)) +
  geom_point() +
  annotate(geom = "text", label = paste("r =", round(cor(df13$X, df13$Y1), digits = 3)), x =5, y = -5)
gg1

df2 <- read_excel("Appendix/計算例/付録_回帰分析編.xlsx", 
    sheet = "図13.1中",
    range = "C7:D27")
gg2 <- ggplot(df2, aes(X, Y2)) +
  geom_point() +
  annotate(geom = "text", label = paste("r =", round(cor(df2$X, df2$Y2), digits = 3)), x =5, y = -5)
gg2

df3 <- read_excel("Appendix/計算例/付録_回帰分析編.xlsx", 
    sheet = "図13.1右",
    range = "C7:D27")
gg3 <- ggplot(df3, aes(X, Y3)) +
  geom_point() +
  annotate(geom = "text", label = paste("r =", round(cor(df3$X, df3$Y3), digits = 3)), x =5, y = 5)
gg3

library(gridExtra)
gridExtra::grid.arrange(gg1, gg2, gg3, nrow = 1)

第14講

library(readxl)
df14.1 <- read_excel("Appendix/計算例/付録_回帰分析編.xlsx", 
    sheet = "§14_15_年齢_体重_成績",
    range = "B8:D38")

library(ggplot2)
gg14.1 <- ggplot(df14.1, aes(年齢, 体重)) +
  geom_point() +
  theme(text = element_text(family = "HiraginoSans-W3"))
gg14.1

library(readxl)
df14.4 <- read_excel("Appendix/計算例/付録_回帰分析編.xlsx", 
    sheet = "図14.4_信頼区間_図",
    range = "B6:C14")

library(ggpubr)

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:huxtable':
## 
##     font

ggscatter(data = df14.4, x = "ｘ", y = "ｙ", add = "reg.line", conf.int = TRUE) +
  theme(text = element_text(family = "HiraginoSans-W3"))

注記: エクセルデータ中、x ではなく全角のｘ、y ではなく全角のｙ。

第16講

library(readxl)
df16 <- read_excel("Appendix/計算例/付録_回帰分析編.xlsx", 
    sheet = "図16_2",
    range = "B6:E36")

library(ggplot2)
gg16 <- ggplot(df16, aes(予測値, 成績)) +
  geom_point() +
  theme(text = element_text(family = "HiraginoSans-W3"))
gg16

cor(df16$予測値, df16$成績)

## [1] 0.8649051

df16.2 <- subset(df16, select = c(年齢, 体重, 成績))

m1 <- lm(成績 ~ 年齢 + 体重, data = df16.2)
summary(m1)

## 
## Call:
## lm(formula = 成績 ~ 年齢 + 体重, data = df16.2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.65054 -0.56227  0.08159  0.59756  2.23442 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.22806    1.02205  -4.137 0.000308 ***
## 年齢         1.15289    0.20043   5.752 4.08e-06 ***
## 体重        -0.05907    0.06110  -0.967 0.342191    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9776 on 27 degrees of freedom
## Multiple R-squared:  0.7481, Adjusted R-squared:  0.7294 
## F-statistic: 40.08 on 2 and 27 DF,  p-value: 8.27e-09

anova(m1)

Df	Sum Sq	Mean Sq	F value	Pr(>F)
1	75.7	75.7	79.2	1.62e-09
1	0.894	0.894	0.935	0.342
27	25.8	0.956

p. 139 基本統計量

summary(df16.2)

##       年齢             体重            成績      
##  Min.   : 6.286   Min.   :12.90   Min.   :1.575  
##  1st Qu.: 7.556   1st Qu.:19.31   1st Qu.:3.420  
##  Median : 8.666   Median :22.96   Median :4.274  
##  Mean   : 8.969   Mean   :22.53   Mean   :4.782  
##  3rd Qu.:10.482   3rd Qu.:26.18   3rd Qu.:6.382  
##  Max.   :11.810   Max.   :31.94   Max.   :7.744

上の結果は、テキストとは縦横が異なる。

res <- cor(df16.2)
round(res, 3)

##       年齢  体重  成績
## 年齢 1.000 0.831 0.860
## 体重 0.831 1.000 0.663
## 成績 0.860 0.663 1.000

第21講

library(readxl)

df21 <- read_excel("Appendix/計算例/付録_Logistic回帰編.xlsx", 
    sheet = "表21.1_ロジスティック回帰のデータ",
    range = "B25:E185")
df21$HyperT <- factor(df21$HyperT, labels = c("高血圧無し", "高血圧有り"))
df21$Age <- factor(df21$Age, labels = c("65未満", "65以上"))
df21$Gender <- factor(df21$Gender, labels = c("男性", "女性"))
df21$AKI <- factor(df21$AKI, labels = c("no AKI", "AKI"))

xtabs(~ Gender + Age + AKI, data = df21)

## , , AKI = no AKI
## 
##       Age
## Gender 65未満 65以上
##   男性     31     25
##   女性     34     28
## 
## , , AKI = AKI
## 
##       Age
## Gender 65未満 65以上
##   男性      9     15
##   女性      6     12

ftable(xtabs(~ HyperT + Age + Gender + AKI, data = df21))

##                          AKI no AKI AKI
## HyperT     Age    Gender               
## 高血圧無し 65未満 男性           17   3
##                   女性           18   2
##            65以上 男性           15   5
##                   女性           16   4
## 高血圧有り 65未満 男性           14   6
##                   女性           16   4
##            65以上 男性           10  10
##                   女性           12   8

xtabs と ftable は、Base R (追加パッケージ不要) ではあるが、出力がシンプルすぎます。

クロス集計を綺麗に出力するパッケージはたくさん有りますが、高血圧・年齢・性別と段階的に表示するのは modelsummary くらいしかなさそうです。crosstable もできますが、表示が独特すぎる。

library(modelsummary)
datasummary(HyperT * Age * Gender ~ AKI, data = df21, fmt = fmt_sprintf("%.0f"))

HyperT	Age	Gender	no AKI	AKI
高血圧無し	65未満	男性	17	3
		女性	18	2
	65以上	男性	15	5
		女性	16	4
高血圧有り	65未満	男性	14	6
		女性	16	4
	65以上	男性	10	10
		女性	12	8

output = “gt”, “kableExtra”, “flextable”, “huxtable” などのオプションがあります。

library(huxtable)
datasummary(HyperT * Age * Gender ~ AKI, data = df21, fmt = fmt_sprintf("%.0f"), output = "huxtable") |>
   set_background_color(2:9, 1:3, "lightgrey")

HyperT	Age	Gender	no AKI	AKI
高血圧無し	65未満	男性	17	3
		女性	18	2
	65以上	男性	15	5
		女性	16	4
高血圧有り	65未満	男性	14	6
		女性	16	4
	65以上	男性	10	10
		女性	12	8

strTemp <- datasummary(HyperT * Age * Gender ~ AKI, data = df21, fmt = fmt_sprintf("%.0f"), output = "markdown")
dfTemp <- read.table(text = strTemp, header = TRUE, sep = '|')
dfTemp <- dfTemp[-1,c(2,3,4,5,6)]
dfTemp$AKI <- as.numeric(dfTemp$AKI)
dfTemp$no.AKI <- as.numeric(dfTemp$no.AKI)
dfTemp$Risk <- dfTemp$AKI / (dfTemp$no.AKI + dfTemp$AKI)
#dfTemp <- data.frame(dfTemp[,"Risk"])

library(huxtable)
hux17 <- hux(dfTemp) |>
   set_background_color(2:9, 1:3, "lightgrey") |>
   set_background_color(1, 4:6, "lightgrey")
rowspan(hux17)[c(2,6), "HyperT"] <- 4
rowspan(hux17)[c(2,4,6,8), "Age"] <- 2
hux17

HyperT	Age	Gender	no.AKI	AKI	Risk
高血圧無し	65未満	男性	17	3	0.15
	65未満	女性	18	2	0.1
	65以上	男性	15	5	0.25
	65以上	女性	16	4	0.2
高血圧有り	65未満	男性	14	6	0.3
	65未満	女性	16	4	0.2
	65以上	男性	10	10	0.5
	65以上	女性	12	8	0.4

さて、ここまで表を作ってきましたが、この表は解析に使うわけではありません。必要なのは最初にエクセルから読み込んだデータフレーム (df21) です。

R でロジスティック回帰分析を実行して、テキスト p. 180 と比較してみましょう。テキストは Stata の出力です。

model21 <- glm(AKI ~ HyperT + Age + Gender, 
                family = binomial(logit), 
                data = df21)
summary(model21)

## 
## Call:
## glm(formula = AKI ~ HyperT + Age + Gender, family = binomial(logit), 
##     data = df21)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -1.8270     0.4144  -4.409 1.04e-05 ***
## HyperT高血圧有り   0.9689     0.3844   2.520   0.0117 *  
## Age65以上          0.8327     0.3816   2.182   0.0291 *  
## Gender女性        -0.4192     0.3761  -1.115   0.2650    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184.21  on 159  degrees of freedom
## Residual deviance: 171.64  on 156  degrees of freedom
## AIC: 179.64
## 
## Number of Fisher Scoring iterations: 4

表示方法が少し (かなり？) 違います。まず、順序としては R では (Intercept) が最初ですが、Stata では一番最後の _cons です。Intercept は切片の英語で、_cons は定数 (constant) を表しています。

列の名称も少し異なります。最初の項目は、R では Estimate ですが、Stata は Coef. となっています。この値は、オッズ比の対数 (底は e) となっています。つまり、高血圧のオッズ比は、以下のように計算することができます。

exp(0.9689)

## [1] 2.635044

これを一気に行いましょう。テキスト p. 184 の図21.2 とほぼ同じ内容になっています。

library("epiDisplay", warn.conflicts=FALSE, quietly=TRUE)
logistic.display(model21)

## 
## Logistic regression predicting AKI : AKI vs no AKI 
##  
##                                  crude OR(95%CI)   adj. OR(95%CI)   
## HyperT: 高血圧有り vs 高血圧無し 2.54 (1.21,5.31)  2.64 (1.24,5.6)  
##                                                                     
## Age: 65以上 vs 65未満            2.21 (1.07,4.57)  2.3 (1.09,4.86)  
##                                                                     
## Gender: 女性 vs 男性             0.68 (0.33,1.38)  0.66 (0.31,1.37) 
##                                                                     
##                                  P(Wald's test) P(LR-test)
## HyperT: 高血圧有り vs 高血圧無し 0.012          0.01      
##                                                           
## Age: 65以上 vs 65未満            0.029          0.026     
##                                                           
## Gender: 女性 vs 男性             0.265          0.263     
##                                                           
## Log-likelihood = -85.8188
## No. of observations = 160
## AIC value = 179.6375

第24講

まず、データがエクセルファイル中にあるので読み込みます。

library(readxl)
df24 <- read_excel("Appendix/計算例/付録_Logistic回帰編.xlsx", 
    sheet = "§24_ROC曲線下面積の比較",
    range = "B11:D61")
df24$diagnosis <- as.factor(df24$diagnosis)

テキストでは、ROC を描く前に AUC を求めています。

library(pROC)

roc(diagnosis ~ reader_1 + reader_2, data = df24)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## $reader_1
## 
## Call:
## roc.formula(formula = diagnosis ~ reader_1, data = df24)
## 
## Data: reader_1 in 22 controls (diagnosis 0) < 28 cases (diagnosis 1).
## Area under the curve: 0.8701
## 
## $reader_2
## 
## Call:
## roc.formula(formula = diagnosis ~ reader_2, data = df24)
## 
## Data: reader_2 in 22 controls (diagnosis 0) < 28 cases (diagnosis 1).
## Area under the curve: 0.9635

それぞれ、0.8701 と 0.9635 となり、テキストと一致しています。

図24.5 は、以下のように作成することができます。見分けられるように reader_2 を点線 (lty = 2) とします。

library(pROC)

plot.roc(df24$diagnosis, df24$reader_1)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

plot.roc(df24$diagnosis, df24$reader_2, lty = 2, add = TRUE)

## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

ggplot2 を使うと、もっと綺麗に出すことができます。ただし、点を描くオプションはない模様。

library(ggplot2)
ggroc(roc(diagnosis ~ reader_1 + reader_2, data = df24),
      aes=c("linetype","color"),
      legacy.axes = TRUE) +
geom_abline(color="dark grey", slope = 1, intercept = 0, size=0.5)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

テキストでは、カットオフポイントの決定方法までは解説されていません。

カットオフの決め方の代表的なものには、以下のようなものがあります。R のパッケージで求めることができます。カッコ内は、順に pROC の print.thres.best.method, OptimalCutpoints での methods, cutpointr の metric での指定方法です。

感度 + (1 - 特異度) が最大の点 (youden, Youden, youden)
左上からの距離が最小の点 (closest.topleft, ROC01, roc01)
感度 × 特異度が最大の点 (…, MaxProdSpSe, prod_sens_spec)
その他多数（CRANパッケージ内の説明を参照）