发帖
雷达卡
用R使用了如下命令：
> y <- rnorm(50)
> x <- matrix(rnorm(100),50)
> fit <- rq(y~x,tau = .4)
Call:
rq(formula = y ~ x, tau = 0.4)
Coefficients:
(Intercept) x1 x2
-0.38203826 0.09384216 0.09347414
Degrees of freedom: 50 total; 47 residual
Call: rq(formula = y ~ x, tau = 0.4)
tau: [1] 0.4
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) -0.38204 0.21700 -1.76056 0.08482
x1 0.09384 0.18203 0.51552 0.60861
x2 0.09347 0.24560 0.38059 0.70522
但是我现在想做的数据结构不是上面的一列数据格式（简称为截面），我现在想做面板（问一句，对于bootstrap，不同数据类型会不同吗，如面板数据和截面数据？），计算系数的程序已经
基本上有了，就是R中quantreg的rq.fit.sfn。但是这个程序只运行出coeff，而t值和p值都没有，所以我想用bootstrap计算（使用block bootstrap即可）
应该怎么实现？
> summary.rq
function (object, se = "nid", covariance = FALSE, hs = TRUE,
...)
{
if (object$method == "lasso")
stop("no inference for lasso'd rq fitting: try rqss (if brave)")
mt <- terms(object)
m <- model.frame(object)
y <- model.response(m)
x <- model.matrix(mt, m, contrasts = object$contrasts)
wt <- model.weights(object$model)
tau <- object$tau
eps <- .Machine$double.eps^(2/3)
coef <- coefficients(object)
if (is.matrix(coef))
coef <- coef[, 1]
vnames <- dimnames(x)[[2]]
resid <- object$residuals
n <- length(resid)
p <- length(coef)
rdf <- n - p
if (!is.null(wt)) {
resid <- resid * wt
x <- x * wt
y <- y * wt
}
if (missing(se)) {
if (n < 1001 & covariance == FALSE)
se <- "rank"
else se <- "nid"
}
if (se == "rank") {
f <- rq.fit.br(x, y, tau = tau, ci = TRUE, ...)
}
if (se == "iid") {
xxinv <- diag(p)
xxinv <- backsolve(qr(x)$qr[1:p, 1:p, drop = FALSE],
xxinv)
xxinv <- xxinv %*% t(xxinv)
pz <- sum(abs(resid) < eps)
h <- max(p + 1, ceiling(n * bandwidth.rq(tau, n, hs = hs)))
ir <- (pz + 1):(h + pz + 1)
ord.resid <- sort(resid[order(abs(resid))][ir])
xt <- ir/(n - p)
sparsity <- rq(ord.resid ~ xt)$coef[2]
cov <- sparsity^2 * xxinv * tau * (1 - tau)
scale <- 1/sparsity
serr <- sqrt(diag(cov))
}
else if (se == "nid") {
h <- bandwidth.rq(tau, n, hs = hs)
if (tau + h > 1)
stop("tau + h > 1: error in summary.rq")
if (tau - h < 0)
stop("tau - h < 0: error in summary.rq")
bhi <- rq.fit.fnb(x, y, tau = tau + h)$coef
blo <- rq.fit.fnb(x, y, tau = tau - h)$coef
dyhat <- x %*% (bhi - blo)
if (any(dyhat <= 0))
warning(paste(sum(dyhat <= 0), "non-positive fis"))
f <- pmax(0, (2 * h)/(dyhat - eps))
fxxinv <- diag(p)
fxxinv <- backsolve(qr(sqrt(f) * x)$qr[1:p, 1:p, drop = FALSE],
fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
fxxinv
scale <- mean(f)
serr <- sqrt(diag(cov))
}
else if (se == "ker") {
h <- bandwidth.rq(tau, n, hs = hs)
if (tau + h > 1)
stop("tau + h > 1: error in summary.rq")
if (tau - h < 0)
stop("tau - h < 0: error in summary.rq")
uhat <- c(y - x %*% coef)
h <- (qnorm(tau + h) - qnorm(tau - h)) * min(sqrt(var(uhat)),
(quantile(uhat, 0.75) - quantile(uhat, 0.25))/1.34)
f <- dnorm(uhat/h)/h
fxxinv <- diag(p)
fxxinv <- backsolve(qr(sqrt(f) * x)$qr[1:p, 1:p, drop = FALSE],
fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
fxxinv
scale <- mean(f)
serr <- sqrt(diag(cov))
}
else if (se == "boot") {
B <- boot.rq(x, y, tau, ...)
cov <- cov(B)
serr <- sqrt(diag(cov))
}
if (se == "rank") {
coef <- f$coef
}
else {
coef <- array(coef, c(p, 4))
dimnames(coef) <- list(vnames, c("Value", "Std. Error",
"t value", "Pr(>|t|)"))
coef[, 2] <- serr
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- if (rdf > 0)
2 * (1 - pt(abs(coef[, 3]), rdf))
else NA
}
object <- object[c("call", "terms")]
if (covariance == TRUE) {
object$cov <- cov
if (se == "iid")
object$scale <- scale
if (se %in% c("nid", "ker")) {
object$Hinv <- fxxinv
object$J <- crossprod(x)
object$scale <- scale
}
else if (se == "boot") {
object$B <- B
}
}
object$coefficients <- coef
object$rdf <- rdf
object$tau <- tau
class(object) <- "summary.rq"
object
}
京公网安备 11010802022788号
