統計解析ソフトRを用いたデータ解析
DAY.1

下川敏雄
和歌山県立医科大学 医学部／附属病院臨床研究センター

シェーマ
日程

内容

1日目(1)

データの要約の方法(1)：連続尺度の場合

1日目(2)

データの要約の方法(2)：離散尺度(カテゴリカル・データ)の場合

1日目(3)

プログラミング(その1)：データの要約の方法

2日目(1)

統計的推測(1)：連続尺度の場合

2日目(2)

統計的推測(2)：離散尺度(カテゴリカル・データ)の場合

2日目(3)

プログラミング(その2)：統計的推測

DAY.1： 記述統計学とグラフの書き方
Section.1： 連続尺度の要約とグラフ表示

CSVファイルの読み込み
ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata1.csvというCSVファイルを読み込む．これは，2019年度の福岡県の市区町村別での交通事故発生件数のデータである．
Input

>

dat <- read.csv("C:/Fukuoka_Seminor/data1.csv",fileEncoding = "cp932")>

head(dat)◼ 関数read.csv()はcsvファイルを読み込むために用いられる．また，このファイルでは市区町村名が日本語で書かれている．
◼ 「<-」は代入を意味しているが，「=」でもかまわない．
◼

において日本語を読み込ませるためには，Shift-JS形式での読み込みを引数「fileEncoding = “cp932”」で指定しなければならない．
◼ head()は，datという変数名で入力されたデータの最初の部分を表示させるための関数である．

◼ なお，データは，データフレーム型という変数形式で代入される(データフレーム型とはExcelのSheetのようなもの)．
Output

City
1
うきは市
2
みやま市
3 鞍手郡鞍手町
4 鞍手郡小竹町
5 遠賀郡芦屋町
6 遠賀郡遠賀町

Accident
7
8
5
3
3
5

◼ ここで，City，Accidentは変数名，左端の1から6までの数字は個体のID番号(行名)を表す．

データフレームの操作
Input
Output>dat$City[1] "うきは市"[7] "遠賀郡岡垣町"[13] "京都郡みやこ町"(以下省略)"みやま市""遠賀郡水巻町""京都郡苅田町""鞍手郡鞍手町""嘉穂郡桂川町""古賀市""鞍手郡小竹町"
"嘉麻市"
"行橋市"

"遠賀郡芦屋町""久留米市""三井郡大刀洗町""遠賀郡遠賀町""宮若市""三潴郡大木町"

Inputで指定されているのは，データフレームdatのなかの変数Cityの中身を表示させることを表している．
なお，dat[,1]と入力しても同じ意味である(データフレームdatの1列目を指定している)．

Input
Output>dat[12,]City Acceident
12 宮若市
17

Inputで指定されているのは，データフレームdatのなかの12番目の個体を指定している．
少し高度な利用方法：交通事故件数が50件以上の市区町村を表示する．

Input
Output>dat[dat$Accident >= 50, ]City Accident
11
久留米市
115
53
飯塚市
60
55
福岡市西区
93
56
福岡市早良区
5457
福岡市中央区98(以下省略)

データフレームの操作
>

Input

colnames(dat)[1] "City"

Output>Input"Accident"

rownames(dat)[1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"[15] "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28"[29] "29" "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42"[43] "43" "44" "45" "46" "47" "48" "49" "50" "51" "52" "53" "54" "55" "56"[57] "57" "58" "59" "60" "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"Output関数colnames()は変数名を抽出するのに用いられ，関数rownames()は個体名を抽出するのに用いられる．
Input

OutputInput

Output>ncol(dat)

[1] 2

>

nrow(dat)[1] 70関数ncol()は変数の数(列数)を抽出するのに用いられ，関数nrow()は個体数(行数)を抽出するのに用いられる．

要約統計量の計算
要約統計量を計算する前に，データフレームdataのなかの変数Accidentの内容をAccに代入する．
Input

Output>Acc <- dat$Accident

>

Acc

[1]
7
[23] 16[45]
6[67] 127

8
5
6
26

5
24
3
35

3
10
7
30

3
17
33

5
12
21

12
44
8

14
38
33

4
14
60

9 115
45 49
32 93

17
1
54

7 29
3
7
98 11221 26
22 49
89 149

4
16
28

6
13
10

41
2
27

37 44 2118 32
6
27 107 117

このように，1個の変数のみで構成されているデータ形式をベクトル型という．また，数値のベクトル型をnumeric(数値)型という．

なお，ここでは，ベクトル型で要約統計量を計算するが，データフレーム型のままでも計算できる．
先ほどのデータフレームとの違いを以下に示す．
データフレームdat

ベクトルAcc

City

Accident

1

2

3

1

うきは市

7

7

8

5

2

みやま市

8

3

鞍手郡鞍手町5

4

鞍手郡小竹町

3

4

・・・
3 ・・・

データフレームは，複数の変数で構成されるため，変数名が存在する．一方で，ベクトルは1変数しか扱えないため，変数名
は存在しない．

・・・

・・・

・・・

> dat[3,2]> Acc[3]それぞれ，5が表示される．

平均値に基づく要約の方法
平均値の計算Input>mean(Acc)

Output [1] 31.7

不偏分散の計算
Input>

var(Acc)

Output [1] 1204.445
不偏標準偏差の計算
Input>

sd(Acc)Output [1] 34.70511

少し高度な利用方法：小数点以下3位を四捨五入して不偏標準偏差を求める

Input

>

round(sd(Acc), 2)

Output [1] 34.71

中央値の計算
Input>

median(Acc)

Output [1] 21
四分位範囲の計算
Input > Q1 <- quantile(Acc, 0.25)

Output>Q3 <- quantile(Acc, 0.75)

>

Q3-Q1

75%
30.75

関数quantile()の引数0.25は第1四分位点(昇順に並べ替えたときの25%に位置する値)，0.75は第3四分位
点(昇順に並べ替えたときの75%に位置する値)を表す．ちなみに，0.50にすると中央値を表す．

少し高度な利用方法：中央値，第1四分位点，第3四分位点を一度に表示させる
Input
Output>quantile(Acc,c(0.5,0.25,0.75))

50%
21.00

25%
75%
7.00 37.75

余談：欠測がある場合の注意点
要約統計量の関数は，いずれも欠測がある場合には，欠測値になってしまう．そのため，オプションとしてna.rm=TRUEを指定して，欠測は除外する必要がある．
ここでは，変数x(ベクトル型)を作成してその内容を理解する (ちなみに，NAは欠測の意味)．
Input
OutputInput>

x <- c(10, 32, NA, 43, 23, 12, 53, 34)>

mean(x)[1] NA>mean(x, na.rm=TRUE)

Output [1] 29.57143
Excelでデータを創る場合のヒント
エクセルでデータを作ったうえで，CSVファイルで保存したものをで読み込ませる場合，欠測は空白にするのではなく，その部分にNAと入力して作る必要がある．空白は，データが記述されていないものと判断され，エラーメッセージが返されるので注意が必要．

ヒストグラム
データの分布をグラフで表す方法がヒストグラムである．

では，hist()という関数を用いる．

hist()
hist(x,breaks,freq,col,border,main,xlab,ylab)：データ(ベクトル形式)：TRUE(密度での記載)，FALSE(度数での記載)：ヒストグラムの線の色の指定：X軸のラベルヒストグラムの棒の数の指定
ヒストグラムを塗りつぶす色の指定
グラフのタイトル
Y軸のラベルHistogram of Acc

0.020
Density

8
6
4

0

0

0.000

2

Frequency

10 12

10 15 20 25 30 35

Histogram of Acc

5Frequency

Histogram of Acc

breaks
col
main
ylab

0.010

x
freq
borderxlab

その他にもオプションはあるが割愛

0

Input
> hist(Acc)

50

100

Acc

150

0

Input

50

100

Acc

> hist(Acc,breaks=30)

150

0

50

100

150

Acc

> hist(Acc,breaks=30, freq=FALSE)

では色の指定方法には４種類ある(1) 色の名前で指定する場合(色の名前は付録参照)(2) 番号で指定する場合(色の名前は付録参照)

Histogram of Acc

0.020
0.000

0.010

Density

0.010
0.000

Density

0.020

Histogram of Acc

0

50

Input

100

150

0

50

100

Input

Acc

> hist(Acc,freq=FALSE,col="yellow",border="red")

150

Acc

> hist(Acc,freq=FALSE,col=3,border=2)
1：黒，2：赤，3：緑，4：青，5：シアン，6：マゼンダ，7：黄色(オレンジ)，8：グレイ
Histogram of Acc(4) 16進数で指定する場合(色の名前は付録参照)

Histogram of Acc

0.000

0.010

Density

0.010
0.000

Density

0.020

0.020

(3) RGBで指定する場合

0

50

Input

100

150

Acc

> hist(Acc,freq=FALSE,col=rgb(1,1,0.3),+ border=rgb(0,0.5,0.2))

0

Input

50

100

150

Acc

> hist(Acc,freq=FALSE,col=3,border=2)

参考：ヒストグラムのパターン

名前：対称分布．
形状：平均値，最頻値および中央値が重なりデータの存在範囲の中心である．平均値を中心に左右対称である．
説明：一般に現れる形．

名前：歯抜け分布．
形状：級の一つおきに度数が少なくなっている．
説明：標本サイズの割に級の数が多い，境界値のとり方が適当でないときに現れる．

名前：左にひずんでいる形，右に裾が長い(L形状) ．
形状：平均値．最頻値および中央値が一致しない場合である．すなわち平均値がデータの存在範囲より左にあり，左右対称ではない．
説明：理論的に下限がおさえられており，下限値以下のデータが存在しない場合に現れる．

名前：右にひずんでいる形，左に裾が長い(Ｊ形状) ．
形状：平均値．最頻値および中央値が一致しない場合である．すなわち平均値がデータの存在範囲より右にあり，左右対称ではない．
説明：理論値または規格値などによって，上限値以上のデータが存在しない場合に現れる．

参考：ヒストグラムのパターン

名前：台形形(高原形) ．
形状：各級に含まれる度数が両端を除いてほぼ等しい．
説明：平均値が異なるデータが混在するときに現れる (統計的評価は困難)．

名前：離れ島形.
形状：対称形のヒストグラムのほかに，それと離れたところに度数がある．
説明：データに外れ値が混入しているときに現れる．

名前：二峰性形．
形状：データの存在する範囲の中心に度数が少なく，両端に山が２つある．
説明：平均値が極端に異なるデータが混在しているときに現れる．

ヒストグラムと位置を表す測度の関係
左に歪んでいる分布(非対称分布)

平均値
中央値

対称分布

右に歪んでいる分布(非対称分布)

平均値
中央値
最頻値

平均値
中央値

最頻値

代表値：中央値(最頻値)

最も一般的な状況
代表値：平均値

最頻値

代表値：中央値(最頻値)

ヒストグラムを完成させる> hist(Acc,col=4,border=3,main="福岡県における交通事故件数", xlab="件数", ylab="度数")10 15 20 25 30 35
0

5

度数

福岡県における交通事故件数

0

50

100

150

件数

左に歪んだ分布形状を示している．すなわち，中央値で要約するほうが良い．
したがって，適切な要約の方法は，中央値 [第1四分位点，第3四分位点]で要約して21.0 [7.0, 37.37]と表すほうことが推奨される．

複数のデータでの要約
ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata2.csvというCSVファイルを読み込む．これは，2017～2020年度の福岡県の市区町村別での交通事故発生件数のデータである．ここでは，Cityという変数に市区町村名が入っているので，それを列名にしたうえで計算する．
Input

Output>dat <- read.csv("C:/Fukuoka_Seminor/data2.csv",fileEncoding = "cp932")>

rownames(dat) <- dat$City>

dat <- dat[,-1]>head(dat)

Y2017 Y2018 Y2019 Y2020
うきは市
9
9
7
3
みやま市
19
18
8
13
鞍手郡鞍手町410
54
鞍手郡小竹町16
30
遠賀郡芦屋町5633
遠賀郡遠賀町1812
56

2017年度～2020年度における年度別での平均値の計算
Input

Output>Ms <- colMeans(dat)

>

Ms

Y2017
Y2018
Y2019
Y2020
41.78873 37.78873 31.61972 27.95775

Inputの1行目は，Msという変数(ベクトル型)に2017～2020年度の交通事故件数の平均値を代入していることを意味する．

2017年度～2020年度における年度別での中央値の計算

ここでは，2017～2020年度における年度別・市区町村別の中央値を計算するとともに，年度別および市区町村別で棒グ
ラフを描く
Input

>

Meds <- apply(dat,2,"median")

>

Meds

Output Y2017 Y2018 Y2019 Y2020
29

25

21

17

関数apply()とは，同じ処理(ここではmedian)を指定するデータフレームに一括して実施するための関数である．ここで，二つ目の引数の「2」は列毎に中央値を求めることを意味する．
25

結果を棒グラフで表す．
>

ここで引数namesは棒グラフのそれぞれの名前をベクトルで指定している．

20
15, col="greenyellow", xlab="年度", ylab="交通事故発生件数中央値")

10+

交通事故発生件数中央値

barplot(Meds, names=c("2017年度","2018年度","2019年度","2020年度")
Input

0

5

交通事故発生件数が減少傾向を示していることがわかる．

2017年度

2018年度

2019年度
年度

2020年度InputOutput>Mean.City <- apply(dat,1,"mean")>

Mean.City

うきは市
7.00
遠賀郡芦屋町4.25(以下省略)

みやま市
14.50
遠賀郡遠賀町
10.25

鞍手郡鞍手町
5.75
遠賀郡岡垣町
10.75

鞍手郡小竹町
2.50
遠賀郡水巻町
15.50

棒グラフで表す． Input
150> par(mar = c(6.5, 3, 1, 1))> barplot(Mean.City, names=names(Mean.City), col="pink", xaxt="n", space=0)> N <- length(Mean.City)

50

100

> axis(side=1, at=1:N-0.5, labels=names(Mean.City), tick=FALSE, las=2, cex.axis=0.7, mgp=c(1,0,1))• 関数par()は，プロットの環境を設定するための関数である． marは，グラフィック・ウィンドウの空白の表しており下側，左側，
上側，右側の空白である．デフォルトは3である．つまり，下側の空白を広くしている(地名なので空白が広く必要なため)．
• 関数barplot()の引数xaxt="n"はX軸の座標を表示しないことを表している．space=0は棒グラフの間の空白を0にするためである．
• 関数axis()は座標を描写するための関数である．side=1はX軸(下側)である．ちなみに，2は左側，3は上側，4は右側である．
atは，座標軸の目盛りを描写する場所である．labelsは目盛りの文字を表している．tickは目盛りに線を引くかどうか(FALSEなので引かない)，lasは目盛り文字の書き方(1は横書き，2は縦書き)，cex.axisは座標文字の大きさ(デフォル
トは1)，mgpは目盛りの位置(今回は上に上げている)．
dev.off()

福岡市博多区
久留米市
北九州市八幡西区
北九州市小倉北区
福岡市東区
福岡市中央区
福岡市南区
北九州市小倉南区
福岡市西区
飯塚市
福岡市早良区
春日市
大牟田市
筑紫野市
北九州市門司区
大野城市
太宰府市
福岡市城南区
糟屋郡粕屋町
糸島市
宗像市
直方市
北九州市八幡東区
北九州市若松区
京都郡苅田町
行橋市
中間市
柳川市
八女市
筑後市
田川市
北九州市戸畑区
糟屋郡志免町
福津市
古賀市
小郡市
糟屋郡新宮町
那珂川市
大川市
朝倉市
遠賀郡水巻町
朝倉郡筑前町
宮若市
田川郡香春町
みやま市
糟屋郡宇美町
糟屋郡須恵町
嘉麻市
豊前市
遠賀郡岡垣町
糟屋郡篠栗町
遠賀郡遠賀町
八女郡広川町
京都郡みやこ町
築上郡築上町
糟屋郡久山町
田川郡糸田町
うきは市
田川郡川崎町
三井郡大刀洗町
鞍手郡鞍手町
田川郡福智町
三潴郡大木町
朝倉郡東峰村
遠賀郡芦屋町
田川郡添田町
嘉穂郡桂川町
築上郡吉富町
鞍手郡小竹町
築上郡上毛町
田川郡大任町0>

関数を作る
先ほどの棒グラフは設定がたくさんあるため，複数のグラフを作る場合には，関数にしておくと便利．

「ファイル」→「新しいスクリプト」を選択 (左側のようなウィンドウ「Rエディッタ」が開く)「Rエディッタ」にプログラムを記述する．

Rエディッターの内容
city.barplot <- function(dat){par(mar = c(6.5, 3, 1, 1))
barplot(dat, col="pink", xaxt="n", space=0)
N <- length(dat)axis(side=1, at=1:N-0.5, labels=names(dat), tick=FALSE, las=2, cex.axis=0.7, mgp=c(1,0,1))}

入力したプログラム全体をドラックして，Ctrl+Rキーを押す(関数が
Input

>

に読み込まれる)．

city.barplot(Mean.City)

先ほどと同じ棒グラフが表示される．
「ファイル」→「別名で保存」でRエディッタの内容が保存できる (このように関数を作ることで，業務が次第に効率化される)．

散布図を作る
ここでは，データフレームdatに含まれている2017年度(Y2017)と2020年度(Y2020)のあいだで散布図を作る
Input>

plot(dat$Y2017, dat$Y2020, xlab="2017年度", ylab="2020年度", pch=16, col="red", xlim=c(0,200), ylim=c(0,200))>

abline(a=0, b=1, col="blue", lty=2)200関数plotの書式 plot(x,y,options)xlab：X軸の名前

ylab：Y軸の名前

100

xlim：X軸の最小，最大(xlim=c(0,200)は0～200まで)
ylim：Y軸の最小，最大(ylim=c(0,200)は0～200まで)

50

pch：データ点の形状 (付録2を参照)
col：データ点の色 (付録1を参照)

0

2020年度150

関数plotのオプション

0

50

100
2017年度

150

200• 関数ablineは，低水準グラフィクスと呼ばれ，プロット(高水準グラフィクス)に切片a，傾きbで直線を引くために用いられる．
• 関数ablineの引数ltyは直線の形状

相関係数の計算において，相関係数を計算する関数はcor()である．

cor()
使い方1： cor(x, y, method)

変数xと変数yの相関係数を求める．

使い方2： cor(M, method)

行列のすべての組み合わせの相関係数を求める．

methodには，"pearson"(デフォルト)，"spearman", "kendall"の3種類がある．
• pearsonはいわゆる相関係数(Pearsonの積率相関係数)である．正規分布に従うことが想定される(分布が歪んでいない)．
• spearmanはSpearmanの順位相関係数，kendallはKendallの順位相関係数であり，ノンパラメトリック相関と呼ばれる．
データが正規分布に従っていない(歪んでいる)場合に用いる．
いいかえれば，平均値で要約する場合はpearson，中央値で要約する場合にはspearmanかkendallを用いる．
2017年度と2020年度のSpearmanの順位相関係数

2017年度と2020年度のPearsonの積率相関係数

Input> cor(dat$Y2017,dat$Y2020)

Output [1] 0.9594974

Input> cor(dat$Y2017,dat$Y2020,method="spearman")

Output [1] 0.8699311

2017年度と2020年度のKendallの順位相関係数

Input> cor(dat$Y2017,dat$Y2020,method="kendall")

Output [1] 0.7262708

Note：欠測がある場合欠測がある場合(NAがある場合)には，オプションとして，
na.rm=TRUEを入れておかないと，NAで返されるので注意．

2017年度から2020年度までのPearsonの積率相関係数

Input
Output> cor(dat[,-1])
Y2017
Y2018
Y2019
Y2020
Y2017 1.0000000 0.9641707 0.9505330 0.9594974
Y2018 0.9641707 1.0000000 0.9756841 0.9764935
Y2019 0.9505330 0.9756841 1.0000000 0.9812387
Y2020 0.9594974 0.9764935 0.9812387 1.0000000

2017年度から2020年度までのSpearmanの順位相関係数

Input
Output> cor(dat[,-1], method="spearman")
Y2017
Y2018
Y2019
Y2020
Y2017 1.0000000 0.8458279 0.8730797 0.8699311
Y2018 0.8458279 1.0000000 0.9495602 0.9529501
Y2019 0.8730797 0.9495602 1.0000000 0.9766983
Y2020 0.8699311 0.9529501 0.9766983 1.0000000

2017年度から2020年度までのKendallの順位相関係数

Input
Output> cor(dat[,-1], method="kendall")
Y2017
Y2018
Y2019
Y2020
Y2017 1.0000000 0.6977526 0.7295368 0.7262708
Y2018 0.6977526 1.0000000 0.8340509 0.8389619
Y2019 0.7295368 0.8340509 1.0000000 0.8884804
Y2020 0.7262708 0.8389619 0.8884804 1.0000000

少し脱線：for文の利用方法
Rエディッタ1 for (i in 1:5){
2

print(i)3 }Output

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

これは，1から5まで1個づつ変数iの値を増やしながら繰り返す繰り返し文といわれるものである．また，print(i)は変数iの内容を表示している．
Rエディッタ1 for (i in 1:5){
2

print(Fukuoka[i,])3 }Output

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

これは，1から5まで1個づつ変数iの値を増やしながら，データフレームFukuokaのi行目を表示させている．

折れ線グラフを作る
ここでは，データフレームdatに含まれている福岡市の7区を抽出したうえで，折れ線グラフを用いて2017年から2020年案
での交通事故発生件数の変化を描写する．

Rエディッターの内容 (CTRL+Rで実行できる)
1 idx <- grep("福岡市", rownames(dat))

2 Fukuoka <- dat[idx,]
3 N <- nrow(Fukuoka)
4 plot(1:4, Fukuoka[1,], type="n", ylim=c(1,200), xaxt="n",xlab="年度", ylab="交通事故発生件数")

5 axis(side=1, at=1:4, labels=c("2017年度","2018年度","2019年度","2020年度"))
6 for (i in 1:N){
7

points(1:4,Fukuoka[i,], pch=i+15, col=i+1)

8lines(1:4,Fukuoka[i,], lty=i, col=i+1)

9 }10 pchs <- 1:N + 15

11 cols <- 1:N + 1
12 ltys <- 1:N
13 legend("topright", legend = rownames(Fukuoka), col = cols, pch = pchs, lty = ltys, cex=0.7)

次ページから，各行の説明を行う

1 idx <- grep("福岡市", rownames(dat))

関数grep(a,b)のなかのaは，検索する文字列，bは検索されるリスト(ベクトル型)を表している．また，関数rownames()は，データ
フレームdatの列名を表す．なお，idxには，列名のなかで「福岡市」を含んでいる番号が返される．
2 Fukuoka <- dat[idx,]

関数grep()で検索した結果，「福岡市」を含むリストの番号idxの列を新しいデータフレームFukuokaとして代入
3 N <- nrow(Fukuoka)

データフレームFukuokaの列数(福岡市の区の交通事故発生件数)をNに代入．
4 plot(1:4, Fukuoka[1,], type="n", ylim=c(1,200), xaxt="n",xlab="年度", ylab="交通事故発生件数")

とりあえず，高水準グラフィクスのplotでアウトラインを作成している．ここでtype="n"とは，プロットを描写しないことを表している．
5 axis(side=1, at=1:4, labels=c("2017年度","2018年度","2019年度","2020年度"))

X軸の描写である．
6 for (i in 1:N){
7

points(1:4,Fukuoka[i,], pch=i+15, col=i+1)

8lines(1:4,Fukuoka[i,], lty=i, col=i+1)

9 }ここでは，for文というものを用いている(次ページで例を挙げて説明)．これは，1からNまで1個づつ値を増やしながらカッコ内を繰り返す繰り返し文といわれるものである．また，関数points，linesは低水準グラフィクスであり，それぞれ，点および線を引くための関数である(pch, lty, colは関数plotと同じ定義である)．

10 pchs <- 1:N + 15

11 cols <- 1:N + 1
12 ltys <- 1:N
13 legend("topright", legend = rownames(Fukuoka), col = cols, pch = pchs, lty = ltys, cex=0.7)

10行目の1:Nは，1からNまで1づつ上昇するベクトルを作成している．それに15を足し合わせている．N=7なので，pchのなかは下記のとおり．
[1] 16 17 18 19 20 21 22

11行目の1:Nは，1からNまで1づつ上昇するベクトルを作成している．それに1を足し合わせている．N=7なので，colsのなかは下記のとおり．
[1] 2 3 4 5 6 7 8200

関数legend()は，凡例を表示させるための低水準グラフィクスである．最初の"topleft"は右上に表示させる関数．legendは凡例の名称，その他のオプションはplotと同じ．

100
50
0

交通事故発生件数150

福岡市城南区
福岡市西区
福岡市早良区
福岡市中央区
福岡市東区
福岡市南区
福岡市博多区2017年度

2018年度

2019年度
年度

2020年度殆どの区で減少傾向を示している．とくに，南区の減少は著しい．

箱ひげ図を作る
ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata3.csvというCSVファイルを読み込む．これは，2020年度の福岡県の市区町村別での交通事故発生件数と2015年度の人口のデータである．ここでは，Cityという変数に
市区町村名が入っているので，それを列名にしたうえで計算する．
Input

Output>dat <- read.csv("C:/Fukuoka_Seminor/data3.csv",fileEncoding = "cp932")>

rownames(dat) <- dat$City>

dat <- dat[,-1]>head(dat)

Y2020 Population
うきは市
3
29509
みやま市
13
38139
鞍手郡鞍手町416007
鞍手郡小竹町
07810
遠賀郡芦屋町314208
遠賀郡遠賀町6
18877

人口10,000人当たりでの交通事故件数を計算する．
Input
Output>

Ac10000 <- dat$Y2020/dat$Population*10000>

Ac10000[1] 1.016639 3.408584 2.498907 0.000000 2.111486 3.178471 2.216593 3.793496[9] 1.481921 2.323000 4.104389 2.845760 3.457985 6.292366 3.623251 3.825121[17] 2.642357 3.527088 3.005960 3.419122 5.244428 3.966680 3.691302 7.294833(以下省略)

人口の中央値を求めて，中央値以上を1，中央値未満を0とする新しい変数groupを作る．
Input

Output>group <- numeric(nrow(dat))>group[dat$Population >= median(dat$Population)] <- 1

>

group

[1] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 1 1 1 1 0 0 1 0 0 0 1 1 0 1 1 0 0 0[37] 1 1 0 0 0 1 1 0 0 0 0 0 0 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1

箱ひげ図を描く

boxplot()hist(formula,col,border,horizontal,main,xlab,ylab)
formulacolhorizontal
xlab

その他にもオプションはあるが割愛：データの形式を記載する．例えば，x~gorupとすれば，group毎にデータxのボックスプロットを描く
border ：箱ひげ図の輪郭を塗りつぶす色の指定：箱ひげ図を塗りつぶす色の指定
main：TRUEにすると水平の箱ひげ図が作成されるグラフのタイトル
ylab
Y軸のラベル：X軸のラベル箱ひげ図の概要

髭以上

髭以下

Q1 − 1.5  IQR

以上の最小値

中央値( Q2)

第3四分位点( Q3)

第1四分位点( Q1)■ 中央値，第1四分位点，第3四分位点に基づいて構成されるグラフ．
■ 外れ値，分布形状，2標本の比較などに用いることができる．

■ 単峰性が仮定されていることに注意

Q3 + 1.5  IQR
以下の最大値

8boxplot(Ac10000~group, names=c("<Median", ">Median"),+col="moccasin", xlab="群", ylab="人口10000人当たりの事故数")

4

>

0

2

人口10000人当たりの事故数

6

箱ひげ図(ボックスプロット)を描く<Median>Median

8

群

ノッチ付き箱ひげ図(ノッチド・ボックスプロット)を描く
boxplot(Ac10000~group, names=c("<Median", ">Median"),notch=TRUE,+

col="moccasin", xlab="群", ylab="人口10000人当たりの事故数")

2

4

ノッチとは，中央値(横線の部分)の信頼度を表しており，ノッチが重なっていない場合には中央値が異なる可能性が高いと判断される．

0

人口10000人当たりの事故数

6><Median>Median群人口10,000万にあたりの交通事故件数ならびに，人口の中央値で分けたデータをデータフレームcompdatに保存するとともに， CドライブのFukuoka_Seminorにファイル名compdata.csvで保存する．

Input

>

compdat <- data.frame(accident=Ac10000, group=group)>

write.csv(compdat, "C:/Fukuoka_Seminor/compdata.csv",row.names=FALSE)

・・・

・・・

DAY.1： 記述統計学とグラフの書き方
Section.2： 離散尺度(カテゴリカルデータ)の要約とグラフ表示

度数の計算ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata4.csvというCSVファイルを読み込む．これは，福岡県内で発生した交通事故に関して，年齢層(高齢者/非高齢者)，時間帯(日中／夜)，天候に関するデータである．
Input
Output>dat <- read.csv("C:/Fukuoka_Seminor/data4.csv",fileEncoding = "cp932")>

head(dat)

Senior Day_Night Weather
1
高齢者
夜
曇
2 非高齢者夜晴
3 非高齢者夜晴
4 非高齢者夜晴
5 非高齢者
昼
晴
6
高齢者
昼
晴

天候別での度数を集計する(このとき，度数が高い順に並べ替える)．
Input

Output>Weather.tbl <- sort(table(dat$Weather),decreasing=TRUE)>

head(Weather.tbl)

晴
1164

曇
雨
738 179

小雨
137

円グラフの作成
円グラフの描写について，

ではpie()という関数を用いる．

pie()
pie(x,labels,radius,clockwise,init.angle,density,col,angle)x：データ(ベクトル形式)
radius：円の大きさを(-1から1，デフォルトは0.8)
init.angle ：始点を確度で指定(デフォルトは90)
col：塗りつぶしの色(ベクトル形式)

その他にもオプションはあるが割愛
labels：各扇状のラベル(文字列ベクトル)
clickwise ：時計回り(TRUE)，反時計回り(FALSE)
deisity：塗分ける際の線の本数
angle：塗分ける際の線の確度
晴

晴

小雨

小雨

雨

雨

Input

曇

Input

曇

>

cols <- 1:length(Weather.tbl)+1>

dens <- (1:length(Weather.tbl)-1)*10>

pie(Weather.tbl, col=cols)>pie(Weather.tbl, density=dens, angle = dens)

関数sprintf()を用いて度数分布表のための関数を作成する．
ここでは，Freq.table()という新しい関数を作成する．
Rエディッターの内容
1 Freq.table <- function(dat){2

tbl <- sort(table(dat),decreasing=TRUE)

各カテゴリの度数を計算

3

N <- sum(tbl, na.rm=T)

データの個数を計算

4pc <- tbl/N*100

各カテゴリの割合を計算

5res <- sprintf("%d (%3.1f)", tbl, pc)

関数sprint()の説明は下記で説明

6

return(res)7 ｝

関数sprintf()は，書式付きの文字列を作る関数である．括弧()内の説明を以下に示す― ""は文字列に含む書式を表している．

― %dの部分には整数が入ることを意味する．
― %3.1fの部分にはすべての桁数が3桁で，そのうち1桁が小数点以下の桁数であることを意味する(なお，マイナスも1文字になるので注意)― ""のなかの括弧()は文字として認識される．
― ""の外側の変数においてtblは%dに入る数値，pcは%3.1fに入る数値．

Input
Output>Freq.table(dat$Weather)

晴"1164 (52.5)"

曇"738 (33.3)"

雨"179 (8.1)"

小雨"137 (6.2)"

クロス集計表の作成
年齢層(Senior)と天候(Weather)のクロス集計表を作る．このとき，高齢者は「高齢者，非高齢者」の順番，天気は，「晴，曇，小雨，雨」の順番にする．
Input
OutputInput

OutputInput

OutputInput

Output>Sn <- factor(dat$Senior)>levels(Sn)[1] "高齢者""非高齢者">

We <- factor(dat$Weather)>levels(We)[1] "雨""小雨" "晴""曇">

levels(We) <- c("晴","曇", "小雨", "雨")>levels(We)[1] "晴""曇""小雨" "雨">dt <- data.frame(Sn, We)>

conj.tbl <- table(dt$Sn, dt$We)>

conj.tbl
晴 曇 小雨 雨
高齢者
46 38 361 275
非高齢者 133 99 803 463

関数sprintf()を用いてクロス集計表のための関数を作成する．
ここでは，Conject.table()という新しい関数を作成する．
Rエディッターの内容
1 Conject.table <- function(X,Y,idx=NULL){2

dat <- data.frame(X,Y)

3tbl <- table(dat$X, dat$Y)4

pct <- prop.table(tbl, idx)*100

5

N <- nrow(tbl)6

M <- ncol(tbl)7res <- matrix(numeric(N*M), ncol=M)

8

for (i in 1:N){

9

for (j in 1:M){

10

res[i,j] <- sprintf("%d (%3.2f)", tbl[i,j], pct[i,j])

11

}

12

}

13

colnames(res) <- colnames(tbl)

14

rownames(res) <- rownames(tbl)

15

return(res)16 ｝

関数prop.table(x, margin)において，xはクロス集計表，marginはパーセントを求める方向を表しており，NULL(デフォルト)は全体パーセント，1は行パーセント(各行で合計すると100%になる)，2は列パーセント(各列で合計すると100%になる)を表している．

全体パーセントを加えたクロス集計表を作る
Input
Output>Conject.table(Sn, We)

晴
曇
小雨
雨
高齢者
"46 (2.07)" "38 (1.71)" "361 (16.28)" "275 (12.40)"
非高齢者 "133 (6.00)" "99 (4.46)" "803 (36.20)" "463 (20.87)"

行パーセントを加えたクロス集計表を作る
Input
Output>Conject.table(Sn, We, 1)

晴
曇
小雨
雨
高齢者
"46 (6.39)" "38 (5.28)" "361 (50.14)" "275 (38.19)"
非高齢者 "133 (8.88)" "99 (6.61)" "803 (53.60)" "463 (30.91)"

列パーセントを加えたクロス集計表を作る
Input
Output>Conject.table(Sn, We, 2)

晴
曇
小雨
雨
高齢者
"46 (25.70)" "38 (27.74)" "361 (31.01)" "275 (37.26)"
非高齢者 "133 (74.30)" "99 (72.26)" "803 (68.99)" "463 (62.74)"

帯グラフの作成
帯グラフは，クロス集計表におけるパーセントと関数barplot()を組み合わせることで描写できる．
t(conj.tbl)
conj.tbl
晴曇小雨
雨

晴 曇 小雨 雨
高齢者
46 38 361 275
非高齢者 133 99 803 463

Input

Input

>

>

barplot(prop.table(t(conj.tbl),2), beside=F, col=c(2,3,4,5))

0.0

0.0

0.2

0.2

0.4

0.4

0.6

0.6

0.8

0.8

1.0

1.0

barplot(prop.table(conj.tbl,2), beside=F, col=c(2,3))

高齢者 非高齢者
46
133
38
99
361
803
275
463

晴

曇

小雨

雨

高齢者

非高齢者

クロス集計表のグラフ表示
クロス集計表をグラフ表示する方法としてモザイク・プロットというものがある．
できる．

では，mosaicplot()でびょうしゃ

mosaicplot()mosaicplot(x,main,xlab,ylab,col)その他にもオプションはあるが割愛：クロス集計表：X軸の名前：塗りつぶしの色(ベクトル形式)

x
xlab
col

main
ylab

：グラフのタイトル：Y軸の名前

年齢層と天候の関係
非高齢者

雨

小雨

曇

晴

高齢者

Input
>

mosaicplot(conj.tbl, col=c(2,3,4,5),main="年齢層と天候の関係")

対応のあるクロス集計表
対応のないクロス集計表(通常のクロス集計表)
いま，ゴミ処理施設を建設することを計画している．施設の周辺環境に対する環境への影響に関して，心配であるか否かを施設建設予定地付近の地域Aと，予定地から離れている地域Bでアンケート調査を行った．

対応のあるクロス集計表
いま，ゴミ処理施設を建設することを計画している．施設の周辺環境に対する環境への影響に関して，心配であるか否かを調査した後で，住民説明会を行った．
その後，不安であるか否かを改めて調査した．

説明後

不安の有無
あり

なし

地域A

(a)

(b)

地域B

(c)

(d)

地域Aの不安度 =(𝐚)𝐚 +(𝐛)

比較

地域Bの不安度 =(𝐜)𝐜 +(𝐝)

説明
前

不安
あり

不安
なし

不安あり

(a)

(b)

不安なし

(c)

(d)

期待すること：(𝐛)

不安あり → 不安なし = 𝐚 + 𝐛 + 𝐜 +(𝐝)あってはならないこと：

不安なし → 不安あり =(𝐜)𝐚 + 𝐛 + 𝐜 +(𝐝)

対応のあるクロス集計表の例
ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata2.csvというCSVファイルを読み込む．このデータを加工する．先ず，2018年度と2017年度を比較して上昇した市区町村にはup，そうでない市区町村にはdownとラベルとつけ，次いで，2019年度と2018年度を比較して，同様のラベルをつける．そのうえで，対応のあるクロス集計表をつくる．
Input

>

H2018 <- H2019 <- rep("Down", nrow(dat))>

H2018[dat$Y2018 > dat$Y2017] <- "Up">H2019[dat$Y2019 > dat$Y2018] <- "Up">

Conject.table(H2018,H2019)

関数rep(a,N)はaをN回繰り返すベクトルを作るための関数である．ここでは，Downという文字をデータフレームdatの列数だけ作ることを意味する．

OutputDown
Up
Down "32 (45.07)" "10 (14.08)"Up"28 (39.44)" "1 (1.41)"

行(縦方向)が2018年度，列(横方向)が2019年度である．両年度ともに前年度比で交通事故件数が減少している
市区町村の割合は45.07%であった．また，2018年度は増加していたが2019年度には減少に転じた市区町村は39.44%であった．一方で，減少していたのに増加した市区町村は14.08%であり，各年度とも上昇傾向を示した
市区町村は1件(1.41%)のみだった．

DAY.1： 記述統計学とグラフの書き方
Section.3： 自作したプログラムを動かしてみる

作成した関数を利用する方法
では，以前に作成した関数を使いまわすことで，業務を大幅に簡略化することができる．ここでは，3種類の関数を用意した．

平均値および標準偏差で要約：Summary.Means()
Summary.Means(dat,grp,dig)dat
dig：データフレーム：表示する桁数(デフォルトは3)grp：グループ変数(オプション)

中央値および第1・3四分位点で要約：Summary.Medians()
Summary.Medians(dat,grp,dig)dat
dig：データフレーム：表示する桁数(デフォルトは3)grp：グループ変数(オプション)

度数分布表を一括して計算：Summary.Freq()
Summary.Freq(dat)dat：データフレーム

自作した関数の読み込み自作した関数は，関数source()を用いて読み込むことができる．ここでは，CドライブのFukuoka_Seminorに関数群が入っていることを想定する．
Input

>

source("C:/Fukuoka_Seminor/Continuous_Summary.R")

Summary.Means()，Summary.Medians()が入っているソースを読み込み

Input

>

source("C:/Fukuoka_Seminor/Categorical_Summary.R")

Summary.Freq() が入っているソースを読み込み

ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata5.csvというCSVファイルを読み込む．これは，2017～2020年度の福岡県の市区町村別での交通事故発生件数のデータである．ここでは，Cityという変数に市区町村名が入っているので，それを列名にしたうえで計算する．
Input

Output>dat.Cont <- read.csv("C:/Fukuoka_Seminor/data5.csv",fileEncoding = "cp932")>

rownames(dat.Cont) <- dat.Cont [,1]>dat.Cont <- dat.Cont[,-1]>head(dat.Cont)Y2017 Y2018 Y2019 Y2020 group
うきは市
9
9
7
3
0
みやま市
19
18
8
13
0
鞍手郡鞍手町410
5
40
鞍手郡小竹町16
300
遠賀郡芦屋町56330
遠賀郡遠賀町18125
60

ファイルdata2.csvにgroupを加えたものである．groupは人口が中央値以上を1，中央値未満を0とした変数である．

dat.Cont
うきは市
みやま市
鞍手郡鞍手町
鞍手郡小竹町
遠賀郡芦屋町
遠賀郡遠賀町[,1] [,2] [,3] [,4] [,5]Y2017 Y2018
9
9
19
18
4
10
1
6
5
6
18
12

Y2019 Y2020
7
3
8
13
5
4
3
0
3
3
5
6

group
0
0
0
0
0
0

各年度の平均および標準偏差を計算する．
Input
Output>Summary.Means(dat.Cont[,1:4])[,1]
Y2017 "41.8(48.8)"
Y2018 "37.8(43.6)"
Y2019 "31.6(35.1)"
Y2020 "28(32.3)"

人口が中央値以上／未満で分けたうえで平均および標準偏差を計算する．
Input

Output>Summary.Means(dat.Cont[,1:4], dat.Cont[,5])0
1Y2017 "13.9(48.8)" "68.9(48.8)"
Y2018 "11.2(43.6)" "63.7(43.6)"
Y2019 "8.46(35.1)" "54.1(35.1)"
Y2020 "7.11(32.3)" "48.2(32.3)"

dat.Cont
うきは市
みやま市
鞍手郡鞍手町
鞍手郡小竹町
遠賀郡芦屋町
遠賀郡遠賀町[,1] [,2] [,3] [,4] [,5]Y2017 Y2018
9
9
19
18
4
10
1
6
5
6
18
12

Y2019 Y2020
7
3
8
13
5
4
3
0
3
3
5
6

group
0
0
0
0
0
0

各年度の中央値および第1・第3四分位点を計算する．
Input
Output>Summary.Medians(dat.Cont[,1:4])[,1]
Y2017 "29 [11.5, 42.5]"
Y2018 "25 [ 9.5, 40.5]"
Y2019 "21 [ 7.0, 38.0]"
Y2020 "17 [ 7.0, 32.5]"

人口が中央値以上／未満で分けたうえで中央値および第1・第3四分位点を計算する．
Input

Output>Summary.Medians(dat.Cont[,1:4], dat.Cont[,5])0
1Y2017 "12 [5.5, 17.5]" "41.5 [33.0, 92.0]"
Y2018 "10 [6.0, 14.0]" "40.5 [32.0, 92.8]"
Y2019 " 7 [4.5, 12.5]" "38.0 [28.0, 69.0]"
Y2020 " 7 [4.0, 9.5]" "32.5 [25.5, 57.2]"

平均値での結果を保存する
Input

>

result <- Summary.Means(dat.Cont[,1:4], dat.Cont[,5])>write.csv(result, "C:/Fukuoka_Seminor/Mean_Summary.csv",fileEncoding = "cp932")

綺麗にExcelのファイルとして保存できる．

各市区町村の4年間の平均および標準偏差，中央値および第1・第3四分位点を計算する．
Input

>

dat.Cont2 <- t(dat.Cont)[-5,]>head(dat.Cont2)

dat.Cont [,1] [,2] [,3] [,4] [,5]
うきは市
みやま市
鞍手郡鞍手町
鞍手郡小竹町
遠賀郡芦屋町
遠賀郡遠賀町

Input

OutputInput
Output>Y2017 Y2018
9
9
19
18
4
10
1
6
5
6
18
12

group
0
0
0
0
0
0

Summary.Means(dat.Cont2)

うきは市
みやま市
鞍手郡鞍手町
鞍手郡小竹町(以下省略)>

Y2019 Y2020
7
3
8
13
5
4
3
0
3
3
5
6

[,1]
"7(2.83)""14.5(5.07)""5.75(2.87)""2.5(2.65)"

Summary.Medians(dat.Cont2)

うきは市
みやま市
鞍手郡鞍手町
鞍手郡小竹町(以下省略)[,1]" 8.0 [ 6.00," 15.5 [ 11.75," 4.5 [ 4.00," 2.0 [ 0.75,

9.00]"
18.25]"
6.25]"
3.75]"

dat.Cont2 [,1]

[,2]

[,3]

[,4]

うきは市
9
9
7
3

みやま市
19
18
8
13

鞍手郡鞍手町410
54

鞍手郡小竹町
1
6
3
0

Y2017
Y2018
Y2019
Y2020

・・・
・・・
・・・
・・・
・・・
・・・

ここでは，CドライブのFukuoka_Seminorというフォルダにあるdata4.csvというCSVファイルを読み込む．これは，福岡県内で発生した交通事故に関して，年齢層(高齢者/非高齢者)，時間帯(日中／夜)，天候に関するデータである．
Input
Output>dat.Categ <- read.csv("C:/Fukuoka_Seminor/data4.csv",fileEncoding = "cp932")>

head(dat.Categ)

Senior Day_Night Weather
1
高齢者
夜
曇
2 非高齢者夜晴
3 非高齢者夜晴
4 非高齢者夜晴
5 非高齢者
昼
晴
6
高齢者
昼
晴

各項目の度数分布表を計算する
Input

Output>result <- Summary.Freq(dat.Categ)>

resultvname
Seniorcategresult1
高齢者 720 (32.5)
2
非高齢者 1498 (67.5)
3 Day_Night
昼 1283 (57.8)
4
夜 935 (42.2)
5
Weather雨
179 (8.1)
6小雨
137 (6.2)
7
晴 1164 (52.5)
8
曇 738 (33.3)

結果を保存する
Input>

write.csv(result, "C:/Fukuoka_Seminor/Freq_Summary.csv",fileEncoding = "cp932",row.names=FALSE)

綺麗にExcelのファイルとして保存できる．

Thank you for your kind attention

shimokaw@wakayama-med.ac.jp
toshibow2000@gmail.com

付録1：色の名前

付録1：色の名前(続き)付録1：色の名前(Hex code)付録2：pchの一覧
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

付録3：平均値と標準偏差を一括して計算するプログラムを作る
1 Summary.Means <- function(dat, grp=NULL, dig=3){
2
P <- ncol(dat)3
if (is.null(grp)==TRUE){4
Ms <- colMeans(dat, na.rm=T)
5
Sds <- apply(dat,2,sd,na.rm=T)
6
res <- matrix(numeric(P),ncol=1)
7
for (i in 1:P){
8
res[i,1] <- sprintf("%s(%s)", format(Ms[i], digit=dig), format(Sds[i], digit=dig))9}10
rownames(res) <- colnames(dat)11}12
if (is.null(grp)==FALSE){
13gname = sort(unique(grp))
14
P.g <- length(gname)15
res <- matrix(numeric(P*P.g), ncol=P.g)
16
for (j in 1:P.g){
17
dt <- dat[grp==gname[j],]18
Ms <- colMeans(dt, na.rm=T)
19
Sds <- apply(dat,2,sd,na.rm=T)20
for (i in 1:P){
21
res[i,j] <- sprintf("%s(%s)", format(Ms[i], digit=dig), format(Sds[i], digit=dig))
22
}
23
}
24
rownames(res) <- colnames(dat)25
colnames(res) <- gname
26}
27
return(res)28 }付録4：中央値と第1,第3四分位点を一括して計算するプログラムを作る
1 Summary.Medians <- function(dat, grp=NULL, dig=3){
2
P <- ncol(dat)3
if (is.null(grp)==TRUE){
4
Meds <- format(apply(dat,2,median, na.rm=T),digit=dig)
5
Q1s <- format(apply(dat,2,quantile,0.25,na.rm=T),digit=dig)
6
Q3s <- format(apply(dat,2,quantile,0.75,na.rm=T),digit=dig)
7
res <- matrix(numeric(P),ncol=1)8
for (i in 1:P){
9
res[i,1] <- sprintf("%s [%s, %s]", Meds[i], Q1s[i], Q3s[i])10}11
rownames(res) <- colnames(dat)12}13
if (is.null(grp)==FALSE){
14gname = sort(unique(grp))
15
P.g <- length(gname)16
res <- matrix(numeric(P*P.g), ncol=P.g)
17
for (j in 1:P.g){
18
dt <- dat[grp==gname[j],]19Meds <- format(apply(dt,2,median, na.rm=T),digit=dig)20
Q1s <- format(apply(dt,2,quantile,0.25,na.rm=T),digit=dig)
21
Q3s <- format(apply(dt,2,quantile,0.75,na.rm=T),digit=dig)
22
for (i in 1:P){
23
res[i,j] <- sprintf("%s [%s, %s]", Meds[i], Q1s[i], Q3s[i])
24
}
25
}
26
rownames(res) <- colnames(dat)27
colnames(res) <- gname28}
29
return(res)30 }付録5：度数分布表を一括して計算するプログラムを作る
1 Summary.Freq <- function(dat){2
P <- ncol(dat)3
Re <- vn <- list()4
for (j in 1:P){
5
tb <- table(dat[,j])6
pr <- prop.table(tb)*100
7
N <- length(tb)8
st <- numeric(N)9vn[[j]] <- rep("",N)10vn[[j]][1] <- colnames(dat)[j]11names(st) <- names(tb)12
for (i in 1:N){
13
st[i] <- sprintf("%d (%3.1f)", tb[i], pr[i])14}15Re[[j]] <- st16}
17
unRe <- unlist(Re)18result <- data.frame(vname=unlist(vn), categ=names(unRe), result=unRe)
19
rownames(result) <- 1:nrow(result)
20
return(result)21 }