lme 线性 b 样条模型的斜率答案

【问题标题】：Slopes for lme linear b-spline modellme 线性 b 样条模型的斜率
【发布时间】：2014-07-15 04:46:27
【问题描述】：

我想知道对于使用线性 b 样条的 lme 模型，如何获得每个段的 SE 和 p 值的斜率估计。我可以使用预测得到斜率估计，但不能使用 SE 和 p 值。

这是一个例子：

rm(list = ls())
library(splines)
library(nlme)

getY <- function(x) ifelse(x < 7, x * 1.3, x * 0.6) + rnorm(length(x))

set.seed(123)

data <- data.frame(Id = numeric(0), X = numeric(0), Y = numeric(0))
for (i in 1:10) {
    X <- sample(1:10, 4)
    Y <- getY(X) + rnorm(1, 0.5)
    Id <- rep(i, 4)
    data <- rbind(data, cbind(Id = Id, X = X, Y = Y))
}

gdata <- groupedData(Y ~ X | Id, data)

mod <- lme(fixed = Y ~ bs(X, degree = 1, knots = 7), data = gdata, random = ~1 | 
Id)

summary(mod)

Linear mixed-effects model fit by REML
 Data: gdata 
    AIC   BIC logLik
  158.2 166.2 -74.09

Random effects:
 Formula: ~1 | Id
        (Intercept) Residual
StdDev:       1.217    1.389

Fixed effects: Y ~ bs(X, degree = 1, knots = 7) 
                              Value Std.Error DF t-value p-value
(Intercept)                   3.098    0.5817 28   5.326   0e+00
bs(X, degree = 1, knots = 7)1 4.031    0.7714 28   5.225   0e+00
bs(X, degree = 1, knots = 7)2 3.253    0.7258 28   4.481   1e-04
 Correlation: 
                              (Intr) b(X,d=1,k=7)1
bs(X, degree = 1, knots = 7)1 -0.597              
bs(X, degree = 1, knots = 7)2 -0.385  0.233       

Standardized Within-Group Residuals:
      Min        Q1       Med        Q3       Max 
-1.469915 -0.628202  0.005586  0.541398  1.748387 

Number of Observations: 40
Number of Groups: 10


plot(augPred(mod))

pred1 <- predict(mod, data.frame(X = 1:2), level = 0)
pred2 <- predict(mod, data.frame(X = 8:9), level = 0)

(slope1 <- diff(pred1))

10.6718

(slope2 <- diff(pred2))

1-0.2594

【问题讨论】：

标签： r

【解决方案1】：

你不会只考虑预测结果的差异吗？

predict(mod, newdata=data.frame(X=1:10, Id=1) )
       1        1        1        1        1        1        1        1        1 
3.449572 4.121362 4.793152 5.464941 6.136731 6.808521 7.480311 7.220928 6.961544 
       1 
6.702161 
attr(,"label")
[1] "Predicted values"

So:

 plot( predict(mod, newdata=data.frame(X=1:10, Id=1) ), ylim=c(-2,8))
 lines( 1:9, diff(predict(mod, newdata=data.frame(X=1:10, Id=1) ), ylim=c(-2,8)) )

【讨论】：

正如您在我的代码中看到的那样，我使用 predict 得到了斜率，但我也想得到标准误差和 p 值。