# ←この記号の後は次の改行までコメントになります。
# データフレームにアクセスして記述統計量を求める（名義変数）。
# まずはデータ空間をクリアしておいてから、
rm(list=ls())
# ここから43行目まではサンプルデータセットの作成です。
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
# Dataset.0というデータフレームが出来上がっているはずです。
# RStudioを使っている方は、右上のウィンドウのEnvironmentタグに表示されています。
# いつものごとく、中身を見てみましょう。
nrow(Dataset.0) #サンプル数
## [1] 800
head(Dataset.0,20) #先頭から20行を見ます。
## Sex Age Severity Comorbidity Intervention Outcome
## 1 Female 77 0 1 Yes No
## 2 Female 64 0 3 Yes Yes
## 3 Female 84 3 0 Yes Yes
## 4 Female 67 3 0 No No
## 5 Female 70 2 4 No No
## 6 Male 87 3 0 Yes No
## 7 Female 82 4 0 Yes No
## 8 Female 79 5 1 No Yes
## 9 Female 64 5 0 No Yes
## 10 Female 80 6 3 Yes Yes
## 11 Male 61 1 3 No No
## 12 Female 76 4 4 Yes No
## 13 Male 72 3 2 No No
## 14 Male 83 3 1 Yes No
## 15 Female 77 2 0 No No
## 16 Female 77 0 3 No No
## 17 Female 78 2 0 Yes No
## 18 Male 86 0 2 No Yes
## 19 Female 77 3 2 Yes Yes
## 20 Female 76 3 0 Yes No
summary(Dataset.0) #データセット全体の要約
## Sex Age Severity Comorbidity Intervention
## Male :254 Min. : 53.00 Min. :0.000 Min. :0.00 No :371
## Female:546 1st Qu.: 70.00 1st Qu.:2.000 1st Qu.:1.00 Yes:429
## Median : 75.00 Median :3.000 Median :1.00
## Mean : 75.02 Mean :3.027 Mean :1.48
## 3rd Qu.: 80.00 3rd Qu.:4.000 3rd Qu.:2.00
## Max. :107.00 Max. :9.000 Max. :8.00
## Outcome
## No :511
## Yes:289
##
##
##
##
# すでにsummary関数でおよその推測統計量は出ていますね。
# 名義変数
# 論文に名義変数の記述統計量を書く場合（例えばTable 1.）、通常はカウントと割合が
# 要求されます。性別（Sex）を対象としてやってみます。
table(Dataset.0$Sex) # カウント数
##
## Male Female
## 254 546
table(Dataset.0$Sex)/sum(table(Dataset.0$Sex)) # sum関数で総数を求めて割合を出した。
##
## Male Female
## 0.3175 0.6825
# 結果を代入しておくと何かと便利です。
tab<-table(Dataset.0$Sex)
tab
##
## Male Female
## 254 546
tab/sum(tab)
##
## Male Female
## 0.3175 0.6825
# 群分けしたカウントデータ（2X2表）を作ってみます。
# 治療介入で群分けし、性別のカウントデータを求めます。
table(Dataset.0$Sex,Dataset.0$Intervention)
##
## No Yes
## Male 115 139
## Female 256 290
# xtabs関数（Cross tablesの意）を使って同じ動作ができます。
# 引数の渡し方がtable関数と少し異なることに注意。
# xtabs(~名義変数1+名義変数2,data=データセット)という書式です。
xtabs(~Sex,data=Dataset.0)
## Sex
## Male Female
## 254 546
xtabs(~Sex+Intervention,data=Dataset.0)
## Intervention
## Sex No Yes
## Male 115 139
## Female 256 290
# そのまま2X2表の検定をすることもできます。
# カイ二乗検定にはchisq.testを、Fisherの正確検定にはfisher.testを使います。
chisq.test(xtabs(~Sex+Intervention,data=Dataset.0))
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: xtabs(~Sex + Intervention, data = Dataset.0)
## X-squared = 0.12191, df = 1, p-value = 0.727
fisher.test(xtabs(~Sex+Intervention,data=Dataset.0))
##
## Fisher's Exact Test for Count Data
##
## data: xtabs(~Sex + Intervention, data = Dataset.0)
## p-value = 0.7036
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.6870372 1.2774363
## sample estimates:
## odds ratio
## 0.9373191
# table関数でも同様のことができます。
chisq.test(table(Dataset.0$Sex,Dataset.0$Intervention))
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(Dataset.0$Sex, Dataset.0$Intervention)
## X-squared = 0.12191, df = 1, p-value = 0.727
fisher.test(table(Dataset.0$Sex,Dataset.0$Intervention))
##
## Fisher's Exact Test for Count Data
##
## data: table(Dataset.0$Sex, Dataset.0$Intervention)
## p-value = 0.7036
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.6870372 1.2774363
## sample estimates:
## odds ratio
## 0.9373191
# xtabs関数とtable関数の動作はほとんど同じです。使いやすい方を憶えておけばよいです。
# Table1の一行を作ってみましょう。ちょっと難しいですが、頑張って読み解いて下さい。
# 治療介入で群分けして、それぞれ女性の割合を求めてみます。
# ~名義変数1+名義変数2 のうち、名義変数2を群分け変数とみなしましょう。
tab<-xtabs(~Sex+Intervention,data=Dataset.0) #クロス集計表の作成
tab #集計表の内容確認
## Intervention
## Sex No Yes
## Male 115 139
## Female 256 290
row.sex<- # 性別のカウントと割合を群分けした1行を作ります。
c("Female sex", # 行の最初の要素は変数のタイトルです
paste( # 文字列を結合するpaste関数を使っています。
tab[2,1], # 2番めの要素は治療なし群の女性のカウントと割合。
" (", # 割合は括弧で囲む
round(tab[2,1]/sum(tab[,1])*100,0),
# 割合を計算しround関数で四捨五入
")", # 割合は括弧で囲む
sep=""), # paste関数のの区切りは、なしを指定
paste( # 3番目の要素は治療あり群の女性のカウントと割合。以下同様。
tab[2,2],
" (",
round(tab[2,2]/sum(tab[,2])*100,0),
")",
sep=""),
format( # 4番目の要素はP-valueです。chisq.testの結果から
round( # 取り出して、小数点以下3桁までを丸めて（round）
chisq.test(tab)$p.value,3),
nsmall=3)) # 小数点以下3桁までの表示を指定しています（format）
# コメント行無しでもう一度やってみます。スクリプトの視認性はこちらが良い
row.sex<-
c("Female sex",
paste(tab[2,1]," (",round(tab[2,1]/sum(tab[,1])*100,0),")",sep=""),
paste(tab[2,2]," (",round(tab[2,2]/sum(tab[,2])*100,0),")",sep=""),
format(round(chisq.test(tab)$p.value,3),nsmall=3))
# タイトル、2群のカウントデータと割合、P-valueが、長さ4のベクターに
# 格納されています。
row.sex
## [1] "Female sex" "256 (69)" "290 (68)" "0.727"
# さらにrbind関数（行を繋げる関数）を使って、表の一番上の部分を作ってみましょう。
Table1<-
rbind(
c("Variables","Control","Active","P-value"),
row.sex)
# オブジェクトTable1は2行4列のマトリックスです。
Table1
## [,1] [,2] [,3] [,4]
## "Variables" "Control" "Active" "P-value"
## row.sex "Female sex" "256 (69)" "290 (68)" "0.727"
# このように1行ずつ集計していくことで記述統計の表を作ることができます。
# 次回は連続変数の集計にチャレンジします。