R：自举混合模型二进制逻辑回归 [英] R: bootstrapped mixed model binary logistic regression
问题描述
非常感谢
我的代码是（我暂时保留复制号码太低，无法调试代码）：
库（启动）
库（lme4）
mixedGLM< - 函数（公式，数据，索引）{
d < - data [indices ，]
（fit< - lmer（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2
+（1 | PatientID）+（0 + Trt | PatientID）
，family = binomial（logit），d））
return（coef（fit））
}
结果< - boot（data = MixedModelData4，statistic = mixedGLM，R = 2，公式= DV〜Demo1 + Demo2 + Demo3 + Demo4 + Trt）
。
。
。
我的错误是：
t.star [r，]< - res [[r]]中的错误：
矩阵上的下标不正确数量
另外：警告消息：
1：在mer_finalize（ans）中：false convergence（8）
2：glm.fit：不收敛
3：glm.fit：拟合概率0或1发生
4：glm.fit：拟合概率0或1发生
5：在mer_finalize（ans）中：假收敛（8）
。
。
。
你还能告诉我如何使boot（）函数给P值呢？它只是给出beta和SE和偏差和CI，但是我也需要P值。
非常感谢。
------------------------------------------------ ---发展故事--------------------------------------------- --------
好的，我很高兴地运行了Henrik的漂亮代码。但代码没有完成运行。首先它给了这个错误：
配件17 lmer（）型号：
[...
错误：pwrssUpdate没有收敛30次迭代
另外：警告消息：
混合（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2 +（1 | PatientID）+：
由于缺少值，将观察数减少到90
>（results2< - mixed（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2
+ results3&混合（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2
然后我删除第一个括号阻止并修改了这个语法：
results3< - mixed（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2
+（0 + Trt | PatientID），
family =二项式（logit），data = MixedModelData4，
method =PB，args.test = list nsim = 2））
这次测试通过了第一步（装配模型），但失败了获得P值，再次给出t他同样的错误和警告：
配件17 lmer（）型号：
[........ .........]
获取16个p值：
[....
错误：pwrssUpdate在30次迭代中没有收敛
另外：警告消息：
1：混合（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2 +（0 + Trt | ：
由于缺少值，将观察值减少到90
2：In（function（fn，par，lower = rep.int（-Inf，n），upper = rep.int（Inf， ：
在10000次评估中无法收敛
3：In（function（fn，par，lower = rep.int（-Inf，n），upper = rep.int（Inf，：
未能收敛于10000个评估
4：In（function（fn，par，lower = rep.int（-Inf，n））upper = rep.int（Inf，：
无法收敛在10000评估
5：In（function（fn，par，lower = rep.int（-Inf，n），upper = rep.int（Inf，：
无法收敛10000个评估
6：In（function（fn，par，lower = rep.int（-Inf，n））upper = rep.int（Inf，：
在10000个评估中无法收敛
我不知道如何调试它，或者问题是我的数据集？我应该补充说，我的数据集是完全以平均为中心的变量），DV只被否定（因为平均对中不允许R工作，否定将对二进制结果做同样的） 。
---------------------------------- ------------------------更新------------------------- ------------------------------------
我将METHOD的PB值改为LRT（作为Henrik推荐），并且完成拟合模型的过程，但没有开始获取P值的过程：
>结果4< - 混合（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2
+ +（0 + Trt | PatientID），
+ family =二项式（logit），data = MixedModelData4 ，
+ method =LRT，args.test = list（nsim = 2））
装配17个lmer（）模型：
[........... ...]
警告信息：
混合（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2 +（0 + Trt |：
由于缺少值，减少的观察数量为90
原来，当LRT不通过自举获取P值因此，结果已经准备好了（虽然没有引导）。
如果你想要p值一个 GLMM 带有一个参数引导，你可以使用函数混合从包 afex 通过 pbkrtest :: PBmodcomp 获取它们：
库（afex）
结果< - 混合（DV〜（Demo1 + Dem o2 + Demo3 + Demo4 + Trt）^ 2
+（1 | PatientID）+（0 + Trt | PatientID），
family =二项式（logit），data = d，
method =PB，args.test = list（nsim = 1000））
您甚至可以先定义本地群集（即使用多个核心）：
cl< - makeCluster（rep（localhost，4））
结果< - 混合（DV〜（Demo1 + Demo2 + Demo3 + Demo4 + Trt）^ 2
+（1 | PatientID）+（0 + Trt | PatientID），
family =二项式（logit），data = d，
method =PB args.test = list（nsim = 1000，cl = cl））
这可能是最好的安装所有三个软件包的开发版本（如当前版本的 pbkrtest 是专为 lme4 1.0而设计的，还没有在起重机上）：
-
afex从r -forge：install.packages（afex，repos =http://R-Forge.R-project.org） -
pbkrtestfrom： http ：//people.math.aau.dk/~sorenh/software/pbkrtest/devel/ -
lme4from github： https://github.com/lme4/lme4
I need to bootstrap my mixed model binary logistic regression. The model itself works fine (and is approved and corrected by an expert friend), but the bootstrapped version is buggy. The bootstrapped version was previously approved by another expert friend (in CrossValidated but later mods removed my post saying it does not belong on CrossValidated). But the same code happened to work for a simple fixed-effects multiple logistic regression (although in that case too there were lots of warnings similar to the warnings here [except this single warning which is for the lmer() function: "In mer_finalize(ans) : false convergence (8)").
Could you please let me know where the error resides and how to debug it?
Many thanks.
My code is (I temporarily kept the replicate numbers too low to debug the code):
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)
. . . My errors are:
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)
. . . Also could you please tell me how to make the boot() function give P values too??! It just gives beta and SE and bias and CI, but I need the P values too.
Many thanks.
--------------------------------------------------- Developing Story -----------------------------------------------------
Ok I gladly ran the nice code of Henrik. But the code did not quite finish running. First it gave this error:
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
Then I removed the first parentheses block and revised the syntax to this one:
results3 <- mixed(DV ~ (Demo1 +Demo2+Demo3 +Demo4 +Trt)^2
+ (0 + Trt | PatientID),
family=binomial(logit), data = MixedModelData4,
method = "PB", args.test = list(nsim = 2))
This time the test passed the first step (fitting the models) but failed at obtaining P values, again giving the same errors and warnings:
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
I have no idea how to debug it, or if the problem is my dataset? I should add that my dataset is fully mean-centered (all variables). The DV is only negated (since mean centering disallowed R to work and negating would do the same for a binary outcome).
---------------------------------------------------------- Update -------------------------------------------------------------
I changed the PB value of METHOD to LRT (as Henrik recommended) and the process of fitting the models finished but the process of obtaining the P values didn't start:
> 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
It turned out the P values are not obtained by bootstrapping when LRT is being used. Therefore, the results were already ready (although non-bootstrapped).
If you want p-values from a GLMM with a parametric bootstrap you can use function mixed from package afex which obtains them via pbkrtest::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))
You could even first define a local cluster (i.e., use multiple cores):
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))
It is probably the best to install the development versions of all three packages (as the current version of pbkrtest is designed for lme4 1.0 which is not yet on cran):
afexfrom r-forge:install.packages("afex", repos="http://R-Forge.R-project.org")pbkrtestfrom: http://people.math.aau.dk/~sorenh/software/pbkrtest/devel/lme4from github: https://github.com/lme4/lme4
这篇关于R：自举混合模型二进制逻辑回归的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
