[R] PROC MIXED vs lme. Split-plot with heterogeneous variances.

[Kevin Wright]

> A curious difference between SAS and R.  I wonder if anyone can explain it.
> 
> Basic idea: Split-plot design (Male = whole plot, Trt = Sub plot).  Rep is random, Rep*Male variance component is 0 and deleted.  Heterogeneous variances - each Trt has different variance.
> 
> The model with only Male and Trt as fixed effects is the same in SAS and R.  
> 
> When I add Male:Trt interaction, the results of the tests of fixed effects are no longer the same (comparing SAS and R) for Male, but are the same for Trt and Male:Trt.
> 
> Is my model specification incorrect?
> 
> Kevin Wright.  Details follow.
> 
> 
> Model with Male, Trt as fixed effects
> 
> proc mixed data=pollen ;
>   class Trt Rep Male;
>   model Yield = Male Trt;
>   random Rep;
>   repeated / group = Trt;
>   lsmeans Trt Male;
> run;
> 
>                                 Type 3 Tests of Fixed Effects
> 
>                                       Num     Den
>                         Effect         DF      DF    F Value    Pr > F
> 
>                         Male            1      47       3.64    0.0624
>                         Trt             9      47       3.80    0.0012
> 
> 
> 
> > pollen.hetero<-lme(Yield~Male+Trt,pollen,random=list(Rep=~1),
> +                    weights=varIdent(form=~1|Trt))
> > 
> > anova(pollen.hetero)
>             numDF denDF   F-value p-value
> (Intercept)     1    47 14.613222  0.0004
> Male            1    47  3.640521  0.0625
> Trt             9    47  3.797328  0.0012
> 
> 
> 
> Now add Male:Trt interaction as a fixed effect.
> 
> proc mixed data=pollen ;
>   class Trt Rep Male;
>   model Yield = Male Trt Male*Trt;
>   random Rep;
>   repeated / group = Trt;
>   lsmeans Trt Male;
> run;
> 
> 
>                                 Type 3 Tests of Fixed Effects
> 
>                                       Num     Den
>                         Effect         DF      DF    F Value    Pr > F
> 
>                         Male            1      38       0.39    0.5384
>                         Trt             9      38       8.40    <.0001
>                         Trt*Male        9      38       2.97    0.0090
> 
> > pollen.hetero<-lme(Yield~Male+Trt+Male:Trt,pollen,random=list(Rep=~1),
> +                    weights=varIdent(form=~1|Trt))
> > 
> > anova(pollen.hetero)
>             numDF denDF   F-value p-value
> (Intercept)     1    38 27.007016  <.0001
> Male            1    38  8.431796  0.0061
> Trt             9    38  8.396913  <.0001
> Male:Trt        9    38  2.964672  0.0090
> 
> 
>