# ←この記号の後は次の改行までコメントになります。
# データフレームにアクセスして記述統計量を求める（連続変数）。
# まずはデータ空間をクリアしておいてから、
rm(list=ls())
# 乱数的に作成した架空のデータセットを使います。
# ここから44行目まではサンプルデータセットの作成です。
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.）
# 平均と標準偏差で表すか、中央値と四分位範囲（25 percentile - 75 percentile）で
# 表されるか、ほとんどの場合はこのどちらかです。
# 年齢（Age）を用いて、どちらもやってみましょう。
mean(Dataset.0$Age) # 平均
## [1] 75.02375
sd(Dataset.0$Age) #標準偏差
## [1] 7.14717
# round関数を使って、小数点以下１桁（digits=1）と有効数字を指定します。
round(c(mean(Dataset.0$Age),sd(Dataset.0$Age)),digits=1)
## [1] 75.0 7.1
# さらに文字列に変えることもできます。この時も小数点以下１桁まで表示桁数を
# 指定する（nsmall=1）ことが必要です。
format(round(c(mean(Dataset.0$Age),sd(Dataset.0$Age)),digits=1),nsmall=1)
## [1] "75.0" " 7.1"
# よく論文にあるような表記にしてみましょう。
# 前回もでてきたpaste関数を使って、文字列を貼り合わせます。
paste(
format(round(mean(Dataset.0$Age),digits=1),nsmall=1), # 平均値、小数点以下１桁
" +/- ", # " +/- "を間に挟む
format(round(sd(Dataset.0$Age),digits=1),nsmall=1), # 標準偏差、小数点以下１桁
sep="") # 区切り記号はなし
## [1] "75.0 +/- 7.1"
# Interventionで群分けしてやってみましょう。
# 群わけのやりかたには、いくつか方法があります。
# subsetを使う
mean(subset(Dataset.0,Intervention=="Yes")$Age) # subsetの中ではDataset.0$Intervention==という表記
## [1] 75.07226
sd(subset(Dataset.0,Intervention=="Yes")$Age) # をただIntervention==と省略できます（Tips）。
## [1] 7.006799
mean(subset(Dataset.0,Intervention=="No")$Age)
## [1] 74.96765
sd(subset(Dataset.0,Intervention=="No")$Age)
## [1] 7.315221
# [...]bracketを使う
mean(Dataset.0[Dataset.0$Intervention=="Yes",]$Age)
## [1] 75.07226
sd(Dataset.0[Dataset.0$Intervention=="Yes",]$Age)
## [1] 7.006799
mean(Dataset.0[Dataset.0$Intervention=="No",]$Age)
## [1] 74.96765
sd(Dataset.0[Dataset.0$Intervention=="No",]$Age)
## [1] 7.315221
# by関数をつかう。データセット（ベクター）、群分けの条件式、関数を
# 引数に指定します。ちょっととっつきにくいですが、慣れるとスマートです。
by(Dataset.0$Age,Dataset.0$Intervention=="Yes",mean)
## Dataset.0$Intervention == "Yes": FALSE
## [1] 74.96765
## --------------------------------------------------------
## Dataset.0$Intervention == "Yes": TRUE
## [1] 75.07226
by(Dataset.0$Age,Dataset.0$Intervention=="Yes",sd)
## Dataset.0$Intervention == "Yes": FALSE
## [1] 7.315221
## --------------------------------------------------------
## Dataset.0$Intervention == "Yes": TRUE
## [1] 7.006799
# 中央値 [四分位範囲]の記載もやってみましょう。
quantile(Dataset.0$Age)
## 0% 25% 50% 75% 100%
## 53 70 75 80 107
quantile(Dataset.0$Age)[c("50%","25%","75%")] # 要素に付けられた名前を元に取り出します。
## 50% 25% 75%
## 75 70 80
quantile(Dataset.0$Age)[c(3,2,4)] # ３番目、２番め、４番目という指定でもOK。
## 50% 25% 75%
## 75 70 80
# これも論文用の表記にしてみましょう。
Qtile<-quantile(Dataset.0$Age)
paste(
format(round(Qtile[3],digits=0),nsmall=0), # 中央値、小数点以下0桁
" [", # [
format(round(Qtile[2],digits=0),nsmall=0), # 25パーセンタイル、小数点以下0桁
", ", # ,
format(round(Qtile[4],digits=0),nsmall=0), # 75パーセンタイル、小数点以下0桁
"]", # ]
sep="") # 区切り記号はなし
## [1] "75 [70, 80]"
# 群分けは同様です。
# subsetを使う
quantile(subset(Dataset.0,Intervention=="Yes")$Age)
## 0% 25% 50% 75% 100%
## 53 70 75 80 96
quantile(subset(Dataset.0,Intervention=="No")$Age)
## 0% 25% 50% 75% 100%
## 57 70 75 80 107
# bracketを使う
quantile(Dataset.0[Dataset.0$Intervention=="Yes",]$Age)
## 0% 25% 50% 75% 100%
## 53 70 75 80 96
quantile(Dataset.0[Dataset.0$Intervention=="No",]$Age)
## 0% 25% 50% 75% 100%
## 57 70 75 80 107
# byを使う
by(Dataset.0$Age,Dataset.0$Intervention=="Yes",quantile)
## Dataset.0$Intervention == "Yes": FALSE
## 0% 25% 50% 75% 100%
## 57 70 75 80 107
## --------------------------------------------------------
## Dataset.0$Intervention == "Yes": TRUE
## 0% 25% 50% 75% 100%
## 53 70 75 80 96
# 検定は2つに分割したベクターをt.test関数に投入します。
t.test(subset(Dataset.0,Intervention=="Yes")$Age,subset(Dataset.0,Intervention=="No")$Age)
##
## Welch Two Sample t-test
##
## data: subset(Dataset.0, Intervention == "Yes")$Age and subset(Dataset.0, Intervention == "No")$Age
## 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:
## -0.8938104 1.1030226
## sample estimates:
## mean of x mean of y
## 75.07226 74.96765
t.test(Dataset.0[Dataset.0$Intervention=="Yes",]$Age,Dataset.0[Dataset.0$Intervention=="No",]$Age)
##
## Welch Two Sample t-test
##
## data: Dataset.0[Dataset.0$Intervention == "Yes", ]$Age and Dataset.0[Dataset.0$Intervention == "No", ]$Age
## 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:
## -0.8938104 1.1030226
## sample estimates:
## mean of x mean of y
## 75.07226 74.96765
t.test(Dataset.0$Age~Dataset.0$Intervention) #こういう表記は便利ですね。
##
## Welch Two Sample t-test
##
## data: Dataset.0$Age by Dataset.0$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
# U-testはwilcox.test関数をつかいます。
wilcox.test(subset(Dataset.0,Intervention=="Yes")$Age,subset(Dataset.0,Intervention=="No")$Age)
##
## Wilcoxon rank sum test with continuity correction
##
## data: subset(Dataset.0, Intervention == "Yes")$Age and subset(Dataset.0, Intervention == "No")$Age
## W = 80301, p-value = 0.8248
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(Dataset.0[Dataset.0$Intervention=="Yes",]$Age,Dataset.0[Dataset.0$Intervention=="No",]$Age)
##
## Wilcoxon rank sum test with continuity correction
##
## data: Dataset.0[Dataset.0$Intervention == "Yes", ]$Age and Dataset.0[Dataset.0$Intervention == "No", ]$Age
## W = 80301, p-value = 0.8248
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(Dataset.0$Age~Dataset.0$Intervention) #こういう表記は便利ですね。
##
## Wilcoxon rank sum test with continuity correction
##
## data: Dataset.0$Age by Dataset.0$Intervention
## W = 78858, p-value = 0.8248
## alternative hypothesis: true location shift is not equal to 0
# 前回作ったTableに平均値をくっつけてみましょう
Table1<-c("Variables","Control","Active","P-value") # Table1の1行目を作ります。
tab<-xtabs(~Sex+Intervention,data=Dataset.0) # 2行目(Sex)のクロス集計表の作成
row.2<- # 2行目の文字列を作る（前回参照）。
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))
tab<- # 2行目(Sex)のクロス集計表の作成
matrix(
c(mean(Dataset.0[Dataset.0$Intervention=="Yes",]$Age), # 必要な平均値と標準偏差を配列に並べます。
mean(Dataset.0[Dataset.0$Intervention=="No",]$Age), # もう少しスマートな方法もあるのですが、、、
sd(Dataset.0[Dataset.0$Intervention=="Yes",]$Age), # いずれ解説します
sd(Dataset.0[Dataset.0$Intervention=="No",]$Age)),
ncol=2)
row.3<- # 2行目(Sex)のクロス集計表の作成
c("Age, year old", # 今回はさらに難しいですね
paste(format(round(tab[1,1],1),nsmall=1)," +/- ",format(round(tab[1,2],1),nsmall=1),sep=""),
paste(format(round(tab[2,1],1),nsmall=1)," +/- ",format(round(tab[2,2],1),nsmall=1),sep=""),
format(round(
t.test(Dataset.0[Dataset.0$Intervention=="Yes",]$Age,
Dataset.0[Dataset.0$Intervention=="No",]$Age)$p.value,3),nsmall=3))
Table1<-
rbind(
Table1,
row.2,
row.3
)
rownames(Table1)<-NULL
Table1
## [,1] [,2] [,3] [,4]
## [1,] "Variables" "Control" "Active" "P-value"
## [2,] "Female sex" "256 (69)" "290 (68)" "0.727"
## [3,] "Age, year old" "75.1 +/- 7.0" "75.0 +/- 7.3" "0.837"
# 演習
# Table1を、性別を"Male sex"に年齢を中央値と四分位での表記にそれぞれ変更して
# 作ってみてください。