R中的LMS（Lambda-Mu-Sigma）方法

Question

我想使用 LMS (Lambda-Mu-Sigma) 方法为我的数据创建百分位曲线。我有以下示例数据。如何使用 LMS 绘制 yvar（y 轴）与年龄（x 轴）的 10th、50th 和 90th 百分位曲线？

age = sample(5:75, 500, replace=T) 
yvar = rnorm(500, age, 20)
mydata = data.frame(age, yvar)
head(mydata)
  age     yvar
1  61 87.16011
2  58 49.73289
3  65 15.60212
4  71 83.32699
5  33 40.89592
6  18 25.04376

plot(age, yvar)

我遇到了 VGAM 包 http://www.inside-r.org/packages/cran/VGAM/docs/lms.bcn 。这是最好的方法吗？我无法真正理解它的示例代码来从上述数据创建简单的百分位曲线。感谢您的帮助。

Answer 1

模拟数据（可重复）：

set.seed(1001)
mydata <- data.frame(
   age = sample(5:75, 500, replace=TRUE))
mydata <- transform(mydata,
   yvar = rnorm(500, age, 20))

由于LMS 方法通常似乎基于Box-Cox transformation 的变体，这需要正值，因此更简单的方法是使用quantile regression。

library("quantreg")
library("ggplot2"); theme_set(theme_bw())
g0 <- ggplot(mydata,aes(x=age,y=yvar))+geom_point()
g0 + geom_smooth(method="rq",tau=c(0.1),se=FALSE,lty=2)+
      geom_smooth(method="rq",tau=c(0.5),se=FALSE)+
      geom_smooth(method="rq",tau=c(0.9),se=FALSE,lty=2)

rq() 本身具有同时拟合所有三个百分位数的能力，但您需要使用this blog post 中建议的策略来更方便地绘制它们：

model.rq <- rq(yvar ~ age, mydata, tau=c(0.1, 0.5, 0.9))
quantile.regressions <- data.frame(t(coef(model.rq)))
colnames(quantile.regressions) <- c("intercept", "slope")
quantile.regressions$quantile <- rownames(quantile.regressions)
g0 + geom_abline(aes(intercept=intercept, slope=slope,
   colour=quantile), show_guide=TRUE, data=quantile.regressions)

或者 可以在 VGAM 中执行此操作，但我不确定这是否是您想要的/结果是否有意义。通过 lms.yjn 进行的 Yeo-Johnson 转换允许您在某些数据值为负数的情况下执行此操作，但您可以查看 ?lms.bcg、?lms.bcn 以了解适用于非负数据的替代方案。

library("VGAM")
fit <- vgam(yvar ~ s(age, df = 4), lms.yjn, data=mydata,
            control=vgam.control(maxit=100),
            trace=FALSE)

我们收到一条警告消息：

## Warning message:
## In vgam.fit(x = x, y = y, w = w, mf = mf, Xm2 = Xm2, Ym2 = Ym2,  :
##  convergence not obtained in 100 iterations

这可能是因为我们使用 4 节样条模型过度拟合数据？

分位数图（跟随example("lms.yjn")）

par(bty = "l", mar = c(5, 4, 4, 3) + 0.1, xpd = TRUE)
qtplot(fit, percentiles = c(10, 50, 90),
        las = 1, ylab = "yvar", lwd = 2, lcol = 4)

这是一个可怕的 hack，但如果您想访问原始值以便自己绘制曲线：

pcurves <- qtplot.lmscreg(fit,show.plot=FALSE,
                          percentiles=c(10,50,90))
vals <- data.frame(age=mydata$age,pcurves$fitted.values)
vals <- vals[order(vals$age),]
matplot(vals$age,vals[,-1],type="l",lty=c(2,1,2),col=1,
        xlab="age",ylab="")