我想在R中拟合随机截距互补对数回归，以检查未观察到的用户异质性。 我搜索了互联网和书籍，在 Stata 中只找到了一个解决方案，也许有人可以将其改编为 R。 在 Stata 中，有 2 个命令可用：
- 两级随机截距的
xtcloglog - 随机系数和更高级别模型的
gllamm
我的数据与人们的活动是否完成并受到阳光的影响有关 -completion是结果变量，sunshine和下面提到的其他变量将是解释变量;这是一个简化版本。
581755 obs. of 68 variables:
$ activity : int 37033 37033 37033 37033 37033 37033 37033 37033 37033 37033 ...
$ people : int 5272 5272 5272 5272 5272 5272 5272 5272 5272 5272 ...
$ completion: num 0 0 0 0 0 0 0 0 0 0 ...
$ active : int 0 2 2 2 2 2 2 2 2 2 ...
$ overdue : int 0 0 0 0 0 0 0 0 0 0 ...
$ wdsp : num 5.7 5.7 7.7 6.4 3.9 5.8 3.5 6.3 4.8 9.4 ...
$ rain : num 0 0 0 0 0 0 0 0 0 0 ...
$ UserCompletionRate: num [1:581755, 1] NaN -1.55 -1.55 -1.55 -1.55 ...
..- attr(*, "scaled:center")= num 0.462
..- attr(*, "scaled:scale")= num 0.298
$ DayofWeekSu : num 0 0 0 0 0 1 0 0 0 0 ...
$ DayofWeekMo : num 0 0 0 0 0 0 1 0 0 0 ...
$ DayofWeekTu : num 1 0 0 0 0 0 0 1 0 0 ...
$ DayofWeekWe : num 0 1 0 0 0 0 0 0 1 0 ...
$ DayofWeekTh : num 0 0 1 0 0 0 0 0 0 1 ...
$ DayofWeekFr : num 0 0 0 1 0 0 0 0 0 0 ...
$ MonthofYearJan : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearFeb : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearMar : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearApr : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearMay : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearJun : num 1 1 1 1 1 1 1 1 1 1 ...
$ MonthofYearJul : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearAug : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearSep : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearOct : num 0 0 0 0 0 0 0 0 0 0 ...
$ MonthofYearNov : num 0 0 0 0 0 0 0 0 0 0 ...
$ cold : num 0 0 0 0 0 0 0 0 0 0 ...
$ hot : num 0 0 0 0 0 0 0 0 0 0 ...
$ overduetask : num 0 0 0 0 0 0 0 0 0 0 ...
原始(简化(数据：
df <- people = c(1,1,1,2,2,3,3,4,4,5,5),
activity = c(1,1,1,2,2,3,4,5,5,6,6),
completion = c(0,0,1,0,1,1,1,0,1,0,1),
sunshine = c(1,2,3,4,5,4,6,2,4,8,4)
到目前为止，我已经将这段代码用于堵塞日志：
model <- as.formula("completion ~ sunshine")
clog_full = glm(model,data=df,family = binomial(link = cloglog))
summary(clog_full)
使用包 glmmML：
model_re <- as.formula("completion ~ sunshine")
clog_re = glmmML(model_re,cluster = people, data= df,
family = binomial(link = cloglog))
summary(clog_re)
使用软件包 lme4：
model_re1<- as.formula("completion ~ (1|people) + sunshine")
clog_re1 <- glmer(model_re1, data=df,
family = binomial(link = cloglog))
summary(clog_re1)
但是，R 不会从中获得任何结果，只是运行它们但永远不会得到结果。我是否必须将人员或活动用作群集？
如果有人对如何使用固定截距运行此模型也有想法，我很高兴知道。
这对我来说很好用(或多或少：请参阅下面的注释(
## added data.frame()
df <- data.frame(people = c(1,1,1,2,2,3,3,4,4,5,5),
activity = c(1,1,1,2,2,3,4,5,5,6,6),
completion = c(0,0,1,0,1,1,1,0,1,0,1),
sunshine = c(1,2,3,4,5,4,6,2,4,8,4)
)
model_re1 <- completion ~ (1|people) + sunshine
clog_re1 <- glmer(model_re1, data=df,
family = binomial(link = cloglog))
- 这很快完成(不到一秒钟(：也许你忘了关闭引号或括号或其他东西...... ？
- 但是，它确实会产生一条消息"边界(奇异(拟合：参见？isSingular"，这是因为您的数据集太小/太嘈杂，以至于人与人之间变异的最佳估计值为零(因为它不能是负数(。
更新：很抱歉地告诉您，混合模型 (GLMM( 的计算密集度明显高于标准 GLM：具有 68 个预测变量的 500K 观测值绝对是一个大问题，您应该预计拟合需要数小时。我有几点建议：
- 您绝对应该尝试数据的子集(包括观测值和预测变量(，以了解计算时间将如何缩放
- 据我所知，
glmmTMB包是 R 中解决此问题的最快选项(lme4会随着大量预测变量而严重扩展(，但可能会遇到内存限制。Mixed Models.jl Julia 包可能会更快。 - 您通常可以打开"详细"或"跟踪"选项，让您知道模型至少正在解决问题，而不是完全卡住(知道需要多长时间才能完成并不可行，但至少你知道正在发生一些事情......
- 如果 Stata 要快得多(我对此表示怀疑，但有可能(，您可以使用它
下面是一个包含 10,000 个观测值和一个预测变量的示例。
n <- 10000
set.seed(101)
dd <- data.frame(people=factor(rep(1:10000,each=3)),
sunshine=sample(1:8,size=n, replace=TRUE))
dd$completion <- simulate(~(1|people)+sunshine,
newdata=dd,
family=binomial(link="cloglog"),
newparams=list(beta=c(0,1),
theta=1))[[1]]
glmer运行 80 秒，然后失败：
system.time(update(clog_re1, data=dd, verbose=100))
另一方面，glmmTMB大约 20 秒内完成此问题(我的计算机上有 8 个内核，glmmTMB使用所有这些内核，因此分配给此作业的 CPU 分配高达 750%;如果内核较少，经过的计算时间将相应增加(。
library(glmmTMB)
system.time(glmmTMB(model_re1,data=dd,family = binomial(link = cloglog),
verbose=TRUE))
如果我再添加四个预测变量(总共 5 个(，计算时间最多为 46 秒(同样使用 8 个内核：所有内核的总计算时间为 320 秒(。 由于预测变量是原来的 13 倍，观测值是原来的 50 倍，您绝对应该期望这是一个具有挑战性的计算。
对异质性的一个非常粗略的评估是拟合齐次模型，并将残差偏差(或皮尔逊残差的平方和(与模型的残余自由度进行比较;如果前者要大得多，那就是某种形式的不拟合(异质性或其他东西(的证据。