MCSE missing for relative bias when any standard errors are NA

[mikesweeting]

Hi Alessandro,

Thanks again for the rsimsum package! I find it very useful.

I think I've discovered a small bug when calculating the relative bias performance measure. The issue arises if any of the standard errors are missing (e.g. if the model returned a point estimate but no standard error).

The relative bias is calculated only on point estimate from simulations where standard errors are non-missing. I think this is intended behaviour of rsimsum? However, the MCSE is not returned for relative bias if there are any missing standard errors, whereas the MCSE is returned if calculating bias.

If the intended behaviour of rsimsum is to exclude any simulations with missing standard errors when calculating bias / relative bias then it would be useful that this is clearly documented. I tend to agree that this behaviour is sensible as a missing SE indicates possible non-convergence of the model.

Here is a worked example. I've compared output to that in the mcstatsim package:

# Worked example

library(rsimsum)
library(mcstatsim)

set.seed(135)
test = data.frame(
  dataset_id = 1:5,
  method = 1, 
  estimate = rnorm(5, mean = 1),
  se = rnorm(5, 1/sqrt(5), 0.05),
  true = 1
)
test
#>   dataset_id method   estimate        se true
#> 1          1      1  0.5549292 0.4594740    1
#> 2          2      1  0.5340348 0.5074508    1
#> 3          3      1  0.8955369 0.4251786    1
#> 4          4      1  2.4083353 0.4683549    1
#> 5          5      1 -0.3160867 0.4954568    1

ss = simsum(
  data = test, estvarname = "estimate", se = "se", true = "true",
  methodvar = "method", x = TRUE, ref = "1")
# Relative Bias
summary(ss, stat = "rbias")
#> Values are:
#>  Point Estimate (Monte Carlo Standard Error)
#> 
#> Relative bias in point estimate:
#>                 1
#>  -0.1847 (0.4456)
# check with mcstatsim
# Note relative bias here is 1+relative bias reported by rsimsum
# rel bias and MCSE agree
crb = calc_relative_bias(test$estimate, true_param = unique(test$true))
crb$rel_bias-1
#> [1] -0.1846501
crb$rel_bias_mcse
#> [1] 0.4455973


# Now lets set one of the SEs to be NA
test2 = test
test2[1,"se"] = NA
test2
#>   dataset_id method   estimate        se true
#> 1          1      1  0.5549292        NA    1
#> 2          2      1  0.5340348 0.5074508    1
#> 3          3      1  0.8955369 0.4251786    1
#> 4          4      1  2.4083353 0.4683549    1
#> 5          5      1 -0.3160867 0.4954568    1

ss2 = simsum(
  data = test2, estvarname = "estimate", se = "se", true = "true",
  methodvar = "method", x = TRUE, ref = "1")
# Relative Bias is returned, but MCSE is NA
summary(ss2, stat = "rbias")
#> Values are:
#>  Point Estimate (Monte Carlo Standard Error)
#> 
#> Relative bias in point estimate:
#>             1
#>  -0.1195 (NA)
# check with mcstatsim
# Note relative bias here is 1+relative bias reported by rsimsum
# rel bias does not agree with rsimsum and MCSE is returned
crb2 = calc_relative_bias(test2$estimate, true_param = unique(test2$true))
crb2$rel_bias-1
#> [1] -0.1846501
crb2$rel_bias_mcse
#> [1] 0.4455973

# If we remove simulation 1 then we get back the relative bias from rsimsum, plus the MCSE
crb3 = calc_relative_bias(test2$estimate[-1], true_param = unique(test2$true))
crb3$rel_bias-1
#> [1] -0.1195449
crb3$rel_bias_mcse
#> [1] 0.5690903

# Bias is also calculated only on the indivduals with non-missing SEs
# But unlike relative bias, the MCSE is returned
summary(ss2, stat = "bias")
#> Values are:
#>  Point Estimate (Monte Carlo Standard Error)
#> 
#> Bias in point estimate:
#>                 1
#>  -0.1195 (0.5691)
calc_bias(test2$estimate[-1], true_param = unique(test2$true))
#> $bias
#> [1] -0.1195449
#> 
#> $bias_mcse
#> [1] 0.5690903

System information:
I have rsimsum_0.13.0 installed.

Many thanks
Mike

[ellessenne]

Hi @mikesweeting, thanks for reporting this!

This is intended behaviour, incomplete (point estimate, SE) pairs are set to NA here:

rsimsum/R/simsum.R

Line 196 in 7ac2502

data <- .na_pair(data = data, estvarname = estvarname, se = se)

Whether that is sensible, that's a different story: I agree with you that incomplete pairs may be a symptom of non-convergence, but SEs are not needed for calculating relative bias (nor its MCSE), so it feels a bit of a strong design choice to require a SE (I may have previously overlooked this aspect). I need to keep thinking about this, any other thoughts on your end?

(also tagging @tpmorris – any thoughts?)

The MCSE should however not be NA, and in fact, there was a bug there. The NA is propagated by forgetting to pass na.rm = control$na.rm when calling sd here:

rsimsum/R/performance.R

Line 95 in 7ac2502

rbias_mcse <- sd(rbias_i) / sqrt(nsim)

Adding that (na.rm) fixes it, and yields numbers that are the same as {mcstatsim}:

# Worked example
library(rsimsum)
library(mcstatsim)
set.seed(135)
test <- data.frame(
  dataset_id = 1:5,
  method = 1,
  estimate = rnorm(5, mean = 1),
  se = rnorm(5, 1 / sqrt(5), 0.05),
  true = 1
)
simsum(
  data = test,
  estvarname = "estimate",
  se = "se",
  true = "true",
  methodvar = "method",
  x = TRUE,
  ref = "1"
) |>
  summary(stat = "rbias")
#> Values are:
#>  Point Estimate (Monte Carlo Standard Error)
#> 
#> Relative bias in point estimate:
#>                 1
#>  -0.1847 (0.4456)
crb <- calc_relative_bias(test$estimate, true_param = unique(test$true))
crb$rel_bias - 1
#> [1] -0.1846501
crb$rel_bias_mcse
#> [1] 0.4455973

test2 <- test
test2[1, "se"] <- NA
simsum(
  data = test2,
  estvarname = "estimate",
  se = "se",
  true = "true",
  methodvar = "method",
  x = TRUE,
  ref = "1"
) |>
  summary(stat = "rbias")
#> Values are:
#>  Point Estimate (Monte Carlo Standard Error)
#> 
#> Relative bias in point estimate:
#>                 1
#>  -0.1195 (0.5691)
crb2 <- calc_relative_bias(
  test2$estimate[-1],
  true_param = unique(test2$true)
)
crb2$rel_bias - 1
#> [1] -0.1195449
crb2$rel_bias_mcse
#> [1] 0.5690903

^{Created on 2026-04-07 with reprex v2.1.1}

I will fix this and submit a patch release to CRAN as soon as possible (since #61 needs to also be fixed urgently to stay on CRAN).

Many thanks!

Alessandro

[ellessenne]

Update: I fixed the missing na.rm in {rsimsum} 0.13.1, now on CRAN.
I'll leave this issue open, as I think it warrants further discussions!