Partly as a learning exercise, I’m attempting to reproduce a published SAS PROC MIXED analysis using R’s lme4 package. A full description of both the application and the statistical method is available in http://www.biomedcentral.com/1472-6750/11/15 . Briefly, the broad purpose is to assess the food safety of a specific genetically-modified crop (GMO). In the experiment, the GMO crop is grown along with a genetically similar but unmodified comparator and with a diverse set of commercially-available reference varieties. The reference varieties are presumed safe and thus provide estimates of the range of acceptable nutrient levels. The crop is harvested and subjected to chemical analyses to measure various constituents. Nutrient levels are then compared between the GMO and the comparator and between the GMO and the set of reference varieties.
The published method defines a variety group variable with 3 levels: GMO, comparator (comp), and reference (ref). An integer index is also defined that singles out the reference varieties. To clarify, assume that we have 4 reference varieties. The variety-related variables are then:
variety vgroup isref comp comp 0 GMO GMO 0 refA ref 1 refB ref 1 refC ref 1 refD ref 1
Focusing only on this part of the statistical model and ignoring the other (important but straightforward) model factors (i.e. locations & blocks), the published SAS PROC MIXED formulation is:
CLASS variety vgroup …; MODEL y=vgroup …; RANDOM isref*variety …;
Modeling vgroup as a fixed factor enables formal tests of (GMO – comp) and (GMO – ref). Note that isref is not designated as a CLASS variable and thus remains a numeric quantity. The SAS model runs successfully. A single variance is estimated for the isref*variety term, indicating that this quantity is a single factor rather than an interaction. The random effect predictions for the isref*variety factor are identically zero for the comp and GMO varieties and are non-zero for the ref varieties. I have yet to replicate these results in R lmer. The following model statement
y ~ vgroup + (0 + factor(isref)|variety) + …
is close, but doesn’t annihilate the isref == 0 variance estimate and the associated random effect predictions. Perhaps as a direct consequence, the vgroup standard errors for the comp and GMO levels are considerably larger than the corresponding SAS results.
So finally, is there a way to specify this exact model in R lmer? If so, how? If not, why?