Интерпретация вывода кусочно-смешанных эффектов
Я пытаюсь понять итоговый вывод из кусочно-смешанной модели эффектов и мог бы использовать некоторые идеи. В частности, я хочу знать, как получить перехваты и наклоны регрессии для линии слева и справа от точки останова. Из того, что я понимаю, пересечение, данное в выходных данных ниже, относится к линии регрессии слева от точки останова, а значение, данное для I(Days * (Days < 6.07)), представляет собой наклон этой линии. Тем не менее, я не думаю, что я (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 (Дни * (Дни <6.07)) + + I (Дни * (Дни>= 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), data = 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
Может ли кто-нибудь сказать мне, как интерпретировать вывод lmer(), когда I() используется в предикторах модели, или подсказать, как запустить модель lmer() без I() в предикторах модели?
Я ценю любые доступные рекомендации, так как я не смог найти ничего на страницах справки R для этого!
Спасибо.
1 ответ
Я думаю, что вы можете получить то, что вы хотите, следующим образом:
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)