線形回帰分析をする時、誤差項が同じ分布に従うという条件があります。分散不均一の場合、一般線形モデル（GLM、General linear model）で分析することが知られています。分散不均一性を検出する方法は、いくつあります。今日は、例を挙げて、Breusch–Pagan検定を紹介します。
Breusch–Pagan検定は、下記URLに参考してください。ここで説明しません。
https://en.wikipedia.org/wiki/Breusch%E2%80%93Pagan_test
簡単に説明します。図Total、A、B、C、D、Eは豆腐の生産量を表しています。生産量を目的関数と、気温、気温差、降水量を説明変数とし、モデル構築します。
図Total、A、B、C、D、Eの真ん中部分の変動は激しい。一方、気温、気温差、降水量の変動はそんなに差はない。直観に、誤差項の分散は同じではないと判断します。では、検定してみましょう。
検定の結果は、分散不均一であることが分かりました。
コードを添付します。
BP.rb
# データの読み取り
tofu<-read.csv("tofu.csv",header=T,sep=",",fill=T)
# 欠損値の処理
tofu<-tofu[-(51:52),]
tofu<-tofu[-46,]
BP<-matrix(0,1,6)
colnames(BP)<-c("total","A","B","C","D","E")
for (i in 1:6) {
# 回帰
R<-lm(tofu[,i]~tofu[,7]+tofu[,8]+tofu[,9])
x<-as.matrix(tofu[,6:9])
x[,1]<-1
xx<-as.matrix(R$coefficients)
# 残差の自乗を計算する
u2<-(tofu[,2]-x%*%xx)^2
# 補助式の回帰
RR<-lm(u2~tofu[,7]+tofu[,8]+tofu[,9])
# BP検定の補助回帰式の決定係数
R2<-summary(RR)$r.squared
BP[1,i]<-nrow(tofu)*R2
}
# True⇒分散均一 False⇒分散不均一
BP<=qchisq(0.05,3,lower.tail=FALSE)