2008年7月アーカイブ

アメリカ各州の犯罪件数のデータ

を利用して主成分分析を行います。まず、データを読み込みます。read.csv中の「row.names=1」はデータの1列目をデータとしてではなく、列の名前として指定するというオプションです。

> crime <- read.csv("crime.csv",row.names=1)

> crime

殺人 強姦 車の窃盗

ALABAMA 14.2 25.2 280.7

ALASKA 10.8 51.6 753.3

ARIZONA 9.5 34.2 439.5

ARKANSAS 8.8 27.6 183.4

CALIFORNIA 11.5 49.4 663.5

COLORADO 6.3 42.0 477.1

CONNECTICUT 4.2 16.8 593.2

DELAWARE 6.0 24.9 467.0

FLORIDA 10.2 39.6 351.4

GEORGIA 11.7 31.1 297.9

HAWAII 7.2 25.5 489.4

主成分分析を行うための関数はprcompです。以下のように利用します。「scale=T」は標準化後のデータについて分析を行うというオプションです。いつでも必ず利用するオプションと考えて下さい。

> result <- prcomp(crime,scale=T)

summaryを使うと、累積寄与率が表示されます。

> summary(result)

Importance of components:

PC1 PC2 PC3

Standard deviation 1.314 0.970 0.577

Proportion of Variance 0.575 0.313 0.111

Cumulative Proportion 0.575 0.889 1.000

第2主成分までの累積寄与率が88.9%であるので、第2主成分までについて分析を進めていきましょう。

結果を代入した変数名$rotationとすると、主成分負荷量が表示されます。

> result$rotation

PC1 PC2 PC3

殺人 -0.6088455 0.49760728 -0.6178140

強姦 -0.6910353 0.04978736 0.7211043

車の窃盗 -0.3895861 -0.86597241 -0.3135514

また、結果を代入した変数名$x とすると、主成分得点が表示されます。

> result$x

PC1 PC2 PC3

ALABAMA -0.83442342 1.300508984 -0.95824549

ALASKA -2.94664437 -1.131057633 0.58807447

ARIZONA -0.99230103 0.026252731 0.13840981

ARKANSAS 0.05770611 1.052383789 0.22314053

CALIFORNIA -2.73467002 -0.649053651 0.47438226

COLORADO -1.06513967 -0.517820358 1.11148111

累積寄与率から第2主成分までに全体の89％の情報が含まれていることが分かりました。ここでは、第1主成分と第2主成分をそれぞれX軸、Y軸にして、各州の主成分得点を散布図にしてみましょう。

> plot(result$x[,1],result$x[,2],pch=" ",xlab="第1主成分",ylab="第2主成分")

> text(result$x[,1],result$x[,2],rownames(result$x))

練習問題

には製薬会社の経営指標が、

には各都道府県の産業別労働従事者比率が記録され得ています。このデータを使って主成分分析を行いなさい。

ファイル

に横浜と川崎を除く神奈川県内の各市における人口とゴミ収集量のデータが保存されています。このデータを使って、

（ゴミ収集量）＝×(人口）

の単回帰分析を行ってみましょう。

データの入力

いつものように read.csvを使ってデータを読み込みます。

> g <- read.csv("garbage.csv")

> g

市 人口ごみ収集量

1 平塚 251991 75938

2 鎌倉 172638 64139

3 藤沢 362088 106898

4 小田原 197460 72415

5 茅ヶ崎 209575 64302

6 逗子 56436 20700

7 相模原 560366 169343

8 三浦 53795 20550

9 秦野 161692 46948

10 厚木 206186 85428

11 大和 202200 76155

12 伊勢原 96436 31048

13 海老名 110638 32124

14 座間 117025 36712

15 南足柄 43314 12939

16 綾瀬 80430 25665

R で回帰分析を行うための関数はlmです。説明変数と被説明変数をそれぞれ、人口とごみ収集量とする回帰分析を行うには以下のように入力します。

> result <- lm(g$ごみ収集量 ~ g$人口)

以上で変数resultに結果が代入されました。結果の要約を表示するにはsummaryという関数を利用します。

> summary(result)

Call:

lm(formula = g$ごみ収集量 ~ g$人口)

Residuals:

Min 1Q Median 3Q Max

-6495 -4545 -3229 1595 18787

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 4.811e+03 3.283e+03 1.466 0.165

g$人口 2.999e-01 1.486e-02 20.175 9.55e-12 ***

---

Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 7598 on 14 degrees of freedom

Multiple R-Squared: 0.9667, Adjusted R-squared: 0.9644

F-statistic: 407 on 1 and 14 DF, p-value: 9.552e-12

まず注目すべき点は、Coefficientsの部分です。ここで、

Estimate ：係数

Std. Error ：標準誤差

Pr(>|t|) ： P値

を意味します。また、

Residual standard error：標準誤差

Multiple R-Squared：決定係数

を表しますので、この回帰分析の結果は

となります。

p値は

帰無仮説：「回帰係数は0である」、 対立仮説：「回帰係数は0ではない」

に関する仮説検定のp値を意味しています。この結果から、5%の有意水準で、「a=0」は受容され、また、「b=0」は棄却されることが分かります。いいかえれば、切片は0であるが、人口はゴミ収集量に影響を与えていると考えられます。

残差を読み取る

次に残差をグラフ化してみましょう。残差はresidualsという関数で得られます。

> residuals(result)

1 2 3 4 5 6 7

-4439.2846 7557.7275 -6494.6658 8390.2205 -3355.7699 -1035.1521 -3508.3591

8 9 10 11 12 13 14

-393.1811 -6350.8365 18786.5078 10708.8111 -2682.1678 -5864.9981 -3192.3023

15 16

-4861.1871 -3265.3622

これをグラフにしてみましょう。

> plot(residuals(result),xlab="データNo.",ylab="残差")

次のようなグラフが表示されます。

もし、回帰分析がうまくいっているなら、残差の点は0を中心にまんべんなくばらついているはずです。ところが、このグラフを見ると、大きくプラスの方向へ外れている点が4点あります。この4点がどの都市に相当するのかを調べてみましょう。

以下のように入力すると、グラフに十字マークが表示されます。マークをプラス方向へ外れている4点に移動させ、マウスで右クリックして下さい。そのデータに対応する都市名が表示されます。表示が終わったらグラフ左上の「stop」メニューから「stop locator」を選択して下さい。

> identify(residuals(result),labels=d$市)

この結果から、これら4点が厚木・大和・小田原・鎌倉市であることが分かりました。これらの市は大都市または観光地です。このことから、ゴミの収集量を説明するには、単に人口だけではなく、都市の規模や、観光地があるかどうかなどを考慮しなければならないことが示唆されます。

練習問題

小売業に関する売上高（百億円）と従業員数（百人）のデータが

に保存されています。このデータを利用して散布図を描画し、単回帰モデルを推定してみましょう。

1. セル E2,に「=C2*D2」と入力します。
2. セル E2からE16までドラッグし、下方向へコピーします。
   

1. 入力Y範囲：B1からB16
2. 入力X範囲：D1からD16
3. 「ラベル」にチェック

1. 入力Y範囲：B1からB16
2. 入力X範囲：D1からE16
3. 「ラベル」にチェック

1. 入力Y範囲：A1からA16
2. 入力X範囲：B1からG16
3. 「ラベル」にチェック

1. セル E2, F2, G2 にそれぞれ「=IF(B2=1,1,0)」「=IF(B2=2,1,0)」「=IF(B2=3,1,0)」と入力します。
2. セル E2からG22までドラッグし、下方向へコピーします。
   

1. 入力Y範囲：C1からC22
2. 入力X範囲：D1からD22
3. 「ラベル」にチェック

1. 入力Y範囲：C1からC22
2. 入力X範囲：D1からG22
3. 「ラベル」にチェック

• 作付面積を横軸に、収穫量を縦軸にとる散布図を作成して下さい。
  

1. 入力Y範囲：B1からB17
2. 入力X範囲：C1からC17
3. 「ラベル」にチェック

1. 入力Y範囲：B1からB17
2. 入力X範囲：C1からD17
3. 「ラベル」にチェック

分散分析

ある同一の部品を製作している3つの班があります。各班の製作効率に差があるかどうかを検証するために、各班が1個の部品を製作するための時間を複数回測定した結果が

に保存されています。

> time <- read.csv("anova.csv")

> time

時間 班

1 8.82 1

2 9.26 1

3 8.70 1

4 8.97 1

5 8.64 1

6 8.29 1

7 9.45 1

8 9.42 1

9 9.25 1

10 8.21 2

11 6.65 2

「時間」の列は製作するためにかかった時間を、「班」の列はどの班が作ったかを表しています。

それでは、分散分析を行います。分散分析にはoneway.testという関数を利用します。

> oneway.test(time$時間 ~ time$班,var.equal=T)

One-way analysis of means

data: time$時間 and time$班

F = 8.0335, num df = 2, denom df = 23, p-value = 0.002260

p-valueが0.00226<0.05なので、5％の有意水準で帰無仮説を棄却します。つまり、1つの部品を製作する平均時間は3班で差があるということになります。

Wilcoxonの順位和検定

以下のようなアンケートを実施しました。

問1 あなたの国籍は何ですか （日本 ・ アメリカ）

問2 あなたは親に反抗しますか？

1. よくする 2.時々する 3.あまりしない 4.ほとんどしない

このアンケートの結果が

に保存されています。

> question1 <- read.csv("wilcoxon1.csv")

> question1

国籍 回答

1 日本 1

2 日本 4

3 日本 2

4 日本 4

5 日本 1

6 日本 4

7 日本 2

8 日本 4

9 日本 3

10 日本 4

11 アメリカ 3

12 アメリカ 1

この結果から、日米で反抗的態度に差があるかどうかをWilcoxonの順位和検定を使って検証しましょう。まず、subsetを使って、日本人とアメリカ人の回答を別々に抽出します。

> q1jp <- subset(question1,国籍 == "日本")

> q1us <- subset(question1,国籍 == "アメリカ")

準備が出来たので、検定を行います。wilcox.testという関数を利用します。

> wilcox.test(x=q1jp$回答,y=q1us$回答,alternative="two.sided")

Wilcoxon rank sum test with continuity correction

data: q1jp$回答 and q1us$回答

W = 78.5, p-value = 0.02698

alternative hypothesis: true mu is not equal to 0

Warning message:

Cannot compute exact p-value with ties in: wilcox.test.default(x = q1jp$回答, y = q1us$回答, alternative = "two.sided")

基本的にオプションの設定の仕方はt.testと同じです。たとえば、片側検定を行いたい場合は、alternative="two.sided"の代わりに、alternative="greater"または、alternative="less"を使います。

p-valueを見ると、0.026<0.05なので、5％の有意水準で帰無仮説が棄却されます。つまり、日米で反抗的態度（の中央値）に差があるということが分かります。

符号検定

ある病院の入院患者を対象に以下のようなアンケートを入院直後と退院直前の2回実施しました。

看護師はあなたの相談を良く聞いてくれましたか？

1. 非常に良く聞いてくれた 2. 良く聞いてくれた 3. 普通

4. あまり聞いてくれなかった 5. ほとんど聞いてくれなかった

このアンケートの結果が

に保存されています。

> q2 <- read.csv("signed.csv")

> q2

回答者番号 入院直後 退院直前

1 1 2 2

2 2 3 2

3 3 3 3

4 4 4 4

5 5 2 2

6 6 1 1

7 7 2 2

8 8 3 2

9 9 4 3

10 10 3 3

符号検定は以下のコマンドを入力することで実行できます。

> x <- q2$入院直後

> y <- q2$退院直前

> binom.test(c(length(x[x>y]),length(x[x<y])),alternative="greater")

Exact binomial test

data: c(length(x[x > y]), length(x[x < y]))

number of successes = 5, number of trials = 7, p-value = 0.2266

alternative hypothesis: true probability of success is greater than 0.5

95 percent confidence interval:

0.3412614 1.0000000

sample estimates:

probability of success

0.7142857

p-valueは0.22(>0.05)より、5%の有意水準で帰無仮説を棄却しません。つまり、入院直後と退院直前で看護に対する印象に差がないと判断されます。「alternative="greater"」は「xがyより大きい」（この場合は入院直後の方が退院直前より大きい）という対立仮説で検定することを意味します。もし、「xよりyが大きい」という対立仮説で検定したい場合、「alternative="less"」とします。

Kruskal-Wallis検定

以下のようなアンケートを実施しました。

問1 あなたの国籍は何ですか （日本 ・ 中国 ・ アメリカ）

問2 あなたは親に反抗しますか？

1. よくする 2.時々する 3.あまりしない 4.ほとんどしない

このアンケートの結果が

に保存されています。

> question3 <- read.csv("wilcoxon3.csv")

検定を行うにはkruskal.testという関数を利用します。

> kruskal.test(question3$回答,question3$国籍)

Kruskal-Wallis rank sum test

data: question3$回答 and question3$国籍

Kruskal-Wallis chi-squared = 9.8041, df = 2, p-value = 0.007431

この結果から、5％の有意水準で帰無仮説を棄却します。つまり、各国で親に対する反抗的態度に差が見られると判断されます。