Comparison to SAS’ PROC MIXED

[W-APrince]

We’re currently evaluating the suitability of the {mmrm} package for use in our biostatistics department, but we are experiencing some issues and are getting different estimates/p-values between the 2 softwares.
The issue seems to be related to the covariate matrix calculations – the unstructured, compound symmetry, and heterogeneous compound symmetry matrix structures seem to be okay, however, we seem to be seeing differences in the others.
For example, when we run the following model on the bcva_data dataset included in the mmrm package:
• R Code:

data("bcva_data", package = "mmrm")
bcva_data$ARMCD <- as.factor(bcva_data$ARMCD)
levels(bcva_data$ARMCD)
bcva_data$ARMCD <- relevel(bcva_data$ARMCD, ref = "CTL")

contrasts(bcva_data$ARMCD) <- contr.treatment(levels(bcva_data$ARMCD), base = which(levels(bcva_data$ARMCD) == "CTL"))
contrasts(bcva_data$AVISIT) <- contr.treatment(levels(bcva_data$AVISIT))

levels(bcva_data$AVISIT)
levels(bcva_data$ARMCD)
write_xpt(bcva_data, "./bcva.XPT", version = 5)

mmrm_model <- mmrm(
  BCVA_CHG ~ ARMCD + AVISIT + ARMCD * AVISIT + BCVA_BL + ar1(AVISIT | USUBJID),
  data = bcva_data,
  reml = TRUE,
  method = "Kenward-Roger",
  vcov = "Kenward-Roger-Linear"
)

• SAS PROC MIXED Code:

<read in XPT file written via R with references corresponding to factor levels>

proc mixed data=bcva  method = reml;
   class ARMCD(ref="1") AVISIT(ref="1") USUBJID;
   model BCVA_CHG = ARMCD AVISIT BCVA_BL ARMCD*AVISIT / solution ddfm=KR;
   repeated AVISIT / subject=USUBJID type=AR(1) r rcorr; /* Auto-regressive (AR1) covariance matrix */
   lsmeans ARMCD*AVISIT/ diff cl alpha=.05;

   /* Output the fixed effects estimates and related statistics */
   ods output SolutionF=fixed_effects; 
run;

We are experiencing the following differences with a 1E-3 tolerance (note, there are more but these are some of the more extreme differences with regard to the p-value):

<style> </style>

TYPE	effect	estimate	StdErr	df	tValue	Probt
R Result	ARMCDTRT:AVISITVIS02	0.191755	0.142338	6073.801	1.347177	0.177973
SAS Result	ARMCDTRT:AVISITVIS02	0.186703	0.173244	7581.078	1.077687	0.281208
Difference	....................	-0.00505	0.030906	1507.278	-0.26949	0.103234
R Result	ARMCDTRT:AVISITVIS03	0.475896	0.164868	8482.299	2.886522	0.003905
SAS Result	ARMCDTRT:AVISITVIS03	0.465548	0.173933	7633.822	2.676596	0.007453
Difference	....................	-0.01035	0.009064	-848.477	-0.20993	0.003548
R Result	AVISITVIS02	0.275391	0.101673	6088.213	2.708602	0.006776
SAS Result	AVISITVIS02	0.276339	0.123635	7587.288	2.235119	0.025439
Difference	....................	0.000948	0.021962	1499.075	-0.47348	0.018663
R Result	AVISITVIS03	0.46004	0.117666	8480.155	3.90971	9.31E-05
SAS Result	AVISITVIS03	0.464771	0.12412	7645.261	3.744529	0.000182
Difference	....................	0.004731	0.006454	-834.894	-0.16518	8.89E-05
R Result	BCVA_BL	-0.00071	0.002757	2227.454	-0.25665	0.797473
SAS Result	BCVA_BL	-0.00073	0.002635	2632.496	-0.27637	0.782282
Difference	....................	-2.1E-05	-0.00012	405.0417	-0.01972	-0.01519

Do you know what is causing these mismatches? Is this something obvious we’re missing, or is there a difference in the way these softwares are calculating the results?

[danielinteractive]

Hi @W-APrince , thanks for reaching out!

Two ideas:

1. With AR1, you don't need vcov = "Kenward-Roger-Linear" to get close to SAS results, therefore I would omit that argument from the mmrm() call.
2. Please double check that the order of the AVISIT levels is the same in R and in SAS. The problem is that SAS puts reference levels last while R puts them first. This would then lead to a different order, and hence results. See e.g. here for SAS code we used for a comparison, putting the last level as reference.

In general, I would say you don't need to worry too much about the matching of SAS results - many many people and companies have meanwhile tested this package extensively (plus we have done all the testing in the integration tests shipped with this package).

Does this help?

[SimonTBate]

Hi @danielinteractive and @W-APrince ,

Thanks for the suggestions on this one. I'm pleased this issue has been opened as I too have been finding AR(1) models tricky to align between Proc Mixed and mmrm.

Daniel, your suggestions make a lot of sense to me and I assumed that re-parameterizing was the answer.

However, when I tried it out I unfortunately found that:

1. The inclusion of the vcov option seemed to improve the alignment between SAS and R (slightly).
2. Re-ordering the order of the "Time" factor (in my example) to align SAS and R made things different (and slightly worse) than the default (different) ordering of the factor levels.

Note I was mainly looking at the overall effects tests using the car and emmeans packages under the contr.treatment option in R:

model <- mmrm(formula = Response ~ Treatment+Day+Treatment*Day + ar1(Day | Subject), data=Data, method = "Kenward-Roger", vcov="Kenward-Roger-Linear")
Anova(model, type = "III", ddf = "Kenward-Roger")
joint_tests(model, cov.reduce = symmint(0))

[danielinteractive]

Hi @SimonTBate , so for Type III test results we are aware that there is some issues at the moment, see also the discussion in #502. We are actively working on this at the moment.

On the other hand, I would expect that all other model results for AR(1) are very close between mmrm and SAS PROC MIXED/GLIMMIX. Can you confirm that?

[SimonTBate]

Hi @danielinteractive , yes I can confirm that all the other results (well at least the ones I've looked at so far) do agree. I've looked at em means, comparisons between the em means and fit statistics (AIC, BIC, log-likelihood)

[danielinteractive]

Hi @SimonTBate , can you please try with the newly released mmrm package version from CRAN if this resolves your issue with the Type III test discrepancies?