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

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

• One thing to bear in mind is that you need to specify an error distribution. In addition, I don't completely follow your "p.values" or "F_value" lists; is the idea that the sampling distributions of F & p be centered on the listed numbers? – gung - Reinstate Monica May 16 '12 at 0:56
• @gung, the distribution is assumed in the test (both here assume normality or I would have had to use something like glm). – Etienne Low-Décarie May 16 '12 at 4:17
• The test certainly does assume a specific distribution, but that doesn't mean that the data were actually generated in that way. You could generate data from any number of distributions for various purposes (e.g., simulation studies, robustness, etc.). – gung - Reinstate Monica May 16 '12 at 4:26
• @gung, I have tried to edit my question to make clear that I want to create a data set that corresponds to the expectations of the test/model. I am not looking for something for a simulation study or for a robustness study, but to provide a dry-run of an experimental design. – Etienne Low-Décarie May 16 '12 at 4:29
• stats.stackexchange.com/questions/11233/… – Etienne Low-Décarie Jun 18 '12 at 16:05

## 3 Answers

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)

• simulate would be great if it did not require "an object representing a fitted model" which itself already requires data or simulated data. – Etienne Low-Décarie May 16 '12 at 4:22
• >getAnywhere("simulate.class") no object named ‘simulate.class’ was found >getS3method("simulate","class") Error in getS3method("simulate", "class") : S3 method 'simulate.class' not found – Etienne Low-Décarie May 16 '12 at 4:31
• One could always construct an initial fitted object with all response values being 0, say, and replace the coefficients in the fit with desired parameter values. Then simulated responses can be obtained using 'simulate'. – NRH May 16 '12 at 18:43
• That is true. The packages in the [experimental design task] (cran.r-project.org/web/views/ExperimentalDesign.html) view has many packages to create the dummy variable setup. Integrating the simulate and experimental design functions would be great! – Etienne Low-Décarie May 25 '12 at 11:14
• Some simulate functions do not create the data only the fit: stats.stackexchange.com/questions/11233/… – Etienne Low-Décarie Jun 18 '12 at 16:06

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)

• I can't seem to find the new TeachingDemos on R-Forge? – Etienne Low-Décarie Jun 18 '12 at 15:57
• install.packages("TeachingDemos", repos="R-Forge.R-project.org") – Etienne Low-Décarie Jun 18 '12 at 16:02
• Try install.packages("TeachingDemos", repos="http://R-Forge.R-project.org"). You need "http://" at the beginning and no ";". – Greg Snow Jun 18 '12 at 19:36

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.

• Thanks, but this seems like @Momo's answer. – Etienne Low-Décarie May 16 '12 at 13:01