Simulating from M-SURGE

[sakrejda]

I am working on running simulations in R for my group based on output from M-SURGE models. One difficulty we have run into is generating design matrices within R which are compatible with the M-SURGE formulation.

Is it possible to call M-SURGE from the command line so that it outputs just the design matrices for a particular data set and GEMACO model? We need to do this since we are dealing with simulations and so the 'data set' is dynamic.

I am also looking at reconstructing the process of generating the matrices solely within R, but I would rather use M-SURGE directly to avoid future incompatibility.

Best regards,

Krzysztof Sakrejda-Leavitt

[CHOQUET]

Hello,
In reply to:

>I am working on running simulations in R for my group based on output >from M-SURGE models. One difficulty we have run into is generating design >matrices within R which are compatible with the M-SURGE formulation.

>Is it possible to call M-SURGE from the command line so that it outputs just >the design matrices for a particular data set and GEMACO model? We need >to do this since we are dealing with simulations and so the 'data set' is >dynamic.

No, there is no possibility to call M-SURGE from R.

>I am also looking at reconstructing the process of generating the matrices >solely within R, but I would rather use M-SURGE directly to avoid future >incompatibility.

From GEMACO, you can save each matrix of constraints in a text file.

Sincerely,

Rémi

[sakrejda]

For predictive population simulations where the population (and therefore design matrix) changes it is not sufficient to generate the design matrix from GEMACO.

It turns out the model matrix machinery in R can generate compatible model matrices. The solution is in using the 'contrasts.arg' argument to the 'model.matrix' function. 'contrasts.arg' expects a named list with the names referring to factors in the data frame. For all factors used in the M-SURGE model, pass the 'contrasts = FALSE' argument. So, e.g.- contrasts.arg = list( stage = list(contrasts = FALSE), species = list(contrasts = FALSE) ) ... once that argument is set, it's pretty straightforward to convert a GEMACO model formula into an R model formula for model.matrix

[sakrejda]

Correction, you have to call

contrasts( column, contrasts = FALSE )

for each factor column you will use as a predictor, then when making the call to model.matrix you can pass the resulting matrices as the 'contrasts.arg' arguments.

data(iris)
sp <- contrasts( iris$Species, contrasts = FALSE )
model.matrix( ~ Species, data = iris, contrast.args = list(Species = sp) )

I had wrapped it all in a function and forgotten about it. Apologies for the confusion.