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

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

(2022年4月17日)統計表を更新,PHP8対応,リファラ・チェック改良。

目次

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

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

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

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

サンプル・プログラム

圧縮ファイルの内容
correlation.phpサンプル・プログラム本体
birthrate.csv統計表
pahooStat.php統計に関わるクラス pahooStat。
使い方は「PHPで相関係数と回帰直線を表示」などを参照。include_path が通ったディレクトリに配置すること。

解説:準備

0035: //jqPlotのあるフォルダ
0036: define('JQPLOT', '../../../../common/jqplot/');
0037: 
0038: //統計表
0039: define('FILE_CSV', 'birthrate.csv');
0040: 
0041: //CSVデリミタ
0042: define('CSV_DELIMITER', ',');
0043: 
0044: //グラフの幅、高さ
0045: define('GRAPH_WIDTH',  600);
0046: define('GRAPH_HEIGHT', 350);
0047: 
0048: //婚姻率、離婚率、合計特殊出生率 推移グラフの名前
0049: define('GRAPH_PROGRESS', 'jqPlot_progress');
0050: 
0051: //合計特殊出生率=婚姻率 相関グラフの名前
0052: define('GRAPH_MARRIAGE', 'jqPlot_marriage');
0053: 
0054: //合計特殊出生率=離婚率 相関グラフの名前
0055: define('GRAPH_DIVORCE', 'jqPlot_divorce');
0056: 
0057: //合計特殊出生率=初婚年齢 相関グラフの名前
0058: define('GRAPH_1STMARRIAGE', 'jqPlot_1stmarriage');
0059: 
0060: //統計に関わるクラス:include_pathが通ったディレクトリに配置
0061: require_once('pahooStat.php');

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

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

毎年の婚姻率、離婚率、合計特殊出生率は、厚生労働省「令和2年 人口動態統計の年間推計について」をExcelに入力し、統計表ファイル FILE_CSV(CSV形式ファイル)として保存してある。また、平均初婚年齢(妻)は、厚生労働省「令和2年 人口動態統計月報年計(概数)の概況」を利用した。
これらのデータを参照し、各種グラフを描いてゆく。

解説:データの検査

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 

解説:散布図

0298: /**
0299:  * jqPlot用のスクリプト:散布図
0300:  * @param   array $x, $y データ配列
0301:  * @param   string $label_x, $label_y 軸ラベル
0302:  * @param   string $name オブジェクト名
0303:  * @return  string スクリプト/FALSE
0304: */
0305: function plotScatter($x$y$label_x$label_y$name) {
0306:     //相関係数、回帰直線
0307:     $pst = new pahooStat();
0308:     $r = $pst->correl($x$y);                //相関係数
0309:     if ($pst->iserror()) return FALSE;
0310:     $st = $pst->correl_strength($r);
0311:     $r = sprintf('%4.3f (%s)', $r$st);
0312:     list($a$b) = $pst->slope($x$y);       //回帰係数、切片
0313:     if ($pst->iserror()) return FALSE;
0314:     $x0 = min($x);
0315:     $y0 = $a * $x0 + $b;
0316:     $x1 = max($x);
0317:     $y1 = $a * $x1 + $b;
0318:     $pst = NULL;
0319: 
0320:     //系列の生成
0321:     $series = '';
0322:     $n = min(count($x), count($y));
0323:     for ($i = 0; $i < $n$i++) {
0324:         $series .= sprintf('[%f, %f], ', $x[$i]$y[$i]);
0325:     }
0326: 
0327:     $width  = GRAPH_WIDTH;
0328:     $height = GRAPH_HEIGHT;
0329: $js =<<< EOT
0330: <script>
0331: jQuery(function() {
0332:     jQuery.jqplot('{$name}',
0333:     [
0334:         [ {$series} ]
0335:     ],
0336:     {
0337:         //タイトル
0338:         title: {
0339:             text: '{$label_x}={$label_y}<br /><span style="font-size:small; color:blue;">相関係数 = {$r}</span>',
0340:             show: true,
0341:             fontSize: '20px',
0342:             textAlign: 'center',
0343:             textColor: 'black'
0344:         },
0345:         //散布図
0346:         seriesDefaults: {
0347:             showLine: false,
0348:             color: 'green',
0349:             markerOptions: { size: 4 }
0350:         },
0351:         //回帰直線
0352:         canvasOverlay: {
0353:             show: true,
0354:             objects: [
0355:             {
0356:                 line: {
0357:                     lineWidth: 3,
0358:                     color: 'blue',
0359:                     shadow: false,
0360:                     lineCap: 'round',
0361:                     start: [ {$x0}, {$y0} ],
0362:                     stop:  [ {$x1}, {$y1} ]
0363:                 }
0364:             }]
0365:         },
0366:         //軸
0367:         axes: {
0368:             xaxis: {
0369:                 label: '{$label_x}',
0370:                 tickOptions: { formatString: '%.1f' }
0371:             },
0372:             yaxis: {
0373:                 label: '{$label_y}',
0374:                 tickOptions: { formatString: '%.1f' }
0375:             }
0376:         }
0377:     }
0378:     );
0379: });
0380: </script>
0381: <div id="{$name}" style="width:{$width}px; height:{$height}px;"></div>
0382: 
0383: EOT;
0384: 
0385:     return $js;
0386: }

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

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

解説:年次推移グラフ

0197: /**
0198:  * jqPlot用のスクリプト:年次推移グラフ
0199:  * @param   array  $items データ配列
0200:  * @param   array  $labels ラベル配列
0201:  * @param   string $name オブジェクト名
0202:  * @param   string $title グラフ・タイトル
0203:  * @return  string スクリプト
0204: */
0205: function plotProgress($items$labels$name$title) {
0206:     //系列の生成
0207:     $series = array();
0208:     $n = count($labels);
0209:     for ($i = 1; $i < $n$i++) {
0210:         $series[$i] = '';
0211:         foreach ($items[$i] as $key=>$val) {
0212:             $series[$i] .= sprintf('[%d, %f], ', $items[0][$key]$val);
0213:         }
0214:     }
0215:     $s1 = '';
0216:     $s2 = '';
0217:     foreach ($series as $key=>$val) {
0218:         $s1 .= "[ {$val} ], ";
0219:         if (max($items[$key]) > 20) {
0220: $s2 .=<<< EOT
0221:                 {
0222:                 label: '{$labels[$key]}',
0223:                 yaxis: 'y2axis'
0224:                 },
0225: 
0226: EOT;
0227:         } else {
0228: $s2 .=<<< EOT
0229:                 {
0230:                 label: '{$labels[$key]}'
0231:                 },
0232: 
0233: EOT;
0234:         }
0235:     }
0236: 
0237:     $width  = GRAPH_WIDTH;
0238:     $height = GRAPH_HEIGHT;
0239: $js =<<< EOT
0240: <script>
0241: jQuery(function() {
0242:     jQuery.jqplot('{$name}',
0243:     [
0244:         {$s1}
0245:     ],
0246:     {
0247:         //タイトル
0248:         title: {
0249:             text: '{$title}',
0250:             show: true,
0251:             fontSize: '20px',
0252:             textAlign: 'center',
0253:             textColor: 'black'
0254:         },
0255:         //ラベル
0256:         series: [
0257:             {$s2}
0258:         ],
0259:         legend: {
0260:             show: true,
0261:             placement: 'inside',
0262:             location: 'ne',
0263:             renderer: $.jqplot.EnhancedLegendRenderer,
0264:             rendererOptions: { numberRows: 1 }
0265:         },
0266:         //グラフ
0267:         seriesDefaults: {
0268:             showLine: true,
0269:             rendererOptions: { smooth: false },
0270:             markerOptions: { size: 0 },
0271:         },
0272:         //軸
0273:         axes: {
0274:             xaxis: {
0275:                 label: '{$labels[0]}',
0276:                 tickOptions: { formatString: '%d' }
0277:             },
0278:             yaxis: {
0279:                 label: '率',
0280:                 tickOptions: { formatString: '%.1f' }
0281:             },
0282:             y2axis: {
0283:                 label: '年齢',
0284:                 tickOptions: { formatString: '%d' }
0285:             }
0286:         }
0287:     }
0288:     );
0289: });
0290: </script>
0291: <div id="{$name}" style="width:{$width}px; height:{$height}px;"></div>
0292: 
0293: EOT;
0294: 
0295:     return $js;
0296: }

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

参考サイト

(この項おわり)
header