www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

乗用車ブランド通称名別順位のデータ分析4 - 日産とホンダの販売台数、順位に違いはあるか?

 

www.crosshyou.info

 の続きです。

tapply関数でメーカー別の販売台数を見てみましょう。

f:id:cross_hyou:20191116152739p:plain

トヨタが204,628台でダントツですね。2位がホンダの63,507台、3位が日産の62,943台です。

メーカーごとの平均順位はどうでしょうか?

f:id:cross_hyou:20191116152939p:plain

トヨタが15.78で1位ですね。2位が日産で24.75、3位がホンダで27.17です。日産とホンダの販売台数を比較してみましょう。

まず、日産とホンダだけのデータフレームを用意します。

f:id:cross_hyou:20191116153658p:plain

Makerの列を見ると、ホンダが18、日産が12となっていますので、ホンダのほうがランクインしている会社が多いということですね。

日産の販売台数のベクトル、ホンダの販売台数のベクトルを作ります。

f:id:cross_hyou:20191116154226p:plain

それぞれの統計値を見てみます。

f:id:cross_hyou:20191116154411p:plain

日産は平均値は5245台、平均値の95%信頼区間は2266台から8223台です。

ホンダはどうでしょうか?

f:id:cross_hyou:20191116154654p:plain

ホンダは平均値は3528台、平均値の95%の信頼区間は2345台から4710台です。

両社の平均値に有意な違いは無さそうです。

データの数が日産は12、ホンダは18と少ないので、Wilcoxon Rank-Sum Testで調べましょう。wilcox.test関数です。

f:id:cross_hyou:20191116155024p:plain

p-value = 0.7321なので両社の販売台数の間に有意な違いはないようです。

グラフで見てみます。

f:id:cross_hyou:20191116155536p:plain

f:id:cross_hyou:20191116155547p:plain

赤いのが日産、青いのがホンダです。日産は上位2種類のクルマが目立ちますね。

両社のランクについても同じようにしてみましょう。

f:id:cross_hyou:20191116160024p:plain

グラフにしてみます。

f:id:cross_hyou:20191116160219p:plain

f:id:cross_hyou:20191116160229p:plain

日産(赤い点)のほうが範囲が大きいですね。統計値を調べます。

f:id:cross_hyou:20191116160445p:plain

日産の順位の平均値の95%信頼区間は12から37、ホンダは20から34ですので両社のあいだに有意な違いは無いですね。wilcox.testで確認します。

f:id:cross_hyou:20191116160735p:plain

p-valueが0.8323なのでやはり両社に違いがあるとは言えないですね。

今回は以上です。

 

乗用車ブランド通称名別順位のデータ分析3 - R言語でヒストグラムや箱ひげ図を作成する。1月と2月の販売台数をt.test関数、wilcox.test関数で検定。

 

www.crosshyou.info

 の続きです。

今回は販売台数をグラフにしてみましょう。

まずは、小さい順グラフ、ヒストグラム、箱ひげ図の3つのグラフを一度に作成する関数を作ります。

f:id:cross_hyou:20191116095130p:plain

gpという名前の関数を作りました。

それでは販売台数をグラフにしてみます。

f:id:cross_hyou:20191116095358p:plain

f:id:cross_hyou:20191116095424p:plain

青い線が中央値、赤い線が平均値なので、平均値のほうが大きいことがわかります。ヒストグラムの形状は右の裾野が広いですね。箱ひげ図では、値が小さいほうには外れ値は無いですが、値が大きいほうに外れ値が多いことがわかります。

1月だけ、2月だけでグラフを作成してみましょう。

f:id:cross_hyou:20191116095855p:plain

f:id:cross_hyou:20191116095906p:plain

f:id:cross_hyou:20191116100046p:plain

f:id:cross_hyou:20191116100058p:plain

1月も2月も分布の形状は同じような形状ですね。

1月だけの販売台数の平均値や信頼区間をみてみます。前回に作成したsummary2関数を使います。

f:id:cross_hyou:20191116100345p:plain

1月の販売台数の平均値は3700台、varianceは7095741、変動係数は0.720、平均値の95%の信頼区間は2943台から4456台です。

2月だけの販売台数はどうでしょうか?

f:id:cross_hyou:20191116101856p:plain

2月の販売台数の平均値は4435台、varianceは9877172、変動係数は0.709, 平均値の95%の信頼区間は3542台から5328台です。

1月と2月の平均値は、信頼区間が重なっていますので、違いがあるとは言えなのではないでしょうか?平均値が同じかどうかの検定をしてみます。

参考図書は STATISTICS : An Introduction using R by Michael J. Crawleyです。

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

 

 まず1月と2月でvarianceが同じかどうかを検定します。

1月のvarianceが7095741で2月のvarianceが9877172です。二つのvarianceの比率がF ratioです。

f:id:cross_hyou:20191116102600p:plain

1月も2月もデータ数は50ですから、自由度は49です。二つのvarianceが有意に違わないという帰無仮説が正しいという確率を求めます。

f:id:cross_hyou:20191116103138p:plain

0.25ですから、2つのvarianceが有意に違うとは言い切れないですね。つまり、2つのvarianceが同じとみなして検定してもいい、ということです。

var.test関数でも確認しましょう。

f:id:cross_hyou:20191116103503p:plain

p-value = 0.2505とpf関数で計算した値と同じになりました。

二つのvarianceが同じとき、平均値が同じかどうかを検定するのはt.test関数を使います。やってみましょう。

f:id:cross_hyou:20191116104016p:plain

p-value = 0.2102なので、0.05よりも大きな値ですから、2つの平均値には有意な違いは無い、という帰無仮説を棄却できません。つまり、2つの平均値には有意な違いが無いということです。

Wilcoxson's Rank-Sum Testもやってみましょう。wilcox.test関数を使います。

f:id:cross_hyou:20191116104506p:plain

p-value = 0.1985と0.05よりも大きいです。この検定でも1月と2月の販売台数の平均値に違いがあるとは言えません。

今回は以上です。

 

 

乗用車ブランド通称名別順位のデータ分析2 - R言語で自作関数を作成。標準偏差、標準誤差、変動係数、信頼区間などを算出。

 

www.crosshyou.info

 の続きです。

今回は標準偏差、標準誤差、変動係数などを計算する関数を作成してみます。

まずは、前回保存したCSVファイルを読み込んでsummary関数で様子を見てみます。

f:id:cross_hyou:20191114195714p:plain

こうなりました。summary関数だと数値データは最小値、第1分位、中央値、平均値、第3分位、最大値、NAの数は表示されますが、標準偏差や変動係数などが表示されません。そこで自作の関数を作ってみました。

f:id:cross_hyou:20191114200828p:plain

こんな感じで作ってみました。

それでは、NumとYoYで試してみます。

まずはNum、販売台数から

f:id:cross_hyou:20191114201122p:plain

r = 2としたら、ちゃんと小数点以下2桁までになりましたね。標準偏差は2921、変動係数は0.72、1標準誤差での平均値の推定は、3775 ~ 4359(n = 100)です。

95%信頼区間での平均値の推定は、34887 ~ 4647(n = 100)です。

前年比(YoY)も見てみます。

f:id:cross_hyou:20191114201508p:plain

こちらもうまくいきました。前年比は8個データが欠損しています。変動係数が1.56ですから販売台数よりもデータのばらつきが大きいことがわかります。

今回は以上です。

 

 

乗用車ブランド通称名別順位のデータ分析1 - R言語のXMLパッケージのreadHTMLTable関数でWebから直接データを読み込む。

今回は、R言語のXMLパッケージのreadHTMLTable関数でWebの表形式のデータを直接読込んでみます。

参考図書は、Rクックブックです。

 

Rクックブック

Rクックブック

 

 まず、XMLパッケージをinstall.packages関数でインストールします。

f:id:cross_hyou:20191109101702p:plain

次にlibrary関数でXMLパッケージを呼び出します。

f:id:cross_hyou:20191109101810p:plain

これで準備完了です。

今回は、

http://www.jada.or.jp/data/month/m-brand-ranking/

www.jada.or.jp

このWebサイトの表データを読込んでみようと思います。

f:id:cross_hyou:20191109102118p:plain

こういう表です。

readHTMLTable関数で読み込みたいWebサイトのURLを指定するだけでいいらしいです。やってみます。

f:id:cross_hyou:20191109102256p:plain

tblsにはWebサイトのすべての表が読み込まれています。何個の表があるかlength関数で確認します。

f:id:cross_hyou:20191109102429p:plain

7個の表があるようです。

一つ目の表がどんなものか見てみましょう。

f:id:cross_hyou:20191109102616p:plain

V1からV5までが2019年1月のデータ、V6が空欄で、V7からV11が2019年2月のデータですね。

str関数でデータ構造を見てみましょう。

f:id:cross_hyou:20191109103050p:plain

観測数が51、変数が11のデータフレーム構造ですね。すべての変数がファクタになっています。

分析しやすいように、2019年1月のデータフレーム、2019年2月のデータフレームに分割しましょう。

f:id:cross_hyou:20191109103545p:plain

V1は順位だからRank, V2は社名だからName, V3は会社名だからMaker, V4は販売台数だからNum, V5は前年比だからYoYに変数名を変更しましょう。

f:id:cross_hyou:20191109104000p:plain

Rankをファクタから数値に、Nameをファクタから文字列に、Numをファクタから数値に、YoYをファクタから数値にデータの型を変換したいですね。

Rankは変換というよりは、1:50に置き換えたほうが簡単ですね。NumとYoYはカンマがあるのでこのままでは変換できないですね。gsub関数でカンマを削除してみます。

f:id:cross_hyou:20191109105909p:plain

うまくカンマが削除できましたね。

それではデータ型を変換してみます。

f:id:cross_hyou:20191109110024p:plain
YoYにNAが4つあります。Webサイトを確認してみます。

f:id:cross_hyou:20191109110150p:plain

37位のインサイト、40位のCR-V、42位のES300H、43位のUX250Hの4つが空白または、(18-10)というデータなので前年比が無いですね。(18-10)というのは2018年の10月に販売開始という意味で、空白は今年販売開始なので前年比が無いということです。

is.na関数で確認してみます。

f:id:cross_hyou:20191109110654p:plain

たしかにそうですね。

これで2019年1月のデータフレームはできましたので、2019年2月のデータフレームも作成しましょう。

f:id:cross_hyou:20191109111611p:plain

いい感じですね。

Jan19にMonthという変数名でJanという値を持つ列を、Feb19にもMonthという変数名でFeb]という値を持つ列を加えます。

f:id:cross_hyou:20191109112123p:plain

それでは、Jan19とFeb19を一つのデータフレームにまとめましょう。rbind関数を使います。

f:id:cross_hyou:20191109112416p:plain

次回以降はこのjcarを分析しようと思います。自販連のWebサイトの構成が変わるといけないので、このjcarをCSVファイルとして書き出します。write.csv関数です。

f:id:cross_hyou:20191109113044p:plain

CSVファイルを読込みます。read.csv関数ですね。

f:id:cross_hyou:20191109113301p:plain

うまくできました!

今回は以上です。

 

都道府県別の凶悪犯認知件数の分析6 - R言語で県内総生産当りの凶悪犯認知件数を人口と可住地面積で回帰分析。

 

www.crosshyou.info

 の続きです。

今回は、県内総生産当りの凶悪犯認知件数を人口と可住地面積で回帰分析しようと思います。

まずは、県内総生産当りの凶悪犯認知件数を算出します。

f:id:cross_hyou:20191107190348p:plain

avgGDPは百万円単位なので、百万をかけて1円当りの件数にしています。富山県が一番少なく5.9件、大阪府が一番多く27件です。

それではグラフで分布の様子を見てみましょう。

f:id:cross_hyou:20191107190904p:plain

f:id:cross_hyou:20191107190915p:plain

右の裾野が広い分布ですね。

度数分布表もみてみましょう。

f:id:cross_hyou:20191107191120p:plain

10件から15件のレンジが一番度数が多いです。21件です。

BadGDPとavgPop, BadGDPとavgAreaの散布図を見てみます。

f:id:cross_hyou:20191107191514p:plain

f:id:cross_hyou:20191107191532p:plain

avgAreaのほうは北海道が大きすぎて散布図がよくわからないですね。

対数にしてみます。

 

f:id:cross_hyou:20191107191937p:plain

f:id:cross_hyou:20191107192004p:plain

可住地面積は対数にしたほうがいいかもしれないですね。

それでは、R言語のlm関数で線形重回帰分析をしましょう。

f:id:cross_hyou:20191107192304p:plain

avgPop:log(avg(Area)^2)はいらないですね。

f:id:cross_hyou:20191107192603p:plain

p値が0.9392なのでmodel1もmodel2も有意な違いはありません。従って単純なmodel2のほうを採用します。

f:id:cross_hyou:20191107192756p:plain

I(log(avgArea)^2)はいらないですね。

f:id:cross_hyou:20191107193308p:plain

p値が0.1605と0.05よりも大きいので、model2とmodel3に統計的な有意な違いはありません。なので、より単純なmodel3をさらに調べます。

f:id:cross_hyou:20191107193543p:plain

すべてが***になりました。p値は2.038e-08なので有意なモデルです。

avgPopの係数が正の値です。

I(avgPop^2)の係数が負の値ですが、係数の値は非常に小さいので、avgPopが大きいほど件数は多くなるということですね。

log(avgArea)の係数が負ですから面積が大きいほど件数は少なくなります。

lm関数で作成したモデルには、fitted.valuesという名前でモデルが予測した値が保存されています。実際の値との差を見てみましょう。

f:id:cross_hyou:20191107194728p:plain

モデルが予測した件数よりも実際の件数が少ない県は神奈川県や愛知県です。その反対にモデルが予測した件数よりも実際の件数が多いのは沖縄県、高知県大阪府、埼玉県などです。

今回は以上です。

 

都道府県別の凶悪犯認知件数の分析5 - R言語で可住地面積当りの凶悪犯認知件数を人口と県内総生産で回帰分析。人口・県内総生産額の大きな県ほど件数は多い

 

www.crosshyou.info

 の続きです。

前回作成した、可住地面積10万ha当りの凶悪犯認知件数を人口と県内総生産で回帰分析してみます。

まずはplot関数で散布図を描いてみます。

f:id:cross_hyou:20191103103352p:plain

f:id:cross_hyou:20191103103403p:plain

右側の対数変換した後の散布図のほうが分析するにはよさそうですね。

県内総生産との散布図も描きましょう。

f:id:cross_hyou:20191103104145p:plain

f:id:cross_hyou:20191103104202p:plain

pch = 21とbg = で色付きのドットにしてみました。県内総生産額も対数変換したほうがいいですね。

それでは、lm関数で重回帰分析をしてみます。

f:id:cross_hyou:20191103104640p:plain

p-valueは3.535e-10なので0.05よりも小さいの有意なモデルです。

Intercept以下全ての変数でPr(>|t|)が0.05以下なのですべて有意ということですね。

log(avgPop):log(avgGDP)の係数の符号がプラスなのでこれが大きいと可住地面積当りの凶悪犯認知件数は増加します。log(avgPop)とloh(avgGDP)の係数の符号はマイナスなのでこれが大きいと可住地面積当りの凶悪犯認知件数は減少します。

avgPopとavgGDPの平均値を調べてみます。mean関数を使います。

f:id:cross_hyou:20191103105437p:plain

avgPopの平均値は271万7604人、avgGDPの平均値は11兆2274億77百万円です。

仮の人口が300万人、県内総生産が10兆円の県(普通県としましょう)と

人口が600万人、県内総生産が20兆円(大型県としましょう)と

人口が150万人、県内総生産が5兆円(小型県としましょう)の3つの県があったと仮定して、それぞれの可住地面積当りの凶悪犯認知件数を推測してみます。

f:id:cross_hyou:20191103110427p:plain

まず、上のように仮定の県のデータのベクトル、estPop, estGDPを作りました。

そして、predict関数でモデルの予測値を計算します。

f:id:cross_hyou:20191103110930p:plain

モデルの反応変数は、log(BadArea)なので、logをはずすため、exp関数を使っています。

小型県は、可住地面積10万haあたり約28件、普通県は約48件、大型県は約120件と予測されました。人口、GDPの規模の大きな県ほど可住地面積当り凶悪犯認知件数は増加するようですね。

念のため、手計算でも予測値を計算してみました。

f:id:cross_hyou:20191103111943p:plain

少し値がpredict関数の計算結果とずれていますが、だいたいあっています。おそらくsummary関数では係数が小数点以下4桁までしか表示されていませんが、predict関数ではもっと細かい桁数まで使っているのでしょう。

今回は以上です。

 

都道府県別の凶悪犯認知件数の分析4 - R言語で可住地面積当りの凶悪犯認知件数を計算。度数分布表を出力する関数を作成。

 

www.crosshyou.info

 の続きです。

今回は可住地面積当りの凶悪犯認知件数を計算してみます。

f:id:cross_hyou:20191102143405p:plain

10万ha当りの件数を算出しました。秋田県は8.4件、大阪府は780.5件です。100倍近い差がありますね。

グラフで分布をみてみます。

f:id:cross_hyou:20191102143729p:plain

f:id:cross_hyou:20191102143714p:plain

 

大きな値に外れ値がいっぱいあることがわかります。

summary関数で平均値などを見てみます。

f:id:cross_hyou:20191102143928p:plain

平均値は35.5で中央値は85です。

今回は、度数分布表を作成する関数を作ろうと思います。

参考図書は

 

現場ですぐに使える!  R言語プログラミング逆引き大全 350の極意

現場ですぐに使える! R言語プログラミング逆引き大全 350の極意

 

 です。

ヒストグラムを作成する関数histはただヒストグラムを描くだけでなくていろいろなデータも同時に計算しています。

f:id:cross_hyou:20191102144726p:plain

f:id:cross_hyou:20191102144738p:plain

$breaksがヒストグラムの階級で、$countが度数です。これを利用します。

f:id:cross_hyou:20191102151625p:plain
こんな感じです。

それでは関数を実行してみます。

f:id:cross_hyou:20191102151646p:plain

 

f:id:cross_hyou:20191102151711p:plain

ヒストグラムも出力するようになっています。

度数分布表を見ると38都道府県、8割の都道府県が100件以下に入っています。

この38件の分布を詳しく見てみましょう。

f:id:cross_hyou:20191102152239p:plain

f:id:cross_hyou:20191102152255p:plain


BadArea[BadArea < 100]とすると、100件以下の都道府県だけに絞り込むことができます。

100件以下の都道府県にすると山型に近い分布になりますね。

今回は以上です。