Fix Standard Error Estimation With Imputed Data

[jackmwolf]

Resolves #30 by estimating the standard errors of the survival and cumulative incidence functions using Rubin's rules.

Reprex showing adjustedcif and adjustedsurv running with updated variance estimation on imputed data:

# Testing updated variance imputation ==========================================
library(adjustedCurves)
library(mice)
#> 
#> Attaching package: 'mice'
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
library(WeightIt)

set.seed(1)
# adjustedsurv() ===============================================================
sim_dat <- sim_confounded_surv(n=50, max_t=1.2)
sim_dat$group <- as.factor(sim_dat$group)

# introduce random missingness in x1 as example
sim_dat$x1 <- ifelse(runif(n=50) < 0.5, sim_dat$x1, NA)

# perform multiple imputation
mids <- mice::mice(data=sim_dat, method="pmm", m=10, printFlag=FALSE)

# IPTW KM using WeightIt on imputed data
adj <- adjustedsurv(data=mids,
                    variable="group",
                    ev_time="time",
                    event="event",
                    method="iptw_km",
                    treatment_model=group ~ x1 + x2 + x5 + x6,
                    weight_method="ps",
                    conf_int = TRUE,
                    force_bounds = TRUE
                    )
head(adj$adj)
#>    group       time      surv       var_w        var_b       var_t         se
#> 1      0 0.00000000 1.0000000 0.000000000 0.000000e+00 0.000000000 0.00000000
#> 3      0 0.03525682 1.0000000 0.000000000 0.000000e+00 0.000000000 0.00000000
#> 5      0 0.05324323 0.9690449 0.001958898 8.799185e-06 0.001968577 0.04436865
#> 7      0 0.08104409 0.9307052 0.004393613 5.481180e-05 0.004453906 0.06673759
#> 9      0 0.09757016 0.8734748 0.007727660 2.056206e-04 0.007953842 0.08918432
#> 11     0 0.14808191 0.8734748 0.007727660 2.056206e-04 0.007953842 0.08918432
#>     ci_lower ci_upper
#> 1  1.0000000 1.000000
#> 3  1.0000000 1.000000
#> 5  0.8820839 1.056006
#> 7  0.7999019 1.061508
#> 9  0.6986767 1.048273
#> 11 0.6986767 1.048273

# adjustedcif() ================================================================
# simulate some data as example
sim_dat <- sim_confounded_crisk(n=100, max_t=1.2)
sim_dat$group <- as.factor(sim_dat$group)

# introduce random missingness in x1 as example
sim_dat$x1 <- ifelse(runif(n=50) < 0.5, sim_dat$x1, NA)

# perform multiple imputation
mids <- mice::mice(data=sim_dat, method="pmm", m=10, printFlag=FALSE)

# IPTW Pseudo using WeightIt on imputed data, for cause = 1
adj <- adjustedcif(data=mids,
                   variable="group",
                   ev_time="time",
                   event="event",
                   method="iptw_pseudo",
                   cause=1,
                   treatment_model=group ~ x1 + x2 + x5 + x6,
                   weight_method="ps",
                   conf_int = TRUE,
                   force_bounds = TRUE
                   )
#> Loading required namespace: prodlim
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
#> Warning in sqrt(surv_sd): NaNs produced
head(adj$adj, 10)
#>    group        time        cif        var_w        var_b        var_t
#> 1      0 0.000000000 0.00000000          NaN 6.731613e-29          NaN
#> 3      0 0.001977503 0.00000000 0.0000000000 0.000000e+00 0.0000000000
#> 5      0 0.002027630 0.00000000 0.0000000000 0.000000e+00 0.0000000000
#> 7      0 0.002422720 0.00000000 0.0000000000 0.000000e+00 0.0000000000
#> 9      0 0.005725169 0.00000000 0.0000000000 0.000000e+00 0.0000000000
#> 11     0 0.011687012 0.00000000 0.0000000000 0.000000e+00 0.0000000000
#> 13     0 0.011867946 0.01206984 0.0001489893 1.748739e-07 0.0001491817
#> 2      0 0.015892799 0.02663836 0.0003652194 1.459200e-06 0.0003668245
#> 4      0 0.031244015 0.02663836 0.0003652194 1.459200e-06 0.0003668245
#> 6      0 0.033138364 0.02663836 0.0003652194 1.459200e-06 0.0003668245
#>            se    ci_lower   ci_upper
#> 1         NaN         NaN        NaN
#> 3  0.00000000  0.00000000 0.00000000
#> 5  0.00000000  0.00000000 0.00000000
#> 7  0.00000000  0.00000000 0.00000000
#> 9  0.00000000  0.00000000 0.00000000
#> 11 0.00000000  0.00000000 0.00000000
#> 13 0.01221399 -0.01186915 0.03600883
#> 2  0.01915266 -0.01090017 0.06417689
#> 4  0.01915266 -0.01090017 0.06417689
#> 6  0.01915266 -0.01090017 0.06417689

^{Created on 2024-07-28 with reprex v2.1.1}

Testthat results:

══ Results ══════════════════════════════════════════════════════════════════════════════════════════════════════
Duration: 236.0 s

[ FAIL 0 | WARN 4 | SKIP 0 | PASS 1862 ]

All warnings do not seem to be related to this PR:

✔ |   4     43 | cif_iptw_pseudo [1.3s]                                                                          
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Warning (test_cif_iptw_pseudo.r:29:3): 2 treatments, with conf_int, no boot
NaNs produced
Backtrace:
    ▆
 1. └─adjustedCurves::adjustedcif(...) at test_cif_iptw_pseudo.r:29:3
 2.   ├─R.utils::doCall(cif_fun, args = args, .ignoreUnusedArgs = FALSE) at adjustedCurves/R/adjustedcif.r:353:5
 3.   ├─R.utils::doCall.default(cif_fun, args = args, .ignoreUnusedArgs = FALSE)
 4.   │ └─base::do.call(.fcn, args = args, envir = envir)
 5.   └─adjustedCurves (local) `<fn>`(...)
 6.     └─adjustedCurves:::method_iptw_pseudo(...) at adjustedCurves/R/method_iptw_pseudo.r:127:3

Warning (test_cif_iptw_pseudo.r:59:3): 2 treatments, with conf_int, with boot
NaNs produced
Backtrace:
    ▆
 1. └─adjustedCurves::adjustedcif(...) at test_cif_iptw_pseudo.r:59:3
 2.   ├─R.utils::doCall(cif_fun, args = args, .ignoreUnusedArgs = FALSE) at adjustedCurves/R/adjustedcif.r:353:5
 3.   ├─R.utils::doCall.default(cif_fun, args = args, .ignoreUnusedArgs = FALSE)
 4.   │ └─base::do.call(.fcn, args = args, envir = envir)
 5.   └─adjustedCurves (local) `<fn>`(...)
 6.     └─adjustedCurves:::method_iptw_pseudo(...) at adjustedCurves/R/method_iptw_pseudo.r:127:3

Warning (test_cif_iptw_pseudo.r:111:3): 3 ways of iptw calculation are equal
NaNs produced
Backtrace:
    ▆
 1. └─adjustedCurves::adjustedcif(...) at test_cif_iptw_pseudo.r:111:3
 2.   ├─R.utils::doCall(cif_fun, args = args, .ignoreUnusedArgs = FALSE) at adjustedCurves/R/adjustedcif.r:353:5
 3.   ├─R.utils::doCall.default(cif_fun, args = args, .ignoreUnusedArgs = FALSE)
 4.   │ └─base::do.call(.fcn, args = args, envir = envir)
 5.   └─adjustedCurves (local) `<fn>`(...)
 6.     └─adjustedCurves:::method_iptw_pseudo(...) at adjustedCurves/R/method_iptw_pseudo.r:127:3

Warning (test_cif_iptw_pseudo.r:121:3): 3 ways of iptw calculation are equal
NaNs produced
Backtrace:
    ▆
 1. └─adjustedCurves::adjustedcif(...) at test_cif_iptw_pseudo.r:121:3
 2.   ├─R.utils::doCall(cif_fun, args = args, .ignoreUnusedArgs = FALSE) at adjustedCurves/R/adjustedcif.r:353:5
 3.   ├─R.utils::doCall.default(cif_fun, args = args, .ignoreUnusedArgs = FALSE)
 4.   │ └─base::do.call(.fcn, args = args, envir = envir)
 5.   └─adjustedCurves (local) `<fn>`(...)
 6.     └─adjustedCurves:::method_iptw_pseudo(...) at adjustedCurves/R/method_iptw_pseudo.r:127:3
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────

[RobinDenz1]

Thank you very much for finding and fixing (!) this error. Will push an update to CRAN soon.