The comparison of output of SAS AND R for mixed model

The comparison of output of SAS AND R for mixed model

Data process

library(nlme)
setwd("C:\\Users\\hed2\\OneDrive - National Institutes of Health\\Mixed model by SAS and R")
head(Orthodont)

## Grouped Data: distance ~ age | Subject
##   distance age Subject  Sex
## 1     26.0   8     M01 Male
## 2     25.0  10     M01 Male
## 3     29.0  12     M01 Male
## 4     31.0  14     M01 Male
## 5     21.5   8     M02 Male
## 6     22.5  10     M02 Male

write.csv(Orthodont,  "Orthodont.csv" ,na="",row.names = FALSE)

Fit the models with unstructured covariance structure- default

# fit the models with unstructured covariance structure of random effect
model1 <- lme(distance ~ age+Sex, data=Orthodont, random= list(Subject = pdSymm(~ 1+age+Sex ))) 
summary(model1)

## Linear mixed-effects model fit by REML
##   Data: Orthodont 
##        AIC   BIC    logLik
##   453.6604 480.2 -216.8302
## 
## Random effects:
##  Formula: ~1 + age + Sex | Subject
##  Structure: General positive-definite
##             StdDev    Corr         
## (Intercept) 3.0339255 (Intr) age   
## age         0.2264246 -0.827       
## SexFemale   2.8221068 -0.984  0.795
## Residual    1.3100420              
## 
## Fixed effects:  distance ~ age + Sex 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 17.710236 0.9091362 80 19.480290  0.0000
## age          0.660185 0.0712530 80  9.265371  0.0000
## SexFemale   -1.660111 0.7888414 25 -2.104493  0.0456
##  Correlation: 
##           (Intr) age   
## age       -0.864       
## SexFemale -0.581  0.335
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.983655835 -0.455860066  0.005090891  0.472520364  3.906689751 
## 
## Number of Observations: 108
## Number of Groups: 27

# https://rpubs.com/samuelkn/CovarianceStructuresInR

model2 <- lme(distance ~ age + Sex  , random = ~1+age+Sex  | Subject, data = Orthodont) 
summary(model2)

## Linear mixed-effects model fit by REML
##   Data: Orthodont 
##        AIC   BIC    logLik
##   453.6604 480.2 -216.8302
## 
## Random effects:
##  Formula: ~1 + age + Sex | Subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr         
## (Intercept) 3.0339748 (Intr) age   
## age         0.2264291 -0.827       
## SexFemale   2.9224528 -0.983  0.768
## Residual    1.3100383              
## 
## Fixed effects:  distance ~ age + Sex 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 17.710238 0.9091410 80 19.480188  0.0000
## age          0.660185 0.0712534 80  9.265319  0.0000
## SexFemale   -1.660143 0.7888386 25 -2.104541  0.0456
##  Correlation: 
##           (Intr) age   
## age       -0.864       
## SexFemale -0.581  0.335
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.983630079 -0.455863516  0.005090832  0.472520465  3.906698506 
## 
## Number of Observations: 108
## Number of Groups: 27

In SAS, the corresponding agument: type = un

# PROC IMPORT DATAFILE="...\Mixed model by SAS and R\orthodont.csv" DBMS=csv
# OUT=orthodont;
# RUN;

# PROC MIXED DATA=orthodont ;
# CLASS sex  subject   ;
# MODEL distance =  sex age /solution;
# RANDOM  sex age/ subject=subject type=un ;
# RUN; QUIT;

Fit the models with compound symmetric covariance structure

model11 <- lme(distance ~ age+Sex, data=Orthodont, random= list(Subject = pdCompSymm(~ 1+age+Sex ))) 
summary(model1)

## Linear mixed-effects model fit by REML
##   Data: Orthodont 
##        AIC   BIC    logLik
##   453.6604 480.2 -216.8302
## 
## Random effects:
##  Formula: ~1 + age + Sex | Subject
##  Structure: General positive-definite
##             StdDev    Corr         
## (Intercept) 3.0339255 (Intr) age   
## age         0.2264246 -0.827       
## SexFemale   2.8221068 -0.984  0.795
## Residual    1.3100420              
## 
## Fixed effects:  distance ~ age + Sex 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 17.710236 0.9091362 80 19.480290  0.0000
## age          0.660185 0.0712530 80  9.265371  0.0000
## SexFemale   -1.660111 0.7888414 25 -2.104493  0.0456
##  Correlation: 
##           (Intr) age   
## age       -0.864       
## SexFemale -0.581  0.335
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.983655835 -0.455860066  0.005090891  0.472520364  3.906689751 
## 
## Number of Observations: 108
## Number of Groups: 27

# https://stats.stackexchange.com/questions/26556/equivalent-mixed-models-yielding-different-results-in-sas

In SAS, the corresponding type = cs

# PROC MIXED DATA=orthodont ;
# CLASS sex  subject   ;
# MODEL distance =  sex age /solution;
# RANDOM  sex age/ subject=subject type=CS ;
# RUN; QUIT;

Other covariance structures

Name	nlme
Compound symmetry	corCompSymm
AR1	corAR1
CAR1	corCAR1
Unstructured	CorSymm

# Many covariance structures are reasonable for these data, such as unstructured (TYPE=UN), autoregressive (TYPE=AR(1)), compound symmetric (TYPE=CS), Toeplitz (TYPE=TOEP), etc. 
# http://support.sas.com/kb/37/107.html

Variance Components (VC) is the default setting in proc mixed in SAS, but unstructured (UN) in R.