# ←この記号の後は次の改行までコメントになります。
# 数値変数の記述統計量をグラフィカルに表現する。
# Lesson 13では名義変数の記述統計量を作図しました。
# 今回は数値変数の記述統計量を作図してみましょう。
# 今まで同様にまず準備から。
# 第3回講習会で使ったデータセットを使います。
# ここから46行目まではサンプルデータセットの作成です。
rm(list=ls())
set.seed(1111)
Dataset.0<-
data.frame(
Sex=factor(rnorm(n=800,mean=0.5,sd=1)>0,labels=c("Male","Female")),
Age=50+rpois(n=800,lambda=25)+round(rnorm(n=800,mean=0,sd=5),0),
Severity=rpois(n=800,lambda=3))
Dataset.0$Comorbidity=
rpois(n=800,lambda=1)+round(Dataset.0$Age*rnorm(n=800,mean=0.03,sd=0.08)^2,0)
Dataset.0$UnknownConfounder=
factor(
with(Dataset.0,
rnorm(n=800,mean=0)+
I(Sex=="Male")*rnorm(n=800,mean=0.02)*0.01+
I(Age-80)*rnorm(n=800,mean=0.002)*0.001+
I(Severity>3)*rnorm(n=800,mean=0.05)*0.01)>1,
labels=c("No","Yes"))
Dataset.0$Intervention<-
factor(
with(Dataset.0,
rnorm(n=800,mean=0)+
I(Sex=="Female")*rnorm(n=800,mean=0.2)+
I(Age-80)*rnorm(n=800,mean=-0.1)*0.1+
I(Severity-2)*rnorm(n=800,mean=0.4)+
I(Comorbidity-1)*rnorm(n=800,mean=-1.2)*0.6+
I(UnknownConfounder=="Yes")*rnorm(n=800,mean=0.2)*0.2)>0,
labels=c("No","Yes"))
Dataset.0$Outcome<-
factor(
with(Dataset.0,
rnorm(n=800,mean=0)*10+
I(Sex=="Male")*rnorm(n=800,mean=0.01)+
I(Age-75)*rnorm(n=800,mean=0.3)*5+
I(Severity-1)*rnorm(n=800,mean=0.7)*6+
I(Comorbidity-1)*rnorm(n=800,mean=0.9)*5+
I(Intervention=="Yes")*rnorm(n=800,mean=-0.6)*4+
I(UnknownConfounder=="Yes")*rnorm(n=800,mean=0.3)*0.7)>20,labels=c("No","Yes"))
Dataset.0$UnknownConfounder<-NULL
# 数値変数の特徴をつかむために最初に使う関数は、summary, hist, boxplotの各関数です。
# このうち、histとboxplotはグラフィックス関数です。
summary(Dataset.0$Age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 53.00 70.00 75.00 75.02 80.00 107.00
hist(Dataset.0$Age)
# 箱ひげ図は数値変数の分布を直感的に表示することができます。
# 箱の真ん中の線は中央値（50%）で、箱の上下端はそれぞれ、第3四分位値（75%）、第1四分位値（25%）です。
# 箱の下からの3つの線がそれぞれ第1、第2、第3の四分位値を表すということです。
# 箱の上下端の幅は四分位範囲 (Interquartile range, IQR) と言われ、中央値とともに連続変数の分布を示す
# 重要な記述統計量としてしばしば記載されます。
# ひげの上端は、箱の上端の上に、四分位範囲の幅の1.5倍さらに離れた位置か最大値のうち小さい方です。
# ひげの下端は、箱の下端の下に、四分位範囲の幅の1.5倍さらに離れた位置か最小値のうち大きい方です。
# ひげの上下端をも外れた値は外れ値とみなされ、丸印でプロットされます。
boxplot(Dataset.0$Age)
# 重症度（Severity）でも作図してみましょう。
# 分布の違いとヒストグラムの表現の違いに注意して見て下さい。
hist(Dataset.0$Severity)
boxplot(Dataset.0$Severity)
# 名義変数（Intervention）で群分けしてみる。
# boxplot関数ではxtabs関数などと同様にdata= でデータセットの指定ができます。
boxplot(Age~Intervention,data=Dataset.0)
boxplot(Severity~Intervention,data=Dataset.0)
# 両群の年齢と重症度にも大きな差は無いように見えます。検定してみましょう。
# boxplotの部分をt.testに変えると検定ができます。
# 正規分布を仮定
t.test(Age~Intervention,data=Dataset.0)
##
## Welch Two Sample t-test
##
## data: Age by Intervention
## t = -0.20567, df = 770.65, p-value = 0.8371
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.1030226 0.8938104
## sample estimates:
## mean in group No mean in group Yes
## 74.96765 75.07226
t.test(Severity~Intervention,data=Dataset.0)
##
## Welch Two Sample t-test
##
## data: Severity by Intervention
## t = -3.9783, df = 789.04, p-value = 7.581e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.7221513 -0.2449571
## sample estimates:
## mean in group No mean in group Yes
## 2.768194 3.251748
# 非正規分布を仮定
wilcox.test(Age~Intervention,data=Dataset.0)
##
## Wilcoxon rank sum test with continuity correction
##
## data: Age by Intervention
## W = 78858, p-value = 0.8248
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(Severity~Intervention,data=Dataset.0)
##
## Wilcoxon rank sum test with continuity correction
##
## data: Severity by Intervention
## W = 65930, p-value = 2.081e-05
## alternative hypothesis: true location shift is not equal to 0
# いずれも重症度は有意に違うようです。ついでに群分けしてsummaryで見てみます。
# 箱ひげ図と見比べて下さい。
# by関数の復習です。
by(Dataset.0$Severity,Dataset.0$Intervention,summary)
## Dataset.0$Intervention: No
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 3.000 2.768 4.000 9.000
## --------------------------------------------------------
## Dataset.0$Intervention: Yes
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 3.000 3.252 4.000 9.000
# グラフに必要な情報を盛り込んで飾りましょう。
boxplot(
Severity~Intervention,
data=Dataset.0,
main="Comparison of disease severity between the groups",
names=c("Intervention (-)","Intervention (+)"),
ylab=("Disease severity"),
ylim=c(0,10)
)
# グラフィックスは数学的に妥当で説得力のある表現をするにはとても重要です。
# Rには様々なグラフィックスが用意されていて、ただの統計ソフトとは思えない
# ほど強力ですが、一方であまりに機能が豊富すぎて学ぶのが大変です。本も幾つ
# かありますが、以前も紹介したRTipsが学びやすくてとても良いです。
# http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html
# この47-58章が非常によくできています。