结合LOESS和Quantreg计算数据的百分位数/分位数

Question

我正在尝试计算相当分散的数据的百分位数或分位数。使用Loess函数可以很好地表示平均值，但是，我无法从该函数中获得百分位数/分位数。

我试着把quantreg和loess结合起来。此图显示的是线性曲线，而不是黄土平滑曲线。

我希望得到类似如下的结果：

​

data(cars)
plot(cars)
lmodel <- loess(cars$dist~cars$speed,span = 0.3, degree = 1)
lpred<-predict(lmodel, newdata= 5:25,se=TRUE)
lines(5:25, lpred$fit,col='#000066',lwd=4)
lines(5:25, lpred$fit - qt(0.975, lpred$df)*lpred$se, lty=2)
lines(5:25, lpred$fit + qt(0.975, lpred$df)*lpred$se, lty=2)


#### combination of quantreg with loess

plot(cars$speed,cars$dist)
xx <- seq(min(cars$speed),max(cars$speed),1)
f <- coef(rq(loess(cars$dist~cars$speed,span = 0.3, degree = 1), tau=c(0.1,0.25,0.5,0.75,0.9)) )
yy <- cbind(1,xx)%*%f
for(i in 1:length(taus)){
  lines(xx,yy[,i],col = "gray")
}

我也尝试了建议的代码，但是，我无法更改平滑的设置。线条显示为波状路径。

library(quantreg)
data(cars)
taus <- c(0.1, 0.25, 0.5, 0.75, 0.9)
lmodel <- loess(dist ~ speed, data = cars, span = 0.9, degree = 1)
rqmodel <- rq(lmodel, tau = taus, data = cars)
f <- coef(rqmodel)
xx <- seq(min(cars$speed), max(cars$speed), length.out = nrow(cars))
yy <- predict(rqmodel)
plot(cars)
matlines(xx, yy, col = "grey",lwd=3)

​

​

Loess函数不像rg那样为分位数提供数据。

然而，黄土函数允许得到一条没有Z字形的曲线。请看代码片段。使用rg函数的tau=0.5的设置是什么，才能产生与Loess函数相同的结果。

data(cars)
lmodel <- loess(dist ~ speed, data = cars, span = 0.9 )
plot(cars)
lines( x=4:25 , y=predict(lmodel, newdata= data.frame(speed=4:25)) ,col="Blue")

​

​

Answer 1

我相信问题中的代码是混合黄土和分位数回归，当它们是不同的方法时，后者不需要前者。

我将尝试将两者都拟合，并绘制出各自的结果。在下面的代码中，我将使用matlines，而不是for循环。

这些代码行很常见。

library(quantreg)

data(cars)

xx <- seq(min(cars$speed), max(cars$speed), length.out = nrow(cars))

首先是黄土模型。

lmodel <- loess(dist ~ speed, data = cars, span = 0.5, degree = 1)
ls_yy <- predict(lmodel, se = TRUE)
ls_yy <- cbind(ls_yy$fit, 
               ls_yy$fit - 2*ls_yy$se.fit, 
               ls_yy$fit + 2*ls_yy$se.fit)

plot(cars)
matlines(xx, ls_yy, col = "darkgrey")

​

​

现在是分位数回归。

taus <- c(0.1, 0.25, 0.5, 0.75, 0.9)
rqmodel <- rq(dist ~ speed, tau = taus, data = cars)

rq_yy <- predict(rqmodel)

plot(cars)
matlines(xx, rq_yy, col = "darkgrey")

​

​

Answer 2

下面的代码(取自“答案”)不正确，不应该包含在正确的解决方案中。这将提供拟合的95%置信区间，以及区间落在真实趋势线上的概率。它不对应于从该移动平均值范围内的数据计算的分位数。建议的基于正态分布的近似需要将ls_yy$se.fit乘以sqrt( ni )，其中ni是特定跨度中的观测值数量。不幸的是，loess不返回ni，因此这不是一个可行的解决方案，除非跨度覆盖整个数据集，并且ni可以设置为n，并且没有异方差。

data(cars)
plot(cars)

lmodel <- loess(dist ~ speed, data = cars, span = 0.5, degree = 1)
ls_yy <- predict(lmodel, se = TRUE)

#wrong - this does not denote quantiles for the input data:
ls_yy <- cbind(ls_yy$fit, 
               ls_yy$fit - 2*ls_yy$se.fit, 
               ls_yy$fit + 2*ls_yy$se.fit)
plot(cars)
matlines(xx, ls_yy, col = "darkgrey")

我们可以使用具有更多观察值的样本数据集来使这一点变得更加明显。除了样本大小(500和1500个观察值)之外，样本1和2是相同的，因此它们应该具有非常相似的分位数。

set.seed(1)
x1 = runif(500,0,10)
y1 = x1 + rnorm(length(x1))

x2 = runif(1500,0,10)
y2 = x1 + rnorm(length(x2))

dfpd = data.frame(x=1:9)

lmodel1 <- loess(y ~ x, data = data.frame(x=x1,y=y1), span = 0.5, degree = 1)
ls_yy1 <- predict(lmodel1, newdata=dfpd, se = TRUE)

lmodel2 <- loess(y ~ x, data = data.frame(x=x2,y=y2), span = 0.5, degree = 1)
ls_yy2 <- predict(lmodel2, newdata=dfpd, se = TRUE)

#the only difference between lmodel1 and lmodel2 is the number of observations
#the quantiles should be very similar, but their se values are a function of sample
#size and are thus quite different
ls_yy1$se
ls_yy2$se


ls_yy1$se / ls_yy2$se

我们可以看到，se值的比例约为60%，这证实了它们不能按原样用于分位数计算。

​

​