計量経済学レポート

Rの本まだ届かず・・・

SPSSを使って作成したレポートを、同じようにRでもやってみるということをしました。


テーマが「公開データを使って、仮説を立てて検証せよ」というものだったのですが
仮説がなかなか立てられなくて夏休み中困っていましたね。

大学で経済学検定なるものを受けなければならなかったのでその勉強もしていたところ
「上級財と下級財」という概念を知り、例に「牛肉は上級財で豚肉は下級財」というものが
あったので、牛肉と豚肉それぞれの購入金額を被説明変数とし可処分所得を説明変数とした
回帰分析を行って提出しました。

今見ると、びっくりするくらい内容薄かったです・・・

ただ当時は大学に入学したての前期で授業にもついていけてない状態だったことを考えると
まぁはじめの一歩としてそれなりにがんばったと思ってもよいかもしれません（笑）


前置きが長くなりましたが・・・

＞＞＞

> mydata=read.csv("上級財と下級財.csv")
> mydata
year shotoku ushi buta
1 2000 474411 26152 21728
2 2001 466003 21157 22396
3 2002 453716 19982 23434
4 2003 440667 21374 21858
5 2004 446288 20918 23362
6 2005 441156 21324 23191
7 2006 441448 20705 23249
8 2007 442504 20868 23923
9 2008 442749 20885 25555
10 2009 427912 20167 24791
11 2010 429967 18964 23957

●●●　牛は上級財？　●●●

> plot(ushi~shotoku)



> summary(lm(ushi~shotoku))

Call:
lm(formula = ushi ~ shotoku)

Residuals:
Min 1Q Median 3Q Max
-1902.19 -428.06 69.94 711.27 2281.40

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.191e+04 1.283e+04 -1.708 0.12185
shotoku 9.651e-02 2.875e-02 3.357 0.00843 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1268 on 9 degrees of freedom
Multiple R-squared: 0.5559, Adjusted R-squared: 0.5066
F-statistic: 11.27 on 1 and 9 DF, p-value: 0.008434

●●●　豚は下級財？　●●●

> plot(buta~shotoku)



> summary(lm(buta~shotoku))

Call:
lm(formula = buta ~ shotoku)

Residuals:
Min 1Q Median 3Q Max
-1824.8 -335.4 -31.0 379.4 1979.6

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.640e+04 9.679e+03 4.794 0.000982 ***
shotoku -5.155e-02 2.169e-02 -2.377 0.041446 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 956.9 on 9 degrees of freedom
Multiple R-squared: 0.3856, Adjusted R-squared: 0.3174
F-statistic: 5.649 on 1 and 9 DF, p-value: 0.04145

＜＜＜

CSVファイル中の数値に桁区切り文字が入っていると散布図がおかしくなりました。
（棒状？になりました）

原因がなかなかわからなくて時間を取られましたが、ひとつ勉強になりました～。

ヒストグラムと正規分布を重ねて描画

できました～ヾ(＠°▽°＠)ﾉ

> hist(nichijishuekiritsu, xlim=c(-3,3), ylim=c(0,25), xlab="" , ylab="" , main="日次収益率")
> par(new =T)
> plot(function(x) 1/sqrt(2*pi)*exp(-x^2/2),-3,3, type = "l", axes = FALSE, ylab = "", col = 'blue')

記述統計量＆各種グラフ

またまた3年前の宿題より。

（1）TOPIX終値の時系列グラフ
> mydata=read.csv("TOPIX.csv")
> mydata
hiduke owarine nichijishuekiritsu
1 2012-01-04 742.99 0.00000000
2 2012-01-05 736.28 -0.90310771
3 2012-01-06 729.60 -0.90726354
4 2012-01-10 731.93 0.31935307
5 2012-01-11 733.47 0.21040263
6 2012-01-12 727.15 -0.86165760
7 2012-01-13 734.60 1.02454789
8 2012-01-16 725.24 -1.27416281
9 2012-01-17 731.53 0.86729910
10 2012-01-18 734.98 0.47161429
　　　　　　　　　・
　　　　　　　　　・
　　　　　　　　　・

> attach(mydata)
> plot(owarine~hiduke)
> lines(owarine~hiduke)




（2）日次収益率の時系列グラフ
> plot(nichijishuekiritsu~hiduke)
> lines(nichijishuekiritsu~hiduke)




（3）日次収益率の平均・分散
> describe(nichijishuekiritsu)
Description of structure(list(x = c(0, -0.903107713, -0.907263541, 0.31935307, 0.210402634, -0.8616576, 1.024547892, -1.27416281, 0.867299101, 0.471614288, 0.775531307, 1.996813739, 0.174725667, 0.080603602, 1.320306311, -0.363565285, -0.45513399, -0.541300435, -0.229851653, 0.356164021, 0.592379545, -0.230834809, 1.204169898, 0.379294668, 1.238402112, 0.274816576, -0.690894721, 0.335014825, 0.654999488, 2.053889171, -0.337501245, 1.274601687, 1.058671109, -0.33454208, 1.116024942, 0.478555852, 0.595647194, 0.115067902, 0.386710566, -0.300543841, -0.528733432, 0.755225245, -0.592012604, -0.661575775, -0.560826736, 1.634840952, 1.500908917, -0.404142758, 0.0059152, 1.393538618, 0.758362404, 0.361274186, 0.186909418, -1.102090171, 0.383101609, -1.106638672, -0.083281527, 2.418351295, -0.915843286, -0.773920387, -0.39522466, 0.198981682, -0.587582501, -1.840144767, -0.333987742, -0.823954743, -1.455716898, -0.031953201, -0.933085822, 0.501340216, 0.691460463, -1.428606465, -0.092059266, 2.014718151, -0.627387796, -0.268998809, -0.295588344, -0.691750871, 0.690350026, 0.075356088, -0.719664239, -1.837691323, 0.428124485, -2.624642123, 0.584151491, -1.383004752, -0.053536686, -0.919756474, -0.22416203, -1.226410107, -1.139951833, 1.120614985, -2.893623856, -0.053753067, 1.128042474, -1.603643653, 0.094238951, -0.01938387, -0.138483057, 0.820956581, -0.469031539, -0.570741549, -1.467706292, -1.892993667, 1.83031157, 1.457133175, 1.696448452, -1.780362641, 1.717892273, -0.780747052, 0.285765562, -0.107372942, 0.125403081, 1.684627772, -0.557653524, 1.721814643, 0.885808334, -0.403204414, -0.759068876, -0.849413596, 0.891878358, 1.788109674, 1.485220279)), .Names = "x", row.names = c(NA, -123L), class = "data.frame")

Numeric
mean median var sd valid.n
x 0.03 -0.02 1.06 1.03 123


（4）日次収益率の歪度・尖度
> install.packages("e1071",rep="http://cran.ism.ac.jp")
> library(e1071)

> skewness(nichijishuekiritsu)
[1] -0.05045585

> kurtosis(nichijishuekiritsu)
[1] -0.1926556


（5）日次収益率のヒストグラム
> hist(nichijishuekiritsu, xlim=c(-3,3), ylim=c(0,25), xlab="" , ylab="" , main="日次収益率")




ヒストグラムに正規分布を重ね合わせるのが色々試してもできなかったのですが
もう明日に回そうと思います。

ジニ係数とローレンツ曲線

Rの本がまだ届かないので、昔の宿題でもやってみることとします。

（１）パッケージのインストール＆呼び出し
> install.packages("ineq",rep="http://cran.ism.ac.jp")
> library(ineq)

（２）ジニ係数を求める
> ineq(c(300,600,800,1100,2200),type="Gini")
[1] 0.344

（３）ローレンツ曲線を描画
> plot(Lc(c(300,600,800,1100,2200)),col="darkred",lwd=2)



＞＞＞

（１）データ読み込み
> mydata=read.csv("H21県民所得.csv")
> mydata
todoufuken H21kenminshotoku
1 東京都 3,907
2 神奈川県 3,086
3 愛知県 2,970
4 滋賀県 2,955
5 静岡県 2,926
6 千葉県 2,917
7 大阪府 2,879
8 埼玉県 2,867
9 栃木県 2,859
10 京都府 2,815
　　　　　・
　　　　　・
　　　　　・

（２）ジニ係数を求める
> ineq(mydata$H21kenminshotoku,type="Gini")
[1] 0.3262411

（３）ローレンツ曲線を描画
> plot(Lc(mydata$H21kenminshotoku),col="darkred",lwd=2)



＜＜＜

3年前のレポートなのですが、最初の問題は一致しているものの
次の問題がRの結果と合っていません。

間違えたかも～～

・・・って、3年後に気づいてもなぁ（汗）

回帰分析（４）

教科書の練習問題をもう1問やってみることにしました。

「大手メーカーの株価収益率（被説明変数）とTOPIXの変化率（説明変数）のデータから
　回帰モデルを推定せよ」


（１）データ読み込み
> shueki<-c(-5.2,-2.3,0.5,4.1,-5.8,-2.7,3.8,-6.8,-2.5,1.9,2.5,3.1,-3.4,-2.4,1.5,2.7,5.1)
> topix<-c(-3.5,-1.1,0.6,2.7,-3.3,-1.5,2.4,-4.6,-1.6,1.2,1.9,2.4,-2.2,-1.4,1.3,1.9,3.4)
> mydata94<-data.frame(shueki,topix)
> mydata94
shueki topix
1 -5.2 -3.5
2 -2.3 -1.1
3 0.5 0.6
4 4.1 2.7
5 -5.8 -3.3
6 -2.7 -1.5
7 3.8 2.4
8 -6.8 -4.6
9 -2.5 -1.6
10 1.9 1.2
11 2.5 1.9
12 3.1 2.4
13 -3.4 -2.2
14 -2.4 -1.4
15 1.5 1.3
16 2.7 1.9
17 5.1 3.4


（２）データ間の関係を見るため、相関行列を求める
> cor(mydata94)
shueki topix
shueki 1.0000000 0.9969249
topix 0.9969249 1.0000000


（３）線形関係を見るため、対散布図を描画する
> pairs(mydata94)




（４）回帰分析
> lm.shueki=lm(mydata94$shueki~mydata94$topix)
> summary(lm.shueki)

Call:
lm(formula = mydata94$shueki ~ mydata94$topix)

Residuals:
Min 1Q Median 3Q Max
-0.56719 -0.20035 0.03701 0.22227 0.40676

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.22201 0.07407 -2.997 0.00902 **
mydata94$topix 1.51842 0.03082 49.271 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3052 on 15 degrees of freedom
Multiple R-squared: 0.9939, Adjusted R-squared: 0.9934
F-statistic: 2428 on 1 and 15 DF, p-value: < 2.2e-16

回帰モデル推定：株価収益率＝－0.22201＋1.51842*TOPIX変化率


一応、教科書の回答とは一致しています。

ただ、簡単すぎて、回答がない問題の結果はまだ自分を信じられないかも（笑）

回帰分析（３）

教科書の例題もやってみました。

「消費支出を被説明変数、所得・預金残高を説明変数とした回帰モデルを推定し
　その結果を吟味せよ。但し誤差項は仮定５つが満たされているものとする」

勉強会の資料は単回帰だったので、ネットで重回帰分析のサンプルを参照しました。


（１）データ読み込み（事前にデータをCSVファイルに入力）
> mydata91=read.csv("ex9.1.csv")
> mydata91
shishutsu shotoku zandaka
1 59 70 271
2 69 74 388
3 50 60 213
4 83 93 450
5 71 81 363
6 57 66 288
7 76 84 425
8 57 63 304
9 63 75 283
10 71 89 304
11 60 75 250
12 61 78 225


（２）データ間の関係を見るため、相関行列を求める
> round(cor(mydata91),4)
shishutsu shotoku zandaka
shishutsu 1.0000 0.9165 0.8856
shotoku 0.9165 1.0000 0.6295
zandaka 0.8856 0.6295 1.0000


（３）線形関係を見るため、対散布図を描画する
> pairs(mydata91)




（４）回帰分析
> mydata91.lm<-lm(shishutsu~.,data=mydata91)
> summary(mydata91.lm)

Call:
lm(formula = shishutsu ~ ., data = mydata91)

Residuals:
Min 1Q Median 3Q Max
-0.79821 -0.34123 0.09471 0.28730 0.53657

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.282561 1.113943 2.947 0.0163 *
shotoku 0.553101 0.018687 29.598 2.80e-10 ***
zandaka 0.062538 0.002457 25.454 1.07e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4852 on 9 degrees of freedom
Multiple R-squared: 0.9978, Adjusted R-squared: 0.9973
F-statistic: 2048 on 2 and 9 DF, p-value: 1.082e-12


（５）回帰診断図描画
> par(mfrow=c(2,2),oma = c(1,1,2,1),mar = c(4, 4, 2, 1))
> plot(mydata91.lm,pch=21,bg=2,col=2,cex=1.5)




有意水準５％では、どの係数も有意であり、係数も教科書と一致したので
使い方は間違っていないものと思われます。

ただ、回帰診断図の見方がよくわからず・・・

回帰分析（２）

男性の場合の単回帰モデル（被説明変数：収縮期血圧／説明変数：肥満度）を
あてはめよとの演習をやってみます。

恐らく何も考えずともできるはずですが、解答がないので
間違えないよう一応慎重になりたいと思います。


（１）散布図を描画（男性のみ）
> man=bp.obese[bp.obese$sex==0,]
> plot(man$bp~man$obese)




（２）manをグローバル化
> attach(man)
The following objects are masked from woman:

bp, obese, sex

womanとかぶっていてダメなようです。

> detach(woman)
> attach(man)
The following objects are masked from man (pos = 3):

bp, obese, sex

最初のattachが効いているということなのでしょうか・・・

> detach(man)
> attach(man)
The following objects are masked from man (pos = 3):

bp, obese, sex

でも同じメッセージがでます。訳わかめ

よくわかりませんが、変数の中身を書き出してみたところ
一致しているようだったのでこのまま進めます。

> man
sex obese bp
1 0 1.31 130
2 0 1.31 148
3 0 1.19 146
4 0 1.11 122
5 0 1.34 140
6 0 1.17 146
7 0 1.56 132
8 0 1.18 110
9 0 1.04 124
10 0 1.03 150
11 0 0.88 120
12 0 1.29 114
13 0 1.26 136
14 0 1.16 118
15 0 1.32 190
16 0 1.37 118
17 0 1.25 130
18 0 1.48 112
19 0 1.58 126
20 0 0.93 162
21 0 1.29 124
22 0 1.06 126
23 0 1.19 134
24 0 0.96 110
25 0 1.13 118
26 0 1.19 110
27 0 0.81 94
28 0 1.11 118
29 0 1.29 140
30 0 1.29 128
31 0 1.28 126
32 0 1.20 140
33 0 1.02 124
34 0 1.09 104
35 0 1.08 134
36 0 1.04 130
37 0 1.14 124
38 0 1.13 110
39 0 1.16 134
40 0 1.57 144
41 0 1.07 116
42 0 1.04 118
43 0 1.37 118
44 0 1.26 132
> bp
[1] 130 148 146 122 140 146 132 110 124 150 120 114 136 118 190 118 130 112 126 162
[21] 124 126 134 110 118 110 94 118 140 128 126 140 124 104 134 130 124 110 134 144
[41] 116 118 118 132
> obese
[1] 1.31 1.31 1.19 1.11 1.34 1.17 1.56 1.18 1.04 1.03 0.88 1.29 1.26 1.16 1.32 1.37
[17] 1.25 1.48 1.58 0.93 1.29 1.06 1.19 0.96 1.13 1.19 0.81 1.11 1.29 1.29 1.28 1.20
[33] 1.02 1.09 1.08 1.04 1.14 1.13 1.16 1.57 1.07 1.04 1.37 1.26


（３）回帰分析（男性のみ）
> lm.bp=lm(bp~obese)
> summary(lm.bp)

Call:
lm(formula = bp ~ obese)

Residuals:
Min 1Q Median 3Q Max
-25.645 -9.533 -1.598 7.060 59.315

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 102.11 17.49 5.839 6.76e-07 ***
obese 21.65 14.50 1.493 0.143
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.36 on 42 degrees of freedom
Multiple R-squared: 0.05038, Adjusted R-squared: 0.02777
F-statistic: 2.228 on 1 and 42 DF, p-value: 0.143

「推定された回帰式」を書きなさいとありますが
obeseが有意じゃなくてもとりあえずは当てはめてみるべきなのか？

よくわからないけど、
とりあえず当てはめたら「血圧＝102.11＋21.65*肥満度」


（５）回帰直線を散布図上に描く
> plot(bp~obese)
> abline(lm.bp,col="red")




womanに比べたら傾向も弱くまとまりもなく外れ値もあるので
こんな感じの結果でよいんでしょうか～なんかあやふやですが（汗）

回帰分析（１）

またまたRの続きをやっています。

いよいよ最終章～というのに気づくのが遅れ、慌てて本を発注～
今度こそは２１日を途切れさせないようにと思ったのに・・・

Rを使って実際に検定をしてみるという作業はなかなか楽しいです。


（１）パッケージのデータセットを使う（id・性別・肥満度・収縮期血圧）
> install.packages("ISwR",rep="http://cran.ism.ac.jp")
Installing package into ‘C:/Users/madoka/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.ism.ac.jp/bin/windows/contrib/3.2/ISwR_2.0-6.zip'
Content type 'application/zip' length 217390 bytes (212 KB)
downloaded 212 KB

package ‘ISwR’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\madoka\AppData\Local\Temp\RtmpcLNG5d\downloaded_packages

> library(ISwR)
> data(bp.obese)
> bp.obese
sex obese bp
1 0 1.31 130
2 0 1.31 148
3 0 1.19 146
4 0 1.11 122
5 0 1.34 140
6 0 1.17 146
　　　　・
　　　　・
　　　　・


（２）散布図を描画（女性のみ）
> woman=bp.obese[bp.obese$sex==1,]
> plot(woman$bp~woman$obese)




（３）womanをグローバル化
> attach(woman)


（４）回帰分析（女性のみ）
> lm.bp=lm(bp~obese)
> summary(lm.bp)

Call:
lm(formula = bp ~ obese)

Residuals:
Min 1Q Median 3Q Max
-23.522 -12.348 -2.203 3.907 71.500

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.517 12.048 6.849 6.14e-09 ***
obese 31.204 8.426 3.703 0.000488 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.56 on 56 degrees of freedom
Multiple R-squared: 0.1967, Adjusted R-squared: 0.1824
F-statistic: 13.72 on 1 and 56 DF, p-value: 0.000488

回帰式：血圧＝82.517＋31.204*肥満度


（５）回帰直線を散布図上に描く
> plot(bp~obese)
> abline(lm.bp,col="red")





演習問題として「男性の場合の単回帰モデルをあてはめなさい」というものが出ているので
続きを別記事にて行いたいと思います。

Rではじめるビジネス統計分析

回帰分析で勉強会資料が終わってしまうので、本を買うことにしました。

Rの入門書とか、Rを使用した統計学の入門書など様々なものがあるようですが
「本書の構成」として紹介されていた内容が実践的に思えたのでこちらを選択しました。

Rではじめるビジネス統計分析/翔泳社

￥3,456
Amazon.co.jp

■Part1:R の基本
・Chapter1:R の基本操作
この章ではR のインストールと基本的な操作方法について解説します。

■Part2:ビジネス統計分析の基本
・Chapter2:データを視覚化する
この章では数字のデータをわかりやすくする「視覚化」について解説します。

・Chapter3:データを要約する
この章では社内に膨大に存在しているさまざまなデータを要約する方法について解説します。
またそれらのデータを分析する前にどのような考え方でのぞめばよいのかについても解説します。

・Chapter4:データの関連性を見る
この章では相関分析を利用してデータの関連性を見る方法について解説します。

・Chapter5:未知のデータを予測する
この章では単回帰分析を利用して、未知のデータを予測する分析手法について解説します。

・Chapter6:データを分類する
この章では主成分分析、因子分析、コレスポンデンス分析、多次元尺度構成法、決定木を利用した分析手法について解説します。

■Part3:本格的なビジネス統計分析
・Chapter7:テキストマイニングを行う
テキストマイニングとは、大量の言語データから、コンピュータを使って有益な情報を取り出す技術のことです。
このような技術を用いることで、アンケートの自由記述文やブログ、Twitter
といったSNS 上のクチコミデータを定量的に分析することが可能になります。

・Chapter8:ログデータからクラスタ分析を行う
この章では、ビジネスの現場で、統計ツールを組み合わせて利用する方法について解説します。
具体的には、ソーシャルゲームのサービス改善にあたり、どのように統計ツールを利用していくのかを紹介します。

独立性の検定

Rの続きです。

２つの質的変数（ここでは「病気」「性別」）が独立か？を調べるというものです。


１）データを変数へ格納（結果のみ表示）
> mydata2
id sex disease
1 1 男 有
2 2 男 有
3 3 男 有
4 4 男 無
5 5 男 無
6 6 男 無
7 7 女 有
8 8 女 有
9 9 女 有
10 10 女 有
11 11 女 無
12 12 女 無
13 13 女 無


２）クロス集計表作成
> table(mydata2$sex,mydata2$disease)

無 有
女 3 4
男 3 3


３）カイ２乗検定
> chisq.test(table(mydata2$sex,mydata2$disease),correct=F)

Pearson's Chi-squared test

data: table(mydata2$sex, mydata2$disease)
X-squared = 0.066327, df = 1, p-value = 0.7968

Warning message:
In chisq.test(table(mydata2$sex, mydata2$disease), correct = F) :
カイ自乗近似は不正確かもしれません

結論：病気に性差はない（但しサンプル数が少なすぎるため警告あり）


４）フィッシャーの正確検定
> fisher.test(table(mydata2$sex,mydata2$disease))

Fisher's Exact Test for Count Data

data: table(mydata2$sex, mydata2$disease)
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.05228419 10.67349535
sample estimates:
odds ratio
0.7668651

結論：病気に性差はない


＞＞＞＞

教科書の例題と回答でも確認してみました。

「ある銀行が融資先企業２００社に対して、その企業の戦略が
　優れているか否かを判定した。優れた企業戦略と企業業績に
　関連があるか有意水準５％で検定しなさい」

　　　　　　業績悪／業績良　計
----------------------------------------
　　戦略悪：４５／４５　　　９０
　　戦略優：３０／８０　　　１１０
----------------------------------------
　　　　　　７５／１２５　　２００


１）クロス集計表を変数へ格納
> table82=matrix(c(45,30,45,80),2,2)
> table82
[,1] [,2]
[1,] 45 45
[2,] 30 80


２）カイ２乗検定
> chisq.test(table82,correct=F)

Pearson's Chi-squared test

data: table82
X-squared = 10.909, df = 1, p-value = 0.0009569

結論；企業戦略と企業戦略には関連がある（＝独立を棄却）


＞＞＞

「３都市で、都市と商品の好みに関連があるかを調べるため
　４商品A・B・C・Dのどれが一番好きかの調査を行った。
　都市と商品の好みに関連があるか有意水準５％で検定しなさい」

　　東京／大阪／福岡　計
------------------------------
　A:３６／１４／１４　６４
　B:１０／１０／１０　３０
　C:２１／１４／１５　５０
　D:１８／２７／１１　５６
------------------------------
　　８５／６５／５０　２００


１）クロス集計表を変数へ格納
> table812=matrix(c(36,10,21,18,14,10,14,27,14,10,15,11),4,3)
> table812
[,1] [,2] [,3]
[1,] 36 14 14
[2,] 10 10 10
[3,] 21 14 15
[4,] 18 27 11

２）カイ２乗検定
> chisq.test(table812,correct=F)

Pearson's Chi-squared test

data: table812
X-squared = 13.879, df = 6, p-value = 0.03102

結論：都市と商品の好みに関連がある（＝独立を棄却）