Contrast ordering - dicrepancy between SAS PROC MIXED and {mmrm}

[kkmann]

Summary

When using KR for df a type 3 ANOVA produces different results between PROC MIXED and {mmrm}. Reordering contrasts manually seems to fix the issue (see below).

Is there a problem / inconsistency with the way contrasts for interaciton terms are sorted?

I am raising this issue on behalf of @InkaRo (thanks again for making me aware of the issue and the excellent diagnosis!).

---

SAS Program (fev data imported via csv)

proc mixed data=ads method=reml order=internal cl;
   class USUBJID RACE(ref='Asian') AVISIT(ref='VIS1');
   model FEV1 = RACE AVISIT RACE*AVISIT FEV1_BL /alpha=0.1 cl ddfm=kr e3 S;
   repeated AVISIT /subject=USUBJID type=un r rcorr; 

   * Type 3 tests of fixed effects, Covariance matrix and LSMeans *;
   ods output tests3 = No1_type3_un 
              r = No1_cov_matrix_r_un
              lsmeans = No1_lsmeans_un; 

SAS output

---

R code to show discrepancy & manual fix.

library(mmrm)
> library(car)
Loading required package: carData
> mod <- mmrm(
+     formula = FEV1 ~  RACE + AVISIT + RACE * AVISIT + FEV1_BL + us(AVISIT|USUBJID),
+     data = fev_data,
+     control = mmrm_control(
+         method = "Kenward-Roger",
+         vcov = "Kenward-Roger-Linear"
+       ),
+     reml = TRUE
+   )
> car::Anova(mod, type = "III") |>
+   print()
Analysis of Fixed Effect Table (Type III F tests)
            Num Df Denom Df F Statistic   Pr(>=F)    
RACE             2   175.66      28.044 2.693e-11 ***
AVISIT           3   202.20     114.386 < 2.2e-16 ***
FEV1_BL          1   197.21      24.404 1.662e-06 ***
RACE:AVISIT      6   198.21       0.916    0.4844    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> tol <- sqrt(.Machine$double.eps)
> # with {mmrm} functions
> mmrm:::h_get_contrast(mod, "AVISIT", "3", tol) |>
+   print()
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]      [,8]      [,9]     [,10]     [,11]     [,12]
[1,]    0    0    0    1    0    0    0 0.3333333 0.0000000 0.0000000 0.3333333 0.0000000
[2,]    0    0    0    0    1    0    0 0.0000000 0.3333333 0.0000000 0.0000000 0.3333333
[3,]    0    0    0    0    0    1    0 0.0000000 0.0000000 0.3333333 0.0000000 0.0000000
         [,13]
[1,] 0.0000000
[2,] 0.0000000
[3,] 0.3333333
> df_md(mod, mmrm:::h_get_contrast(mod, "AVISIT", "3", tol)) |>
+   print()
$num_df
[1] 3

$denom_df
[1] 202.1992

$f_stat
[1] 114.3857

$p_val
[1] 2.472796e-43

> # manually change contrasts for interaction terms
> ctr_changed <- mmrm:::h_get_contrast(mod, "AVISIT", "3", tol)
> ctr_changed[, 10] <- ctr_changed[, 9]
> ctr_changed[, 11] <- ctr_changed[, 10]
> ctr_changed[, 9] <- ctr_changed[, 8]
> ctr_changed[, 12] <- ctr_changed[, 13]
> print(ctr_changed)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]      [,8]      [,9]     [,10]     [,11]     [,12]
[1,]    0    0    0    1    0    0    0 0.3333333 0.3333333 0.0000000 0.0000000 0.0000000
[2,]    0    0    0    0    1    0    0 0.0000000 0.0000000 0.3333333 0.3333333 0.0000000
[3,]    0    0    0    0    0    1    0 0.0000000 0.0000000 0.0000000 0.0000000 0.3333333
         [,13]
[1,] 0.0000000
[2,] 0.0000000
[3,] 0.3333333
> print(df_md(mod, ctr_changed))
$num_df
[1] 3

$denom_df
[1] 154.8794

$f_stat
[1] 130.6001

$p_val
[1] 3.243576e-42

R session info

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Dublin
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] car_3.1-2        carData_3.0-5    mmrm_0.3.12.9000

loaded via a namespace (and not attached):
 [1] utf8_1.2.4        generics_0.1.3    stringi_1.8.3     lattice_0.21-9   
 [5] digest_0.6.34     magrittr_2.0.3    grid_4.3.2        pkgload_1.3.3    
 [9] fastmap_1.1.1     Matrix_1.6-1.1    pkgbuild_1.4.3    sessioninfo_1.2.2
[13] backports_1.4.1   urlchecker_1.0.1  promises_1.2.1    purrr_1.0.2      
[17] fansi_1.0.6       abind_1.4-5       Rdpack_2.6        cli_3.6.2        
[21] shiny_1.8.0       rlang_1.1.3       rbibutils_2.2.16  ellipsis_0.3.2   
[25] remotes_2.4.2.1   cachem_1.0.8      devtools_2.4.5    tools_4.3.2      
[29] parallel_4.3.2    memoise_2.0.1     checkmate_2.3.1   httpuv_1.6.13    
[33] vctrs_0.6.5       R6_2.5.1          mime_0.12         lifecycle_1.0.4  
[37] stringr_1.5.1     fs_1.6.3          htmlwidgets_1.6.4 usethis_2.2.2    
[41] miniUI_0.1.1.1    pkgconfig_2.0.3   pillar_1.9.0      later_1.3.2      
[45] glue_1.7.0        profvis_0.3.8     Rcpp_1.0.12       tibble_3.2.1     
[49] rstudioapi_0.15.0 xtable_1.8-4      htmltools_0.5.7   nlme_3.1-163     
[53] TMB_1.9.10        compiler_4.3.2   

[danielinteractive]

Thanks @kkmann , did you check the vignette details? I remember @clarkliming described known differences between our implementation and SAS: https://openpharma.github.io/mmrm/latest-tag/articles/hypothesis_testing.html#hypothesis-testing-in-sas

[clarkliming]

this is a bug that the order of the interaction term will design matrix. using the following model it is correct

mod2 <- mmrm(
  formula = FEV1 ~  AVISIT + RACE + AVISIT * RACE + FEV1_BL + us(AVISIT|USUBJID),
  data = fev_data,
  control = mmrm_control(
      method = "Kenward-Roger",
      vcov = "Kenward-Roger-Linear"
    ),
  reml = TRUE
)

apply(mmrm:::h_get_contrast(mod2, "AVISIT", "3", 1e-8), 1, \(x) names(coef(mod2))[x != 0])

I will fix this issue

[kkmann]

Hey, sry - just coming back from holidays :)

Thanks for looking into it @clarkliming!

[clarkliming]

Hi @kkmann I have complete my fix in #460, would you please also have a look? (maybe there are corner cases that I did not cover)

[kkmann]

will take a look asap, maybe we can stall merging until next week?

[clarkliming]

will take a look asap, maybe we can stall merging until next week?

thank you, we can wait until the checks are finished