# ←この記号の後は次の改行までコメントになります。
# 名義変数の記述統計量をグラフィカルに表現する。
# Lesson 10-12では記述統計量を作表しました。
# 今回は作図してみましょう。
# 今まで同様にまず準備から。
# 第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
# 基本はplot関数です。
# plot関数はスカラーまたはベクターが与えられて動作します。
# Rでは関数がオブジェクトの種類によって動作を少しづつ変えることにより、
# 複数のオブジェクトに対して動作が変わりながら対応可能な関数があります。
# plot関数もそのひとつです。plotは数値ベクターと数式に対応して動作する
# だけでなく、様々な関数の出力結果の可視化を受け持っています。
plot(1) # スカラーの指定
plot(c(1,3,2,4)) # ベクターの指定
plot(c(1,3,5)~c(2,4,6)) # y~xで指定
x<-c(0:30/10) # xに代入 0.0, 0.1, 0.2, ... ,3.0
y<-x^2-x-2 # yはxの二次式
plot(y~x) # xとyの関係（放物線）を作図
plot(y~x,type="l") # 折れ線で結ぶ(ほぼ曲線に見えますが)
grid() # グリッドを追加
# plotからコマンドラインで直接ファイルに保存することもできますが、
# 最初は右下のウィンドウ → Export から保存したほうが良いと思います。
# 名義変数の統計量を可視化する
barplot(table(Dataset.0$Sex)) # シンプルな棒グラフ
# クロス集計表を作っておいて、
tab<-xtabs(~Sex+Intervention,data=Dataset.0)
tab
## Intervention
## Sex No Yes
## Male 115 139
## Female 256 290
# 群分けした性別のカウント数を表示
barplot(tab)
tab.perc<- # 割合に計算し直してみます
tab/c(c(sum(tab[,1]),sum(tab[,1]),sum(tab[,2]),sum(tab[,2])))
barplot(tab.perc) # 割合を積み上げグラフで表示
barplot(tab.perc,horiz=T) # 横向きにしてみる
# 円グラフはあまり使いませんが…
pie(tab[,1])
pie(tab[,2])
# 3水準以上だとこういうグラフが特に意味を持ちます。
tab<-xtabs(~Comorbidity+Intervention,data=Dataset.0)
tab.perc<-
cbind(tab[,1]/sum(tab[,1]),tab[,2]/sum(tab[,2]))
barplot(tab.perc,horiz=T)
# 飾ってみましょう
barplot(tab.perc,
# 横向き
horiz=T,
# グラフそれぞれのタイトル
names.arg=c("Yes","No"),
# 積み上げグラフの色の指定
col=rainbow(9),
# 指定した色の脚注
legend.text=c(0:8),
# 積み上げグラフの区切り線なし
border=NA,
# 図全体のタイトル
main="Figure. Distribution of comorbidity index in groups
dichotomized by presence of intervention"
)
# 作図後はExportボタンから、eps (encripted postscript)形式に保存し、
# Adobe Illustratorなどで加工しています。高価なアプリケーションですが、
# 研究チームでひとつ共有しておいたほうが良いものです。