はてなから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)
ちなみに、JMPで実行した場合の推定値は
- 予測値: 19.09091
- 下側95%:18.09069
- 上側95%:19.99694
であるので、ほぼ一致している。