[程式] [R] Rfit的robust信賴區間
---
雖然McKean and Hettmansperger (2015, p.10) 書中提供兩樣本比較的CI，
但其中的 R 範例無法運作...
我認為應該是作者寫錯了，
矩陣 mat 和 vc 沒有定義？
根據p.6的公式，既然迴歸係數
beta_hat approximately ~ 多元常態 MVN(beta_0, tau^{2}(X'X)^{-1})
其中 tau 是 scale parameter
那直接拿 rfit() 輸出結果的 Estimate + t臨界值*Std. Error 不就得了？
R 範例程式碼：
library("Rfit")
data(serumLH)
# full model fit
fitmod <- rfit(serum~factor(light.regime)+
factor(LRF.dose)+ factor(light.regime)*
factor(LRF.dose), data=serumLH)
# hvec picks the contrast
hvec <- rep(0,60); hvec[27] <- -1; hvec[57] <- 1
# estimate of contrast
contr <- hvec%*%fitmod$fitted.values
# error in CI
se2 <- t(hvec)%*%mat%*%vc%*%t(mat)%*%hvec
# error term in the confidence interval.
err <- qt(.975,50)*sqrt(se2)
#(McKean and Hettmansperger, 2015, p.10)
Ref:
McKean, J. W., & Hettmansperger, T. P. (2015, April). Rank-Based Analysis of
Linear Models and Beyond: A Review. Paper presented at the International
Conference on Robust Rank-Based and Nonparametric Methods, Kalamazoo,
Michigan, U.S.
感謝撥冗回答~
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 218.164.23.158
※ 文章網址: https://www.ptt.cc/bbs/Statistics/M.1526644187.A.28A.html
※ 編輯: MTIS (218.164.23.158), 05/18/2018 19:52:25
------
補充：
McKean回信並提供了code：
> library(Rfit)
> head(serumLH)
serum light.regime LRF.dose
1 72 Constant 0
2 64 Constant 0
3 78 Constant 0
4 20 Constant 0
5 56 Constant 0
6 70 Constant 0
> attach(serumLH)
The following objects are masked from serumLH (pos = 3):
light.regime, LRF.dose, serum
> fitmod <-
rfit(serum~factor(light.regime)+factor(LRF.dose)+factor(light.regime)*factor(LRF.dose),data=serumLH)
> hvec=rep(0,60);hvec[27]=-1;hvec[57]=1
> contr = hvec%*%fitmod$fitted.values
> contr
[,1]
[1,] 200.3648
> mat <- fitmod$x
> vc <- vcov(fitmod)
> se2 <- (t(hvec)%*%mat%*%vc%*%t(mat)%*%hvec)^.5
> err<- qt(.975,50)*se2
> err
[,1]
[1,] 65.35394
※ 編輯: MTIS (59.126.233.13), 05/30/2018 00:25:09