I am working with physical activity data and follow-up pain data. I have a large dataset but for the shake of this example I have created a small one with the variables of my interest.
As my physical activity data are compositional in nature I am using compositional data analysis before using those variables as predictors in my mixed-effects model. My goal is to use the predict() function to predict some new data that I have created but I am receiving the folowing:
Error in rep(0, nobs) : invalid 'times' argument
I have googled it and I saw a post that was posted a few year ago but the answer did not work for mine.
Below is the dataset and my code:
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)))))
Any help would be greatly appreciated.
Thanks
lpa_w
,mvpa_w
,lpa_l
,mvpa_l
come from/how are they defined?comp_w
orcomp_l
. I get two columns returned byacomp()
, which can't get thrown into a single column, so I ended up with an empty column namedcomp_w
and nocomp_l
. Are you getting something different? Did you inspect what you expected?