Applying boostrapping to MANCOVA in R (boot and MANCOVA functions)

Question

I have been attempting implement a MANCOVA within the boot function and keep running into the following error. Any help would be appreciated!

library(boot)
library (jmv)
DATA <- read.csv("DATA.csv", header = T)
STAT <- function(DATA) {
    mancova(data = DATA,
        deps = cbind('FAC1', 'FAC2', 'FAC3', 'FAC4', 'FAC5'),
        factors = 'subtype',
        covs = 'age',
        multivar = 'pillai',
        boxM = FALSE,
        shapiro = FALSE,
        qqPlot = FALSE)
}
BS.OUT = boot(DATA, statistic = STAT, R = 1000, sim = 'parametric')

Error in t.star[r, ] <- res[[r]] : number of items to replace is not a multiple of replacement length

Not quite sure what to do to fix this. Any advice would be much appreciated!

Thanks, A

Welcome to StackOverflow. See [how to make a great reproducible example](https://stackoverflow.com/a/5963610/2359523), namely by supplying example data. — Anonymous coward, Aug 27 '18 at 22:25

Anonymous coward · Answer 1 · 2018-08-29T16:22:51.037

The problem you are running into is because you're getting different shaped data. You can specify the shape of the data in your function. See if these links help at all.

I've supplied some random data as an example.

Var1 = runif(10, 3, 30)
Var2 = runif(10, 2, 5)
Y <- cbind(Var1, Var2)
Age <- runif(10, 20, 30)
Sex <- gl(2, 5, labels = c("male", "female"))
df <- as.data.frame(cbind(Y, Age, Sex))
df$Sex <- as.factor(df$Sex)

stat_test_func <- function(x, indices) {
  x$Var1 <- x$Var1[indices]
  model <- glm(Var1 + Var2 ~ Sex*Age, data = x)
  return(coefs = coef(model))
}

bs.out <- boot(df, statistic = stat_test_func, R = 10, sim = "parametric")

> bs.out

PARAMETRIC BOOTSTRAP

Call:
boot(data = df, statistic = stat_test_func, R = 10, sim = "parametric")

Bootstrap Statistics :
       original  bias    std. error
t1*  131.650088       0           0
t2* -132.423883       0           0
t3*   -4.659926       0           0
t4*    5.804188       0           0

Edit

Upon reading ?boot some more, it looks like for parametric bootstraps you need to supply a random generation function, ran.gen. But I'm having some trouble with the same error. It may be that each bootstrap in unequal, i.e. the factors aren't equal, but I'm not sure. Removing age didn't change it.

stat_test_func <- function(x) {
  mancova(x, deps = c("Var1", "Var2"), factors = "Sex", multivar = "pillai")  
    }

stat_test_rg <- function(x, mle) {
  out <- x
  out$Var1 <- rexp(nrow(out), 1/mle)
  out$Var2 <- rexp(nrow(out), 1/mle)
  out
  print(out)
}

bs.out <- boot(df, statistic = stat_test_func, R = 10, sim = "parametric", ran.gen = stat_test_rg, mle = mean(df$Var1))

Error in t.star[r, ] <- res[[r]] : 
  number of items to replace is not a multiple of replacement length

If you treat Age as a factor, I get a different error. So this answer may not be so useful anymore.

Var1 = runif(10, 3, 30)
Var2 = runif(10, 2, 5)
Age <- gl(2, 5, labels = c("young", "old"))
# Age <- runif(10, 20, 30)
Sex <- gl(2, 5, labels = c("male", "female"))
Y <- cbind(Var1, Var2)
df <- as.data.frame(cbind(Y, Age, Sex))
df$Sex <- as.factor(df$Sex)
df$Age <- as.factor(df$Age)

Error: Table$setRow(): value 'stat[pillai]' is not atomic

I don't understand the `indices` part of your code. Is that used in your example? — DanY, Aug 27 '18 at 23:05
@DanY That looks like it was supposed to be used in the `boot` function, but removing it doesn't change anything, so I'm not sure what the author intended it to do. The `statistic` argument requires `data` and a vector of indices. And actually, I just noticed that if `parametric` is used, it wants `ran.gen`, which generates a random subset of the original data. I'll edit to provide an example using `ran.gen`. — Anonymous coward, Aug 29 '18 at 15:05

Applying boostrapping to MANCOVA in R (boot and MANCOVA functions)

1 Answers1

Edit