5

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

5
  • Could you please share all your data? Now your code is not reproducible.
    – Quinten
    Commented Sep 26, 2022 at 15:53
  • Specifically: is the structure you present supposed to be assigned to a variable? WHere do lpa_w, mvpa_w, lpa_l, mvpa_l come from/how are they defined?
    – Ben Bolker
    Commented Sep 26, 2022 at 20:45
  • 1
    Hey, apologies. Now the code is corrected. lpa_w, mvpa_w, lpa_l and mvpa_l become two compositions for leisure and work hence w = work and l = leisure. lpa = light physical activity and mvpa = moderate-to-vigorous physical activity. Commented Sep 27, 2022 at 6:23
  • It doesn't look like you've actually created comp_w or comp_l. I get two columns returned by acomp(), which can't get thrown into a single column, so I ended up with an empty column named comp_w and no comp_l. Are you getting something different? Did you inspect what you expected?
    – Kat
    Commented Sep 28, 2022 at 4:24
  • Hi Kat, yes I get 4 extra columns in my dataset which are called comp_w[,1], comp_w[,2], comp_l[,1] and comp_l[,2]. If you type colnames(dataset) you will only get comp_w and comp_l. Commented Sep 28, 2022 at 6:23

1 Answer 1

6
+25

Assigning the results of acomp to an element of a data frame gives a weird data structure that messes things up downstream.

Constructing this data set (without messing up the original 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]))

The result is so weird that str(dataset_weird), the usual way of investigating the structure of an R object, fails with

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

If we run sapply(dataset_weird, class) we see that these elements have class acomp. (They also appear to have an odd print() method: when we print(dataset_weird$comp_w) the results are a matrix of strings, but if we unclass(dataset_weird$comp_w) we can see that the underlying object is numeric [!])

This whole problem is kind of tricky since you're working with n-column matrices that are getting converted to special acomp() objects that are then getting converted to (n-1)-dimensional matrices (isometric log-ratio-transformed compositional data), the columns of which are then getting used as predictors. The basic point is that lme4's machinery will get confused if you have elements in your data frame that are not simple one-dimensional vectors. So you have to do the work of creating data frame columns yourself.

Here's what I came up with, with one missing piece (described below):

## 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")))

This seems to work.

The reason I did not include the _l composition/irl in the model or the predictions is that I couldn't understand what this statement was doing:

comp_l = do.call("rbind", replicate(n=nrow(griddata), mean(acomp(dataset[,12])), simplify = FALSE))
2
  • Thanks Ben, that's exactly what I did as well and it worked. Thank you very much for all the help and the explanation. Much appreciated! That statement was calculating the mean composition of all 5 activities during leisure time. Commented Oct 3, 2022 at 13:35
  • if the answer solved your problem, you are encouraged to click the check-mark to accept it ...
    – Ben Bolker
    Commented Oct 3, 2022 at 14:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.