鉱工業出荷内訳表の分析４ - ２つのデータの分布の違いを検定する。2標本のK-S検定(ks.test関数)

今回は、鉱工業出荷内訳表のデータを使って、2つのデータの分布の違いを検定してみたいと思います。具体的には、前回、前々回でも取り上げた、平均値が最大のデータのはん用.国内と平均値が最小のデータの情報通.国内を比較します。もう一組は、標準偏差が最大のデータの情報通.輸出と標準偏差が最小の化学工.出荷を比較します。

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

4つのデータの平均値と標準偏差を一度に計算したいので、4つのデータだけのマトリックスを作成します。matrix関数です。

matrix関数でデータのマトリックスを作成し、colnamesでそれぞれの列に名前を付けています。

summary関数でデータの要約を表示します。

はん用.国内の平均値は111.1です。情報通.国内の平均値は66.0です。

apply関数とsd関数で、標準偏差を計算します。

あ！、情報通.国内のほうが情報通.出荷より大きな標準偏差ですね。すみません、間違えてしまいました。このブログの趣旨はR言語の操作と簡単な統計分析の練習なので、とりあえずこのまま進めます。

分布を比較する2つのデータのヒストグラムを作成してみましょう。hist関数ですね。同じ画面に2つのヒストグラムを作図したいので、par(mfrow=c(2,1))というコマンドを加えます。

あきらかに分布が違いますね。col=で色を指定、main=でタイトルを指定、breaks=でヒストグラムの区間を指定しています。

同じように、情報通.出荷と化学工.出荷もヒストグラムにしてみます。

こちらもあきらかに分布の形状は違いますね。

このような、2つのデータの分布の違いを検定するには、2標本のK-S検定をします。

帰無仮説H0は、データの分布は、はん用.国内と情報通.国内で差がない。

対立仮設H1は、データの分布は、はん用.国内と情報通.国内で異なる。

検定の有意水準はα=0.05とします。

ks.test関数を実行します。

p-value < 2.2e-16 < 0.05 なので、帰無仮説は棄却されます。はん用.国内と情報通.国内の分布の形状は異なります。

次は、情報通.出荷と化学工.出荷の分布の形状の検定です。

帰無仮説は、情報通.出荷と化学工.出荷のデータの分布は差がない。

対立仮設は、情報通.出荷と化学工.出荷のデータの分布は異なる。

検定の推移水準は、α=0.05 とします。

ks.test関数を実行します。

p-value < 2.2e-16 < 0.05 ですから帰無仮説は棄却されます。情報通.出荷と化学工.出荷の分布は異なります。

せっかくですので、情報通.国内と情報通.出荷のペアで検定してみましょう。

まずは、ヒストグラムを並べてみましょう。

似ていますね。それでは、検定しましょう。

帰無仮説は、情報通.国内と情報通.出荷のデータの分布には差がない。

対立仮設は、情報通.国内と情報通.出荷のデータの分布は異なる。

検定の有意水準は、α=0.05 とします。

ks.test関数でK-S検定をします。

p-value = 0.1497 > 0.05 ですから、帰無仮説を棄却できません。つまりこの2つのデータの分布に差があるとはいえません。

鉱工業出荷内訳表の分析３ - データのバラツキをヒストグラムにして見える化する。(hist関数)

今回は、鉱工業出荷内訳表のデータのヒストグラムを作成しようと思います。

48のデータの種類全部をヒストグラムにするのは骨が折れるので、前回と同じく平均値が一番大きいデータのはん用.国内、平均値が一番小さいデータの情報通.国内、標準偏差の一番大きいデータの情報通.輸出、標準偏差の一番小さいデータの化学工.出荷の4つのデータについてヒストグラムを作成しようと思います。

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

各データの平均値を計算して、平均値の大きい順に並び替えましょう。

apply関数とmean関数を使って各データの平均値を一括して計算します。そして、order関数で並び替えます。

はん用.国内が111.0736で一番大きく、情報通.国内が66.0024で一番小さいですね。

apply関数とsd関数で各データの標準偏差を一括して計算します。そして、order関数で並び替えます。

情報通.輸出の標準偏差が22.913790で一番大きく、化学工.出荷の標準偏差が4.024985で一番小さいことがわかります。

それでは、ヒストグラムを作成しましょう。hist関数を使います。

breaks = c(30,40,50,60,70,80,90,100,110,120,130,140,150,160)というオプション引数がヒストグラムの区間の区切りを指定しています。

次は、情報通.国内のヒストグラムです。

はん用.国内と違って、左に大きく偏っていますね。

次は、情報通.輸出です。

標準偏差が一番大きい、情報通.輸出ですから、幅が広いのがわかります。

標準偏差が一番小さい、化学工.出荷はこうなります。

化学工.出荷のデータは80から110の範囲しかないですね。

こうしてヒストグラムにすると、各データのバラツキがよくわかります。

最後に、4つのヒストグラムを1つの画面に配置してみましょう。

par(mfrow=c(4,1))という関数を作動させてから、4つのヒストグラムを作成します。

par(mfrow=c(4,1))という関数はグラフを、4行1列の配置で1つの画面に、という関数です。そして、col="red"などで色を付けました。

鉱工業出荷内訳表の分析２ - データの推移をグラフにして観る。(plot関数)

今回は、鉱工業出荷内訳表のデータをグラフにしようと思います。とはいっても、データの種類が48もあるので、前回の分析で判明した、平均値が一番大きかった、はん用.国内、平均値が一番小さかった、情報通.国内、標準偏差が一番大きかった、情報通.輸出、標準偏差が一番小さかった、化学工.出荷の4つのデータ系列についてグラフを作成しようと思います。

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

それではグラフを作成します。plot関数を使います。type = "l" とオプション引数を追加してライングラフにします。

まずは、平均値が一番大きいデータ、はん用.国内です。

ylim=c(0,150)というオプション引数でY軸の範囲を0から150に設定しています。はじめのほうに急落して、だんだんと回復という動きです。

次は、平均値が一番小さいデータ、情報通.国内です。

はじめのころに大きく落ち込み、大きく回復した後、減少傾向が続いています。

次は、標準偏差が一番大きいデータ、情報通.輸出です。

はじめのころに大きく落ち込んだ後、回復しないで減少傾向です。

次は、標準偏差が一番小さいデータ、化学工.出荷です。

標準偏差が一番小さいだけあって、動きが今までの３つとは違ってかなり小さいですね。

それでは最後に、この4つのデータを一つの画面にグラフ化します。

まず、はん用.国内のグラフを表示します。このとき、main=""でグラフのタイトルを空白に、col="black"でグラフの色を黒に、ylab=""でY軸のラベルを空白にしています。

そして、そのあとに、par(new=TRUE)でグラフを重ねるように指示しています。

このように、データによってはグラフの形状がかなり違うことがわかります。

鉱工業出荷内訳表の分析１ - 基本統計値の算出(summary関数、apply関数とsd関数\max関数\min関数\mean関数\median関数)

今回は、鉱工業出荷内訳表のデータを分析してみましょう。e-Stat(政府統計の総合窓口)のサイトを訪問したら、今日の新着データで、「鉱工業出荷内訳表」というデータがありました。

鉱工業生産は聞いたことがありましたが、鉱工業出荷内訳表というのは、お恥ずかしながら初耳でした。どんなデータなんでしょうか？クリックしてみました。

鉱工業出荷内訳表は、国内の鉱工業製品の出荷が国内と海外輸出のどちらに向けられたのかを表すもので、鉱工業出荷指数と財務省の貿易統計（輸出）を元に、業種別、財別の国内向け出荷指数、輸出向け出荷指数を作成しています。とのことです。工業製品が国内向けか、国外向けかということを表すデータのようです。クリックしてみます。

この平成22年(2010年)基準というのをクリックしてみます。

ファイルが２つあるようです。月次のExcelファイルをダウンロードしてみました。

各業種別に出荷、輸出、国内の３つの区分があって、2008年01月からの月次のデータです。このファイルを分析しやすいように、下のように加工しました。

年と月を別々にして、データ番号として通し番号を加えました。

read.csv関数でR言語に読込み、summary関数でデータを要約してみます。

データ番号の最大値が、125ですから、データの行数は125だとわかります。1年が12か月で12行ですから、約10年半分のデータです。

データ項目の数は何個でしょうか？dim関数で表示します。

125が行数で、51が列数です。データ番号、年、月の3つを引いて、48列がデータの数です。1つの業種に出荷、輸出、国内の3つのデータがありますから、業種の数は、48/3=16です。

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

apply関数は、apply(データフレーム, 2, 計算する関数) で列ごとのデータ計算ができます。2を1にすると、行ごとのデータ計算になります。

48もデータがあると、データを一目で把握するのは難しいですね。

とりあえず、各データの最大値をだして、大きい順で並べてみましょう。

情報通信.輸出が153.8で一番ですね。

データの最小値をだして、小さい順に並べてみます。

情報通信の国内が39.7で最小です。

標準偏差をだして大きい順に並べてみます。

情報通信の輸出が一番標準偏差が大きいです。

平均値を出してみましょう。

標準偏差の大きな情報通信の平均値が低いです。情報通信のデータの変動がすごく大きいということですね。

最後にデータの中央値を計算しましょう。

はん用.国内が平均値と同じく、一番大きいですね。

花き産業振興総合調査の分析５ - 花き産業は衰退傾向と言えるのか？回帰分析で確認(lm関数)

今回は、花き産業振興総合調査のデータで、花き産業は衰退しているのかを見てみたいと思います。

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

出荷額の推移を年ごとにみてみましょう。tapply関数とsum関数で年ごとの出荷額の合計値を出します。

2009年の出荷額合計が93億6900万円で、2016年が64億4700万円です。barplot関数でグラフにしてみます。

全体としては減少傾向ですが、2013年と2014年は増加しています。

この出荷額の傾向を回帰分析で分析してみたいと思います。lm関数を使います。

はじめに、年 <- c(2009:2016)で2009年から2016年のベクトルを作っています。

次に、回帰分析結果 <- lm(出荷額合計 ~ 年) という式で回帰分析を実行しています。これは、

出荷額合計 = a * 年 + 誤差項

とう回帰式の回帰分析で、その分析結果を「回帰分析結果」という名前で格納しています。そして、sumary(回帰分析結果)で分析結果の要約を表示しています。

この要約のEstimateが回帰式の係数です。つまり、

出荷額合計 = -346.61 * 年 + 705598.25

ということです。1年経過するたびに、出荷額が3億4661万円減少するっていうことです。やっぱり、衰退傾向なんですね。

一番下の p-value = 0.005943 が、この回帰式が意味があるかどうかを表しています。p-value = 0.005943 < 0.05 ですから有意です。

そして、

この部分の赤線のところが、年の係数、-346.61が有意かどうかを表していて、0.00594 < 0.05 ですから有意とわかります。全体としての花き産業の出荷額は、年々減少していると言えます。

また、

Multiple R-squared: 0.7426 とあります。これがいわゆる決定係数(R2)ですね。出荷額と年の相関係数の2乗が、0.7426ということです。

そういえば、出荷額が増加していましたよね。モミジの出荷額で同じように回帰分析してみましょう。

まずは、モミジだけのデータフレームを作ります。

出荷額の棒グラフを作成してみます。

barplot関数が棒グラフを作成する関数で、main="~~"でグラフのタイトルｗ設定し、ylim=c(0,1000)でY軸の範囲を0から1000にして、names.arg=c(2009:2016)でX軸のラベルを2009から2016にしています。

減少傾向ではないですね。lm関数で回帰分析します。

推定された回帰式は、

モミジの出荷額 = 13.27 * 年 - 25932.42

です。1年経過すると出荷額が1327万円増える！と思いますが、一番下のp-valueを見てください。p-value = 0.2586 > 0.05 ですから、この回帰式は有意ではない、ということです。モミジの出荷額が増加傾向にある、とは言えないです。

今回は、lm関数で出荷額を年で回帰分析しました。その結果、出荷額は年々減少していることが統計的に有意だとわかりました。

花き産業振興総合調査の分析４ - 各データの散布図(plot関数)と相関を見る(cor関数)

今回は、花き産業振興総合調査のデータを相関関係から見てみたいと思います。

まずは、csvファイルに保存してあるデータをread.csv関数で読み込みます。

相関係数を求める関数は、cor関数です。cor関数はデータフレームが数値データだけだと、そのデータフレームの相関マトリックスを表示してくれます。この「花きデータ」は、1列目は種類、2列目は年なので、1列目と2列目を除外してcor関数を使います。

cor(花きデータ[ , c(-1,-2)])と入力します。

一番相関の強い組合せは、作付面積と出荷数量です。一番相関の弱い組合せは、栽培農家数と出荷額です。

plot関数で散布図を作成します。plot関数はデータフレームが数値データだけだと、散布図のマトリックスを作成します。

こうして散布図を見ると、栽培農家数と出荷額が一番バラバラになっているとわかりますね。そして、それぞれの散布図のデータですが、いくつかのカタマリになっているように見えます。これは、花きの種類ごとにカタマリになっていると推測できます。

種類ごとの相関マトリックスと散布図マトリックスを作成しましょう。まずは、それぞれの種類ごとのデータフレームを作成しました。

サツキの相関マトリックスはこうなりました。

全体では、0.38と低い値だった栽培農家数と出荷額の相関係数が0.84と高い値になりました。散布図マトリックスはこうなりました。

つぎは、ツツジです。

つぎは、カイヅカイブキです。

カイヅカイブキは相関係数が高いですね。

お次は、タマイブキです。

ツバキはこうです。

次はモミジです。

モミジは前回の分析でもそうでしたが、他の種類の花きとは傾向が違いますね。

出荷額が他のデータを逆相関になっています。

ヒバ類はどうでしょうか？

最後は、ツゲ類です。

ツゲ類の栽培農家数と出荷額もマイナスの相関係数でした。

花き産業振興総合調査の分析３ - 種類ごとの経年変化を見る。棒グラフで見える化(barplot関数)

今回は、ツツジ、サツキなど種類ごとの経年変化を見たいと思います。

まずは、csvファイルに保存してあるデータをR言語のread.csv関数で読込み、summary関数でデータの要約をします。

サツキだけのデータフレーム、ツツジだけのデータフレーム、など種類ごとのデータフレームを作りましょう。

サツキのデータフレームです。dfサツキと名付けました。

ツツジのデータフレームです。dfツツジと名付けました。

カイヅカイブキのデータフレームです。dfカイヅカイブキと名付けました。

タマイブキだけのデータフレームです。dfタマイブキと名付けました。

ツバキだけのデータフレームです。dfツバキと名付けました。

モミジだけのデータフレームです。dfモミジと名付けました。

ヒバ類だけのデータフレームです。dfヒバ類と名付けました。

ツゲ類だけのデータフレームです。dfツゲ類と名付けました。

こうして、種類ごとのデータを見ると、どの種類も作付面積、出荷数量、出荷額、栽培農家数が低下傾向にあるようです。例外は、モミジの出荷額ぐらいですね。

出荷額のデータを棒グラフでみてみます。barplot関数ですね。

main="サツキの出荷額" でグラフのタイトルを指定し、

ylim=c(0,2000) でY軸の範囲を、0から2000に指定し、

names.arg=c(2009:2106) でそれぞれのバーの名前を2009から2016に指定しています。

同じように他のものも棒グラフ化します。

ツツジは2013年、2014年と出荷額は増えていましたね。

モミジの出荷額の推移は、他の種類とは形が違いますね。

ヒバ類は2014年までは、横ばいという感じでしたが、2015年、2016年と減少しています。

ツゲ類の出荷額は、2013年と2014年に急増していますが、15年、16年に減少しています。