Is there a general method for simulating data from a formula or analysis available?

Question

De novo simulation of data from an experimental design data frame.
With a focus on R (though other language solution would be great).

In designing an experiment or a survey, simulating data and conducting an analysis on this simulated data can provide terrific insight into advantages and weaknesses of the design.

Such an approach can also be essential to the understanding and proper use of statistical tests.

However, this process tends to be somewhat tedious and many are led to skip past this important step in an experiment or survey.

Statistical models and test contain most of the information required to simulate the data (including an assumption or an explicit statement of the distribution).

Given an an analysis model (and its associated assumptions eg. normality and balance), the levels of a factor and a measure of significance (such as p-value), I would like to obtain simulated data (ideally with a generalized function akin to print(), predict(), simulate()).

Is such a generalized simulation framework possible?

If so, is such a framework currently available?

Example, I would like a function, such as:

 sim(aov(response~factor1+factor2*factor3),
          p.values=list(factor1=0.05,
                        factor2=0.05,
                        factor3=0.50,
                        factor2:factor3=0.05),
          levels=list(factor1=1:10,
                      factor2=c("A", "B", "C"),
                      factor3=c("A", "B", "C")))

ie, a generalized version of:

sim.lm<-function(){
library(DoE.base)
design<-fac.design(nlevels=c(10,3,3),
                   factor.names=c("factor1", "factor2", "factor3"),
                   replications=3,
                   randomize=F)

response<-with(design, as.numeric(factor1)+
                      as.numeric(factor2)+
                      as.numeric(factor3)+
                      as.numeric(factor2)*as.numeric(factor3)+
                      rnorm(length(factor1)))

simulation<-data.frame(design, response)}

OR

sim(glm(response~factor1+factor2*factor3, family=poisson),
         p.values=list(factor1=0.05,
                       factor2=0.05,
                       factor3=0.50,
                       factor2:factor3=0.05),
         levels=list(factor1=1:10,
                     factor2=c("A", "B", "C"),
                     factor3=c("A", "B", "C")))

OR

  library(lme4)
  sim(lmer(response~factor1+factor2 + (factor2|factor3)),
           F_value=list(factor1=50,
                        factor2=50),
           levels=list(factor1=1:10,
                       factor2=c("A", "B", "C"),
                       factor3=c("A", "B", "C")))

that would create a complete corresponding data.frame

potential examples of specific functions (please edit at will)
- arima.sim

function exist to create a data.frame of the factor levels, without the modelled response:
eg. conf.design
http://cran.r-project.org/web/views/ExperimentalDesign.html

Answer 1

There actually is an S3 generic simulate that even returns the data frame (or other list) you want. Type

?simulate  

It has methods for classes lm (works also for glm or for your aov example) and glm.nb (in MASS) already. You can now write S3 simulate methods for other classes of objects, e.g. for objects from lme4. You can check for which classes there are methods by typing

getAnywhere("simulate.class"), getAnywhere("simulate")  

or

getS3method("simulate","class"), methods(simulate) 

Answer 2

There is a new function called simfun in the TeachingDemos package for R (it is currently only in the development version on R-forge, it will be a while before it is on CRAN). It is intended to help with creating functions to do simulations.

One of the intended uses is that a teacher would use the simfun function to create a function and distribute it to students (possibly have a web interface at some point as an alternative). The students would then create a data frame of the factors representing an experimental design, pass this data frame to the created function and get returned the data frame with an additional column of the response simulated according to the parameters and error distribution set up by the teacher, the student can then analyze the data. This allows the teacher to set up a "True" relationship, but allow students to try many different experimental designs to explore ways to get at the "Truth" in much less time than doing an actual experiment and requiring less work for the teacher than to create or find a bunch of different sample datasets representing the different possible designs.

The simfun function is designed to be flexible so the teacher/creator can base the simulations on a fitted regression model, parameters supplied by the teacher/creator, or parameters suplied by the student/user.

The created function can also be easily used in simulations (with the replicate command) to explore power, sample size, effect size, etc. Though the resulting simulation may be slower than hand crafting the simulation process.

This looks like what you describe with the exception that it does not take p-values to create the data, but use of power. functions or pwr. from the pwr package could be incorporated to create simulations based on specifying power and alpha rather than means and differences.

Here is an example from the help page (there are several more examples) that assumes that you are measuring heights of subjects (male and female) that are nested in cities which are nested in states, there is a random effect for state with a SD of 1 and a random effect for city (within state) with a SD of 0.5 then the "error" SD is 3, females have a simulation mean of 64 inches and males have a mean of 69 inches (the error SD and means are realistic, the random effects are a bit contrived). The simfun function is used to create a new function called simheight, then a data frame is created with state ID's, city ID's and a column for sex of the subject (the experimental design or sampling design), this is passed to simheight and the result in a new data frame with the simulated heights (in addition to the other variables) that could then be analyzed using appropriate tools.

# simulate a nested mixed effects model
simheight <- simfun({
  n.city <- length(unique(city))
  n.state <- length(unique(state))
  n <- length(city)
  height <- h[sex] + rnorm(n.state,0,sig.state)[state] + 
    rnorm(n.city,0,sig.city)[city] + rnorm(n,0,sig.e)
}, sig.state=1, sig.city=0.5, sig.e=3, h=c(64,69),
  drop=c('sig.state','sig.city','sig.e','h','n.city','n.state','n'))

tmpdat <- data.frame(state=gl(5,20), city=gl(10,10), 
  sex=gl(2,5,length=100, labels=c('F','M')))
heightdat <- simheight(tmpdat)

Answer 3

Typing methods(simulate) or getAnywhere("simulate") should work. The former gives you a few methods, if package lme4 is loaded:

[1] simulate.lm* simulate.merMod* simulate.negbin* simulate.polr*

Lm-objects are used for both lm and glm models.