我正在处理身体活动数据和后续疼痛数据。我有一个很大的数据集，但为了这个例子的摇动，我用我感兴趣的变量创建了一个小数据集。
由于我的身体活动数据本质上是组合的，因此我在使用这些变量作为混合效应模型中的预测变量之前使用组合数据分析。我的目标是使用 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)))))
任何帮助将不胜感激。
谢谢
将结果分配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]))
结果非常奇怪str(dataset_weird)，以至于研究 R 对象结构的常用方法失败了
$ comp_w : unclass(x)[i, , drop = drop] 中的错误：（下标）逻辑下标太长
如果我们运行，sapply(dataset_weird, class)我们会看到这些元素具有 class acomp。（他们似乎还有一个奇怪的print()方法：当我们的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在模型或预测中包含 Composition/irl 的原因是我无法理解这个语句在做什么：
comp_l = do.call("rbind", replicate(n=nrow(griddata), mean(acomp(dataset[,12])), simplify = FALSE))