統計学 対数変換後、元のスケールに戻して信頼区間を算出する方法

IPAが公開しているソフトウエア開発データ白書で、
信頼区間　５０％、９５％の具体的な計算方法がわかりません。
計算で信頼区間を算出し、散布図に描画したのですが、公開されているデータ、図とぜんぜん合いません。


対象データ：
https://www.ipa.go.jp/ikc/publish/whitepaper_dl. …

対象ファイル：
graph_data_06.xlsx
「6-2-1」シート

対応する図：
ソフトウェア開発データ白書 2018-2019_本編.pdf
P136 図表6-2-2


■やったこと：
信頼区間の算出については、以下のサイトを参考にしました。
https://cyclo-commuter.hatenablog.jp/entry/2018/ …

エクセルで、データ白書のデータについて対数変換してから、
信頼区間をプロットしてみたのですが、IPAの図表の信頼区間と全く合致しません。。
（データと近似曲線（累乗近似）は、エクセルでIPAどおりに描画できています）

統計については、エクセルやR言語、EZRで最近勉強はじめまた初心者です。
エクセル、R言語、どちらの方法でも結構ですので、対数スケールに変換してから通常スケールに戻す場合の信頼区間の計算方法についてご教示ください。

■補足１
X軸:実績工数（人時）
Y軸:実績月数

データ白書の説明（p28-29）によれば、プロジェクトのデータは正規分布していない事が多いが、
対数に変換するとほぼ正規分布と見なせる。このため、対数スケールに変換してから通常スケールに
戻して信頼区間付き散布図を作成したとあります。

■補足２
データ白書のP29に「※ 信頼区間の算出は、
田中豊、脇本和昌、｢多変量統計解析法｣、現代数学社、1983
涌井良幸、涌井貞美、｢図解でわかる回帰分析、日本実業出版社、2002 を参考にした。」
とあり、まず「図解でわかる回帰分析」を読んでいるところですが、道は長いです。
非線形回帰の話が出てくるのですが、単回帰で計算できるはずなので、深いところまでは
追っていません。統計は奥が深いですね。

No.7 ベストアンサー 回答者： kamiyasiro 回答日時： 2021/02/06 11:49

#1です。

ちなみに、yを一旦対数変換し、等分散仮定で、ｙ～log(x)のモデルで近似式を作り、全体をexpで戻してプロットした図を添付します。

ご質問者様も、こんなグラフになったのではないでしょうか？

この回答へのお礼

ありがとうございます。そうです、こういう図になりました。
IPAが公開されている図の信頼区間とずれていますが、
計算結果はこうなのですね。

Rスクリプト、および、図の解説ありがとうございます！

＞あの図6-2-1の信頼区間は絶対におかしいです。あれが再現できても意味ないですよ。いいですか。信頼区間っていうのは、いわば回帰線の存在範囲ですよ。にもかかわらず、あの図の95％信頼区間は、まるで予測区間（データの存在範囲）のような位置に来ています。

私も予測区間で算出し、プロットしてみたのですが、
ご指摘のとおり９５％の信頼区間と似ていると感じました。

大変勉強になりました。一人では、ここまで理解できなかったです。
この結果をどう評価するのか・・・・はおいておいて、
ひとまず、自分で試してみるという理解・道筋への大きな
一歩となりそうです。

大変ありがとうございます。

お礼日時：2021/02/06 13:35

No.8 回答者： kamiyasiro 回答日時： 2021/02/06 12:09

#1です。

信頼限界幅の係数ですが、参考にされたＨＰの式は係数がルートの中に入っているので、F値で良いです。

この回答へのお礼

ご教示いただいたRスクリプトで
一般化線形モデル(GLM）の初歩的なやり方を勉強させていただきました。
https://www1.doshisha.ac.jp/~mjin/R/Chap_16/16.h …
GLMは「正規分布を含んだ分布族にデータを対応させ、
非線形の現象を線形モデルの場合と同じく簡単に扱え
かつ不自然な尺度で解釈しないよう
工夫したデータ解析方法」とあり、なるほど。

近似値の計算のところで、どちらが良いか調べている部分が、理解がついていかなかったのですが、以下の文献で概要がわかりました。
https://www.agr.nagoya-u.ac.jp/~seitai/document/ …
⇒3.3 予測値とその 95% 信頼区間の推定

この文献では、dのことを「最尤推定値」といっているようでした。
また、その後の処理で、一般化線形モデルによる予測(predict）
を使っている部分も、こうやるのかと勉強になりました。

あとは、そのデータを使って
信頼区間９５％や５０％を算出してプロットすればよいのですね。
当然ですが、R言語はエクセルで手計算でやるより簡単ですね。

Rスクリプトから沢山のことが分かりました。
改めて、お礼申し上げます。

お礼日時：2021/02/10 19:06

No.6 回答者： kamiyasiro 回答日時： 2021/02/06 09:49

#1です。

rm(list = ls())

x <- read.csv("6-2-1.csv")
colnames(x) <- c("実績工数（プロジェクト全体）［人時］",
"実績月数（プロジェクト全体）［月］")


par(mar=c(10.1,4.1,8.1,2.1))
plot(x, pch = 4, col = 4, cex = 0.6, cex.lab = 0.8, las = 1,
xaxs = "i", yaxs = "i", xlim = c(0, 300000), ylim = c(0, 40),
main = "図表 6-2-1 ● プロジェクト全体の工数と工期（新規開発）（信頼区間 50%、95% 付き）",
cex.main = 0.8)
axis(side = 2, tck = 1.0, labels = FALSE)

n <- nrow(x)
mtext(paste("N =",n), line = 0, at = 280000)
phi <- n - 2 # 原点を通しているから１かもしれない

# -------------------------------------------------------
# GLM

colnames(x) <- c("x", "y")
result <- glm(y ~ log(x), data = x, family = "Gamma")
summary(result)

# 近似線を求める計算点を準備
xpred <- data.frame(x = seq(0, 300000, by=1000))

# 近似値の計算と描画

d <- 1/summary(result)$dispersion
d

library("MASS")
d <- gamma.shape(result)$alpha
d # こちらの方が良いらしい

preds <- predict(result, xpred, type = "response", dispersion = d, se.fit=TRUE)
lines(cbind(xpred, preds$fit))

# 95%信頼区間
critp <- qt(0.975, phi)
lines(cbind(x = xpred, y = preds$fit + critp * preds$se.fit), col = 2)
lines(cbind(x = xpred, y = preds$fit - critp * preds$se.fit), col = 2)

# 50%信頼区間
critp <- qt(0.75, phi)
lines(cbind(x = xpred, y = preds$fit + critp * preds$se.fit), col = 6)
lines(cbind(x = xpred, y = preds$fit - critp * preds$se.fit), col = 6)



# 凡例を描く（省略）

par(mar=c(5.1,4.1,4.1,2.1)) #デフォルト

No.5 回答者： kamiyasiro 回答日時： 2021/02/06 09:46

#1です。

対数を取ってから平均を取るために足し算すると、対数の足し算は掛け算だから幾何平均を取っていることになるのです。

回帰に関してはガウスマルコフの定理っていうのがあり、それを満たさなければなりません。そのひとつにばらつきの期待値E(ε)は０、つまり上下の平均を通るというのがあるのですが、幾何平均を通るとガウスマルコフの定理に反し、不偏推定量にならないのです。

あの図6-2-1の信頼区間は絶対におかしいです。あれが再現できても意味ないですよ。いいですか。信頼区間っていうのは、いわば回帰線の存在範囲ですよ。にもかかわらず、あの図の95％信頼区間は、まるで予測区間（データの存在範囲）のような位置に来ています。

手順に則って描くと添付図のようになります。誤差モデルはGamma（変動係数一定）にしています。

続いて、Rスクリプトを投稿します。

この回答へのお礼

IPAに本件を質問し、以下の回答をいただきました。
---------
信頼区間については、以下のSECジャーナルVol. 8, No.1（通巻28号）(2012年）に
掲載の「エンタプライズ系ソフトウェアプロジェクトにおける層別生産性とその信頼
区間」で紹介しております。
https://www.ipa.go.jp/files/000024515.pdf

この論文のp.19にあるように「予測値Yの50%及び95％の信頼区間」
（＝予測区間）になっています。

実測値から求めた平均値（推測値：Xに対応する回帰直線上の値）の信頼区間と
誤解する表現になっており、申し訳ありませんでした。

後はご参考になりますが、上記の区間の50％を求める式はp.19の(1)式の　
Yhat ± t0.5(n-2)×sqrt({1 + 1/n + (X-Xbar)^2/nsXX}Ve)
但しVeは回帰直線からの誤差分散になります。
(1)式の平方根(sqrt)の中の1/n＋(X-Xbar)^2/nsXXの部分が平均値（推測値）の
信頼区間部分（らっぱ形状はここから出ていて、標本数が大きくなると、
この部分は小さくなります）で、平方根の中の1の部分がデータの母集団の
ばらつきに対応する部分となります。
---------
#1様のご指摘どおり、予測区間でした。

お礼日時：2021/02/18 06:30

No.4 回答者： kamiyasiro 回答日時： 2021/02/06 01:44

#1です。

普通はyがリニアスケールのときは、回帰線は平均を通っていきますよね。

ｙを対数変換したときに平均を通るってことは、expで戻したときは幾何平均を通るってことですよね。そんなの回帰線の定義から外れますよね。

だから、あのグラフ、回帰線をあえて入れてないのですかね。

この回答へのお礼

ありがとうございます。勉強になります。

＞普通はyがリニアスケールのときは、回帰線は平均を通っていきますよね。
ここは分かります。

＞ｙを対数変換したときに平均を通るってことは、expで戻したときは幾何平均を通るってことですよね。そんなの回帰線の定義から外れますよね。
こちらは、幾何平均を通るというのが、勉強不足でわかっていないです。
調べてみます。

■やったことの補足
対数への変換ですが、x、ｙともにlogで対数をとって
logx、logyの値をもとに単回帰、信頼区間９５％の上下限を
算出しています。

※データ白書(P29)の図表3-3-7の左側みたいにプロットはできていました。（作った散布図の見た目は似てる）

※信頼区間５０％の算出は、自信ないのです。
https://cyclo-commuter.hatenablog.jp/entry/2018/ …
の信頼区間の幅で、F(1,個体の数-2;0.50)が５０％に該当するだろうと思い
計算しています。

問題は、元のスケールにどうやって戻すかわかっておらず、logをとったので、元のスケールに戻すためにexpを使って信頼区間下限、上限の値を
もどしました。No.3のご回答にも「expで全体を戻す」とあるので
考え方はあってると思うのですが。。。

この結果が、図表3-3-7の右の信頼区間のようになりません。
（累乗近似の線に近いところにひかれてしまいます。）

お礼日時：2021/02/06 08:16

No.3 回答者： kamiyasiro 回答日時： 2021/02/06 00:39

#1です。

『データ白書の説明（p28-29）によれば、プロジェクトのデータは正規分布していない事が多いが、対数に変換するとほぼ正規分布と見なせる。このため、対数スケールに変換してから通常スケールに戻して信頼区間付き散布図を作成したとあります。』

これって、ｙの誤差の話？

だったら、上下非対称というのも頷けます。

ｙを対数変換した空間で回帰し信頼区間を求めたうえで、expで全体を戻すんですかね。
そのときも、ｙ～A＋x^B の式で良いのかしらん？

誤差モデルとしては、なんだか変な話ですよね。

面白い問題ですね。

No.2 回答者： kamiyasiro 回答日時： 2021/02/06 00:12

#1です。

いやー、あの図の信頼区間はおかしいですよ。なぜなら、ｙはリニアスケールなのに、上下で非対称！

poissonでやってもGAMMAでやっても合わないので、2時間くらい悩んで無駄な時間を費やしてしまいました。

なぜ、最初に気づかなかったのだろう。ちゃんとした解析者がやったと思って結果を信じたから？

ご質問者様が合わないのも無理ないですよ。

あとで私が作成したRスクリプトを提供します。

この回答へのお礼

ご回答ありがとうございます。
図は6-2-1でした。失礼しました。


信頼区間がおかしいというご指摘に、ものすごいショックを受けております。

統計をあまり知らない私からすると、あの図の信頼区間は
なんとなく納得感があったのですが・・・・。

IPAのデータ白書は、毎年刊行され、数年の歴史があります。
図の信頼区間もほぼ同じです。（なにか意図があるのかなぁ・・・）

有識者の方から、アドバイスをいただき、あるべき正しい信頼区間
を求めたいと思いますので、
ご教示のほど、よろしくお願いします。

お礼日時：2021/02/06 00:56

No.1 回答者： kamiyasiro 回答日時： 2021/02/05 21:56

企業で統計を推進する立場の者です。

グラフ見ました。

一般化線形モデルでやっているとみて、今、Rスクリプト書き始めたんですが、プロットが違います。対応する図は6-2-2ではなく6-2-1ですね。

ご質問様のやられている数値変換して線形回帰するのとは全く異なる方法です。誤差の仮定が違います。

一致する図ができたら、Rスクリプトを提供します。