PHPで相関係数と回帰直線を表示

(1/1)
PHP で人口を推計する」では、日本の将来推計人口をシミュレーションした。現実には、合計特殊出生率が、なかなか上昇せず、人口減少が加速している。
なぜ出生率は上昇しないのか――婚姻率の低下と離婚率の上昇が関係しているのではないかと考え、それを確かめるため、高校数学で習った相関係数を求めてみることにする。
あわせて、データ集合の散布図回帰曲線を描いてみる。

サンプル・プログラムの実行例

PHPで相関係数と回帰直線を表示
当初の予想通り、婚姻率と出生率の間には強い正の相関関係が、離婚率と出生率の間には負の相関関係があることが検証できた。

ここで注意してほしいことがある――婚姻率や離婚率と出生率の間に相関関係があるからといっても、それは、婚姻は離婚が原因で出生率が低下しているという結果(因果関係になっているわけではない。つまり、相関関係と因果関係は別物であるということだ。

2 つのデータ集合 X, Y があり、X が変化すれば Y も変化することを相関関係という。
これに対し、X(原因)→Y(結果)という一方通行の関係が因果関係である。相関関係は因果関係の必要条件になり得るが、十分条件ではないのである。

サンプル・プログラム

解説:準備

0021: //jqPlotのあるフォルダ
0022: define('JQPLOT', '../../../../common/jqplot/');
0023: 
0024: //元データ
0025: define('FILE_CSV', 'birthrate.csv');
0026: 
0027: //CSVデリミタ
0028: define('CSV_DELIMITER', ',');
0029: 
0030: //グラフの幅、高さ
0031: define('GRAPH_WIDTH',  600);
0032: define('GRAPH_HEIGHT', 350);
0033: 
0034: //婚姻率、離婚率、合計特殊出生率 推移グラフの名前
0035: define('GRAPH_PROGRESS', 'jqPlot_progress');
0036: 
0037: //合計特殊出生率=婚姻率 相関グラフの名前
0038: define('GRAPH_MARRIAGE', 'jqPlot_marriage');
0039: 
0040: //合計特殊出生率=離婚率 相関グラフの名前
0041: define('GRAPH_DIVORCE', 'jqPlot_divorce');
0042: 
0043: //合計特殊出生率=初婚年齢 相関グラフの名前
0044: define('GRAPH_1STMARRIAGE', 'jqPlot_1stmarriage');
0045: 
0046: //統計に関わるクラス:include_pathが通ったディレクトリに配置
0047: require_once('pahooStat.php');

年次推移や散布図を描くのに、jQuery プラグイン「jqPlot」を利用する。
jqPlot の入手方法、インストール方法については、「PHP で NHK 政治意識月例調査をグラフ表示」を参照してほしい。

相関係数PECL の統計ライブラリを導入すれば計算できるが、今回は、自前のクラスを用意することにした。
クラスについては「PHP でクラスを使ってテキストの読みやすさを調べる」で解説している。
クラスファイル "pahooStat.php" は、組み込み関数  require_once  を使って読めるディレクトリに配置する。ディレクトリは、設定ファイル php.ini に記述されているオプション設定 include_path に設定しておく。
pahooStat クラスの内容については後述する。

毎年の婚姻率、離婚率、合計特殊出生率は、厚生労働省「平成 29 年度 人口動態統計の年間推計」を Excel に入力し、CSV ファイル FILE_CSV として保存してある。また、平均初婚年齢(妻)は、政府統計の総合窓口から「人口動態調査 年次別平均婚姻年齢及び夫妻の年齢差」を利用した。
これらのデータを参照し、各種グラフを描いてゆく。

解説:データの検査

0061: /**
0062:  * 配列かどうか、要素が1つ以上あるかどうか検査する
0063:  * @param array $x データ配列
0064:  * @return bool TRUE:配列である/FALSE:配列でない
0065: */
0066: function is_data($x) {
0067:     $res = TRUE;
0068:     if (! is_array($x)) {
0069:         $this->error = TRUE;
0070:         $this->errmsg = '配列でない';
0071:         $res = FALSE;
0072:     } else if (count($x) < 1) {
0073:         $this->error = TRUE;
0074:         $this->errmsg = '配列要素がない';
0075:         $res = FALSE;
0076:     }
0077:     return $res;
0078: }

クラス pahooStat では、配列を計算することが多い。
そこで、引数として渡された変数が配列であり、かつ要素が 1 つ以上あることを検査するメソッド is_data を用意した。

解説:相加平均

0080: /**
0081:  * 相加平均
0082:  * @param array $x データ配列
0083:  * @return double 平均/NULL
0084: */
0085: function arithmetic_mean($x) {
0086:     if (! $this->is_data($x))     return NULL;
0087: 
0088:     return array_sum($x) / count($x);
0089: }

日常、よく使う平均は「相加平均(算術平均)」と呼ばれ、次の計算式によって導かれる。
Excel の AVERAGE関数に相当する。
 mimetex 

解説:分散

0091: /**
0092:  * 分散
0093:  * @param array $x データ配列
0094:  * @return double 分散
0095: */
0096: function variance($x) {
0097:     if (! $this->is_data($x))     return NULL;
0098: 
0099:     $xm = $this->arithmetic_mean($x);
0100:     $xs = 0;
0101:     $n = count($x);
0102:     $xs = 0;
0103:     for ($i = 0; $i < $n$i++) {
0104:         $xs += ($x[$i] - $xm) * ($x[$i] - $xm);
0105:     }
0106: 
0107:     return $xs / $n;
0108: }

分散は、データの散らばり度合いを表す値で、偏差(各々の値と相加平均の差)を二乗し、その相加平均を求める。
Excel の VARP関数に相当する。
 mimetex 

解説:標準偏差

0130: /**
0131:  * 標準偏差
0132:  * @param array $x データ配列
0133:  * @return double 標準偏差
0134: */
0135: function standard_deviation($x) {
0136:     if (! $this->is_data($x))     return NULL;
0137: 
0138:     return sqrt($this->variance($x));
0139: }

標準偏差は、データの散らばり度合いを表す値で、分散の平方根を求める。
Excel の STDEV.P関数に相当する。
 mimetex 

解説:共分散

0110: /**
0111:  * 共分散
0112:  * @param array $x, $y データ配列
0113:  * @return double 共分散
0114: */
0115: function covariance($x$y) {
0116:     if (! $this->is_data($x))     return NULL;
0117:     if (! $this->is_data($y))     return NULL;
0118: 
0119:     $xm = $this->arithmetic_mean($x);
0120:     $ym = $this->arithmetic_mean($y);
0121:     $n = min(count($x), count($y));
0122:     $xs = 0;
0123:     for ($i = 0; $i < $n$i++) {
0124:         $xs += ($x[$i] - $xm) * ($y[$i] - $ym);
0125:     }
0126: 
0127:     return $xs / $n;
0128: }

共分散は、2種類のデータ集合の関係を示す指標で、2 つのデータ集合の偏差の積の平均を求める。
Excel の COVAR関数に相当する。
 mimetex 

解説:相関係数

0141: /**
0142:  * 相関係数
0143:  * @param array $x, $y データ配列
0144:  * @return double 相関係数
0145: */
0146: function correl($x$y) {
0147:     if (! $this->is_data($x))     return NULL;
0148:     if (! $this->is_data($y))     return NULL;
0149: 
0150:     return $this->covariance($x$y) / ($this->standard_deviation($x) * $this->standard_deviation($y));
0151: }

相関係数は、2種類のデータ集合の関係を示す指標で、単位の影響を受けずに関係を示すことができる。
Excel の CORREL関数に相当する。
 mimetex 

解説:相関の強さ

0153: /**
0154:  * 相関の強さ
0155:  * @param double $r 相関係数
0156:  * @return string 相関の強さ(テキスト)
0157: */
0158: function correl_strength($r) {
0159:     $table = array(
0160: array(-0.7, '強い負の相関'),
0161: array(-0.4, '負の相関'),
0162: array(-0.2, '弱い負の相関'),
0163: array(+0.2, 'ほとんど相関がない'),
0164: array(+0.4, '弱い正の相関'),
0165: array(+0.7, '正の相関'),
0166: array(+1.1, '強い正の相関')
0167: );
0168: 
0169:     foreach ($table as $item) {
0170:         if ($r < $item[0])      return $item[1];
0171:     }
0172:     return FALSE;
0173: }

相関係数の値と、相関の強さの基準は下表の通りである。
相関係数 r 相関の強さ
-1.0 ≦ r < -0.7強い負の相関
-0.7 ≦ r < -0.4負の相関
-0.4 ≦ r < -0.2弱い負の相関
-0.2 ≦ r < +0.2ほとんど相関がない
+0.2 ≦ r < +0.4弱い正の相関
+0.4 ≦ r < +0.7正の相関
+0.7 ≦ r ≦ +1.0強い正の相関

解説:回帰直線

0175: /**
0176:  * 回帰係数、切片
0177:  * @param array $x, $y データ配列
0178:  * @return array(回帰係数,切片)
0179: */
0180: function slope($x$y) {
0181:     if (! $this->is_data($x))     return NULL;
0182:     if (! $this->is_data($y))     return NULL;
0183: 
0184:     $a = $this->covariance($x$y) / $this->variance($x);
0185:     $b = $this->arithmetic_mean($y) - $a * $this->arithmetic_mean($x);
0186: 
0187:     return array($a$b);
0188: }

2種類のデータ集合を散布図にプロットしたとき、相関があれば、直線上にデータが並ぶ。この直線のことを回帰直線と呼び、最小自乗法(最小二乗法)によって求めることができる。
Excel の SLOPE関数INTERCEPT関数に相当する。

回帰直線を  mimetex  とするとき、a, b は次の式によって求める。
 mimetex 

 mimetex 

解説:散布図

0245: /**
0246:  * jqPlot用のスクリプト:散布図
0247:  * @param array $x, $y データ配列
0248:  * @param string $label_x, $label_y 軸ラベル
0249:  * @param string $name オブジェクト名
0250:  * @return string スクリプト/FALSE
0251: */
0252: function plotScatter($x$y$label_x$label_y$name) {
0253:     //相関係数、回帰直線
0254:     $pst = new pahooStat();
0255:     $r = $pst->correl($x$y);                //相関係数
0256:     if ($pst->iserror()) return FALSE;
0257:     $st = $pst->correl_strength($r);
0258:     $r = sprintf('%4.3f (%s)', $r$st);
0259:     list($a$b) = $pst->slope($x$y);       //回帰係数、切片
0260:     if ($pst->iserror()) return FALSE;
0261:     $x0 = min($x);
0262:     $y0 = $a * $x0 + $b;
0263:     $x1 = max($x);
0264:     $y1 = $a * $x1 + $b;
0265:     $pst = NULL;
0266: 
0267:     //系列の生成
0268:     $series = '';
0269:     $n = min(count($x), count($y));
0270:     for ($i = 0; $i < $n$i++) {
0271:         $series .= sprintf('[%f, %f], ', $x[$i]$y[$i]);
0272:     }
0273: 
0274:     $width  = GRAPH_WIDTH;
0275:     $height = GRAPH_HEIGHT;
0276: $js =<<< EOD
0277: <script>
0278: jQuery(function() {
0279:     jQuery.jqplot('{$name}',
0280:     [
0281:         [ {$series} ]
0282:     ],
0283:     {
0284:         //タイトル
0285:         title: {
0286:             text: '{$label_x}={$label_y}<br /><span style="font-size:small; color:blue;">相関係数 = {$r}</span>',
0287:             show: true,
0288:             fontSize: '20px',
0289:             textAlign: 'center',
0290:             textColor: 'black'
0291:         },
0292:         //散布図
0293:         seriesDefaults: {
0294:             showLine: false,
0295:             color: 'green',
0296:             markerOptions: { size: 4 }
0297:         },
0298:         //回帰直線
0299:         canvasOverlay: {
0300:             show: true,
0301:             objects: [
0302:             {
0303:                 line: {
0304:                     lineWidth: 3,
0305:                     color: 'blue',
0306:                     shadow: false,
0307:                     lineCap: 'round',
0308:                     start: [ {$x0}, {$y0} ],
0309:                     stop:  [ {$x1}, {$y1} ]
0310:                 }
0311:             }]
0312:         },
0313:         //軸
0314:         axes: {
0315:             xaxis: {
0316:                 label: '{$label_x}',
0317:                 tickOptions: { formatString: '%.1f' }
0318:             },
0319:             yaxis: {
0320:                 label: '{$label_y}',
0321:                 tickOptions: { formatString: '%.1f' }
0322:             }
0323:         }
0324:     }
0325:     );
0326: });
0327: </script>
0328: <div id="{$name}" style="width:{$width}px; height:{$height}px;"></div>
0329: 
0330: EOD;
0331: 
0332:     return $js;
0333: }

配列に格納したデータを、jqPlot スクリプトへ展開し散布図とする処理は、ユーザー関数 plotScatter で行う。

回帰直線は、X軸の最小値を始点とし、同じく最大値を終点とする直線で、jqPlot プラグインの「jqplot.dateAxisRenderer.js」(グラフに描き加える線)で行う。

解説:年次推移グラフ

0144: /**
0145:  * jqPlot用のスクリプト:年次推移グラフ
0146:  * @param array  $items データ配列
0147:  * @param array  $labels ラベル配列
0148:  * @param string $name オブジェクト名
0149:  * @param string $title グラフ・タイトル
0150:  * @return string スクリプト
0151: */
0152: function plotProgress($items$labels$name$title) {
0153:     //系列の生成
0154:     $series = array();
0155:     $n = count($labels);
0156:     for ($i = 1; $i < $n$i++) {
0157:         $series[$i] = '';
0158:         foreach ($items[$i] as $key=>$val) {
0159:             $series[$i] .= sprintf('[%d, %f], ', $items[0][$key]$val);
0160:         }
0161:     }
0162:     $s1 = '';
0163:     $s2 = '';
0164:     foreach ($series as $key=>$val) {
0165:         $s1 .= "[ {$val} ], ";
0166:         if (max($items[$key]) > 20) {
0167: $s2 .=<<< EOD
0168:                 {
0169:                 label: '{$labels[$key]}',
0170:                 yaxis: 'y2axis'
0171:                 },
0172: 
0173: EOD;
0174:         } else {
0175: $s2 .=<<< EOD
0176:                 {
0177:                 label: '{$labels[$key]}'
0178:                 },
0179: 
0180: EOD;
0181:         }
0182:     }
0183: 
0184:     $width  = GRAPH_WIDTH;
0185:     $height = GRAPH_HEIGHT;
0186: $js =<<< EOD
0187: <script>
0188: jQuery(function() {
0189:     jQuery.jqplot('{$name}',
0190:     [
0191:         {$s1}
0192:     ],
0193:     {
0194:         //タイトル
0195:         title: {
0196:             text: '{$title}',
0197:             show: true,
0198:             fontSize: '20px',
0199:             textAlign: 'center',
0200:             textColor: 'black'
0201:         },
0202:         //ラベル
0203:         series: [
0204:             {$s2}
0205:         ],
0206:         legend: {
0207:             show: true,
0208:             placement: 'inside',
0209:             location: 'ne',
0210:             renderer: $.jqplot.EnhancedLegendRenderer,
0211:             rendererOptions: { numberRows: 1 }
0212:         },
0213:         //グラフ
0214:         seriesDefaults: {
0215:             showLine: true,
0216:             rendererOptions: { smooth: false },
0217:             markerOptions: { size: 0 },
0218:         },
0219:         //軸
0220:         axes: {
0221:             xaxis: {
0222:                 label: '{$labels[0]}',
0223:                 tickOptions: { formatString: '%d' }
0224:             },
0225:             yaxis: {
0226:                 label: '率',
0227:                 tickOptions: { formatString: '%.1f' }
0228:             },
0229:             y2axis: {
0230:                 label: '年齢',
0231:                 tickOptions: { formatString: '%d' }
0232:             }
0233:         }
0234:     }
0235:     );
0236: });
0237: </script>
0238: <div id="{$name}" style="width:{$width}px; height:{$height}px;"></div>
0239: 
0240: EOD;
0241: 
0242:     return $js;
0243: }

配列に格納したデータから、合計特殊出生率、婚姻率、離婚率、平均初婚年齢(妻)の年次推移を jqPlot スクリプトへ展開し、折れ線グラフとする処理は、ユーザー関数 plotProgress で行う。

参考サイト

(この項おわり)
header