在qtreg框架之外的R中引导分位数回归的CI

Question

我有一个对象，model1，由分位数回归产生。在model1中，我有3列和99行，步长为1百分位，如下所示:

> model1
tau intercept  julian_day
1  0.01  17.25584 -0.02565709
2  0.02  19.11576 -0.04716498
3  0.03  20.92364 -0.07115870
4  0.04  22.74398 -0.09706853
5  0.05  24.54071 -0.12248334
6  0.06  26.19349 -0.14133451
7  0.07  27.83602 -0.15839657
8  0.08  29.64554 -0.18364227

我想计算并绘制截距值和变量值的置信区间，使用自举方法。

我没有使用quantreg包进行分位数回归，而是使用以下方法:

#Functions
minimize.logcosh <- function(par, X, y, tau) {
diff <- y-(X %*% par)
check <- (tau-0.5)*diff+(0.5/0.7)*logcosh(0.7*diff)+0.4
return(sum(check))
}

smrq <- function(X, y, tau){
p <- ncol(X)
op.result <- optim(
rep(0, p),
fn = minimize.logcosh,
method = 'BFGS',
X = X,
y = y,
tau = tau
)
beta <- op.result$par
return(beta)
}

run_smrq <- function(data, fml, response) {
x <- model.matrix(fml, data) #modify
y <- data[[response]]
#X <- cbind(x, rep(1,nrow(x)))
X <- x

n <- 99
betas <- sapply(1:n, function(i) smrq(X, y, tau=i/(n+1)))
return(betas) 
}

#Callers
smrq_models <- df %>%
group_by(lat_grouped) %>%
group_map(~ run_smrq(data=., fml=julian_day~year_centered, response="julian_day"))

df是这样的:

> dput(head(df, 20))
structure(list(lat = c("59", "59", "55", "59", "59", "63", "59", 
"59", "59", "59", "63", "59", "59", "59", "57", "56", "56", "59", 
"63", "63"), long = c(18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 
18, 18, 18, 18, 18, 18, 18, 18, 18, 18), date = c("1951-03-22", 
"1951-04-08", "1952-02-03", "1952-03-08", "1953-02-22", "1953-03-12", 
"1954-01-16", "1954-02-06", "1954-03-14", "1954-03-28", "1954-04-02", 
"1955-01-23", "1955-03-06", "1955-03-13", "1955-04-08", "1955-04-11", 
"1955-04-12", "1956-03-25", "1956-04-01", "1956-04-02"), julian_day = c(81, 
98, 34, 68, 53, 71, 16, 37, 73, 87, 92, 23, 65, 72, 98, 101, 
102, 85, 92, 93), year = c(1951L, 1951L, 1952L, 1952L, 1953L, 
1953L, 1954L, 1954L, 1954L, 1954L, 1954L, 1955L, 1955L, 1955L, 
1955L, 1955L, 1955L, 1956L, 1956L, 1956L), decade = c("1950-1959", 
"1950-1959", "1950-1959", "1950-1959", "1950-1959", "1950-1959", 
"1950-1959", "1950-1959", "1950-1959", "1950-1959", "1950-1959", 
"1950-1959", "1950-1959", "1950-1959", "1950-1959", "1950-1959", 
"1950-1959", "1950-1959", "1950-1959", "1950-1959"), time = c(10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L), lat_grouped = c("1", "1", "1", 
"1", "1", "2", "1", "1", "1", "1", "2", "1", "1", "1", "1", "1", 
"1", "1", "2", "2"), year_centered = structure(c(-36, -36, -35, 
-35, -34, -34, -33, -33, -33, -33, -33, -32, -32, -32, -32, -32, 
-32, -31, -31, -31), class = "AsIs")), row.names = 24:43, class = "data.frame")

你有什么建议如何实现我的目标吗?我已经尝试了许多在互联网上找到的启动功能，但到目前为止还没有能够使一些工作。

非常感谢您的帮助，如果我错过了什么，请不要犹豫，我可以更新我的帖子!

Answer 1

首先，问题的代码中缺少一个函数。

logcosh <- function(x) log(cosh(x))

对于这个问题，编写一个函数对data.frame的行进行替换采样，并像问题中一样运行分位数回归。没有必要将它们保存在对象smrq_models中，因为结果将立即返回给调用者。

注意，我已经从函数run_smrq中删除了常数n = 99，现在它是一个参数。

然后，在循环中调用boot函数。

最后，为了获得自举估计，必须首先将分位数回归代码返回的列表制作成数组，以便apply计算平均值。赋行名，我们就完成了。

library(tidyverse)
boot_fun <- function(data, n) {
i <- sample(nrow(data), nrow(data), replace = TRUE)
df <- data[i, ]
df %>%
group_by(lat_grouped) %>%
group_map(~ run_smrq(data=., fml=julian_day~year_centered, response="julian_day", n=n))
}
set.seed(2022)
n <- 99L
R <- 10L
boot_smrq_models <- vector("list", length = R)
for(i in seq.int(R)) {
boot_smrq_models[[i]] <- boot_fun(df, n)
}
l <- length(boot_smrq_models[[1]])
smrq_models_all <- vector("list", length = l)
smrq_models_int <- vector("list", length = l)
for(i in seq.int(l)) {
tmp <- array(dim = c(R, dim(boot_smrq_models[[1]][[i]])))
for(j in seq.int(R)) {
tmp[j, , ] <- boot_smrq_models[[j]][[i]]
}
smrq_models_all[[i]] <- t(apply(tmp, 2:3, mean))
smrq_models_int[[i]] <- apply(tmp, 2:3, quantile, probs = c(0.025, 0.975))
rownames(smrq_models_all[[i]]) <- sprintf("tau_%03.02f", (1:99)/(99+1))
}
lapply(smrq_models_all, head)
#> [[1]]
#>                [,1]       [,2]
#> tau_0.01 -23.173185 -1.4624225
#> tau_0.02   3.509690 -0.6640458
#> tau_0.03 -13.748280 -1.2172158
#> tau_0.04 -38.373045 -1.9780045
#> tau_0.05  -8.166044 -1.0890235
#> tau_0.06  51.701066  0.6258957
#> 
#> [[2]]
#>               [,1]      [,2]
#> tau_0.01  1.387311 -2.315114
#> tau_0.02 36.060786 -1.311070
#> tau_0.03  9.177731 -2.111396
#> tau_0.04 22.320567 -1.731685
#> tau_0.05 40.143578 -1.212108
#> tau_0.06 19.900275 -1.812859

^{在2022-05-22由reprex包(v2.0.1)创建}