Dr. Bates and Luca,

Thank you so much for sending this over.  This is a huge time saver and I really appreciate the help from both of you.  Also, thank you so much for submitting this as a batch program.  After adjusting the variables to be categorical I was running into a huge roadblock with running out of memory and looking out how to submit this to our local cluster.  However after seeing this email your analysis was exactly what I needed.

Thanks again for all the assistance,
Josh

From: dmbates@gmail.com [mailto:dmbates@gmail.com] On Behalf Of Douglas Bates
Sent: Wednesday, May 29, 2013 11:08 AM
To: lborger
Cc: Price, Josh; R-mixed models mailing list
Subject: Re: [R-sig-ME] Duplicating SAS Proc HPMixed results in lme 4 and lmerTest

I enclose an R source file and the output from running this code in batch mode using

R CMD BATCH --vanilla cutsim.R

As Luca mentioned, it is important to convert the categorical variables Animal, Farm and Species to factors in R before fitting models.  In general it is advisable to check the structure of a data frame, using the str() function, and a summary of the columns, using summary(), before fitting models.  This is especially important when dealing with large data sets.

One additional check that I did was to check whether Animal is nested within Species, which one would expect to be true, and also nested within Farm.  By extracting the unique rows of the subset of the data frame formed by taking only the Animal, Species and Farm columns and looking at the row count we see that there are 10000 rows which is the number of distinct Animals.  Thus Species and Farm are nested within Animal. (If they weren't there would be more than 10000 rows.)

The first model that I fit has a fixed-effect for Species and random effects for Animal and Farm:Species.  I know you wanted to model Farm:Species as a fixed effect but with 500 levels I would generally use a random effect for that factor.  The next model is the one you wanted to fit.  As suggested by Luca, it will take a long time and use a lot of memory.  This is because there are fixed-effects for over 500 coefficients.  The model matrix for the fixed-effects is very sparse (over 99% of the elements are zeros) but it is being treated as a dense matrix.  This is why fitting the model takes over 10 times as long as the previous model fit.

To test the significance of the Farm:Species interaction I fit the model without the term and compare it to the model with the term.  In this case the evidence of the significance of the term is overwhelming, although that could be seen from the first model fit too.  The reason I use maximum likelihood fits for the models is to allow the likelihood ratio test.


	[[alternative HTML version deleted]]

