医学情報処理演習：2008年度第3回課題の解答例

課題

前回も使ったp02-2008.txtを用いて，世界各国の平均寿命（変数LIFEEXP）について，購買力平価ベースの一人当たりGDPが1万米ドル以上の国（変数RPの値が"RICH"）と1万米ドル未満の国（変数RPの値が"POOR"）に分けて，別々に記述統計量（平均値，不偏標準偏差，中央値，四分位偏差）を計算しなさい。学籍番号・氏名とともに，下のフォームにRのコードと結果を貼り付けて送信すること。

解答例

1行目は前回と同じくデータの読み込みである。2行目はtapply( )関数を使って，dat$RPごとにdat$LIFEEXPの平均値を計算し，MEANという変数に保存している。変数MEANは名前付きリストという型になる。3行目と4行目は同様に不偏標準偏差SD，中央値MEDIANを計算している。5行目と6行目で四分位偏差を計算する。5行目で五数要約値をFIVENUMという変数に保存し，6行目でFIVENUMの1番目のリスト（RPの値が"POOR"な国）の4番目の要素（FIVENUM[[1]][4]）として第３四分位を取り出し，2番目の要素である（FIVENUM[[1]][2]）第１四分位との差を取って2で割ることによりFIVENUMの1番目のリストの四分位偏差，2番目のリスト（RPの値が"RICH"な国）について同様な計算によりFIVENUMの2番目のリストの四分位偏差が得られる。7行目で4つの名前付きリストをcbind( )関数でカラムごとに結合すると，結果が表の形で表示される。

dat <- read.delim("http://phi.med.gunma-u.ac.jp/medstat/p02-2008.txt")
MEAN <- tapply(dat$LIFEEXP,dat$RP,mean)
SD <- tapply(dat$LIFEEXP,dat$RP,sd)
MEDIAN <- tapply(dat$LIFEEXP,dat$RP,median)
FIVENUM <- tapply(dat$LIFEEXP,dat$RP,fivenum)
SIQR <- c((FIVENUM[[1]][4]-FIVENUM[[1]][2])/2,(FIVENUM[[2]][4]-FIVENUM[[2]][2])/2)
cbind(MEAN,SD,MEDIAN,SIQR)

（注１）SIQRを計算している6行目は，以下でも同じ。

SIQR <- c((FIVENUM$POOR[4]-FIVENUM$POOR[2])/2,(FIVENUM$RICH[4]-FIVENUM$RICH[2])/2)

（注２）第１回でちらっと触れた関数定義を使えば，四分位偏差を計算する関数を定義することができる。つまり，５行目と６行目は，以下で置き換え可能である。思想的にはこちらの方が美しい。

siqr <- function(X) { (fivenum(X)[4]-fivenum(X)[2])/2 }
SIQR <- tapply(dat$LIFEEXP,dat$RP,siqr)

結果の表は下の通り得られる。

         MEAN        SD MEDIAN SIQR
POOR 63.08364 11.534186  67.22 9.63
RICH 75.98348  7.172874  77.93 2.32

別解（１）

どうせ関数定義するなら，4つの記述統計量を計算する関数を先に定義しておく手もある。

descstat <- function(X) {
 siqr <- function(Y) { (fivenum(Y)[4]-fivenum(Y)[2])/2 }
 list(mean(X),sd(X),median(X),siqr(X))
}

dat <- read.delim("http://phi.med.gunma-u.ac.jp/medstat/p02-2008.txt")
as.table(tapply(dat$LIFEEXP,dat$RP,descstat))

結果は以下の通り。

                                    POOR 
   63.08364, 11.53419, 67.22000, 9.63000 
                                    RICH 
75.983478, 7.172874, 77.930000, 2.320000

別解（２）

以下のように[ ]を使って地道にやる手もある。ただし1行ずつ結果が出力されるので読みにくい。

dat <- read.delim("http://phi.med.gunma-u.ac.jp/medstat/p02-2008.txt")
attach(dat)
print("Descriptive statistics of LIFEEXP for POOR countries")
mean(LIFEEXP[RP=="POOR"])
median(LIFEEXP[RP=="POOR"])
sd(LIFEEXP[RP=="POOR"])
Q <- fivenum(LIFEEXP[RP=="POOR"])
(Q[4]-Q[2])/2
print("Descriptive statistics of LIFEEXP for RICH countries")
mean(LIFEEXP[RP=="RICH"])
median(LIFEEXP[RP=="RICH"])
sd(LIFEEXP[RP=="RICH"])
Q <- fivenum(LIFEEXP[RP=="RICH"])
(Q[4]-Q[2])/2
detach(dat)

結果は以下の通り（この出力は，上のコードの実行前にsink("clipboard")を実行して結果の出力先をクリップボードに変え，実行後にsink()でRコンソールに戻し，クリップボードに溜まった出力内容を貼り付けたものである）。

[1] "Descriptive statistics of LIFEEXP for POOR countries"
[1] 63.08364
[1] 67.22
[1] 11.53419
[1] 9.63
[1] "Descriptive statistics of LIFEEXP for RICH countries"
[1] 75.98348
[1] 77.93
[1] 7.172874
[1] 2.32

リンクと引用について