predictInterval fails when specifying 'newdata'

[aaraujo-lih]

I am getting the following error from predictInterval:

Error in [.data.frame(tmp, alllvl) : undefined columns selected

This is the code:

# Load lme4 package
library(lme4);

# Load merTools package
library(merTools);

# Hardcode data
ndata <- data.frame(
  "Patient"=as.factor(
    c(
      1, 1, 1, 1, 1, 1,
      2, 2, 2, 2, 2, 2,
      3, 3, 3, 3, 3, 3,
      4, 4, 4, 4, 4, 4,
      5, 5, 5, 5, 5, 5,
      6, 6, 6, 6, 6, 6,
      7, 7, 7, 7, 7, 7,
      8, 8, 8, 8, 8, 8,
      9, 9, 9, 9, 9, 9,
      10, 10, 10, 10, 10, 10,
      11, 11, 11, 11, 11, 11,
      12, 12, 12, 12, 12, 12
    )
  ),
  "Treatment"=as.factor(
    c(
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B",
      "A", "B", "A", "B", "A", "B"
    )
  ),
  "Cycle"=as.factor(
    c(
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3,
      1, 1, 2, 2, 3, 3
    )
  ),
  "Y"=c(
    2394, 2686, 2515, 2675, 2583, 2802,
    2746, 2726, 2592, 2867, 2743, 2742,
    2668, 2560, 2542, 2584, 2491, 2737,
    2397, 2696, 2411, 2895, 2499, 2760,
    3179, 3221, 2952, 3096, 2600, 3192,
    2643, 2496, 2759, 2847, 2651, 2860,
    2678, 2843, 2492, 2763, 2801, 2890,
    2887, 2862, 2875, 3083, 2689, 2967,
    2490, 2841, 2648, 3044, 2688, 2914,
    2268, 2576, 2413, 2493, 2344, 2699,
    2617, 2923, 2629, 2832, 2732, 2866,
    2627, 2759, 2712, 2698, 2572, 2826
  )
);

# Mixed model using cycle and treatment
# within patients as block structure
full0.lmm <- lmer(
  formula=Y~
    Treatment+(1|Patient)+
    (1|Patient:Cycle)+
    (1|Patient:Treatment),
  data=ndata,
  REML=TRUE
);

# Individual treatment effects
pi.full0.lmm <- predictInterval(
  merMod=full0.lmm,
  newdata=expand.grid(
    "Patient"=levels(ndata$Patient),
    "Treatment"=levels(ndata$Treatment),
    "Cycle"=levels(ndata$Cycle)[1] # constant Cycle
  ),
  level=0.95,
  n.sims=10000,
  stat="mean",
  type="linear.prediction",
  include.resid.var=TRUE,
  returnSims=TRUE
);

Either I am using 'predictInterval' the wrong way or this is a bug. This works with 'predict' from 'lme4' package. If I input the full data.frame to the 'newdata' argument it works. But that is not what I intend.

[carlbfrederick]

I located the offending piece of code in predictInterval and my guess is that it results from the combination of using the colon syntax for interactions (1|Patient:Cycle) and not including all of the factor levels in one of the interaction variables.

Should be fixable, just need to fix the way we have been assembling terms. I will let you know what I find.

[aaraujo-lih]

Dear carlbfrederick,

Many thanks.

That makes sense. If I supply the full dataset to the 'newdata' argument it works. However in my case I am interested in the individual difference between treatments while keeping all the other variables constant. This is the same regardless of the Cycle I choose. Using all the factor levels will be a waste of computational time and memory resources.

Wish you good luck on finding a quick solution.

[carlbfrederick]

I have patched this issue on a separate branch for now. The reproducible example you supplied works as intended. You can install it using

devtools::install_github("jknowles/merTools", ref="issue53")

Let us know if you run into any issues due to this patch. I will leave this issue open until it has been throughly tested and is merged back onto the main branch.

Cheers!

[aaraujo-lih]

Dear carlbfrederick,

Thanks a lot.

It seems to be working now. The full code I am using is attached.

merTools_example.zip

It is also working when I use all the levels of the factors to compute the predictions.

[carlbfrederick]

closed by branch Issue53