研究・教育の最近のブログ記事

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

を利用して主成分分析を行います。まず、データを読み込みます。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. 「ラベル」にチェック

あなたの関心ある事について、2つの母集団からデータを抽出し、その平均に差があるかどうかを有意水準5%で検定しなさい。もし、そのデータに対応がない場合ははじめに等分散の検定を行い、その結果を考慮した上で、平均の差に関する仮説検定を行いなさい。分析ツールの出力結果と、分析に使ったデータを添付すること。

2つの母集団の分散が等しいかどうかを検定します。母集団1と母集団2の分散をそれぞれとすると、帰無仮説と対立仮説は、

となります。

この検定方法に相当するものとして、Excelの分析ツールに「F検定：2標本を使った分散の検定」というのがあります。しかし、この方法は誤った結果を表示するので使えません。

このファイルには、測定装置A、測定装置Bによってある原料の有効成分含有率を測定した結果が保存されています。

ファイルを開くと次のようなデータが表示されます。

「ツール(T)」-「分析ツール（D)」を選択し、「Ｆ検定：2標本を使った分散の検定」を選択します。

新しいダイアログボックスが表示されるので、次のように設定します。

「OK」ボタンを押すと次のような結果が表示されるはずです。

「観測された分散比」が検定統計量を表し、「Ｆ境界値 片側」が片側検定における臨界値を表しています。

両側検定の臨界値が計算されていないので、付け加えましょう。まず、セルA24に以下のように入力します。

次に以下の式を入力します。

ここで、FINV(x,m,n)は自由度(m,n)のF分布の右裾の面積がxになるような値を計算する関数です。

これで、棄却域が0.23以下と、4.53以上ということが分かりました。検定統計量は0.17なので、棄却域に入ります。このことから有意水準5%で2つの母集団の分散は異なると考えられます。

2つの母集団の平均の差の検定を行います。ここでは、2つの母集団の分散が等しくないと仮定する場合を説明します。

以下のデータは、生の落花生と、炒った落花生をそれぞれ10匹のネズミに与えたときのタンパク質の摂取量です。生の落花生といった落花生でタンパク質の平均摂取量が異なるかどうか検定します。

「ツール(T)」-「分析ツール（D)」を選択し、「ｔ検定：分散が等しくないと仮定した2標本による検定」を選択します。

新しいダイアログボックスが表示されるので、次のように設定します。

αは有意水準を表します。この場合5%です。

「OK」ボタンを押すと次のような結果が表示されるはずです。

B9は検定統計値、B11とB13はそれぞれ片側検定、両側検定の時の臨界値を示します。この場合、検定統計値が臨界値より小さい値なので、帰無仮説は受容されます。つまり、生の落花生と、炒った落花生のどちらを与えてもタンパク質平均摂取量に違いはないということになります。

B10,B12はそれぞれ片側検定、両側検定の時のP値を示します。両者とも0.05より大きいので、片側検定、両側検定どちらにおいても、帰無仮説を受容することがわかります。

2つの母集団の平均の差の検定を行います。ここでは、2つの母集団の分散は未知だが、等しいと仮定する場合を説明します。

まず、例として次のようなデータを考えます。これは1978年と1988年の新生児の体重について、10人ずつ標本抽出したものです。

「ツール(T)」-「分析ツール（D)」を選択し、「ｔ検定：等分散を仮定した2標本による検定」を選択します。

新しいダイアログボックスが表示されるので、次のように設定します。

仮説平均との差異は、検定したい2つの母平均の差を入力します。今回は母平均に差はないという仮説を検定するので0とします。

αは有意水準を表します。この場合、5%です。

以下のような結果が表示されるはずです。

B10に検定統計値が、B12,B14に臨界値が表示されています。この場合、2.22は臨界値を超えているので、2つの母集団の平均に差がないという帰無仮説は棄却されます。

2つの母集団の平均の差の検定を行います。ここでは、2つの母集団の分散は未知だが、等しく、また、データに対応がある場合を考えます。

次のようなデータを入力しましょう。これは、10人の被験者の体温を電子式体温計と水銀式体温計で測った結果です。

「ツール(T)」-「分析ツール（D)」を選択し、「ｔ検定：一対の標本による平均の検定」を選択します。

次のようにダイアログボックスに入力します。

「OK」を押すと、次のような結果が得られます。

B10に表示されているのが検定統計値です。片側検定と両側検定の臨界値はB12、B14に表示されています。

この場合、検定統計値は棄却域に入るので、電子式体温計と水銀式体温計の測定結果には差があると考えられます。