2016-12-07 11 views
0

lmer 모델에서 '용어', 특히 ns 스플라인을 예측하려고합니다. mtcars 데이터 세트 (기술적으로 좋지 않은 예제이지만 문제를 해결하기 위해 노력)로 문제를 재현했습니다.R - lmer 모델에서 ns 스플라인 오브젝트를 추출하고 새 데이터를 예측합니다.

data(mtcars) 
mtcarsmodel <- lm(wt ~ ns(drat,2) + hp + as.factor(gear), data= mtcars) 
summary(mtcarsmodel) 
coef(mtcarsmodel) 
test <- predict(mtcarsmodel, type = "terms") 

완벽한을 : 여기

내가 선형 모델을 할 노력하고있어 것입니다. 그러나 lmer predict (unresolved issue here)에 해당하는 'terms'옵션은 없습니다. 동등한 존재하지 감안할 때

mtcarsmodellmer <- lmer(wt ~ ns(drat,2) + (hp|as.factor(gear)), data= mtcars) 
summary(mtcarsmodellmer) 
coef(mtcarsmodellmer) 
ranef(mtcarsmodellmer) 

'예측은, 용어'기능, 나는 위의 고정 및 임의 계수를 추출하고 mtcars 데이터에 계수를 적용하려고하지만, NS 스플라인을 추출하는 방법에 대한 아무 생각이 없었다 모델에서 객체를 가져 와서 새로운 데이터로 '예측'합니다. 같은 '폴리'변환 된 변수 예를 들어갑니다. 폴리 (drat, 2) - 이것을 얻을 수 있다면 특별한 명성.

답변

1

직접하는 것은 어렵지 않습니다.

library(lme4) 
library(splines) 
X <- with(mtcars, ns(drat, 2)) ## design matrix for splines (without intercept) 
## head(X) 
#    1   2 
#[1,] 0.5778474 -0.1560021 
#[2,] 0.5778474 -0.1560021 
#[3,] 0.5738625 -0.1792162 
#[4,] 0.2334329 -0.1440232 
#[5,] 0.2808520 -0.1704002 
#[6,] 0.0000000 0.0000000 

## str(X) 
# ns [1:32, 1:2] 0.578 0.578 0.574 0.233 0.281 ... 
# - attr(*, "dimnames")=List of 2 
# ..$ : NULL 
# ..$ : chr [1:2] "1" "2" 
# - attr(*, "degree")= int 3 
# - attr(*, "knots")= Named num 3.7 
# ..- attr(*, "names")= chr "50%" 
# - attr(*, "Boundary.knots")= num [1:2] 2.76 4.93 
# - attr(*, "intercept")= logi FALSE 
# - attr(*, "class")= chr [1:3] "ns" "basis" "matrix" 

fit <- lmer(wt ~ X + (hp|gear), data= mtcars) 

beta <- coef(fit) 
#$gear 
#   hp (Intercept)  X1   X2 
#3 0.010614406 2.455403 -2.167337 -0.9246454 
#4 0.014601363 2.455403 -2.167337 -0.9246454 
#5 0.006342761 2.455403 -2.167337 -0.9246454 
# 
#attr(,"class") 
#[1] "coef.mer" 

우리가 ns 기간을 예측하려면, 단지

## use `predict.ns`; read `?predict.ns` 
x0 <- seq(1, 5, by = 0.2) ## example `newx` 
Xp <- predict(X, newx = x0) ## prediction matrix 
b <- with(beta$gear, c(X1[1], X2[1])) ## coefficients for spline 
y <- Xp %*% b ## predicted mean 

plot(x0, y, type = "l") 

enter image description here