logo标签列表

在使用lmer和predict时出现无效的“times”参数

rlme4predict
5

我正在处理身体活动数据和随访疼痛数据。我有一个大的数据集，但为了举例，我创建了一个只包含我感兴趣的变量的小型数据集。

由于我的身体活动数据具有组成性质，因此在将这些变量用作混合效应模型的预测变量之前，我使用组成数据分析。我的目标是使用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))

#create compositions of physical activity

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]))

#Make a grid to use for predictions for composition of lpa_w and mvpa_w
mygrid=rbind(
expand.grid(lpa_w = seq(min(2), max(98),5),
          mvpa_w = seq(min(12), max(87), 5)))

griddata <- acomp(mygrid)

#run the model
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)))))

任何帮助都将不胜感激。
谢谢。
请问您能分享所有数据吗？现在您的代码无法再现。 - Quinten
具体来说，您所呈现的结构是否应该分配给变量？lpa_w、mvpa_w、lpa_l、mvpa_l是从哪里来的/如何定义的？ - Ben Bolker
1
抱歉，现在代码已经被更正了。lpa_w、mvpa_w、lpa_l和mvpa_l成为了休闲和工作的两个组合，因此w代表工作，l代表休闲。lpa代表轻度体育活动，mvpa代表中至高强度体育活动。 - user13069688
看起来你实际上没有创建 comp_w 或 comp_l。我通过 acomp() 得到了两列返回值，这些值无法放入单个列中，所以我最终得到了一个名为 comp_w 的空列和没有 comp_l。你得到了不同的结果吗？你检查了你期望得到的东西吗？ - Kat
嗨，Kat，是的，我的数据集中有4列额外的列，它们分别被称为comp_w [，1]，comp_w [，2]，comp_l [，1]和comp_l [，2]。如果您键入colnames（dataset），您只会得到comp_w和comp_l。 - user13069688
1个回答
6

将acomp的结果分配给数据框的一个元素会产生奇怪的数据结构，会在下游造成混乱。

构建这个数据集（不会弄乱原始的dataset）：

dataset_weird <- dataset
dataset_weird$comp_w <- acomp(cbind(lpa_w = dataset[,7], 
                          mvpa_w = dataset[,8]))
dataset_weird$comp_l <- acomp(cbind(lpa_l = dataset[,9], 
                                    mvpa_l = dataset[,10]))

结果非常奇怪，以往查看R语言对象的结构的方式str(dataset_weird)会失败并报错如下：

$ comp_w :Error in unclass(x)[i, , drop = drop] : (subscript) logical subscript too long

如果我们运行sapply(dataset_weird, class)，可以看到这些元素的类是acomp（它们还似乎有一种奇怪的print()方法：当我们打印dataset_weird$comp_w时，结果是字符串矩阵，但如果我们运行unclass(dataset_weird$comp_w)，就可以看到底层对象是数值型[!]）。
整个问题非常棘手，因为你正在处理n列矩阵，这些矩阵正在被转换为特殊的acomp()对象，然后又被转换为(n-1)维矩阵（等比对数比例变换组成数据），其中的列被用作预测变量。基本上，lme4的机制会因为数据框中包含非简单一维向量的元素而变得混乱。因此，您必须自己创建数据框列。
下面是我的解决方案，只有一个缺失的部分（描述如下）：
## utility function: *either* uses a matrix argument (`comp_data`)
## *or* extracts relevant columns from a data frame (`data`):
## returns ilr-transformed values as a matrix, with appropriate column names
ilr_dat <- function(data, suffix = NULL, comp_data = NULL) {
    if (!is.null(suffix) && is.null(comp_data)) {
        comp_data <- as.matrix(data[grep(paste0(suffix,"$"), names(data))])
    }
    ilrmat <- ilr(acomp(comp_data))
    colnames(ilrmat) <- paste0("ilr", suffix, ".", 1:ncol(ilrmat))
    return(ilrmat)
}

## augment original data set (without weird compositional elements)
## using data.frame() rather than $<- or rbind() collapses matrix arguments
## to data frame rows in a way that R expects
dataset2 <- data.frame(dataset, ilr_dat(dataset, "_l"))
dataset2 <- data.frame(dataset2, ilr_dat(dataset, "_w"))

mygrid <- rbind(
    expand.grid(lpa_w = seq(min(2), max(98),5),
                mvpa_w = seq(min(12), max(87), 5)))

## generate ilr data for prediction
griddata <- as.data.frame(ilr_dat(comp_data=mygrid, suffix="_w"))

#run the model: ilr(comp_l) **not** included, see below
model <- lmer(pain ~ ilr_w.1 + age + sex  + ## ilr(comp_l) +
                  (1 | work / department / worker),
          data = dataset2)

## utility function for replication
xfun <- function(s) rep(dataset[[s]], nrow(griddata))
predict(model, newdata = data.frame(griddata,
                                    age = mean(dataset$age, na.rm=TRUE),
                                    sex = "1",
                                    work = xfun("work"),
                                    department = xfun("department"),
                                    worker = xfun("worker")))

这似乎有效。

我之所以没有在模型或预测中包含_l组合/irl，是因为我无法理解此语句的作用：

comp_l = do.call("rbind", replicate(n=nrow(griddata), mean(acomp(dataset[,12])), simplify = FALSE))
谢谢Ben，我也是这样做的，而且它起作用了。非常感谢您提供的所有帮助和解释。非常感激！该语句计算了休闲时间内所有5项活动的平均组成。 - user13069688
如果这个答案解决了你的问题，我们鼓励你点击勾选标记来接受它。 - Ben Bolker
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接