使用lmer与predict函数时出现'invalid times argument'错误求助

处理体力活动与随访疼痛数据时，构建了包含核心变量的小型数据集。由于体力活动数据具有组成型特征，先进行组成数据分析，再将其作为预测变量构建混合效应模型，但使用predict()函数预测新创建的数据时，出现错误：Error in rep(0, nobs) : invalid 'times' argument，旧帖解决方案无效。以下是数据集及原始代码：

library("tidyverse")
library("compositions")
library("robCompositions")
library("lme4")

dataset <- structure(list(work = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 
                              3L, 3L, 4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"), 
           department = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 
                                    3L, 4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"), 
           worker = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
                                4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"), 
           age = c(45, 43, 65, 45, 76, 34, 65, 23, 23, 45, 32, 76), 
           sex = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
                             2L, 2L), .Label = c("1", "2"), class = "factor"), pain = c(4, 
                                                                                        5, 3, 2, 0, 7, 8, 10, 1, 4, 5, 4), lpa_w = c(45, 65, 43, 
                                                                                                                                     76, 98, 65, 34, 56, 2, 3, 12, 34), mvpa_w = c(12, 54, 76, 
                                                                                                                                                                                   87, 45, 23, 65, 23, 54, 76, 23, 54), lpa_l = c(54, 65, 34, 
                                                                                                                                                                                                                                  665, 76, 87, 12, 34, 54, 12, 45, 12), mvpa_l = c(12, 43, 
                                                                                                                                                                                                                                                                                   56, 87, 12, 54, 76, 87, 98, 34, 56, 23)), class = "data.frame", row.names = c(NA, 
                                                                                                                                                                                                                                                                                                                                                                 -12L))

# 创建体力活动的组成数据
dataset$comp_w <- acomp(cbind(lpa_w = dataset[,7], 
                          mvpa_w = dataset[,8]))

dataset$comp_l <- acomp(cbind(lpa_l = dataset[,9], 
                          mvpa_l = dataset[,10]))

# 构建用于预测的lpa_w和mvpa_w网格
mygrid <- expand.grid(lpa_w = seq(min(2), max(98),5),
                      mvpa_w = seq(min(12), max(87), 5))

griddata <- acomp(mygrid)

# 拟合混合效应模型
model <- lmer(pain ~ ilr(comp_w) + age + sex + ilr(comp_l) +
            (1 | work / department / worker),
          data = dataset)

# 尝试预测（报错代码）
(prediction = predict(model, newdata = list(comp_w = griddata,
                                        age = rep(mean(dataset$age, na.rm=TRUE),nrow(griddata)), 
                                        sex = rep("1", nrow(griddata)),
                                        comp_l = do.call("rbind", replicate(n=nrow(griddata), mean(acomp(dataset[,12])), simplify = FALSE)),
                                        work = rep(dataset$work, nrow(griddata)),
                                        department = rep(dataset$department, nrow(griddata)),
                                        worker = rep(dataset$worker, nrow(griddata))))

错误原因 #

• 分组变量长度不匹配：rep(dataset$work, nrow(griddata))生成的向量长度是12 * 320 = 3840，而其他变量（如age）的长度是320，导致模型计算观测数时出错，触发invalid 'times' argument。
• 组成数据构造错误：comp_l的构造方式有误，mean(acomp(dataset[,12]))返回单个组成数据对象，使用do.call("rbind", ...)无法正确生成与griddata行数匹配的组成数据列。
• newdata结构问题：使用list作为newdata时，lme4对组成数据的处理兼容性不如data.frame。

解决方案 #

方案1：仅预测固定效应（忽略随机效应） #

如果不需要预测随机效应，设置re.form=NA，无需提供分组变量：

# 计算comp_l的均值组成
comp_l_mean <- mean(dataset$comp_l)
# 构造匹配行数的comp_l数据
comp_l_pred <- acomp(matrix(rep(comp_l_mean@x, nrow(griddata)), ncol=2, byrow=TRUE))

# 构造newdata为data.frame
new_data <- data.frame(
  comp_w = griddata,
  age = rep(mean(dataset$age, na.rm=TRUE), nrow(griddata)),
  sex = factor(rep("1", nrow(griddata)), levels = levels(dataset$sex)),
  comp_l = comp_l_pred
)

# 预测固定效应
prediction <- predict(model, newdata = new_data, re.form = NA)
head(prediction)

方案2：包含随机效应预测 #

若要包含随机效应，需确保分组变量长度与griddata匹配，这里选择重复原数据集的分组组合：

# 计算comp_l的均值组成
comp_l_mean <- mean(dataset$comp_l)
comp_l_pred <- acomp(matrix(rep(comp_l_mean@x, nrow(griddata)), ncol=2, byrow=TRUE))

# 重复分组变量到匹配griddata行数
group_rep <- rep(1:4, each = nrow(griddata)/4) # 原数据集有4个work/worker组合
new_data <- data.frame(
  comp_w = griddata,
  age = rep(mean(dataset$age, na.rm=TRUE), nrow(griddata)),
  sex = factor(rep("1", nrow(griddata)), levels = levels(dataset$sex)),
  comp_l = comp_l_pred,
  work = factor(group_rep, levels = levels(dataset$work)),
  department = factor(group_rep, levels = levels(dataset$department)),
  worker = factor(group_rep, levels = levels(dataset$worker))
)

# 预测（包含随机效应）
prediction <- predict(model, newdata = new_data)
head(prediction)

内容的提问来源于stack exchange，提问作者user13069688