1

混合モデルのバイナリ ロジスティック回帰をブートストラップする必要があります。モデル自体は正常に動作します (専門家の友人によって承認および修正されています) が、ブートストラップ バージョンにはバグがあります。ブートストラップされたバージョンは、以前に別の専門家の友人によって承認されていました (CrossValidated でしたが、後で mod が CrossValidated に属していないと言って私の投稿を削除しました)。しかし、同じコードがたまたま単純な固定効果多重ロジスティック回帰で機能しました (ただし、その場合も、ここでの警告と同様の警告がたくさんありました [lmer() 関数に対するこの単一の警告を除いて: "mer_finalize( ans) : 偽収束 (8)")。

エラーの場所とデバッグ方法を教えてください。

どうもありがとう。

私のコードは次のとおりです(一時的にレプリケート数を低くしてコードをデバッグできませんでした):

library(boot)
library(lme4)

mixedGLM <- function(formula, data, indices) {
        d <- data[indices, ]
        (fit <- lmer(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                     + (1 | PatientID) + (0 + Trt | PatientID)
                     , family=binomial(logit), d))
        return(coef(fit))
      }

results <- boot(data=MixedModelData4 , statistic = mixedGLM, R= 2, formula= DV~Demo1 +Demo2 +Demo3 +Demo4 +Trt)

. . . 私のエラーは次のとおりです。

Error in t.star[r, ] <- res[[r]] : 
  incorrect number of subscripts on matrix
In addition: Warning messages:
1: In mer_finalize(ans) : false convergence (8)
2: glm.fit: algorithm did not converge 
3: glm.fit: fitted probabilities numerically 0 or 1 occurred 
4: glm.fit: fitted probabilities numerically 0 or 1 occurred 
5: In mer_finalize(ans) : false convergence (8) 

. . . また、boot() 関数に P 値も与える方法を教えてください??! ベータと SE、バイアスと CI が得られるだけですが、P 値も必要です。

どうもありがとう。

-------------------------------------------------- - 進行中の話 - - - - - - - - - - - - - - - - - - - - - - - - ------

わかりました Henrik の素敵なコードを喜んで実行しました。しかし、コードの実行はまだ完了していません。最初にこのエラーが発生しました:

Fitting 17 lmer() models:
[...
Error: pwrssUpdate did not converge in 30 iterations
In addition: Warning message:
In mixed(DV ~ (Demo1 + Demo2 + Demo3 + Demo4 + Trt)^2 + (1 | PatientID) +  :
  Due to missing values, reduced number of observations to 90
> (results2 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2
+ results3 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 

次に、最初の括弧ブロックを削除し、構文を次のように修正しました。

results3 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                 + (0 + Trt | PatientID),
                 family=binomial(logit), data = MixedModelData4,
                 method = "PB", args.test = list(nsim = 2))

今回は、テストは最初のステップ (モデルのフィッティング) に合格しましたが、P 値の取得に失敗し、同じエラーと警告が再び表示されました。

Fitting 17 lmer() models:
[.................]
Obtaining 16 p-values:
[....
Error: pwrssUpdate did not converge in 30 iterations
In addition: Warning messages:
1: In mixed(DV ~ (Demo1 + Demo2 + Demo3 + Demo4 + Trt)^2 + (0 + Trt |  :
  Due to missing values, reduced number of observations to 90
2: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
3: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
4: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
5: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
6: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations

デバッグ方法がわからない、または問題がデータセットにあるのか? 私のデータセットは完全に平均中心 (すべての変数) であることを付け加えておきます。DV は否定されるだけです (平均センタリングは R の動作を許可せず、否定はバイナリ結果に対して同じことを行うため)。

-------------------------------------------------- - - - - アップデート - - - - - - - - - - - - - - - - - - - - - --------------------

METHOD の PB 値を LRT に変更し (Henrik が推奨したように)、モデルを適合させるプロセスは終了しましたが、P 値を取得するプロセスは開始されませんでした。

> results4 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
+                   + (0 + Trt | PatientID),
+                   family=binomial(logit), data = MixedModelData4,
+                   method = "LRT", args.test = list(nsim = 2))
Fitting 17 lmer() models:
[.................]
Warning message:
In mixed(DV ~ (Demo1 + Demo2 + Demo3 + Demo4 + Trt)^2 + (0 + Trt |  :
  Due to missing values, reduced number of observations to 90

LRT を使用している場合、ブートストラップでは P 値が得られないことがわかりました。したがって、結果はすでに準備ができていました (ただし、ブートストラップされていません)。

4

1 に答える 1

4

パラメトリック ブートストラップを使用してから p 値が必要な場合は、次の方法で取得するパッケージGLMMの関数を使用できます。mixedafexpbkrtest::PBmodcomp

library(afex)
results <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                     + (1 | PatientID) + (0 + Trt | PatientID),
                     family=binomial(logit), data = d,
                     method = "PB", args.test = list(nsim = 1000))

最初にローカル クラスターを定義することもできます (つまり、複数のコアを使用します)。

cl <- makeCluster(rep("localhost", 4))
results <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2 
                     + (1 | PatientID) + (0 + Trt | PatientID),
                     family=binomial(logit), data = d,
                     method = "PB", args.test = list(nsim = 1000, cl = cl))

おそらく、3 つすべてのパッケージの開発バージョンをインストールするのが最善の方法です (現在のバージョンのpbkrtestは 1.0 用に設計されlme4ており、まだクランには含まれていません)。

于 2013-08-25T12:12:18.010 に答える