在 confint 选项中使用配置文件和引导方法，使用 glmer 模型
using profile and boot method within confint option, with glmer model
我将 glmer 与 logit link 一起用于高斯误差模型。
当我尝试使用配置文件或带有 confint 选项的引导方法获取置信区间时，我收到了使用配置文件可能性和引导程序的错误：
> Profile: Computing profile confidence intervals ... Error in
> names(opt) <- profnames(fm, signames) : 'names' attribute [2] must
> be the same length as the vector [1]
>
> Boot: Error in if (const(t, min(1e-08, mean(t, na.rm = TRUE)/1e+06)))
> { : missing value where TRUE/FALSE needed In addition: Warning
> message: In bootMer(object, FUN = FUN, nsim = nsim, ...) : some
> bootstrap runs failed (9999/10000)
我查看了关于如何通过安装 lme4 开发工具来克服配置文件问题的在线建议，并且我还从数据集中消除了所有 NA。但是，在这两种情况下，我仍然收到相同的两条错误消息。
是否有人能够就这是 lme4 问题还是更基本的问题提供任何帮助。
这是生成前 20 个观测值和模型格式的代码：
df2 <- data.frame(
prop1 = c(0.46, 0.471, 0.458, 0.764, 0.742, 0.746, 0.569, 0.45, 0.491, 0.467, 0.464,
0.556, 0.584, 0.478, 0.456, 0.46, 0.493, 0.704, 0.693, 0.651),
prop2 = c(0.458, 0.438, 0.453, 0.731, 0.738, 0.722, 0.613, 0.498, 0.452, 0.451, 0.447,
0.48, 0.576, 0.484, 0.473, 0.467, 0.467, 0.722, 0.707, 0.709),
site = c(1,1,2,3,3,3,4,4,4,4,4,5,5,5,6,6,7,8,8,8)
)
df2$site <- factor(df2$site)
model <- glmer(prop2 ~ prop1 + (1|site),
data=df2, family=gaussian(link="logit"))
summary(model)
响应是比例 (0,1)，协变量也是。我使用 logit link 将响应的预期值保持在 (0,1) 之间，而不是不受正常值的约束，尽管该数据相当适合 LMM。
我还将分析两个比例之间的差异，为此我预计一些差异（和拟合值）会违反 (0,1) 边界 - 因此我将使用身份 link 具有高斯误差，或 logit link 在比例差异上（将其强制为 (0,1)）。
另外，鉴于每个站点可能只有1-5条记录，我自然会考虑运行线性回归或GLM（比例），并在建模过程中将站点视为固定效应，如果随机效应方差的估计值为零（或非常小）。
对于您的示例，它适用于 bootstrap:
confint(model, method = "boot")
# 2.5 % 97.5 %
# .sig01 12.02914066 44.71708844
# .sigma 0.03356588 0.07344978
# (Intercept) -5.26207985 1.28669024
# prop1 1.01574201 6.99804555
考虑到在您提出的模型下，虽然您的估计值始终在 0 和 1 之间，但预计会观察到小于 0 和大于 1 的值。
您已确定当前版本 lme4 的错误，现已在 Github 上部分修复（它适用于 method="boot"）。 （devtools::install_github("lme4/lme4") 如果安装了开发工具，应该可以安装它。）
library(lme4)
fit_glmer <- glmer(prop2 ~ prop1 + (1|site),
data=df2, family=gaussian(link="logit"))
Profiling/profile 置信区间仍然不起作用，但有一个更有意义的错误：
try(profile(fit_glmer))
## Error in profile.merMod(fit_glmer) :
## can't (yet) profile GLMMs with non-fixed scale parameters
引导确实有效。有很多警告，很多改装尝试都失败了，但我希望那是因为提供的数据集很小。
bci <- suppressWarnings(confint(fit_glmer,method="boot",nsim=200))
我想推荐几个其他选项。您可以使用 glmmADMB or glmmTMB，并且这些平台还允许您使用 Beta 分布对比例进行建模。我会考虑通过 "melting" 数据对比例类型之间的差异进行建模（以便有一个 prop 列和一个 prop_type 列）并包括 prop_type 作为预测变量（可能具有个体水平的随机效应，确定哪些比例是在同一个人身上测量的）...
library(glmmADMB)
fit_ADMB <- glmmadmb(prop2 ~ prop1 + (1|site),
data=df2, family="gaussian",link="logit")
## warning message about 'sd.est' not defined -- only a problem
## for computing standard errors of predictions, I think?
library(glmmTMB)
fit_TMB <- glmmTMB(prop2 ~ prop1 + (1|site),
data=df2, family=list(family="gaussian",link="logit"))
听起来您的数据可能更适合 Beta 模型？
fit_ADMB_beta <- glmmadmb(prop2 ~ prop1 + (1|site),
data=df2, family="beta",link="logit")
fit_TMB_beta <- glmmTMB(prop2 ~ prop1 + (1|site),
data=df2,
family=list(family="beta",link="logit"))
比较结果（比需要的更漂亮）
library(broom)
library(plyr)
library(dplyr)
library(dwplot)
tab <- ldply(lme4:::namedList(fit_TMB_beta,
fit_ADMB_beta,
fit_ADMB,
fit_TMB,
fit_glmer),
tidy,.id="model") %>%
filter(term != "(Intercept)")
dwplot(tab) +
facet_wrap(~term,scale="free",
ncol=1)+theme_set(theme_bw())
我将 glmer 与 logit link 一起用于高斯误差模型。
当我尝试使用配置文件或带有 confint 选项的引导方法获取置信区间时，我收到了使用配置文件可能性和引导程序的错误：
> Profile: Computing profile confidence intervals ... Error in
> names(opt) <- profnames(fm, signames) : 'names' attribute [2] must
> be the same length as the vector [1]
>
> Boot: Error in if (const(t, min(1e-08, mean(t, na.rm = TRUE)/1e+06)))
> { : missing value where TRUE/FALSE needed In addition: Warning
> message: In bootMer(object, FUN = FUN, nsim = nsim, ...) : some
> bootstrap runs failed (9999/10000)
我查看了关于如何通过安装 lme4 开发工具来克服配置文件问题的在线建议，并且我还从数据集中消除了所有 NA。但是，在这两种情况下，我仍然收到相同的两条错误消息。
是否有人能够就这是 lme4 问题还是更基本的问题提供任何帮助。
这是生成前 20 个观测值和模型格式的代码：
df2 <- data.frame(
prop1 = c(0.46, 0.471, 0.458, 0.764, 0.742, 0.746, 0.569, 0.45, 0.491, 0.467, 0.464,
0.556, 0.584, 0.478, 0.456, 0.46, 0.493, 0.704, 0.693, 0.651),
prop2 = c(0.458, 0.438, 0.453, 0.731, 0.738, 0.722, 0.613, 0.498, 0.452, 0.451, 0.447,
0.48, 0.576, 0.484, 0.473, 0.467, 0.467, 0.722, 0.707, 0.709),
site = c(1,1,2,3,3,3,4,4,4,4,4,5,5,5,6,6,7,8,8,8)
)
df2$site <- factor(df2$site)
model <- glmer(prop2 ~ prop1 + (1|site),
data=df2, family=gaussian(link="logit"))
summary(model)
响应是比例 (0,1)，协变量也是。我使用 logit link 将响应的预期值保持在 (0,1) 之间，而不是不受正常值的约束，尽管该数据相当适合 LMM。
我还将分析两个比例之间的差异，为此我预计一些差异（和拟合值）会违反 (0,1) 边界 - 因此我将使用身份 link 具有高斯误差，或 logit link 在比例差异上（将其强制为 (0,1)）。
另外，鉴于每个站点可能只有1-5条记录，我自然会考虑运行线性回归或GLM（比例），并在建模过程中将站点视为固定效应，如果随机效应方差的估计值为零（或非常小）。
对于您的示例，它适用于 bootstrap:
confint(model, method = "boot")
# 2.5 % 97.5 %
# .sig01 12.02914066 44.71708844
# .sigma 0.03356588 0.07344978
# (Intercept) -5.26207985 1.28669024
# prop1 1.01574201 6.99804555
考虑到在您提出的模型下，虽然您的估计值始终在 0 和 1 之间，但预计会观察到小于 0 和大于 1 的值。
您已确定当前版本 lme4 的错误，现已在 Github 上部分修复（它适用于 method="boot"）。 （devtools::install_github("lme4/lme4") 如果安装了开发工具，应该可以安装它。）
library(lme4)
fit_glmer <- glmer(prop2 ~ prop1 + (1|site),
data=df2, family=gaussian(link="logit"))
Profiling/profile 置信区间仍然不起作用，但有一个更有意义的错误：
try(profile(fit_glmer))
## Error in profile.merMod(fit_glmer) :
## can't (yet) profile GLMMs with non-fixed scale parameters
引导确实有效。有很多警告，很多改装尝试都失败了，但我希望那是因为提供的数据集很小。
bci <- suppressWarnings(confint(fit_glmer,method="boot",nsim=200))
我想推荐几个其他选项。您可以使用 glmmADMB or glmmTMB，并且这些平台还允许您使用 Beta 分布对比例进行建模。我会考虑通过 "melting" 数据对比例类型之间的差异进行建模（以便有一个 prop 列和一个 prop_type 列）并包括 prop_type 作为预测变量（可能具有个体水平的随机效应，确定哪些比例是在同一个人身上测量的）...
library(glmmADMB)
fit_ADMB <- glmmadmb(prop2 ~ prop1 + (1|site),
data=df2, family="gaussian",link="logit")
## warning message about 'sd.est' not defined -- only a problem
## for computing standard errors of predictions, I think?
library(glmmTMB)
fit_TMB <- glmmTMB(prop2 ~ prop1 + (1|site),
data=df2, family=list(family="gaussian",link="logit"))
听起来您的数据可能更适合 Beta 模型？
fit_ADMB_beta <- glmmadmb(prop2 ~ prop1 + (1|site),
data=df2, family="beta",link="logit")
fit_TMB_beta <- glmmTMB(prop2 ~ prop1 + (1|site),
data=df2,
family=list(family="beta",link="logit"))
比较结果（比需要的更漂亮）
library(broom)
library(plyr)
library(dplyr)
library(dwplot)
tab <- ldply(lme4:::namedList(fit_TMB_beta,
fit_ADMB_beta,
fit_ADMB,
fit_TMB,
fit_glmer),
tidy,.id="model") %>%
filter(term != "(Intercept)")
dwplot(tab) +
facet_wrap(~term,scale="free",
ncol=1)+theme_set(theme_bw())