StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

アッシェンフェルターのワイン価格予測式をトレースしてみた

1990年頃にアッシェンフェルターさん(Orley Ashenfelter)がビンテージワインの価格予測式(回帰式)を構築しました。精度はかなりよかったものの、当然ワイン評論家からはフルボッコにされました。という話を以前tokyo.RのLTで聞きました。

LTで時間が足りないせいもあってか(?)、式の各項の単位がよくわからなかったり、質は線形で増えるが価格はexpで増えるという説明があったりでちょっと気になって自分でその後調べてみました。ググってみるとこんな記事こんな記事がヒットします。背景の物語には詳しくなりましたが、肝心の式についてはますます謎は深まるばかり。

その後どうもしっくりこなかったので自分で原論文をもとにトレースしました。ソースはここにあり、データもここにありました。

Rで読み込んでちゃっちゃと重回帰すると

df <- read.table("winedata.txt", header=T, sep="\t")
ans <- lm(LPRICE2 ~ WRAIN + DEGREES + HRAIN + TIME_SV, df)
options(scipen=10)
summary(ans)

Call:
lm(formula = LPRICE2 ~ WRAIN + DEGREES + HRAIN + TIME_SV, data = df)

Residuals:
     Min       1Q   Median       3Q      Max
-0.46027 -0.23864  0.01347  0.18600  0.53446

Coefficients:
               Estimate  Std. Error t value    Pr(>|t|)
(Intercept) -12.1453336   1.6881027  -7.195 0.000000328 ***
WRAIN         0.0011668   0.0004820   2.421     0.02420 *
DEGREES       0.6163924   0.0951755   6.476 0.000001626 ***
HRAIN        -0.0038606   0.0008075  -4.781 0.000089726 ***
TIME_SV       0.0238474   0.0071667   3.328     0.00306 **
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2865 on 22 degrees of freedom
  (11 observations deleted due to missingness)
Multiple R-squared: 0.8275,     Adjusted R-squared: 0.7962
F-statistic: 26.39 on 4 and 22 DF,  p-value: 0.00000004058

となりまして、回帰の係数やRMSEは有効数字4~5桁まで原論文の結果と完全一致します。なのでこれでトレースは正しいのだと思います。結論を日本語で少し詳細に書きますと以下のようになります。

f:id:StatModeling:20201114163113p:plain

日本語で巷にあふれている解説と比較して気づくことを列挙しておきます。

  • 切片の符号が違う…
  • TIME_SVの項が落ちてる…(古いほうが高くなる)
  • 左辺は「ワインの質」ではない。質∝値段 と考えれば質もexpで増えることになる。

また原論文では降水量の部分の単位が「ML」になっているのですが、一般的には世界でも「mm」や「inch」などの高さで測られており、またここなどを見るとどうもここは「mm」が正しいように思います。詳しい方いましたら教えてください。いやー適当なコピペは罪ですよ、ほんとに…。