読者です 読者をやめる 読者になる 読者になる

Rで逆推定

R

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

f:id:Rion778:20150805235054p:plain
ここで、統計ソフトのJMPで同じようなことをしたとすると、この後に逆推定というやつができる。
すなわち、上記の例で言えば、dayが特定の値(例えば30)になるようなtempの値を求めるという操作である。このとき、tempの信頼区間も同時に出力される。
ちなみに、上記データに基づいてJMPでtemp=30について信頼水準=0.95で逆推定を行うと、

  • 予測値: 19.09091
  • 下側95%:18.09069
  • 上側95%:19.99694

という推定値が出力される。
どうもこれがRではそんなに簡単にはできないらしい。関数やパッケージの調査不足かもしれないので、良いやり方をご存じの方があれば教えて頂きたい*1

回帰式を変形させる

y = a + bx
という回帰式を、
y = A + b(x-c)
という風に変形する。このとき、Aは目的とするyの値(上記例で言えば30)であるとすると、x-c=0のときにy=Aとなるので、cが求めるべきxの値に等しいとういことになる。あとは、bcについて、非線形回帰を用いて求め、cの信頼区間を求める。
(参考)

非線形回帰および信頼区間

非線形回帰は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といういかにもな名前の関数が入っているが、それではイマイチ上手く行かなかった。使い方が間違っているのかもしれない。