Possible improvement to D1 function

[fwcavalcante]

I would like to propose a small change to the D1 function in mice.
In specific situations it can produce wrong results that can go undetected:

Here is a reproducible example explaining the source of the issue:

library(mice)

data("airquality")
summary(airquality)

#data = airquality[!is.na(airquality$Ozone), ]
data = airquality
fit1 = glm(data = data, Ozone ~ Wind * Temp) 
fit0 = glm(data = data, Ozone ~ Wind + Temp) 

# LRT test p-value interaction in the unimputed data
anova(fit1, fit0) # PR(>F) = 6.569e-05 ***  Degrees of Freedom Df = 112, 113 

# now we impute (naive imputation)
pred  = quickpred(data, minpuc = 0.5, include = c('Ozone', 'Wind', 'Temp'))
mimp = mice(data, predictorMatrix = pred, m = 5, seed = 12345, printFlag = F)

fit1imp = with(mimp, glm(Ozone ~ Wind * Temp))
fit0imp = with(mimp, glm(Ozone ~ Wind + Temp))

# If we use the imputed dataset and D1 we get a much  higher p-value
mice::D1(fit1imp, fit0imp) # PR(>F) = 0.0127296 /// F-score = 18.4, Degrees of Freedom Df(!) = 1, 4 , dfcom=149

#changing the dfcom does not affect the results (it will be clear later why)
mice::D1(fit1imp, fit0imp, dfcom = NULL)
mice::D1(fit1imp, fit0imp, dfcom = 149)

# if we use D2 we get results much closer to the unimputed data
mice::D2(fit1imp, fit0imp). #p-value = 3.137262e-05

if we look at each imputed dataset individually, each of the 5 imputation sets produce similar p-values to the original data...so the pooled p-value with D1 doesnt look right in this case

  for(x in 1:mimp$m){
    temp = complete(mimp, action = x)
    t1 = with(temp, glm(Ozone ~ Wind * Temp))
    t0 = with(temp, glm(Ozone ~ Wind + Temp))
    res = anova(t1, t0)
    print(res$`Pr(>F)`[2])
  
  }

the issue appears to be the degrees of freedom df1 and df2 getting much lower (1 and 4) than they should be. this is caused by the mitmt::d1() called by mice to compute D1. MITML::D1 uses the Reiter Vf Degree of freedom correction for small samples by default when a dfcom is passed, and mice::d1 always passes a dfcom to mitml....(even if we set dfcom to null)
There is nothing wrong with using the Reiter formula, except that if we have:
-5 imputed sets and
-only 1 extra parameter between the nested models (interaction between two continuous, one binary and one continuous and two binary variables) the Reiter formula cant be used (results in divisions by zero)...

this happens whenever t =4 , (where t = k(m-1) , k = number of new parameters in fit1, m = number of imputations)
if t = 4 the c0 parameter in the Reiter formula will include divisions by zero (c0 = 1 / (t-4)). c0 is then used in the Reiter formula, producing wrong 'corrected degrees of freedom'.*

currently, there is no way to prevent Reiter formula from being used using mice::D1. I would like to suggest a small change to D1 conde to include a warning (either advsing to increase the number of imputations or giving the options to skip the Reiter correction).

using mitml itself it is possible to prevent the Reiter vf formula from being used:

mitml1 = as.mitml.result(fit1imp)
mitml0 = as.mitml.result(fit0imp)
mitml::testModels(mitml1, mitml0, method = 'D1') ### here we get results that appear more consistent with the original data 

mitml::testModels(mitml1, mitml0, method = 'D1', df.com = 149) # here we get the much higher p valuees and incorrect degree of freedom (df2 changed from 225 to 4)

we can confirm this is the case by repeating the analysis but using 6 instead of 5 imputations

mimp6 = mice(data, predictorMatrix = pred, m = 6, seed = 12345, printFlag = F)
fit1imp6 = with(mimp6, glm(Ozone ~ Wind * Temp))
fit0imp6 = with(mimp6, glm(Ozone ~ Wind + Temp))

mice::D1(fit1imp6, fit0imp6) # now the corrected reiter VF df2 jumps to a more reasonable 68 instead of 4 and the results appear consistent with the other analyses
mice::D1(fit1imp6, fit0imp6, dfcom = NULL)

ROOT CAUSE

# the way mice D1 operates ALWAYS passess a df.com to mitml: https://rdrr.io/cran/mice/src/R/D1.R
# and the way mitml testmodel() and internal .D1() function operates to ALWAYS use the Reiter vf degrees of freedom if the df.com is not NULL: https://github.com/cran/mitml/blob/master/R/testModels.R
# below i copied parts of mitml source code and anotated were the issue happens:


# MITML D1 function (and auxiliary functions),


#mitml test model will use the model and null model and dfcom (if given) to compute D1
model = as.mitml.result(fit1imp6)
null.model =  as.mitml.result(fit0imp6)


# first it detects which terms are different between the models ('par.diff'), in our case it will only be the interaction termm, it uses the extract paramter function to do that, i included it here just so you can run the code directly if you want
.extractParameters <- function(model, diagonal = FALSE, include.extra.pars = FALSE){
  
  # number of imputations
  m <- length(model)
  
  # extract parameter estimates and variance-covariance matrices
  Qhat <- lapply(model, coef, include.extra.pars = include.extra.pars) # i made a small change here just to avoid uncessary lengthy code
  Uhat <- lapply(model, coef, include.extra.pars = include.extra.pars) # i made a small change here just to avoid uncessary lengthy code
  p <- length(Qhat[[1]])
  nms <- names(Qhat[[1]])
  
  # preserve parameter labels (if any)
  attr(nms, "par.labels") <- attr(Qhat[[1]], "par.labels")
  
  # ensure proper dimensions
  stopifnot(all(p == dim(Uhat[[1]])))
  Qhat <- matrix(unlist(Qhat), nrow = p, ncol = m)
  Uhat <- array(unlist(Uhat), dim = c(p, p, m))
  
  # extract diagonal
  if(diagonal){
    Uhat <- apply(Uhat, 3, diag)
    if(is.null(dim(Uhat))) dim(Uhat) <- dim(Qhat)
  }
  
  out <- list(Qhat = Qhat, Uhat = Uhat, nms = nms)
  return(out)
  
}

# here mitml will get the parameters from fit 1 and fit -
est <- .extractParameters(model, diagonal = FALSE)
est.null <- .extractParameters(null.model, diagonal = FALSE)

# and identify the number of 'extra' parameters (this will be K for the formula below)
par.diff <- est$nms[!(est$nms %in% est.null$nms)]
par.ind <- match(par.diff, est$nms)
if(length(par.diff) == 0L) stop("The 'model' and 'null.model' appear not to be nested or include the same set of parameters.")

k <- length(par.diff)

# now mitmt uses the estimates and variances of the new terms (just the interaction terms in our example) to compute D1
Qhat <- est$Qhat[par.ind,, drop = FALSE]
Uhat <- est$Uhat[par.ind, par.ind,, drop = FALSE]


# here is were the issue happens:
# Please check annotations inside the function, they explain the issue....
.D1 <- function(Qhat, Uhat, df.com){
  # pooling for multidimensional estimands (D1, Li et al., 1991; Reiter, 2007)
  
  k <- dim(Qhat)[1]
  m <- dim(Qhat)[2]
  
  # D1
  Qbar <- apply(Qhat, 1, mean)
  Ubar <- apply(Uhat, c(1, 2), mean)
  
  B <- cov(t(Qhat))
  r <- (1+m^(-1))*sum(diag(B%*%solve(Ubar)))/k
  Ttilde <- (1 + r)*Ubar
  
  val <- t(Qbar) %*% solve(Ttilde) %*% Qbar / k
  
  # compute degrees of freedom (df2)
  t <- k*(m-1)                  #(!!!) <------ IF WE ONLY HAVE 1 PARAMETER DIFFERENT BETWEEN THE MODELS (ONE EXTRA INTERACTION TERM) AND M == 5, THIS VALUE WILL BE EQUAL TO 4
  if(!is.null(df.com)){         #(!!!) <------ MICE WILL ALWAYS FOLLOW THIS LINE AS DF.COM WILL ALWAYS BE PASSED (EVEN IF DF.COM IS SET TO NULL IN MICE::D1() CALL)
    
    # small-sample degrees of freedom (Reiter, 2007; Eq. 1-2)
    a <- r*t/(t-2)
    vstar <- ( (df.com+1) / (df.com+3) ) * df.com
    
    c0 <- 1 / (t-4)            #(!!!) <------- AS T == 4,  C0 WILL BE 1 DIVIDED BY ZERO...THAT IS WHAT MAKES THE 'REITER CORRECTED DEGREE OF FREEDOM' WRONG
    c1 <- vstar - 2 * (1+a)
    c2 <- vstar - 4 * (1+a)
    
    z <- 1 / c2 +
      c0 * (a^2 * c1 / ((1+a)^2 * c2)) +
      c0 * (8*a^2 * c1 / ((1+a) * c2^2) + 4*a^2 / ((1+a) * c2)) +
      c0 * (4*a^2 / (c2 * c1) + 16*a^2 * c1 / c2^3) +
      c0 * (8*a^2 / c2^2)
    
    v <- 4 + 1/z
    
  }else{
    
    if (t > 4){
      v <- 4 + (t-4) * (1 + (1 - 2*t^(-1)) * (r^(-1)))^2
    }else{
      v <- t * (1 + k^(-1)) * ((1 + r^(-1))^2) / 2
    }
    
  }
  
  return(list(F = val, k = k, v = v, r = r))
  
}

Suggestion

I think the most logical ways to avoid this issue would be to print a warning whenever mice::d1 is called with a 't' value equal to 4...
The warning can instruct users to either increase the number of imputations, if they want to procced Reiter correction cannot be used... (then D1 code has to change to no pass df.com to mitml::testmodel()

it would be easy to calculate T within the mice::D1 call and adjust the following line in the mice::D1() call to only pass dfcom if t>4:
tmr <- mitml::testModels(fit1, fit0, method = "D1", df.com = dfcom)

final note

as it is probably clear by now, i am not in any way a statistician, possibly I made some mistakes in my explanation and if something is really stupid please accept my apologies and feel free to correct me

[stefvanbuuren]

Thanks a lot for your deep analysis and suggested action.

It is both surprising and odd that mice reports df2 = 4. By default, mice calculates dfcom from the residuals of the more complicated model, and send this calculated dfcom down to mitml::testModels() for further processing. For the interaction model fitted to the airquality data, dfcom = nrow(airquality) - 4 = 149.

Your suggestion to bypass dfcom calculation would optically solve the problem. However, note that the solution produces a df2 estimate of 225.6, a value higher than nrow(airquality). The point of the Reiter (2007) paper was to derive df2 such that it would be lower than dfcom. For cases where this does not hold, Reiter wrote:

This can result in a larger proportion of p-values below desired significance levels than would be expected under the null hypothesis for a test with valid frequentist properties.

In the case that t = k * (m - 1) = 4 and df.com is not specified, mitml::testModels() makes a large-sample assumption and reports this by the message Unadjusted hypothesis test as appropriate in larger samples. This is a reasonable solution if we do not know df.com.

However, if we do know df.com - and mice always calculates it - a better solution would be to limit df2 to df.com. There is still a downward bias in the $P$-value, but it is smaller than assuming infinite df.com.

A slighly better approach is to temporarily set tstar <- 5 if t == 4, and use tstar instead of t in the block that follows if(!is.null(df.com)). This mimmicks your m = 6 solution and will yield a $P$-value very close to the appropriate one.

In cases where t < 4, it is probably wisest to report NaN(which is what testModels() already does), so that the user sees the need to increase the number of imputations m.

Such changes might best go into mitml, so I include @simongrund1 here.

[simongrund1]

Thanks for tagging. The t==4 issue with Reiter's df was recently brought up to me via email as well, and testModels indeed does not handle some edge cases very well. It will be addressed in the next dev version of mitml. Meanwhile, the t==4 case is one with very few imputations (fewer than one probably should use).

[stefvanbuuren]

Great. Thanks for addressing. Agree that m = 5 maybe too low for this application.

[fwcavalcante]

Dear Simon, and Stef,

many thanks for your response.
indeed i contacted Simon about it via email before.

I dont have the background to comment on the best solution but i agree the best would be to warn the user to use more imputations (which is what we decided to do in our project)

many thanks!
Fabiano