2016-12-07 10 views
0

lmerモデルから 'terms'、特に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スプラインを抽出する方法については考えているましたモデルからオブジェクトを抽出し、新しいデータに「予測」します。同じことが、例えばポリのような変換された変数にも適用されます。 poly(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

を行います
関連する問題