混合效应逻辑回归中每个集群/组的可变样本量答案

【问题标题】：Variable sample size per cluster/group in mixed effects logistic regression混合效应逻辑回归中每个集群/组的可变样本量
【发布时间】：2013-11-08 00:00:23
【问题描述】：

我正在尝试运行混合效应逻辑回归模型，但我担心每个集群/组中的可变样本大小，以及某些模型中的“成功”数量非常少。

我有大约 700 棵树，分布在 163 个田地（即集群/组）中，从 2004 年到 11 年每年访问一次。我正在为研究的每一年拟合单独的混合效应逻辑回归模型（以下称为 GLMM），以将此输出与来自共享脆弱模型的推论（即随机效应的生存分析）进行比较。

每个地块的树木数量从 1 到 22 不等。此外，有些年份的“成功”数量非常少（即病树）。例如，在 2011 年，694 次“失败”（即健康的树木）中只有 4 次成功。

我的问题是：（1）当推理的重点仅在于估计 GLMM 中的固定效应时，是否存在理想样本数的一般规则，以及（2）当存在这样的成功与失败的比例差异很大。

感谢您对来源的任何意见或建议。

-莎拉

【问题讨论】：

因为这个（好）问题的编程内容很少，最好在stats.stackexchange.com 上提问。请务必在移动之前将其删除，以免最终被交叉发布。
非常同意乔希的观点。我确实认为在您重新发布时提出一些暂定的 R 代码仍然会使这更加“具体”。 (4/ 694) 对于泊松模型来说，事件发生率并不是特别低。

标签： r logistic-regression

【解决方案1】：

（嗨，莎拉，抱歉我之前没有通过电子邮件回复...）

通常很难回答这些问题——你被困住了用你的数据，对吧？所以这不是功率分析的问题。如果你想确保你的结果是合理的可靠，可能最好的办法是运行一些模拟。我将展示lme4 的一个相当新的功能（在开发版本1.1-1，在Github上），这是为了模拟给定公式和一组参数的 GLMM 数据。

首先我必须模拟预测变量（你不会必须这样做，因为你已经有了数据——虽然您可能想尝试改变地块数量的范围，每块地的树等）。

set.seed(101)
## simulate number of trees per plot
## want mean of 700/163=4.3 trees, range=1-22
## by trial and error this is about right
r1 <- rnbinom(163,mu=3.3,size=2)+1
## generate plots and trees within plots
d <- data.frame(plot=factor(rep(1:163,r1)),
            tree=factor(unlist(lapply(r1,seq))))
## expand by year
library(plyr)
d2 <- ddply(d,c("plot","tree"),
        transform,year=factor(2004:2011))

现在设置参数：我将假设年份是固定的效果和总体疾病发病率是plogis(-2)=0.12，除了 2011 年的时候是plogis(-2-3)=0.0067。区间标准差为 1（在 logit 尺度上），与图中树中的标准一样偏差：

beta <- c(-2,0,0,0,0,0,0,-3)
theta <- c(1,1)  ## sd by plot and plot:tree

现在模拟：年份为固定效果，情节和树中情节为随机效应

library(lme4)
s1 <- simulate(~year+(1|plot/tree),family=binomial,
     newdata=d2,newparams=list(beta=beta,theta=theta))
d2$diseased <- s1[[1]]

总结/检查：

d2sum <- ddply(d2,c("year","plot"),
           summarise,
           n=length(tree),
           nDis=sum(diseased),
           propDis=nDis/n)
library(ggplot2)
library(Hmisc)  ## for mean_cl_boot
theme_set(theme_bw())
ggplot(d2sum,aes(x=year,y=propDis))+geom_point(aes(size=n),alpha=0.3)+
    stat_summary(fun.data=mean_cl_boot,colour="red")

现在拟合模型：

g1 <- glmer(diseased~year+(1|plot/tree),family=binomial,
        data=d2)
fixef(g1)

你可以多试几次，看看结果多可靠...

【讨论】：

【解决方案2】：

正如 Josh 所说，这是 CrossValidated 的一个更好的问题。

逻辑回归没有硬性规定，但一条经验法则是设计中每个单元（在本例中为集群）需要 10 次成功和 10 次失败乘以模型中连续变量的数量。

在您的情况下，我认为该模型如果收敛，将是不稳定的。您可以通过引导固定效应估计的误差来检查这一点。

【讨论】：