1

我正在尝试了解分段混合效果模型的摘要输出，并且可以使用一些见解。具体来说，我想知道如何获得断点左右线的回归截距和斜率。据我了解，下面输出中给出的截距是断点左侧的回归线，而 I(Days * (Days < 6.07)) 给出的值是该线的斜率。但是，我不认为 I(Days * (Days >= 6.07)) 是断点右侧线的斜率，也不是两个斜率的差异。

library(lme4)
sleepstudy<-as.data.frame(sleepstudy)

我从上一个线程中提取了断点：https ://stats.stackexchange.com/questions/19772/estimating-the-break-point-in-a-broken-stick-piecewise-linear-model-with-rando

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ I(Days * (Days < 6.07)) + I(Days * (Days >= 6.07)) +      (1 | Subject) 
   Data: sleepstudy 

REML criterion at convergence: 1784.369 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1377.6   37.12   
 Residual              965.7   31.08   
Number of obs: 180, groups: Subject, 18

Fixed effects:
                         Estimate Std. Error t value
(Intercept)              252.2663    10.0545  25.090
I(Days * (Days < 6.07))   10.0754     1.3774   7.315
I(Days * (Days >= 6.07))  10.4513     0.8077  12.940

Correlation of Fixed Effects:
            (Intr) I(*(<6
I(D*(D<6.07 -0.409       
I(D*(D>=6.0 -0.374  0.630

我试图通过消除随机效应来简化：当 I() 包含在 lm 模型中时，斜率/截距与上面的混合模型非常相似，我仍然感到困惑。

mod_lm<-lm(反应 ~ I(Days*(Days < 6.07)) + I(Days*(Days>= 6.07)), data = sleepstudy) 摘要(mod_lm)

Call:
lm(formula = Reaction ~ I(Days * (Days < 6.07)) + I(Days * (Days >= 
    6.07)), data = sleepstudy)

Residuals:
     Min       1Q   Median       3Q      Max 
-111.581  -27.632    1.614   26.994  141.443 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               252.266      7.629  33.066  < 2e-16 ***
I(Days * (Days < 6.07))    10.075      2.121   4.751 4.17e-06 ***
I(Days * (Days >= 6.07))   10.451      1.243   8.405 1.37e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 47.84 on 177 degrees of freedom
Multiple R-squared:  0.2867,    Adjusted R-squared:  0.2786 
F-statistic: 35.57 on 2 and 177 DF,  p-value: 1.037e-13

但是，当 I() 从 lm 公式中删除时，我理解输出，并且结果是有意义的。

mod_lm<-lm（反应〜天*（天< 6.07）+天*（天> = 6.07），数据= sleepstudy）摘要（mod_lm）

Call:
lm(formula = Reaction ~ Days * (Days < 6.07) + Days * (Days >= 
    6.07), data = sleepstudy)

Residuals:
     Min       1Q   Median       3Q      Max 
-114.214  -27.833    0.603   27.254  141.693 

Coefficients: (2 not defined because of singularities)
                      Estimate Std. Error t value Pr(>|t|)   
(Intercept)            207.008     64.211   3.224  0.00151 **
Days                    16.050      7.985   2.010  0.04595 * 
Days < 6.07TRUE         45.908     64.671   0.710  0.47872   
Days >= 6.07TRUE            NA         NA      NA       NA   
Days:Days < 6.07TRUE    -6.125      8.265  -0.741  0.45965   
Days:Days >= 6.07TRUE       NA         NA      NA       NA   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 47.91 on 176 degrees of freedom
Multiple R-squared:  0.2887,    Adjusted R-squared:  0.2766 
F-statistic: 23.81 on 3 and 176 DF,  p-value: 5.526e-13

当 I() 项从 lmer 公式中删除时， lmer 将不会运行。

mod1<-lmer(Reaction ~ Days*(Days < 6.07) + Days*(Days>= 6.07) + (1|Subject), data = sleepstudy)
Error in lme4::lFormula(formula = Reaction ~ Days * (Days < 6.07) + Days *  : 
  rank of X = 4 < ncol(X) = 6

有人可以告诉我在模型预测器上使用 I() 时如何解释 lmer() 输出，或者告诉我如何在没有模型预测器上的 I() 的情况下运行 lmer() 模型？

我很感激任何可用的指导，因为我无法在 R 帮助页面上找到任何关于这个的指导！

谢谢你。

4

1 回答 1

5

我认为你可以得到你想要的如下：

library(lme4)
sleepstudy <- transform(sleepstudy,period=(Days<6.5))
(m0 <- lmer(Reaction ~ Days+ (1 | Subject), sleepstudy))
(m2 <- lmer(Reaction ~ Days*period+ (1 | Subject), sleepstudy))
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days * period + (1 | Subject) 
##    Data: sleepstudy 
## REML criterion at convergence: 1773.86 
## Random effects:
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 37.12   
##  Residual             31.06   
## Number of obs: 180, groups: Subject, 18
## Fixed Effects:
##     (Intercept)             Days       periodTRUE  Days:periodTRUE  
##         207.008           16.050           45.908           -6.125

您的结果I()是构造数字变量而不是分类变量（转换为虚拟变量）。也许您感到困惑的主要原因是您的第一组模型不允许按周期进行单独的截距，只允许单独的斜率......

lmer对你的第二组模型不起作用的原因是它lmer不像现在那样容忍过度参数化（多重共线预测变量）lm，尽管开发版本（在 Github 上可用，很快就会发布）是：如果你运行mod1它将拟合模型并打印一条消息“固定效应模型矩阵秩不足，因此删除 2 列/系数”（与 不同lm，它不保留删除的列和NA系数，只是完全删除它们）。

更新：

sleepstudy <- transform(sleepstudy,cDays=Days-6.5)
m3 <- lmer(Reaction ~ cDays:period+ (1 | Subject), sleepstudy)
library(ggplot2); theme_set(theme_bw())    
library(reshape2)
g0 <- ggplot(sleepstudy,aes(Days,Reaction,group=Subject))+geom_line()
pframe <- data.frame(Days=seq(0,8,length=101))
pframe <- transform(pframe,cDays=Days-6.5,period=Days>6.5)
## next line assumes latest version of lme4 -- you may need REform instead
pframe$Reaction <- predict(m3,newdata=pframe,re.form=NA)
pframe$Reaction2 <- predict(m0,newdata=pframe,re.form=NA)

很难看出斜率的差异——非常微妙。

g0 + geom_line(data=pframe,colour=2,aes(group=NA))+
     geom_line(data=pframe,colour=2,lty=2,
         aes(y=Reaction2,group=NA))+
     geom_vline(xintercept=6.5,lty=2)
于 2014-03-11T02:47:44.243 回答