R mice code for delta sensitivity analysis for a binary variable : apparently not working

Question

I am doing multiple imputation with R mice before my primary analysis which is a binary logistic regression with functional decline as outcome and a number of variables as predictors (continuous and binary). I am using pmm for continuous variables and logreg for binary variables.

I need to perform a sensitivity analysis with a delta-adjustment for the outcome functional decline assuming that the missing are more likely to be decliners by an amount that I have calculated to be equal to log=4.3 according to Rezvan et al.

I'm using this code, from a paper by Leacy et al. to customize the mice.impute.logreg function:

library(mice)
mice.impute.logreg.sens <- function (y, ry, x, delta, ...) {
  x <- cbind(1, as.matrix(x))
  expr <- expression(glm.fit(x[ry, ], y[ry], family = binomial(link = logit)))
  fit <- suppressWarnings(eval(expr))
  fit.sum <- summary.glm(fit)
  beta <- coef(fit)
  beta[1] <- beta[1] + delta
  rv <- t(chol(fit.sum$cov.unscaled))
  beta.star <- beta + rv % % rnorm(ncol(rv))
  p <- 1 / (1 + exp(-(x[!ry, ] % % beta.star)))
  vec <- (runif(nrow(p)) <= p)
  vec[vec] <- 1
  if (is.factor(y)) {
    vec <- factor(vec, c(0, 1), levels(y))
  }
  return (vec)
}

and I apply it to my dataset like so:

myvars <- c("FunctionalDecline", "Age", "Sex", "GDS4")
realdataset <- dataset[myvars]
realdataset$FunctionalDecline <- factor(realdataset$FunctionalDecline)
realdataset$Sex <- factor(realdataset$Sex)
ini <- mice(realdataset, maxit = 0)
meth <- ini$meth
meth["FunctionalDecline"] = "logreg.sens"
imputeddataset <- mice(realdataset, meth = meth, seed = 3, m = 5, maxit = 10, delta = log(4.3), print = F)
summary(pool(with(imputeddataset, glm(FunctionalDecline~Age + Sex + GDS4, family = binomial))))

When I run the code, I think the results are wrong, because:

1. With delta=log(1.0) (i.e. no adjustment) the results of the logistic regression are different than those of the logistic regression on the dataset with no delta-adjustment

2. After the delta-adjustment with delta=log(4.3) all predictors in the binary logistic regression lose significance. If, however, I do a "rougher" worst-scenario sensitivity analysis by putting all imputed functional decline values equal to 1 (i.e. decline; post-imputation or within the imputation itself using the mice post function) most predictors retain significance.

3. If, to approximate the worst scenario analysis, I assume 99% of decliners among the missing and use the corresponding delta=log(280) in the code, the results are totally inconsistent with the "rougher" worst-scenario analysis described in 2.

If I increase the number of iterations and of imputations or change the visiting sequence, the results don't change substantially.

Why does the code give the wrong results?

Answer 1

I have done further researching and I think I now know “where” but still not “why” the code does not work.

The problem could be in the way the mice.impute.logreg function is customized in the 2017 paper.

In fact, looking at the lates documentation on the on the mice.impute.logreg function here (https://github.com/amices/mice/blob/master/R/mice.impute.logreg.R) the code is different and it could well be that the function is not working properly with the “old” code.

I therefore customized the function using the same offset delta for the linear predictor values as done in the 2017 paper.

This is the new customized code:

    library(mice)

    mice.impute.logreg.sens <- function(y, ry, x, wy = NULL, delta, ...) {
        if (is.null(wy)) wy <- !ry
        aug <- augment(y, ry, x, wy)
        x <- aug$x
        y <- aug$y
        ry <- aug$ry
        wy <- aug$wy
        w <- aug$w

        x <- cbind(1, as.matrix(x))
        expr <- expression(glm.fit(
            x = x[ry, , drop = FALSE],
            y = y[ry],
            family = quasibinomial(link = logit),
            weights = w[ry]
        ))
        fit <- eval(expr)
        fit.sum <- summary.glm(fit)
        beta <- coef(fit)
        beta[1] <- beta[1] + delta
        rv <- t(chol(sym(fit.sum$cov.unscaled)))
        beta.star <- beta + rv %*% rnorm(ncol(rv))

        p <- 1 / (1 + exp(-(x[wy, , drop = FALSE] %*% beta.star)))
        vec <- (runif(nrow(p)) <= p)
        vec[vec] <- 1
        if (is.factor(y)) {
            vec <- factor(vec, c(0, 1), levels(y))
        }
        vec
    }

I then applied it to my dataset:

    myvars <- c("FunctionalDecline", "Age", "Sex", "GDS4")
    realdataset <- dataset[myvars]
    realdataset$Functionaldecline <- factor(realdataset$Functionaldecline)
    realdataset$Sex <- factor(realdataset$Sex)
    ini <- mice(realdataset, maxit = 0)
    meth <- ini$meth
    meth["FunctionalDecline"] = "logreg.sens"
    imputeddataset <- mice(realdataset, meth = meth, seed = 3, m = 5, maxit = 10, delta = log(4.3), print = F)
    summary(pool(with(imputeddataset, glm(FunctionalDecline ~ Age + Sex + GDS4, family = binomial))))

However, the results I get are the same as if I run the regression on the imputed dataset without the delta-adjustment.

So it all boils done to the fact that the “old” code gives the wrong results while mice does not “recognize” the new customized impute.logreg.sens function.

I think this is one step further towards the solution, if someone else can take from it here…