JMPには逆推定というコマンドがある
次のようなデータを用意して、線形回帰を行うとする。
temp <- rep(c(15, 20, 25), c(5, 5, 5)) day <- c(34, 33, 36, 37, 35, 30, 28, 29, 25, 28, 23, 25, 22, 24, 26)
こんなかんじで適当に。
> result <- lm(day~temp) > plot(day~temp) > abline(result) > summary(result) Call: lm(formula = day ~ temp) Residuals: Min 1Q Median 3Q Max -4.00 -1.00 0.00 1.25 2.50 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 51.000 2.307 22.110 1.07e-11 *** temp -1.100 0.113 -9.734 2.46e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.787 on 13 degrees of freedom Multiple R-squared: 0.8794, Adjusted R-squared: 0.8701 F-statistic: 94.76 on 1 and 13 DF, p-value: 2.457e-07
ここで、統計ソフトのJMPで同じようなことをしたとすると、この後に逆推定というやつができる。
すなわち、上記の例で言えば、dayが特定の値（例えば30）になるようなtempの値を求めるという操作である。このとき、tempの信頼区間も同時に出力される。
ちなみに、上記データに基づいてJMPでtemp=30について信頼水準=0.95で逆推定を行うと、
- 予測値: 19.09091
- 下側95%:18.09069
- 上側95%:19.99694
という推定値が出力される。
どうもこれがRではそんなに簡単にはできないらしい。関数やパッケージの調査不足かもしれないので、良いやり方をご存じの方があれば教えて頂きたい*1。
回帰式を変形させる
という回帰式を、
という風に変形する。このとき、は目的とする
の値（上記例で言えば30）であるとすると、
のときに
となるので、
が求めるべき
の値に等しいとういことになる。あとは、
と
について、非線形回帰を用いて求め、
の信頼区間を求める。
（参考）
非線形回帰および信頼区間
非線形回帰はnlsを用いて行う。startの値はどう選ぶのが適切かはよく分からないが、選び方がまずいと動かない時がある。
result2 <- nls(day~30+b*(temp-c), start=c(b=0.1,c=0))
予測値はnlsオブジェクトの中身を見れば書いてある（これは回帰直線上の話だからわざわざ非線形回帰をしなくても分かるが）。
> result2 Nonlinear regression model model: day ~ 30 + b * (temp - c) data: parent.frame() b c -1.10 19.09 residual sum-of-squares: 41.5 Number of iterations to convergence: 8 Achieved convergence tolerance: 1.162e-09
パラメーターの信頼区間は、confintで求める。
> confint(result2) Waiting for profiling to be done... 2.5% 97.5% b -1.344124 -0.8558761 c 18.090557 19.9969627
JMPの結果とほんの少し違うのが気になるところ。
*1:Rではlmオブジェクトに対してpredictという関数が用意されているが、これは例で言えばtempが特定の値のときのdayの推定値と信頼区間を求めるというもので、逆推定はできない。また、chemCalというパッケージの中にinverse.predictといういかにもな名前の関数が入っているが、それではイマイチ上手く行かなかった。使い方が間違っているのかもしれない。
