openpharma / mmrm

Mixed Models for Repeated Measures (MMRM) in R.
https://openpharma.github.io/mmrm/
Other
135 stars 23 forks source link

fix anova type 3 #460

Closed clarkliming closed 2 months ago

clarkliming commented 3 months ago

close #459

clarkliming commented 3 months ago

after #455 merged we may further improve the code

github-actions[bot] commented 3 months ago

badge

Code Coverage Summary

Filename                    Stmts    Miss  Cover    Missing
------------------------  -------  ------  -------  ----------------------------
R/between-within.R             59       0  100.00%
R/component.R                  69       0  100.00%
R/cov_struct.R                 97       1  98.97%   407
R/empirical.R                   7       0  100.00%
R/fit.R                       229       3  98.69%   420, 481, 511
R/interop-car.R               130       1  99.23%   9
R/interop-emmeans.R            51       0  100.00%
R/interop-parsnip.R            59       1  98.31%   12
R/kenwardroger.R               92       2  97.83%   41, 63
R/mmrm-methods.R              140       0  100.00%
R/residual.R                    8       0  100.00%
R/satterthwaite.R             116      12  89.66%   238-249
R/skipping.R                    8       0  100.00%
R/testing.R                    64       4  93.75%   29, 31, 80, 82
R/tidiers.R                    72       2  97.22%   46-47
R/tmb-methods.R               283       3  98.94%   279-280, 336
R/tmb.R                       295       0  100.00%
R/utils-formula.R              27       0  100.00%
R/utils-nse.R                  16       0  100.00%
R/utils.R                     199      12  93.97%   279-289, 437, 466
R/zzz.R                        70      24  65.71%   7-22, 55-60, 90, 118, 138
src/chol_cache.h               63       0  100.00%
src/covariance.h              101       1  99.01%   177
src/derivatives.h             126       0  100.00%
src/empirical.cpp              72       0  100.00%
src/exports.cpp                47       0  100.00%
src/jacobian.cpp               47       1  97.87%   54
src/kr_comp.cpp                56       0  100.00%
src/mmrm.cpp                   76       0  100.00%
src/predict.cpp                93       0  100.00%
src/test-chol_cache.cpp        58       5  91.38%   9, 18, 26, 55, 62
src/test-covariance.cpp       123       5  95.93%   9, 29, 40, 61, 72
src/test-derivatives.cpp      108       7  93.52%   44, 53, 62, 85, 94, 106, 124
src/test-utils.cpp            195       7  96.41%   9, 16, 24, 34, 44, 57, 119
src/testthat-helpers.h         15       5  66.67%   36-37, 41, 50, 53
src/utils.h                    78       0  100.00%
TOTAL                        3349      96  97.13%

Diff against main

Filename           Stmts    Miss  Cover
---------------  -------  ------  -------
R/interop-car.R      +58      -2  +3.40%
TOTAL                +58      -2  +0.11%

Results for commit: 1bda9308b7eab0eee73425abea7165055193e6ad

Minimum allowed coverage is 80%

:recycle: This comment has been updated with latest results

github-actions[bot] commented 3 months ago

Unit Tests Summary

    1 files     46 suites   27s :stopwatch:   505 tests   465 :white_check_mark: 40 :zzz: 0 :x: 1 954 runs  1 908 :white_check_mark: 46 :zzz: 0 :x:

Results for commit 1bda9308.

:recycle: This comment has been updated with latest results.

github-actions[bot] commented 3 months ago

Unit Test Performance Difference

Additional test case details | Test Suite | $Status$ | Time on `main` | $±Time$ | Test Case | |:-----|:----:|:----:|:----:|:-----| | car | 👶 | | $+0.01$ | h_first_contain_categorical_works_as_expected | | car | 👶 | | $+0.26$ | h_get_contrast_works_even_if_only_interaction_term_exists | | car | 👶 | | $+0.18$ | h_get_contrast_works_even_if_the_interaction_term_order_changes | | car | 👶 | | $+0.18$ | h_get_contrast_works_for_higher_order_interaction | | car | 👶 | | $+0.06$ | h_get_contrast_works_if_intercept_is_not_given | | car | 👶 | | $+0.01$ | h_get_index_works_as_expected | | car | 👶 | | $+0.01$ | h_obtain_lvls_works_as_expected |

Results for commit 61873e583be98a5e891019fa78134fd5dc8d981c

♻️ This comment has been updated with latest results.

InkaRo commented 3 months ago

Thank you for the quick fix! At first glance, it works fine now. I wanted to run some more thorough checks and also have a look at 3-way interactions. Thanks again!

kkmann commented 3 months ago

It seems that for the 3 way case, there is still a discrepancy.

> library(mmrm)
> library(car)
> 
> mod <- mmrm(
+     formula = FEV1 ~ ARMCD + RACE + AVISIT + RACE * AVISIT * ARMCD + 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)    
ARMCD                  1   169.24      32.726 4.712e-08 ***
RACE                   2   168.43      24.655 4.058e-10 ***
AVISIT                 3   155.86     129.283 < 2.2e-16 ***
FEV1_BL                1   187.76      31.571 6.856e-08 ***
RACE:AVISIT            6   213.76       1.754    0.1099    
ARMCD:RACE             2   137.39       0.199    0.8201    
ARMCD:AVISIT           3   156.17       1.685    0.1725    
ARMCD:RACE:AVISIT      6   198.34       1.269    0.2732    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SAS reports

image

The difference is in the two-way interactions now, 3-way and marginals are fine.

I must admit that I am struggling a bit with the code in h_get_contrast() without any comments it is quite dense. Maybe it would help to flesh out the contrast definition in the documentation first?

clarkliming commented 3 months ago

Hi @kkmann I updated again, hopefully this cover all cases

kkmann commented 3 months ago

@InkaRo could you share the current discrepancy case directly in this PR so that we can discuss it here?

@clarkliming it seems that the contrast matrices are stil la bit off in the three way case. I am sorry, but it is quite difficult for Inka and me to reverse engineer the h_* functions without more specs / comments.

Could we maybe document hiw we expect the the h_get_contrast() function to behave and test it directly for higher order interactions?

https://github.com/openpharma/mmrm/blob/875322ede9d44446a397374d074116b39329350a/R/interop-car.R#L17-L28

clarkliming commented 2 months ago

hi @kkmann the details can be found at https://openpharma.github.io/mmrm/latest-tag/articles/hypothesis_testing.html; and this vignettes is just following the SAS. So in what cases do you observe a difference in the contrast?

danielinteractive commented 2 months ago

@kkmann @InkaRo @clarkliming From the sidelines: Maybe it would be easier to have a short meeting to clarify on this topic?

InkaRo commented 2 months ago

A meeting is a good idea. Do you want to set sth up @danielinteractive?

danielinteractive commented 2 months ago

@InkaRo we actually have one today at 2.30 pm CET if you have time? You could let me know your email via daniel@rconis.com and I can add you

kkmann commented 2 months ago

Hi @clarkliming, @InkaRo and I just double checked on the latest published commit 875322ede9d44446a397374d074116b39329350a and the problem persists. The problem seems to be that the contrasts still need to be sorted in the order of appearance in the formula statement.

InkaRo commented 2 months ago

@clarkliming @danielinteractive sorry, my bad, I did not use load_all(), so I still had the old functions. It works fine now, contrasts look good.

kkmann commented 2 months ago

Apologies from my end as well. Missed it to. Thanks @clarkliming for seeing this through - hope we get that fix to CRAN soon x)