[R-sig-ME] PROC MIXED RANDOM equivalence in R nlme

[Ben Bolker]

Or make them equivalent by using pdDiag() to enforce a diagonal
covariance matrix.

FitTHC <- lme(ln_thc ~ rv + t5 + t9 + ar + ol + ox + su + bz,
          data = emiss,
          random = list(study_vehicle=pdDiag(~1 + rv + t5 + t9 + ar + ol
+ ox + su + bz)))

In lme4 you can do this with the double-bar operator.

On 17-08-11 05:27 PM, Jake Westfall wrote:
> Dennis,
> 
> The two models are not equivalent because the R model uses an unstructured
> covariance for the random effects (estimates all possible correlations
> among random effects) while the SAS model uses a diagonal covariance matrix
> (assumes no correlations among random effects). You can make them
> equivalent by adding type=un to the end of the RANDOM statement in the SAS
> model and you'll find that the SAS model takes a lot longer to run then :)
> 
> Jake
> 
> On Fri, Aug 11, 2017 at 4:03 PM, Dennis F. Kahlbaum <kbomb at umich.edu> wrote:
> 
>> I am trying to reproduce some old SAS PROC MIXED code using R and nlme.
>>
>> As background, the data consist of 300+ vehicle emission test results
>> (grams/mile of THC over a given driving cycle) versus the gasoline fuel
>> properties, which include: Read Vapor Pressure (rv), Temperature for 50%
>> Evaporation (t5), Temperature for 90% Evaporation (t9), Aromatics (ar),
>> Olefins, (ol), Sulfur (su), Oxygenate (ox), Benzene (bz), and a code
>> identifying the study and vehicle (study_vehicle). All variables are real
>> numbers except "study_vehicle", which is character. Unfortunately, since
>> the data are confidential, I'm unable to provide an example.
>>
>> The goal is to determine the relationship between the emissions and the
>> fuel properties.
>>
>> The original SAS v6.12 code is provided below:
>>
>> ------------------------------------------------------------------
>> /* READ DATA */
>> DATA emiss;
>>    INFILE 'data.tab' LRECL=8000 FIRSTOBS=2 DLM='09'X MISSOVER DSD;
>>    INPUT study_vehicle $ thc rv t5 t9 ar ol ox su bz;
>>
>> /* CREATE NEW VARIABLE */
>> ln_thc = log (thc);
>>
>> /* PERFORM ANALYSIS RELATING LOG EMISSIONS TO FUEL PROPERTIES */
>> PROC MIXED DATA=emiss MAXITER=1000 CONVH=1E-8 METHOD=REML NOCLPRINT
>> NOITPRINT;
>> CLASS study_vehicle;
>>
>> MODEL ln_thc = rv t5 t9 ar ol ox su bz
>>                 /S DDFM=RES;
>>
>> RANDOM         int rv t5 t9 ar ol ox su bz
>>                /SUB=study_vehicle;
>> RUN;
>> ------------------------------------------------------------------
>>
>> The R code I've devised for reading/processing the data and to replicate
>> the PROC MIXED statement is shown below:
>>
>> ------------------------------------------------------------------
>> # LOAD REQUIRED LIBRARY
>> library (“lmne”)
>>
>> # FORCE R TO USE THE SAME "CONTRAST" METHOD AS SAS
>> options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))
>>
>> # READ DATA (study_vehicle, thc, rv, t5, t9, ar, ol, ox, su, bz)
>> emiss <- read.table("data.csv", header=TRUE, sep=",") # NA for missing
>> data (default)
>>
>> # CREATE NEW VARIABLE
>> ln_thc <- log(thc)
>>
>> # PERFORM ANALYSIS RELATING LOG EMISSIONS TO FUEL PROPERTIES
>> FitTHC <- lme(ln_thc ~ rv + t5 + t9 + ar + ol + ox + su + bz,
>>           data = emiss,
>>           random = ??????? | study_vehicle )
>>
>> # OUTPUT THE RESULTS
>> summary(FitTHC)
>> ------------------------------------------------------------------
>>
>> A review of the dataframe "emiss" shows that all variables, except
>> "study_vehicle", are of type "num". As expected, the character variable
>> "study_vehicle" is of type "factor".  These properties match those from SAS.
>>
>> As indicated by the questions marks, the problem I'm having is in
>> constructing the equivalent R random code for the SAS RANDOM and any
>> remaining settings (e.g. "/S DDFM=RES"). I've tried
>>
>> random = ~1 + rv + t5 + t9 + ar + ol + ox + su + bz | study_vehicle)
>>
>> but R gets caught in a processing loop that produces no errors or
>> warnings, and no results, even after 15+ minutes of execution time. As I
>> recall, the SAS code took less than 15 seconds back in 1999.
>>
>> I have also tried
>>
>> random = ~1 | study_vehicle)
>>
>> and R immediately produces output.  However, although the regression
>> coefficients are close to those from SAS, the degrees of freedom from SAS
>> are higher than those from R.
>>
>> Therefore, what is the equivalent R random code for the original SAS
>> RANDOM?
>>
>> Thanks!
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>