2

指定された時間に予測平均を取得したい経時データがあります。モデルには、2 つの項、それらの相互作用、および時間変数のスプライン項が含まれます。予測平均を取得しようとすると、「mm %*% fixef(m4) のエラー: 非適合引数」が表示されます。

私の問題を説明するために、lmer の睡眠データセットを使用しました。まず、データをインポートし、インタラクション用の変数「age」を作成します

sleep <- as.data.frame(sleepstudy)  #get the sleep data
# create fake variable for age with 3 levels
set.seed(1234567)
sleep$age <- as.factor(sample(1:3,length(sleep),rep=TRUE))

次に、lmer モデルを実行します

library(lme4)
library(splines)
m4 <- lmer(Reaction ~ Days + ns(Days, df=4) + age + Days:age + (Days | Subject), sleep)

最後に、予測平均を取得するために必要なデータと行列を作成します

#new data frame for predicted means
d <- c(0:9)  # make a vector of days = 0 to 9 to obtain predictions for each day
newdat <- as.data.frame(cbind(Days=d, age=rep(c(1:3),length(d))))
newdat$Days <- as.numeric(as.character(newdat$Days))
newdat$age <- as.factor(newdat$age)

# create a matrix 
mm<-model.matrix(~Days + ns(Days, df=4) + age + Days:age, newdat)  
newdat$pred<-mm%*%fixef(m4)

エラーが発生するのはこの時点です: Error in mm %*% fixef(m4) : non-conformable arguments

予測を使用して平均を取得できます

newdat$pred <- predict(m4, newdata=newdat, re.form=NA)

これは正常に機能しますが、信頼区間を計算できるようにしたいので、適合する行列が必要です。

問題は lmer がエイリアスを作成することである可能性があることをどこかで読みました（その投稿が見つかりません）。このコメントは、同様のタスクに effect() を使用できないことに関して行われました。この問題を克服する方法がよくわかりませんでした。さらに、投稿が少し古かったことを思い出し、エイリアスの問題がもはや関係ないことを願っていました。

私が間違っている可能性があることについて誰かが提案をしている場合は、フィードバックをいただければ幸いです。ありがとう。

4

2 に答える 2

1

ここにはいくつかのことがあります。

  • 列をドロップして、実際にフィッティングされた固定効果ベクトルに対応するモデル マトリックスを作成する必要があります (つまり、共線列をドロップした後、実際にフィッティングに使用されたモデル マトリックスに対応します)。
  • 追加の混乱のために、たまたまサンプル年齢が 2 歳と 3 歳でした (可能性のある {1,2,3} のうち)。

コードを少し整理しました...

library("lme4")
library("splines")
sleep <- sleepstudy  #get the sleep data
set.seed(1234567)
## next line happens to sample only 2 and 3 ...
sleep$age <- as.factor(sample(1:3,length(sleep),rep=TRUE))
length(levels(sleep$age))  ## 2

フィットモデル:

m4 <- lmer(Reaction ~ Days + ns(Days, df=4) +
    age + Days:age + (Days | Subject), sleep)
## message; fixed-effect model matrix is 
##    rank deficient so dropping 1 column / coefficient

固定効果を確認する:

f1 <- fixef(m4)
length(f1)  ## 7
f2 <- fixef(m4,add.dropped=TRUE)
length(f2)  ## 8

この拡張バージョンの固定効果 (値が含まれる) を使用することもできますが、これは計算によって値をNA伝播することで混乱を招くだけです...NA

モデル マトリックスを確認します。

X <- getME(m4,"X")
ncol(X)  ## 7
(which.dropped <- attr(getME(m4,"X"),"col.dropped"))
## ns(Days, df = 4)4 
##             6

予測平均の新しいデータ フレーム

d <- 0:9  
## best to use data.frame() directly, avoid cbind()
##   generate age based on *actual* levels in data
newdat <- data.frame(Days=d,
   age=factor(rep(levels(sleep$age),length(d))))

マトリックスを作成します。

mm <- model.matrix(formula(m4,fixed.only=TRUE)[-2], newdat)
mm <- mm[,-which.dropped]   ## drop redundant columns
## newdat$pred <- mm%*%fixef(m4)    ## works now

sianagh によって追加されました: 信頼区間を取得してデータをプロットするためのコード:

predFun <- function(x) predict(x,newdata=newdat,re.form=NA)
newdat$pred <- predFun(m4)
bb <- bootMer(m4,
   FUN=predFun,
    nsim=200)  
## nb. this produces an error message on its first run, 
## but not on subsequent runs (using the development version of lme4)
bb_ci <- as.data.frame(t(apply(bb$t,2,quantile,c(0.025,0.975))))
names(bb_ci) <- c("lwr","upr")
newdat <- cbind(newdat,bb_ci)

プロット：

plot(Reaction~Days,sleep)
with(newdat,
    matlines(Days,cbind(pred,lwr,upr),
            col=c("red","green","green"),
            lty=2,
            lwd=c(3,2,2)))
于 2015-12-12T02:59:22.217 に答える