random model formula object

Question

I want to put formula in random model, but I think following error is due to wrong formula object (?), but could not fix it.

set.seed(1234)
mydata <- data.frame (A = rep(1:3, each = 20), B = rep(1:2, each = 30),
                      C = rnorm(60, 10, 5))
mydata$A <- as.factor(mydata$A)
mydata$B <- as.factor(mydata$B)
myfunction <- function (mydata, yvars, genovar, replication) { 
              require("lme4")
              formula = paste ("yvars" ~ 1|"genovar" + 1|"replication")
              model1 <- lmer(formula, data = dataframe, REML = TRUE)
              return(ranef(model2))
}

myfunction(mydata=dataf, yvars = "C", genovar = "A", replication = "B")

Error: length(formula <- as.formula(formula)) == 3 is not TRUE

Answer 1

There were several wonky things in here, but this is I think close to what you want.

set.seed(1234)
mydata <- data.frame (A = factor(rep(1:3, each = 20)),
                      B = factor(rep(1:2, each = 30)),
                      C = rnorm(60, 10, 5))
require("lme4")
myfunction <- function (mydata, yvars, genovar, replication) { 
  formula <- paste (yvars,"~ (1|",genovar,") + (1|",replication,")")
  model1 <- lmer(as.formula(formula), data = mydata, REML = TRUE)
  return(ranef(model1))
}
myfunction(mydata=mydata, yvars = "C", genovar = "A", replication = "B")

Beware, however, that lmer doesn't work the way that classical random-effects ANOVA does -- it may perform very badly with such small numbers of replicates. (In the example I tried it set the variance of A to zero, which is at least not unreasonable.) The GLMM FAQ has some discussion of this issue. (Random-effects ANOVA would have exceedingly low power in that case but might not be quite as bad.) If you really want to do random-effects models on such small samples you might want to consider reconstructing the classical method-of-moments approach (as I recall there is/was a raov formula in S-PLUS that did random-effects ANOVA, but I don't know if it was ever implemented in R).

Finally, for future questions along these lines you may do better on the r-sig-mixed-models@r-project.org mailing list -- Stack Overflow is nice but there is more R/mixed-model expertise over there.