一般には説明変数にも誤差を含む

　これまでの線形回帰のパラメータの推定においては、残差平方和を最小化するようなパラメータを求めていた。残差平方和とは、下図のように、回帰直線と実測値yiの差の二乗和のことを指す。

残差二乗和、すなわち回帰直線と実測値のずれが最小の回帰直線こそ、最も当てはまりの良い直線であると判断するわけだ。ただし、ここにはある仮定が存在し、それは「被説明変数だけに誤差がある」点である。この仮定は、説明変数を人間がコントロールしている場合、妥当である。このような被説明変数だけに誤差を仮定する回帰をI型回帰Type I regressionと呼ぶ。I型回帰は説明変数から被説明変数を予測するときに用いられる。

　しかし、野外データにおいて、複数の変数間の関係を考える場合はそうもいかない。「説明変数にも被説明変数にも誤差が存在する」のが普通である。そこで、説明変数側の誤差も考慮した回帰直線と実測値との距離の取り方として、下記のようなものが考えられる。

このような、説明変数と被説明変数の誤差を両方考慮する回帰をII型回帰Type II regressionと呼ぶ。II型回帰は予測よりも変数間の量的関係を明らかにするときに用いられる。このII型回帰に当たる解析には、主軸回帰Major Axis regression（MA回帰、主成分回帰、最大軸回帰）が存在する。さらに、その特殊な状況として、標準化主軸回帰Standarized Major Axis regression（SMA回帰、縮小主軸 (Reduced major axis, RMA) 回帰，幾何平均 (Geometric mean) 回帰）がある。

　ここからはまず、MA回帰について紹介し、その後、MA回帰の特殊な条件の解析であるSMA回帰について紹介する。

MA回帰の係数の推定

　では、上記の「距離平方和R」を最小にするαとβを偏微分を使って推定しよう。αで偏微分すると以下の通り。この式が0となるαを求める。

上記のように、αの推定値は線形回帰の時と同じになることがわかる。

　続いて、βで偏微分するが、その前に上記のαの推定値を使ってRの式変形を終えておこう。

上記の式をβで偏微分する。偏微分した式が0となる値を求める。

ここで、Varは標本分散、Covは標本共分散である。最終的なβの推定値は、かなりけったいな式になってしまった。βの値は2次方程式の解なので、2つ候補が存在するが、このうちの一つがRの最小値を与え、もう一つが最大値を与える。下記のようにCov(x,y) > 0なら大きい方のβ、Cov(x,y) < 0なら小さい方のβが、Rの最小値を与える。Cov(x,y)の符号に注意すると、いずれの場合でも正の平方根のほうがRの最小値を与えるβであることがわかる。

この傾きの推定値、^βは主成分分析の固有ベクトルと一致する。2変数の主成分分析において、以下の関係を満たすsinφ、cosφを使うと、固有ベクトルと傾きは以下のようにあらわせるのだった。さらに、傾きを計算していくと、以下のようになり、x1 = x、x2 = yとすれば、上記と一致する。

以上のように、MA回帰がやっていることは、主成分分析と変わりがなく、それゆえに主成分回帰とも呼ばれるのである。

SMA回帰

　MA回帰は生のデータをそのまま解析するが、SMA回帰では説明変数と被説明変数をともに平均0、標準偏差1に標準化Standarizationしてから、MA回帰を行う。この考え方は、主成分分析と同様で、桁数の違いがある変数間の比較において妥当な変換である。具体的には、変数は下記のような変換を行う。

標準化後にMA回帰を行うと、Var(x') = Var(y') = 1であるから、^β' = ±1となる。また、E(x') = E(y') = 0でもあるので、^α' = 0である。ゆえに、標準化後のMA回帰直線は以下のようにあらわせる。これを変数変換前に戻すことで、SMA回帰の結果を得られる。

すなわち、SMA回帰の結果は、^αは線形回帰やMA回帰と同様の結果であるが、^βの値が変わって、Var(y)/Var(x)の平方根となる。MA回帰と異なり、ずいぶんと簡潔な表記である。

・SMA回帰で最小化しているもの

　MA回帰とSMA回帰では類似したことを行っているが、^βは異なる値となった。上記のように、最小化問題の解のひとつとして得られたということは、SMA回帰でも何かの最小化する問題に帰着できることは想像に難くない。実際、SMA回帰の結果は、変換前の変数において、下記の図に示した回帰直線と実測値が作る三角形の面積の最小化問題の解として得られる。

　β > 0のときは

　β < 0のときは

　以上をまとめると以下のような式になり、β > 0のときについて、面積和をαおよびβで偏微分すると、

　β < 0のときも同様に

以上の情報をもとに増減表を書くと下記。

βの正負、いずれの解が最小値をとるかは、Covの正負に依存しており、下の議論から、Covが正なら正の解、負なら負の解が最小値を与える。

以上から、SMA回帰で求めるパラメータは、回帰直線と実測値がなす三角形の面積の最小化問題の解である。

Rで行うMA回帰/SMA回帰

　では、RでMA回帰およびSMA回帰を行ってみよう。これらの解析を行えるパッケージはlmodel2パッケージやsmatrパッケージがある。今回はsmatrパッケージを使ってみよう。

------------------------------------------------------

library(plotn)

library(smatr)

------------------------------------------------------

では、下記のようにサンプルサイズ50のデータを生成し、xとする。同じデータをyにも格納しよう。すると、当然のことながら、これらのデータの対応関係は一直線上に乗る。

------------------------------------------------------

n <- 50

x <- runif(n, -2, 2)

y <- x

plotn(y ~ x)#図1の図示

------------------------------------------------------

ではここに誤差をのせることを考える。まずは、yだけに誤差をのせてみよう。

------------------------------------------------------

ye <- y + rnorm(n)

plotn(ye ~ x)#図2の図示

------------------------------------------------------

このとき、I型回帰の仮定に沿っている。ゆえに、通常の線形回帰を行う。この時の傾きは、以前紹介したように共分散と分散から計算することができる。

------------------------------------------------------

res1 <- lm(ye ~ x)

summary(res1)

## 

## Call:

## lm(formula = ye ~ x)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -2.43124 -0.69412  0.00517  0.69923  1.91451 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -0.08887    0.13607  -0.653    0.517    

## x            1.07841    0.11618   9.282 2.74e-12 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.9612 on 48 degrees of freedom

## Multiple R-squared:  0.6422, Adjusted R-squared:  0.6347 

## F-statistic: 86.15 on 1 and 48 DF,  p-value: 2.741e-12

cov(x, ye)/var(x)

## [1] 1.078409

------------------------------------------------------

もともと、y = xであるから、傾きは1、切片は0と推定されるはずであり、実際、正しく推定されている。これらの結果から、回帰直線は以下のように引くことができる。

------------------------------------------------------

plotn(ye ~ x)#図3の図示

xx <- seq(min(x), max(x), length = 200)

yy <- coef(res1)[1] + coef(res1)[2]*xx

overdraw(points(xx, yy, type = "l", col = "red"))

------------------------------------------------------

では、本題として、yだけでなくxにも誤差をのせてみる。

------------------------------------------------------

xe <- x + rnorm(n)

ye <- y + rnorm(n)

plotn(ye ~ xe)#図4の図示

------------------------------------------------------

このとき、仮定としてはII型回帰が妥当であるが、通常の線形回帰とともに解析をしてみよう。MA回帰やSMA回帰はともに、smatrパッケージのsma関数を使えば実行できる。method引数をMAとすればMA回帰、SMAとすればSMA回帰である。デフォルトではSMA回帰になる。

------------------------------------------------------

res1 <- lm(ye ~ xe)

summary(res1)

## 

## Call:

## lm(formula = ye ~ xe)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -2.5903 -0.6584 -0.1094  0.7660  2.9100 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  0.04651    0.15839   0.294     0.77    

## xe           0.46001    0.08691   5.293 2.96e-06 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 1.119 on 48 degrees of freedom

## Multiple R-squared:  0.3685, Adjusted R-squared:  0.3554 

## F-statistic: 28.01 on 1 and 48 DF,  p-value: 2.956e-06

------------------------------------------------------

通常の線形回帰では、0.46と傾きが随分と小さく推定されている。では、MA回帰ではどうなるだろうか。

------------------------------------------------------

res2 <- sma(ye ~ xe, method = "MA")

res2

## Call: sma(formula = ye ~ xe, method = "MA") 

## 

## Fit using Major Axis 

## 

## ------------------------------------------------------------

## Coefficients:

##               elevation     slope

## estimate     0.03576341 0.6390932

## lower limit -0.29662223 0.4141873

## upper limit  0.36814904 0.9216762

## 

## H0 : variables uncorrelated

## R-squared : 0.3685172 

## P-value : 2.9558e-06

------------------------------------------------------

MA回帰では、傾きが0.639となり、かなりマシな推定になる。そしてSMA回帰も行ってみよう。

------------------------------------------------------

res3 <- sma(ye ~ xe, method = "SMA")

res3

## Call: sma(formula = ye ~ xe, method = "SMA") 

## 

## Fit using Standardized Major Axis 

## 

## ------------------------------------------------------------

## Coefficients:

##               elevation     slope

## estimate     0.02863951 0.7577649

## lower limit -0.32658994 0.6029003

## upper limit  0.38386897 0.9524089

## 

## H0 : variables uncorrelated

## R-squared : 0.3685172 

## P-value : 2.9558e-06

------------------------------------------------------

すると傾きの推定は0.758となり、やはり線形回帰よりもマシな推定となる。SMA回帰は変数間の桁数の違いの影響を受けなくなるので、基本的な選択肢としてはSMA回帰を行うことが多いだろう。

　上記の傾きはすべて以下のように計算可能である。

------------------------------------------------------

#線形回帰の傾き

cov(xe, ye)/var(xe)

## [1] 0.4600059

#MA回帰の傾き

Vx <- sum((xe-mean(xe))^2)/length(xe)

Vy <- sum((ye-mean(ye))^2)/length(ye)

C <- sum((xe-mean(xe))*(ye-mean(ye)))/length(xe)

((Vy-Vx)+sqrt((Vx-Vy)^2+4*C^2))/(2*C)

## [1] 0.6390932

#SMA回帰の傾き

sign(cov(xe,ye))*sqrt(var(ye)/var(xe))

## [1] 0.7577649

------------------------------------------------------

また、MA回帰の傾きは、主成分分析の固有ベクトルがなす傾きと一致する。

------------------------------------------------------

tmp <- prcomp(cbind(xe, ye))

tmp$rotation[2,1]/tmp$rotation[1,1]

## [1] 0.6390932

------------------------------------------------------

では上記の結果を使って、回帰直線を描いてみる。線形回帰およびSMA回帰の結果を用いよう。

------------------------------------------------------

plotn(ye ~ xe)#図5の描画

xx <- seq(min(xe), max(xe), length = 200)

yy1 <- coef(res1)[1] + coef(res1)[2]*xx

yy2 <- coef(res3)[1] + coef(res3)[2]*xx

overdraw(points(xx, yy1, type = "l", col = "red"),

         points(xx, yy2, type = "l", col = "blue"))

------------------------------------------------------

では、実際、SMA回帰がどれほど妥当な解析になるか、確かめてみよう。以下のように、10000回、xとyに誤差を加えて、線形回帰およびSMA回帰を行い、その時の係数を記録しておく。

------------------------------------------------------

b0_1 <- NULL

b1_1 <- NULL

b0_2 <- NULL

b1_2 <- NULL

for(i in 1:10000){

  xe <- x + rnorm(n)

  ye <- y + rnorm(n)

  res1 <- lm(ye ~ xe)

  res2 <- sma(ye ~ xe, method = "SMA")

  b0_1 <- c(b0_1, coef(res1)[1])#線形回帰の切片

  b1_1 <- c(b1_1, coef(res1)[2])#線形回帰の傾き

  b0_2 <- c(b0_2, coef(res2)[1])#SMA回帰の切片

  b1_2 <- c(b1_2, coef(res2)[2])#SMA回帰の傾き

}

------------------------------------------------------

では各係数のヒストグラムを描く。

------------------------------------------------------

histn(b0_1)#図6の描画

histn(b1_1)#図7の描画

histn(b0_2)#図8の描画

histn(b1_2)#図9の描画

------------------------------------------------------

各パラメータの平均値は以下のようになる。

------------------------------------------------------

mean(b0_1)

## [1] -0.02102343

mean(b1_1)

## [1] 0.5893767

mean(b0_2)

## [1] 0.001234758

mean(b1_2)

## [1] 1.006381

------------------------------------------------------

確かに、線形回帰における傾きは過少推定になっているのに対して、SMA回帰では平均的に見れば正しい値の1に近い推定となっている。