我需要引导我的混合模型二元逻辑回归。该模型本身运行良好（并且得到了专家朋友的批准和纠正），但自举版本有问题。引导版本之前已被另一位专家朋友批准（在 CrossValidated 中，但后来的 mods 删除了我的帖子，说它不属于 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 值吗？？！它只给出 beta 和 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 值不是通过引导获得的。因此，结果已经准备好（尽管不是自举的）。