作者MTIS ( )
看板Statistics
标题[程式] [R] Rfit的robust信赖区间
时间Fri May 18 19:49:45 2018
---
虽然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://webptt.com/cn.aspx?n=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