Dear list members,

I am trying to duplicate SAS results from the SAS procedure HPMixed in R's 'lme4' and' lmerTest' and having a little bit of trouble.  The data is named "cutsim" and can be found in the following location: https://dl.dropboxusercontent.com/u/54480631/Cutsim.csv

The data is a reduced version of the "All the cows in Canada" study presented in the following paper: http://support.sas.com/resources/papers/proceedings09/256-2009.pdf. The data contains 400000 observations and includes 1 random effect (animal).  The model is specified in an identical fashion to what is presented in the paper using both software packages with an additional statement to estimate the fixed effects estimates.  Furthermore, the author of the paper specifies the fixed effects as Species and Farm*Species but does not include Farm by itself as a main effect

As shown below, the variance estimates are nearly identical in both software packages. However, the numerator degrees of freedom and the F-values are not matching.  I believe that the source of the problem might be that I am incorrectly representing the number of numerator degrees of freedom in the data set but I am not 100 percent certain that this is the problem. In addition, the lmertest anova includes Farm as a fixed effect even though it is not specified as such in the lmer statement.  Correcting this may also help better reflect the results presenting in the sas model.

Do you happen to know what I can do to get the results to match up a little better?



I have specified the model in R as follows:

>cutsim <- read.csv(file="c:/Users/Jump1180/Dropbox/Public/cutsim.csv", header=T)

>library("lme4")
>library("lmerTest")

>model2=lmer(Yield ~ Species + Species*Farm + (1|Animal), data=cutsim, REML=TRUE, family=gaussian)

>model2
>fixef(model2)
>anova(model2, ddf="lme4")


The R results are shown below:

> model2
Linear mixed model fit by REML
Formula: Yield ~ Species + Species * Farm + (1 | Animal)
   Data: cutsim
   AIC   BIC  logLik deviance REMLdev
2e+06 2e+06 -999829    2e+06   2e+06
Random effects:
Groups   Name        Variance Std.Dev.
Animal   (Intercept) 4.0139   2.0035
 Residual             8.0479   2.8369
Number of obs: 400000, groups: Animal, 10000

Fixed effects:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.8318926  0.0989447     8.4   <2e-16 ***
Species       1.0490039  0.0295851    35.5   <2e-16 ***
Farm          1.0022088  0.0016881   593.7   <2e-16 ***
Species:Farm -0.0005319  0.0005026    -1.1     0.29
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Specis Farm
Species     -0.907
Farm        -0.870  0.788
Species:Frm  0.791 -0.871 -0.907
> fixef(model2)
(Intercept)      Species         Farm Species:Farm
 0.831892610  1.049003854  1.002208792 -0.000531927
> anova(model2, ddf="lme4")
Analysis of Variance Table
             Df   Sum Sq  Mean Sq    F value
Species       1    56495    56495 7.0198e+03
Farm          1 15796032 15796032 1.9627e+06
Species:Farm  1        9        9 1.0852e+00





The SAS model is specified as follows:

proc hpmixed data=CutSim;
      class Species Farm Animal;
      model Yield = Species Farm*Species /ddfm=residual;
      random animal /cl ;
test Species Farm*Species;
run;


Here is a portion of the SAS results:


Covariance Parameter
Estimates

Cov Parm    Estimate

Animal        4.0297
Residual      8.0479

Type III Tests of Fixed Effects

Num      Den
Effect            DF       DF    F Value    Pr > F

Species            4      4E5    1158.56    <.0001
Species*Farm     495      4E5    3969.97    <.0001


Research Computing Support
University of Tennessee -  Office of Information Technology
524 Greve Hall,  821 Volunteer Blvd
Knoxville, Tennessee, 37996-3395

Email: Josh.Price@utk.edu     Phone: 865-974-1595     Fax:865-974-3960


	[[alternative HTML version deleted]]

