The mmrm from openpharma

Add function `h_record_all_output()`

To do:

See design document here
Add Craig and Alessandro as authors in the DESCRIPTION file (see https://github.com/insightsengineering/rbmi for emails etc)
Add function h_record_all_output() in file R/output.R
Add unit tests in file testthat/test-output.R
- If needed please improve the function code
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Add function `h_cov_estimate()`

To do:

See design document here
Add function h_cov_estimate() in file R/covariance.R
Add unit tests in file testthat/test-covariance.R
- If needed please improve the function code
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Add `assert_data()` function including helpers

To do:

Please see design doc:

mmrm/design/design_fit.Rmd

Line 132 in a9b4972

assert_data <- function(vars, data) {
Add assert_data() function including helper functions in R/assert_data.R
Add roxygen2 documentation chunks for each function
Add unit tests in testtthat/test-assert_data.R
Add to pkgdown.yaml the function names

few things around `component()`

To do:

remove AIC, BIC, deviance, logLik since we have direct accessors already
remove method
- don't print it in print method
- I guess otherwise not needed
Add @details section with md list explaining the possible return values

Design the flow for fitting MMRM

To do:

Add `print` method for `mmrm_tmb`

This is nice for handling the results and will make testing of the package easier.

To do:

Add `summary()` method

This is important since we get the fixed effect estimates table from that.

To do:

Add homogeneous Toeplitz covariance structure

To do:

Robust sandwich estimator

To do:

Just fyi: Here is some simple example of gls() with clubSandwich - might be nice to play around with

library(nlme)
m <- gls(
  value ~ time*method,
  correlation =corCompSymm(form=~1|patient_id), 
  data = x, 
  na.action = na.omit, 
  method = "REML"
)
 
emmeans::emmeans(
  m, 
  specs = ~time*method,
  vcov = matrix(clubSandwich::vcovCR(m, type = "CR3"), nrow=6)
)

Add possibility to set `weights` to obtain weighted MMRM

To do:

Fix typo, improve error message for using wrong covariance structure string

Problem:

To do:

Search for arh1 and replace with ar1h across the package
Improve error message for using wrong covariance structure string
add test for this

Add vignette explaining model definition and fitting algorithm

Model definition
Maximum likelihood estimation
REML estimation

Find reference for corrected AIC formula

To do:

See AIC.mmrm_tmb method definition / code
Find a literature reference for the corrected AIC formula:
Cite that in the method
- import Rdpack (see https://github.com/insightsengineering/hermes/blob/b5f50a1c1788527d289e1f10c4219a9c227d4449/DESCRIPTION#L65 as example)
- specify Rdmacros (see https://github.com/insightsengineering/hermes/blob/b5f50a1c1788527d289e1f10c4219a9c227d4449/DESCRIPTION#L84 as example)
- see https://github.com/insightsengineering/hermes/blob/55f3c08a1b58ee32666326841df0f2ac40e40761/R/differential.R#L18
- include references in https://github.com/openpharma/mmrm/blob/main/inst/REFERENCES.bib

Add function `h_build_formula()`

To do:

See design document here
Add function h_build_formula() in file R/formula.R
Add unit tests in file testthat/test-formula.R
- If needed please improve the function code
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Add `h_free_cores()` function

To do:

See design document here
Add function h_free_cores() in file R/parallel.R
Add simple unit test in file testthat/test-parallel.R
- In this case e.g. just have one test that checks if the function runs silently without messages e.g.
- If needed please improve the function code
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Design for Type II and Type III tests a la `car::Anova()`

To do:

Read ?car::Anova
Go through code material:
- glmmTMB implementation: https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/R/Anova.R
- car package code: https://github.com/cran/car/tree/master/R
Write a design doc first (so no production package code yet)
- strategy (and several options for it if it makes sense)
- prototype code

Note (from ?car::Anova? Details):
The designations "type-II" and "type-III" are borrowed from SAS, but the definitions used here do not correspond precisely to those employed by SAS.

Type-II tests are calculated according to the principle of marginality, testing each term after all others, except ignoring the term's higher-order relatives; This definition of Type-II tests corresponds to the tests produced by SAS for analysis-of-variance models, where all of the predictors are factors, but not more generally (i.e., when there are quantitative predictors).
so-called type-III tests violate marginality, testing each term in the model after all of the others. Be very careful in formulating the model for type-III tests, or the hypotheses tested will not make sense.

Add `refit_multiple_optimizers()`

depends on #25

To do:

Add function `mmrm_tmb()` that fits an MMRM with `TMB`

Idea: Add a function that replaces glmmTMB::glmmTMB() call in

mmrm/R/fit.R

Line 37 in b0a1ff3

glmmTMB::glmmTMB(

.
(This can potentially also replace fit_single_optimizer, let's see)

To do:

Add antedependence correlation

To do:

Add vignette for covariance matrix models

To do:

Document details of each covariance matrix structure included
Resolve reference in the algorithm vignette

Avoid or improve `enum` mapping

Motivation: Currently we are mapping the covariance structures to their C++ enum integer codes here

mmrm/R/tmb.R

Line 169 in b1b6122

cov_type <- as.integer(switch(formula_parts$cov_type,

Ideal outcome: Can we avoid this completely and just pass the cov_type string (e.g. ar1) directly to C++ and then decide based on that?

Otherwise: Somehow do the enum conversion in a more robust way.

Check `NEWS` entry etc. for internal pre-release

Please re-read the NEWS.md file and make it human readable, i.e. re-organise, change order, combine entries, group them if applicable, full sentence, no new entries on the old version (!) etc.

Replace copyright holder

Use F. Hoffmann-La Roche Ltd., 2022 instead of openpharma. Can discuss again when somebody else is contributing code.

Add Heterogeneous 1st order autoregressive correlation

To do:

Comparison with other implementations / packages

For marketing/usability, we should provide users with a description as to how/why mmrm would be more appropriate than the other available packages.

In particular: lme4, nlme, and glmmTMB.

This documentation could fit well in a separate vignette "Comparison with other packages", so that potential users can immediately understand where/when we think our development can help them compared to packages with higher levels of adoption.

Design for Kenward-Roger d.f.

Material:

Start from implementation notes: https://github.com/runehaubo/lmerTestR/blob/master/pkg_notes/implementation.Rmd
Also look at part here: https://drive.google.com/file/d/1imAEKUPCbSieDMoozXjexkPBsl1YNkcV/view?usp=sharing
Look at code: https://github.com/runehaubo/lmerTestR/blob/35dc5885205d709cdc395b369b08ca2b7273cb78/R/contest.R#L282
How does it work together with pbkrtest: https://people.math.aau.dk/~sorenh/software/pbkrtest/
explore autodiff:jacobian in TMB
explore TMB plus Rcpp integration
create prototypes of deriving P, Q, R
match results with SAS

Allow `character` or `factor` for ID variable

Otherwise it is too cumbersome for users. Required for insightsengineering/tern.mmrm#79

To do:

If it is a character, then convert to factor inside.
If there are non-convergence related errors from the single optimizer fit, fail with those already, otherwise the error message is not useful.

Add more functions towards Satterthwaite and LS means functionality

To do:

Note that we don't need:

h_general_jac_list() (since no other jacobian calculations besides from our TMB fit)
vars() (since we can stick with the formula interface for this package)
h_vcov_theta() (since we have that explicitly as part of the object already)
separate mmrm() as another wrapper of fit_model() (since we don't need vars interface and we can just calculate Jacobian in fit_model() directly)
diagnostics() (this is more a wrapper for tern.mmrm)

Design a hex sticker for `mmrm`

To do:

Have fun!
Create svg and png versions
Upload to repo
Use in README and vignette
- see https://github.com/insightsengineering/hermes as example for how to do this

Add homogeneous ante-dependence covariance structure

To do:

Add introduction to introduction vignette

To do:

Add the main part "introduction" to the introduction vignette
To cover:
- simple model fit (similar like in README)
- changing to REML, changing optimizer, changing covariance structure
- extracting coefficients table from summary result, refer to mmrm_tmb_methods for others
- explain lower level function h_mmrm_tmb() - mainly useful if you don't need Satterthwaite d.f. or try other optimizers, "barebones" function
- Show how to use emmeans package, refer to emmeans_support

Support models where original design matrix is not full rank

Background:
It can happen that the design matrix that is created from the model formula does not have full rank. This can be due to some visit and arm combinations not being present in the data set, or too many (categorical) covariates that lead to exact collinearities between columns.

Reprex:

dat <- fev_data[11:25, ]
fit <- h_mmrm_tmb(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = dat
)

gives error

Error in h_mmrm_tmb_assert_start(tmb_object) : 
negative log-likelihood is NaN at starting parameter values

Comparing this with the lm result:

linmod <- lm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT,
  data = dat
)
summary(linmod)

we can see that the coefficients for ARMCDTRT:AVISITVIS2 and ARMCDTRT:AVISITVIS3 are not defined due to singularities. We also note that lm() has the corresponding argument singular.ok = TRUE controlling whether this works or gives an error.

To do:

Add `fit_model()`

To do:

See design document here
Add function fit_model() in file R/fit.R
Add unit tests in file testthat/test-fit.R
- If needed please improve the function code
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Add `h_summarize_all_fits()`

To do:

See design document here
Add function h_summarize_all_fits() in file R/fit.R
- Add unit tests in file testthat/test-fit.R
Add logLik method in file R/tmb-methods.R
- Add corresponding unit test
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Add function `h_labels()` with helpers

To do:

See design document here
Add function h_labels() in file R/labels.R
Add unit tests in file testthat/test-labels.R
- If needed please improve the function code
Add roxygen2 documentation chunks
- Run Build > Document, Build > Load All, and check that documentation looks good
Add function name to pkgdown.yaml
Run Build > Check to ensure that checks pass

Address warnings and notes (as much as possible) from `R CMD check`

On enableR:
The following notes and warnings are generated when R CMD check is run against this package.

* checking compilation flags used ... NOTE
Compilation used the following non-portable flag(s):
  ‘-Werror=format-security’ ‘-Wp,-D_FORTIFY_SOURCE=2’
  ‘-Wp,-D_GLIBCXX_ASSERTIONS’
* checking compiled code ... WARNING
File ‘mmrm/libs/mmrm.so’:
  Found ‘abort’, possibly from ‘abort’ (C)
    Objects: ‘mmrm.o’, ‘test-covariance.o’, ‘test-runner.o’,
      ‘test-utils.o’, ‘tmb.o’
  Found ‘printf’, possibly from ‘printf’ (C)
    Objects: ‘mmrm.o’, ‘test-covariance.o’, ‘test-runner.o’,
      ‘test-utils.o’, ‘tmb.o’
File ‘mmrm/libs/mmrm.so’:
  Found no calls to: ‘R_registerRoutines’, ‘R_useDynamicSymbols’

Compiled code should not call entry points which might terminate R nor
write to stdout/stderr instead of to the console, nor use Fortran I/O
nor system RNGs.
It is good practice to register native routines and to disable symbol
search.

See ‘Writing portable packages’ in the ‘Writing R Extensions’ manual.

To do:

be able to reproduce this - on enableR?...
fix it if possible
- it is not clear where the non-portable flags come from. Possibly from TMB. However this warning does not show up without CRAN flag or on R 4.2 so should be ok.
- it is not clear where the abort and printf symbols come from. I checked the C++ code and we don't have any of this or similar. Maximum we have is Rf_error which is standard R interface for an error so cannot be the problem. Note that this warning does not occur on R 4.2

On Daniel's Macbook with R CMD check --as-cran on the built package tarball:

* checking compiled code ... NOTE
Datei ‘mmrm/libs/mmrm.so’:
  Found no calls to: ‘R_registerRoutines’, ‘R_useDynamicSymbols’

It is good practice to register native routines and to disable symbol
search.

See ‘Writing portable packages’ in the ‘Writing R Extensions’ manual.

To do:

See https://stackoverflow.com/questions/42313373/r-cmd-check-note-found-no-calls-to-r-registerroutines-r-usedynamicsymbols how to fix this.
fix it
ongoing work, as it is more complicated than initially thought:
- see https://github.com/kaskr/adcomp/blob/00046bb67271fa3db266cabd5bd815b58fa0045d/TMB/inst/include/tmb_core.hpp
- see https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/init.h

Moreover parallel processes are not allowed in CRAN:

 ✖ | 2      47 | fit [0.7s]                                                      
  ────────────────────────────────────────────────────────────────────────────────
  Error (test-fit.R:175:3): mmrm falls back to other optimizers if default does not work
  Error in `.check_ncores(cores)`: 3 simultaneous processes spawned
  Backtrace:
    1. testthat::expect_silent(mmrm(formula, data_small))
         at test-fit.R:175:2
    9. mmrm::mmrm(formula, data_small)
   10. mmrm::refit_multiple_optimizers(fit, n_cores = n_cores, accept_singular = accept_singular)
   13. parallel::mclapply(...)
   14. parallel:::.check_ncores(cores)
  
  Error (test-fit.R:185:3): mmrm fails if no optimizer works
  Error in `.check_ncores(cores)`: 3 simultaneous processes spawned
  Backtrace:
    1. testthat::expect_error(...)
         at test-fit.R:185:2
    7. mmrm::mmrm(formula, data_small, reml = FALSE)
    8. mmrm::refit_multiple_optimizers(fit, n_cores = n_cores, accept_singular = accept_singular)
   11. parallel::mclapply(...)
   12. parallel:::.check_ncores(cores)
  ────────────────────────────────────────────────────────────────────────────────

To do:

disable parallel processing in the 2 tests

`_pkgdown.yaml` should only contain exported objects

Add compound symmetry covariance structures

To do:

Add a first spatial covariance structure: spatial exponential

To do:

Error when `AVISIT` column has SCREENING variable

please refer to the discussion on tern#628. This error should be more informative or fixed before the call.

To do:

define reprex for mmrm
include as unit tests
fix problem
confirm all checks pass

Add 1st order autoregressive correlation

To do:

Hex Sticker renders twice on index.html

The hex sticker renders twice on index.html.

This is because we need to add
<p align="center"> <img src="man/figures/logo.svg" align="right" alt="mmrm-logo" style="width: 150px"> </p>
into README.Rmd so that the .md file in GitHub will show the sticker.

Then, in the website, pkgdown builder adds a header with logo.svg for all of the pages rendered to HTML. In this way, index.html has two hex stickers.

I thought that maybe we could follow some sort of approach with conditional code depending on output format, and we'd want the .md to render with the above code, and the .html to render without.

Since .Rd is the first step in the .Rmd -> .md -> html pipeline, the knitr::opts_knit are identical for both .md and .html output, se we couldn't use this to differentiate.

If we care about this appearance issue, we would probably need to remove the code from the .Rmd and inject it into the .md file. Doing this by hand would be quite unsustainable.

Add the functionality to have separate covariances matrices per group / treatment arm

Note: In that case we need covariance structure formula part like this:
us(visit | group / subjid)

To do:

Clean up exports

To do:

decide which objects do not need to be exported
add appropriate badge experimental, stable, deprecated to all exported objects
_pkgdown.yaml should only contain exported objects
unexported objects must have documentation but they need @keywords internal tag and must not have examples (since those won't run)
- when removing examples from those, please make sure that the same thing is covered in existing unit tests, or otherwise add the example as a new unit test - so that we don't lose covering that case basically.
Use this script to easily detect what is left to do.

#' packages having keyword internal and matching specific type
rd_index_installed <- function(pkg, type = NULL) {
  db <- tools::Rd_db(pkg)
  elo <- tools:::.build_Rd_index(tools:::Rd_contents(db), type = type)
  elo$Name
}

man_files <- function(path) {
  list.files(file.path(path, "man"), full.names = TRUE, pattern = ".Rd")
}

rd_index <- function(path = ".") {
  all_docs <- man_files(path)
  res <- vapply(all_docs, FUN.VALUE = logical(1), FUN = function(x) {
    lines <- readLines(x)
    any(grepl("keyword\\{.*internal\\}", lines))
  })
 sort(gsub(".Rd$", "", basename(names(res[!res]))))
}

badge <- function(path = ".") {
  all_docs <- man_files(path)
  res <- vapply(all_docs, FUN.VALUE = logical(1), FUN = function(x) {
    lines <- readLines(x)
    any(grepl("figure\\{lifecycle", lines))
  })
  sort(gsub(".Rd$", "", basename(names(res[res]))))
}


path <- "."
pkg <- "mmrm"


# has a badge but it's internal
setdiff(badge(path), rd_index(path))

# is in documentation index but has no badge
setdiff(rd_index(path), badge(path))

# unexported functions in the doc index
setdiff(rd_index(path), sort(getNamespaceExports(pkg)))
setdiff(rd_index_installed(pkg), sort(getNamespaceExports(pkg)))

Try to use TMB directly to fit an MMRM and obtain Satterthwaite d.f.

Motivation: We still have trouble fitting a true MMRM (i.e. without residual variance, and obtaining correct Satterthwaite) with glmmTMB, see https://github.com/glmmTMB/glmmTMB/blob/satterthwaite_df/glmmTMB/vignettes/satterthwaite_unstructured_example2.Rmd

So one idea is to directly use TMB since then we have more freedom how to define the model. Plus TMB supports Hessian and Jacobian calculations (with latest 1.9.0 version) which should help us with the Satterthwaite d.f. calculations to avoid numDeriv::jacobian().

Material:

Wiki
- Tutorial
The comprehensive TMB documentation
- Matrix operations
Use TMB:::setupRStudio() to setup R Studio (once)

To do:

fit_single_optimizer <- function(formula,
                                 data,
                                 start = NULL,
                                 optimizer = c("L-BFGS-B", "BFGS", "CG", "nlminb"))

and then inside the function:

control <- glmmTMB::glmmTMBControl(
    optimizer = if (optimizer == "nlminb") stats::nlminb else stats::optim,
    optArgs = if (optimizer == "nlminb") list() else list(method = optimizer),
    parallel = 1L
  )

Add a getter function `component()` for `mmrm_tmb` objects

Idea: Instead of replicating code across the package to e.g. look up number of subjects etc. it is better to have a getter function that can be used e.g. as component(object, "n_subjects").

To do:

Add component(object, name) function in https://github.com/openpharma/mmrm/blob/main/R/tmb-methods.R (although it is not an S3 method makes sense here I guess)
- compare getS3method("getME", class = "glmmTMB") but no need to do it exactly like that
Look through the mmrm package and see where we are directly accessing mmrm_tmb contents.
- e.g. object$beta_est
Add corresponding name option to the component function
- e.g. name == "beta" and that will return object$beta_est
replace direct access code with component() calls

openpharma / mmrm Goto Github PK

mmrm's People

Contributors

Stargazers

Watchers

Forkers

mmrm's Issues

Recommend Projects

Recommend Topics

Recommend Org