crosshyou

主にクロス表(分割表)分析をしようかなと思います。

牛乳乳製品の生産動向分析3 - 主要乳製品を因子分析(factanal関数)

今回は前回に引き続き、牛乳乳製品生産動向のデータの主要乳製品の生産量を因子分析したいと思います。

前回までのR言語のコマンドを書くと以下のようになります。


# 牛乳乳製品生産動向 主要乳製品生産量
milk <- read.csv("milk.csv")
head(milk)

f:id:cross_hyou:20180815152317j:plain

 


# 生産量だけのデータフレーム
生産量 <- milk[ , c(2,4,6,8,10,12,14,16)]
head(生産量)

 

f:id:cross_hyou:20180815152337j:plain


# 生産量どうしの相関係数
round(cor(na.omit(生産量)),2)

f:id:cross_hyou:20180815152358j:plain

この相関係数マトリックスを眺めていたら、なんだか、

第1グループ:バター、脱脂粉乳

第2グループ:れん乳類、全紛乳、調整紛乳

第3グループ:チーズ、クリーム

その他:アイスクリーム

という4つのグループにわけられそうだな~~と思ったのでした。

今回は因子分析で検証してみようと思います。

因子分析はfactanal関数です。

f:id:cross_hyou:20180815153128j:plain

NAが入っていると計算できないので、na.omit(生産量)としてNAのデータを削除しています。出力するときに、小さい値のものも出力するようにprint関数でcutoff=0としています。

p-value is 0.129 > 0.05 ですから、因子数は4で十分(またはそれ以上)であることを示しています。

f:id:cross_hyou:20180815153648j:plain

この画像の一番下の列を見ると、Cumulative Var とあります。Factor1で0.544です。これはFactor1で説明される各データの分散の割合が54.4%であることを示しています。

Factor4までで97%の分散が説明できるということになります。

Loadingsの部分のマトリックスを作成してFactorごとに並び替えでみましょう。

f:id:cross_hyou:20180815154848j:plain

まずは、Factor1の大きい順に並び替えます。

f:id:cross_hyou:20180815155558j:plain

全紛乳、練乳、調整紛乳の3つが0.9以上ですね。Factor1はこれらのグループの因子という位置づけです。

Factor2で並び替えます。

f:id:cross_hyou:20180815155841j:plain

バターと脱脂粉乳が値が0.8以上ですね。Factor2はこれらのグループの因子ですね。

Factor3で並び替えます。

f:id:cross_hyou:20180815160108j:plain

アイスクリームが0.925と高いですね。Factor4はアイスクリーム因子ですね。

Factor4で並び替えます。

f:id:cross_hyou:20180815160358j:plain

チーズが一番大きく、クリームが一番小さいです。解釈が難しいですね。因子分析、まだ勉強し始めたばかりなので難しいです。

今回は、因子分析に挑戦してみました。R言語のfactanal関数で簡単にできます。

 

 

 

牛乳乳製品の生産動向分析2 - 主要乳製品の相関関係を見る(cor関数)

今回は主要乳製品の相関関係を見てみます。

まずは、csvファイルにあるデータをread.csv関数でR言語に読込みます。

f:id:cross_hyou:20180815144338j:plain

このように、各製品の生産量(トン、アイスクリームはキロリットル)と前年伸び率(パーセント)のデータです。

バター、脱脂粉乳、れん乳類、全紛乳、チーズ、調整紛乳、クリーム、アイスクリームの8種類があります。これらの相関係数をみてみましょう。

まずは、生産量だけ、伸び率だけのデータフレームを作成します。

f:id:cross_hyou:20180815144932j:plain

f:id:cross_hyou:20180815145152j:plain

うまく作成できました。

相関係数を算出するには、cor関数を使います。

f:id:cross_hyou:20180815145540j:plain

cor関数は、NAがあると計算できないので、まず、na.omit(生産量)としてやってNAのデータを削除します。それをcor関数で相関係数を計算し、round関数で小数点以下2桁までで表示しています。

バターと脱脂粉乳は0.96で高い相関です。絶対値で0.9以上の高相関の組み合わせを上げると、正の相関の組み合わせは、「バターと脱脂粉乳」「れん乳類と全紛乳」「全紛乳と調整紛乳」の3つです。逆相関の組み合わせは、「れん乳類とチーズ」「全紛乳とチーズ」「全紛乳とクリーム」「調整紛乳とクリーム」の組み合わせです。

伸び率の相関も計算してみましょう。

f:id:cross_hyou:20180815150527j:plain

伸び率どうしの相関係数は、生産量どうしの相関係数よりも低いです。

こうしてみると、8つの主要製品は、
第1グループ(バター、脱脂粉乳)

第2グループ(れん乳類、全紛乳、調整紛乳)

第3グループ(チーズ、クリーム)

その他(アイスクリーム)

という感じで分類できるような気がします。

次回の因子分析で確かめてみましょう。

 

牛乳乳製品の生産動向分析1 - 主要乳製品生産量の基本統計量(summary関数, apply関数とsd関数)

政府統計の総合窓口、e-Statのサイトに訪問をしたら、牛乳乳製品の生産動向という統計が新着でありました。

f:id:cross_hyou:20180815092725j:plain

乳製品好きとしては、チェックしないではいられません。早速ファイルを見てみます。

f:id:cross_hyou:20180815092900j:plain

本統計では、全国の生乳生産量、用途別処理量、地域別生乳生産量、主要乳製品生産量、バター、脱脂粉乳の期末在庫等を毎月提供しています。とのことです。

ファイルを開いてみましたところ、主要乳製品の生産量はこんな感じのファイルでした。

f:id:cross_hyou:20180815093147j:plain

このデータを使って、R言語の練習をしましょう。まずは、このデータを簡単にしたものをcsvファイルに保存しました。

f:id:cross_hyou:20180815094210j:plain

単位は、トンです。(アイスクリームは千キロリットル)。伸び率はパーセント表示です。

R言語のread.csv関数でデータを読込みましょう。

f:id:cross_hyou:20180815095515j:plain

head関数で始めの6行を表示します。ちなみに、head(milk, 3)のように、数字をオプションで加えると、その数字の行数分だけデータを表示します。

データの基本統計量をだしてみましょう。summary関数を使うと最小値、第1分位値、中央値、平均値、第3分位値、最大値をだしてくれます。

f:id:cross_hyou:20180815100247j:plain

年度は数値データではないので、それぞれの値の度数が表示されています。

生産量MAXなのは、脱脂粉乳ですね。201,997トンもあります。チーズが145,632トンですからそれよりも多いですね。

summary関数では標準偏差は表示されませんので、apply関数とsd関数を組み合わせて算出します。

f:id:cross_hyou:20180815100748j:plain

あ、クリームやアイスクリームは欠損値(NA)が含まれていたので標準偏差が計算されません。年度は数値データではないのでNAでもいいですが、アイスクリームとクリームの標準偏差も出したいですね。こういうときは、na.rm=TRUEを加えて計算します。

f:id:cross_hyou:20180815101213j:plain

これでクリーム、アイスクリームの標準偏差も計算できました。

せっかくなので、apply関数とその他の関数を組みわわせて統計量を計算しましょう。

apply関数とmin関数で最小値です。

f:id:cross_hyou:20180815101529j:plain

apply関数とmax関数で最大値です。

f:id:cross_hyou:20180815101806j:plain

apply関数とmedian関数で中央値です。

f:id:cross_hyou:20180815101956j:plain

apply関数とmean関数で平均値です。

f:id:cross_hyou:20180815102946j:plain

mean関数は少しテクニックが必要です。まず、年度のデータが数値ではないので、年度の列も入ったままだとNAで計算できなくなります。milk[ , c(-1)]として1列目を削除しています。そして、round関数を使って小数点以下1桁にしています。round関数を使わないと、

f:id:cross_hyou:20180815103423j:plain

このように、小数点以下6桁表示になって、e+04やe-01などがくっついた表記になって、よくわからなくなります。

 

 

鉱工業出荷内訳表の分析8 - 月によって、輸出 / 国内 の比率に違いはあるか?

前回は、鉱工業出荷内訳表のデータを使って、輸出 / 国内の比率の傾向を調べました。2015年以降は輸出のほうが多くなっていることがわかりました。

今回は、月によって、輸出 / 国内 の比率に傾向があるかどうかを調べたいと思います。

まずは、read.csv関数でデータを読込み、summary関数でデータの要約を表示します。

f:id:cross_hyou:20180815071957j:plain

こんな感じです。

次に、輸出だけ、国内だけ、のデータフレームを作成します。

f:id:cross_hyou:20180815072147j:plain

データフレーム[ , c(必要な列1, 必要な列2, ...)]というように、データフレームと[ ]で必要な列を指定して、輸出だけ、国内だけのデータフレームを作成しています。作成したデータフレームを割り算すると、輸出 / 国内 の比率になります。

f:id:cross_hyou:20180815072603j:plain

head(データフレーム, 表示行数)で、初めの2行だけを表示しています。

列名が、鉱工業.輸出のように、.輸出がついていますので、新しく列名を付け直します。

f:id:cross_hyou:20180815072808j:plain

summary関数で要約します。

f:id:cross_hyou:20180815072909j:plain

比率の最小値は、石油が記録した0.5660で、最大値は情報通の1.6620ですね。

cut関数で1未満を国内超、1より大きいを輸出超とファクタに変換します。

f:id:cross_hyou:20180815073227j:plain

こうして作成したファクタのデータ群ともとのデータフレームの年、月のデータを結合させます。cbind関数を使います。

f:id:cross_hyou:20180815073523j:plain


年と月が数値データとして取り込まれていますので、factor関数を使って、ファクタに変換します。

f:id:cross_hyou:20180815073808j:plain

このように、ファクタに変換されました。

月と各業種のクロス表を作成します。table関数を使います。

f:id:cross_hyou:20180815074722j:plain

こうして、作成したクロス表を合計します。単純に足し算するとOKです。

f:id:cross_hyou:20180815075340j:plain

1月と2月は輸出超が多い感じですが、有意な違いといえるのかな?カイ自乗検定で確かめてみましょう。

帰無仮説:月と輸出超・国内超の度数は関連がない。

対立仮設:月と輸出超・国内超の度数は関連がある。

です。

f:id:cross_hyou:20180815080045j:plain

p-vakue = 0.801 > 0.05 ですから、帰無仮説は棄却されません。

月によって、輸出超・国内超の度数に違いはないようです。

 

 

 

鉱工業出荷内訳表の分析7 - 輸出 / 国内 の比率をファクタに変換して、クロス表分析(cut関数, table関数, chisq.test関数)

今回は、鉱工業集荷内訳表のデータを加工した、輸出 / 国内 の比率を1.0より大きいなら輸出超、1.0以下なら国内超という2つのファクタに変換して、クロス表分析をしたいと思います。このブログの名前がCross_Hyouなのに最近、全然クロス表分析をしていませんでしたからね。

まずは、前回までのおさらいで、データを 輸出 / 国内 の比率にするところまでをやります。


 # 鉱工業出荷内訳表データの読込み
 内訳表 <- read.csv("shukka.csv")


# 輸出だけのデータ
輸出だけ <- 内訳表[ , c(5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50)]
 
 
# 国内だけのデータ
国内だけ <- 内訳表[ , c(6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51)]


# 輸出 / 国内
比率 <- 輸出だけ / 国内だけ
 
 
# 比率のデータフレームの列名の変更
colnames(比率) <- c("鉱工業","鉄鋼業","非鉄金","金属製","はん用","電子部","電気機","情報通", "輸送機","窯業","化学工","化学除","石油","プラス","パルプ","繊維工")

これで、「比率」という名前のデータフレームができました。

ためしに、鉱工業の比率のデータを表示してみます。

f:id:cross_hyou:20180814120525j:plain

このような数値データをファクタに変換するには、cut関数を使います。

f:id:cross_hyou:20180814121142j:plain

まず、kubun <- c(0,1,2)というコマンドで、区分を表すベクトルを指定してから、cut関数を使います。labels=で2つのファクタに名前をつけます。ちゃんと2つのファクタになっているか、確認してみます。

f:id:cross_hyou:20180814121448j:plain

はい、うまくいきました。cut関数では、~より大きく~以下、というように区切られますので、一番初めのデータの1.000000は国内超に区分されます。

この調子で残りの業種もファクタに変換します。

f:id:cross_hyou:20180814122852j:plain

こうして作成したファクタのデータを一番はじめの内訳表のデータのデータ番号、年、月と合体させます。cbind関数です。

f:id:cross_hyou:20180814124053j:plain

summary関数で表示すると、年と月は数値データなので、平均値等が表示されています。これをfactor関数でファクタに変換します。

f:id:cross_hyou:20180814124516j:plain

このように年、月の度数が表示されています。数値データがファクタに変わりました。

それでは、年 と鉱工業ファクタのクロス表を作ります。table関数です。

f:id:cross_hyou:20180814125110j:plain

おお~~2015年以降は全て、輸出超ですね。こんな感じでクロス表ができます。

全ての業種のクロス表を作ります。

f:id:cross_hyou:20180814131126j:plain

こうして、全てのクロス表を作ったら、合計します。

f:id:cross_hyou:20180814131821j:plain

2008年、2009年、2012年などは国内超が多く、2015年以降は輸出超が多いですね。では、この年によって、国内超、輸出超の比率が違うのは統計的に有意かどうかをカイ自乗検定で確認します。chisq.test関数です。

f:id:cross_hyou:20180814132212j:plain

p-value < 2.2e-16 < 0.05ですから、年と国内超、輸出超の度数は関連があることがわかります。調整済み残差を表示します。

f:id:cross_hyou:20180814132702j:plain

絶対値が1.96以上の部分は有意な違いがる部分です。こうしてみると、2010年、2013年、2014年は有意な差ではないですが、それ以外の年は国内超、輸出超のどちらかに偏っていて、それは有意な偏りであることがわかります。2015年からは輸出超が続いています。ここ数年の鉱工業生産の好調さは、輸出の恩恵ということですね。

 

鉱工業出荷内訳表の分析6 - 輸出 / 国内 比率の推移をグラフにして見る。(plot関数)

今回は、鉱工業出荷内訳表のデータを使って、前回のブログで作成した、輸出 / 国内 比率の推移をグラフにしてみようと思います。

まずは、前回までのおさらいで、R言語のコマンドを実行します。実行したコマンドは以下のようなものです。

# 鉱工業出荷内訳表データの読込みます。
内訳表 <- read.csv("shukka.csv")

# 輸出だけのデータを作成します。
輸出だけ <- 内訳表[ , c(5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50)]

# 国内だけのデータを作成します。
国内だけ <- 内訳表[ , c(6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51)]

# 輸出 / 国内 比率を作成します。
比率 <- 輸出だけ / 国内だけ

# 比率のデータフレームの列名の変更を変更します。
colnames(比率) <- c("鉱工業","鉄鋼業","非鉄金","金属製","はん用","電子部","電気機","情報通", "輸送機","窯業","化学工","化学除","石油","プラス","パルプ","繊維工")

こうして「比率」というデータフレームが作成できました。これの基本統計量を表示しましょう。summary関数で平均値、中央値、最小値、最大値、第1分位値、第3分位値を表示、sd関数で標準偏差を計算します。

f:id:cross_hyou:20180813115722j:plain

最大値を記録した業種は、情報通で1.6620でした。

最小値を記録した業種は、石油で0.5660でした。

標準偏差の最大値を記録した業種は、情報通で0.19094568でした。

標準偏差の最小値を記録した業種は、鉱工業で0.05909062でした。

今回は、情報通、石油、鉱工業の3つの業種をグラフにしましょう。plot関数を使います。

f:id:cross_hyou:20180813120510j:plain

f:id:cross_hyou:20180813120702j:plain

f:id:cross_hyou:20180813121015j:plain

f:id:cross_hyou:20180813121027j:plain

f:id:cross_hyou:20180813121335j:plain

f:id:cross_hyou:20180813121349j:plain

このようになりました。鉱工業は標準偏差が一番小さいだけあって、変動があまりないですね。

plot関数は、main="~~"でグラフのタイトル、ylim=c(xx,yy)でY軸の範囲、type="~"でグラフのタイプを指定します。l(エル)だと、線グラフになります。

この3つのグラフを重ねて作図してみましょう。

par(new=TRUE)というコマンドを使います。

f:id:cross_hyou:20180813122431j:plain

f:id:cross_hyou:20180813122445j:plain

ylab="~"でY軸のラベルを指定しています。

 

鉱工業出荷内訳表の分析5 - 輸出と国内の比率を計算し、基本統計量を計算する。(summary関数, apply関数とsd関数)

今回は、鉱工業出荷内訳表のデータから輸出と国内の比率を計算して、その比率の基本統計量を算出してみようと思います。出荷内訳、輸出と国内の内訳ですから、その比率が重要かな?と思います。

まずは、csvファイルに保存してあるデータをread.csv関数で読込み、head関数で始めの2行を表示してみます。

f:id:cross_hyou:20180811151747j:plain

f:id:cross_hyou:20180811151802j:plain

このようにデータが格納されています。輸出の列は、5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50番目の列ですね。

そして国内の列は、6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51ですね。

まずは、輸出だけ、国内だけのデータフレームを作りましょう。

f:id:cross_hyou:20180811152505j:plain

同じように、国内だけのデータフレームを作りましょう。

f:id:cross_hyou:20180811152801j:plain

輸出 / 国内 を計算します。

f:id:cross_hyou:20180811153032j:plain

ほんとうに正しく比率が計算されているかチェックしてみましょう。鉄鋼業の1行目は0.7864711です。鉄鋼業.輸出は96.5で、国内が122.7ですから、96.5 / 122.7を計算してみると、0.7864710676となりますので合っていますね。このデータフレームの列名を付け直しましょう。

f:id:cross_hyou:20180811153741j:plain

colnames属性に新しい列名を入れました。それでは、summary関数で基本統計量を算出しましょう。

f:id:cross_hyou:20180811154108j:plain

列名の名前もすっきりしましたね。比率の最小値は、石油の0.5660です。そして最大値は電気機の1,5760です。輸出 / 国内 の比率ですから、石油は輸出するよりも国内のほうが多いでしょうし、電気機は輸出が多いでしょうからこの結果は納得感があります。

標準偏差も計算しましょう。apply関数とsd関数です。

f:id:cross_hyou:20180811154639j:plain

標準偏差の一番大きいのは、情報通の0.19094568です。一番小さいのは鉱工業の0.05909062です。