Rで逆推定(2)

はてなから1年前のブログとかいうメールが来てこんな記事を書いていたのを思い出した。
rion778.hatenablog.com
再び少し調べたところ、これはinvestrというパッケージで実行できるらしい。

データとモデルオブジェクトの作成

データは前回と同じものを使うが、データフレームにしておく。これは、モデルオブジェクト中の$dataがNULLだと逆推定を行う際にデータの指定が必要となるためだ。

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)
d <- data.frame(temp, day)
result <- lm(day ~ temp, data = d)

逆推定

逆推定はinvest()関数で行う。さらに、plotFit()関数で推定区間と共にプロットできる。

library(investr)
(inv <- invest(result, y0 = 30, interval = "inversion", mean.response = TRUE))
# estimate    lower    upper 
# 19.09091 18.09070 19.99693 
plotFit(result, interval = "confidence", shade = TRUE)
abline(h = 30, v = c(inv$lower, inv$estimate, inv$upper), lty = 2)

f:id:Rion778:20160806220620p:plain
ちなみに、JMPで実行した場合の推定値は

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

であるので、ほぼ一致している。