rms | 如何绘制模型带置信区间的预测曲线

最新推荐文章于 2024-07-18 00:30:00 发布

R语言学堂

最新推荐文章于 2024-07-18 00:30:00 发布

阅读量3.9k

点赞数 1

文章标签： python java 机器学习 javascript ajax

本文链接：https://blog.csdn.net/weixin_54000907/article/details/122389342

版权

前段时间，有读者咨询如何绘制带置信区间的预测曲线。本篇介绍的rms工具包可以很简便地解决这个问题。

该包是Regression Modeling Strategies的附加学习资源，小编从网上找到了这本书的目录：https://link.springer.com/content/pdf/bfm%3A978-3-319-19425-7%2F1.pdf。

1 `lm()`/`ols()`

在rms工具包中，使用Predict()函数预测的模型结果可以直接使用plot()或ggplot()函数进行绘制，并自带置信区间，但是建模函数必须使用rms工具包中的相关函数。比如对线性模型而言，基础包stats中的函数是lm()，而rms工具包中是ols()函数。

ols()函数的语法结构如下：

ols(formula, data=environment(formula),
    weights, subset, na.action=na.delete, 
    method="qr", model=FALSE,
    x=FALSE, y=FALSE, se.fit=FALSE, linear.predictors=TRUE,
    penalty=0, penalty.matrix, tol=1e-7, sigma,
    var.penalty=c('simple','sandwich'), ...)

可以看出，ols()函数和lm()函数的用法基本上是一致的。

library(rms)
model.1 <- lm(mpg ~ wt, data = mtcars)
model.2 <- ols(mpg ~ wt, data = mtcars)

summary(model.1)
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10
print(model.2)
## Linear Regression Model
##  
##  ols(formula = mpg ~ wt, data = mtcars)
##  
##                 Model Likelihood     Discrimination    
##                    Ratio Test           Indexes        
##  Obs      32    LR chi2     44.73    R2       0.753    
##  sigma3.0459    d.f.            1    R2 adj   0.745    
##  d.f.     30    Pr(> chi2) 0.0000    g        5.822    
##  
##  Residuals
##  
##      Min      1Q  Median      3Q     Max 
##  -4.5432 -2.3647 -0.1252  1.4096  6.8727 
##  
##  
##            Coef    S.E.   t     Pr(>|t|)
##  Intercept 37.2851 1.8776 19.86 <0.0001 
##  wt        -5.3445 0.5591 -9.56 <0.0001

两个函数的模型参数估算结果也是一样的。

2 `predict()`/`Predict()`

在stats包中，与模型预测相关的有两个函数：fitted()和predict()。前者是根据模型结果对原数据中的样本的响应变量（因变量）进行预测：

head(fitted(model.1))
##         Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive 
##          23.28261          21.91977          24.88595          20.10265 
## Hornet Sportabout           Valiant 
##          18.90014          18.79325

predict()函数则可以对新数据的样本进行预测：

## S3 method for class 'lm'
predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf,
        interval = c("none", "confidence", "prediction"),
        level = 0.95, type = c("response", "terms"),
        terms = NULL, na.action = na.pass,
        pred.var = res.var/weights, weights = 1, ...)

object：模型对象；
newdata：新数据；必须为数据框结构；若省略，则默认为原数据；
interval：区间类型，有两种类型：置信区间（confidence）和预测区间（prediction）；
level：置信水平；默认为95%。

关于置信区间和预测区间的区别，小编不是很清楚，有兴趣的可以参见这个链接的解释：https://blog.sciencenet.cn/blog-54276-288414.html。

predict(model.1, data.frame(wt = seq(1,6,0.5)),
        interval = "prediction")
##          fit       lwr      upr
## 1  31.940655 25.135231 38.74608
## 2  29.268419 22.654123 35.88271
## 3  26.596183 20.128114 33.06425
## 4  23.923947 17.554110 30.29378
## 5  21.251711 14.929874 27.57355
## 6  18.579476 12.254262 24.90469
## 7  15.907240  9.527355 22.28712
## 8  13.235004  6.750452 19.71956
## 9  10.562768  3.925916 17.19962
## 10  7.890533  1.056933 14.72413
## 11  5.218297 -1.852790 12.28938

predict(model.1, data.frame(wt = seq(1,6,0.5)),
        interval = "confidence")
##          fit      lwr       upr
## 1  31.940655 29.18042 34.700892
## 2  29.268419 27.02030 31.516535
## 3  26.596183 24.82389 28.368481
## 4  23.923947 22.55284 25.295059
## 5  21.251711 20.12444 22.378987
## 6  18.579476 17.43342 19.725534
## 7  15.907240 14.49018 17.324295
## 8  13.235004 11.40347 15.066543
## 9  10.562768  8.24913 12.876406
## 10  7.890533  5.06154 10.719526
## 11  5.218297  1.85595  8.580644

而在rms工具包中，模型预测对应的是Predict()函数（首字母大写）：

Predict(object, ..., fun=NULL, funint=TRUE,
        type = c("predictions", "model.frame", "x"),
        np = 200, conf.int = 0.95,
        conf.type = c("mean", "individual","simultaneous"),
        usebootcoef=TRUE, boot.type=c("percentile", "bca", "basic"),
        posterior.summary=c('mean', 'median', 'mode'),
        adj.zero = FALSE, ref.zero = FALSE,
        kint=NULL, ycut=NULL, time = NULL, loglog = FALSE,
        digits=4, name,
        factors=NULL, offset=NULL)

conf.type：区间类型；置信区间（mean）、预测区间（individual）。

Predict()函数的“新数据”使用...表示，直接写出参与预测的变量即可：

predict.2 <- Predict(model.2, wt = seq(1,6,0.5)) 
predict.2
##     wt      yhat    lower     upper
## 1  1.0 31.940655 29.18042 34.700892
## 2  1.5 29.268419 27.02030 31.516535
## 3  2.0 26.596183 24.82389 28.368481
## 4  2.5 23.923947 22.55284 25.295059
## 5  3.0 21.251711 20.12444 22.378987
## 6  3.5 18.579476 17.43342 19.725534
## 7  4.0 15.907240 14.49018 17.324295
## 8  4.5 13.235004 11.40347 15.066543
## 9  5.0 10.562768  8.24913 12.876406
## 10 5.5  7.890533  5.06154 10.719526
## 11 6.0  5.218297  1.85595  8.580644
## 
## Response variable (y): mpg 
## 
## Limits are 0.95 confidence limits

可以看出，predict.2与predict()函数预测的置信区间一致。

也可以不指定参与预测变量的数值，但需要提前定义变量的范围和等分的个数。指定范围的方法如下：

ddist <- datadist(mtcars)
options(datadist = "ddist")

等分数目由np参数确定，默认值为200：

Predict(model.2, wt, np = 10) 
##          wt      yhat     lower    upper
## 1  1.736000 28.007124 25.989733 30.02451
## 2  2.131194 25.895018 24.237593 27.55244
## 3  2.526389 23.782913 22.429583 25.13624
## 4  2.921583 21.670807 20.520507 22.82111
## 5  3.316778 19.558702 18.453202 20.66420
## 6  3.711972 17.446596 16.210345 18.68285
## 7  4.107167 15.334491 13.837242 16.83174
## 8  4.502361 13.222385 11.388690 15.05608
## 9  4.897556 11.110280  8.898861 13.32170
## 10 5.292750  8.998174  6.385598 11.61075
## 
## Response variable (y): mpg 
## 
## Limits are 0.95 confidence limits

3 `plot()`/`ggplot()`

Predict()函数预测的结果可以直接使用plot()函数和ggplot()函数进行绘制。在加载rms工具包时会自动加载ggplot2工具包。

需要注意的是，rms包中的plot()函数采用的是lattice绘图系统，而非基础绘图系统（graphics），因此不能叠加graphics包中的函数。

plot(predict.2)

以下代码会报错：

plot(predict.2)
points(mtcars$wt, mtcars$mpg)

Error in plot.xy(xy.coords(x, y), type = type, ...) : plot.new has not been called yet

在ggplot2绘图系统中，不需要额外使用几何图形函数：

ggplot(predict.2)

还可以使用图层叠加功能：

ggplot(predict.2) +
  geom_point(data = mtcars, aes(wt, mpg))

4 其他情况的简要说明

如果要向上述模型中添加二次项，在lm()函数中可以使用I(wt^2)或poly(wt,2)的形式，但若想使用Predict()函数对模型进行预测，就必须使用rms工具包中的函数进行建模，该包添加二次项的函数为pol()。

以下代码会报错：

ols(mpg ~ wt + I(wt^2), data = mtcars)

Error in if (!length(fname) || !any(fname == zname)) { : 需要TRUE/FALSE值的地方不可以用缺少值

以下为正确的写法：

model.3 <- ols(mpg ~ pol(wt,2), data = mtcars)
predict.3 <- Predict(model.3, wt = seq(1,6,0.5)) 
predict.3

ggplot(predict.3) +
  geom_point(data = mtcars, aes(wt, mpg))

如果模型中有多个自变量，可以使用全部变量进行预测，也可以只使用个别变量进行预测。

使用全部变量进行预测：

ddist <- datadist(mtcars)
options(datadist = "ddist")
model.4 <- ols(mpg ~ wt + drat, data = mtcars)
predict.41 <- Predict(model.4, np = 10)
plot(predict.41)

使用单个变量进行预测:

predict.42 <- Predict(model.4, wt = seq(1,6,0.5))
plot(predict.42)

未参与预测的变量实际上取的是中位数：

median(mtcars$drat)
## [1] 3.695

也可以指定未参与预测变量的取值：

predict.43 <- Predict(model.4, wt = seq(1,6,0.5),
                      drat = 4,
                      np = 10)
plot(predict.43)

rms包中的模型函数基本上是自成体系的，如logistic回归用lrm()函数，以及一系列的样条曲线函数，如rcs()函数（linear tail-restricted cubic spline function）。

ddist <- datadist(mtcars)
options(datadist = "ddist")
model.5 <- ols(mpg ~ rcs(wt) + drat, data = mtcars)
predict.5 <- Predict(model.5, wt)
plot(predict.5)