Break between v0.4.62 and v0.5.1?

[BarkleyBG]

With inferference_0.5.1,

full_fit <-  try(inferference::interference(
    formula = my_formula,
    data = my_dataset,
    allocations = c(.3, .6),
    randomization = 2/3,
    ),
    runSilent = FALSE
))
[1] "Calculating matrix of IP weights..."
[1] "Calculating array of IP weight derivatives..."
[1] "Calculating matrix of scores..."
Error in grad.default(X = structure(list(`(Intercept)` = c(1, 1, 1, 1,  : 
  function returns NA at 7.75248344440477e-059.44252957302433e-063.03965082818232e-060.0001242369543700136.31106479280438e-050.00012825501761554 distance from x.

However, for inferference_0.4.62,

full_fit <-  try(inferference::interference(
    formula = my_formula,
    data = my_dataset,
    allocations = c(.3, .6),
    randomization = 2/3,
    ),
    runSilent = FALSE
))
[1] "Calculating matrix of IP weights..."
[1] "Calculating array of IP weight derivatives..."
[1] "Calculating matrix of scores..."
[1] "Computing effect estimates..."
[1] "Interference complete"

full_fit
 --------------------------------------------------------------------------
                               Model Summary                    
 --------------------------------------------------------------------------

I can provide you with the dataset.

[BarkleyBG]

Minimal working example from the inferference_intro.Rmd vignette:

Richardson method on the vaccinesim data runs fine

example1 <- interference(
    formula = Y | A | B ~ X1 + X2 + (1|group) | group, 
    allocations = c(.3, .45,  .6), 
    data = vaccinesim, 
    randomization = 2/3,
    # method = 'simple', ##comment this out
    runSilent = FALSE)

But, create a very large group, and the Richardson derivative method throws an error

vaccinesim2 <- vaccinesim
vaccinesim2$group[1:2500] <- 1

example1.2 <- interference(
  formula = Y | A | B ~ X1 + X2 + (1|group) | group, 
  allocations = c(.3, .45,  .6), 
  data = vaccinesim2,  ##large group
  randomization = 2/3,
  # method = 'simple', ## Richardson method throws error
  runSilent = FALSE)

Note that the error does not reproduce with the simple derivative method

example1.2 <- interference(
  formula = Y | A | B ~ X1 + X2 + (1|group) | group, 
  allocations = c(.3, .45,  .6), 
  data = vaccinesim2,  ##large group
  randomization = 2/3,
  method = 'simple', ## simple derivative does not throw error
  runSilent = FALSE)

I'm not sure if this is a dividing by NA error for large group propensity scores, or if maybe the Richardson method with large groups somehow attempts to evaluate a negative value for the random effect's variance parameter (which is by default positive).

Let me know how I can help

[bsaul]

Brian, The issue isn't with the propensity scores, but with evaluating the scores of the mixed effects model. The loglihood is -Inf for the large group in vaccinesim2:

mod <- lme4::glmer(B ~ X1 + X2 + (1|group), data = vaccinesim2, family = binomial(link = 'logit'))
theta_mod <- unlist(lme4::getME(mod, c('beta', 'theta')))
X1.2 <- model.matrix(B ~ X1 + X2, data = subset(vaccinesim2, group == 1))
B1.2 <- subset(vaccinesim2, group == 1)$B

# "by hand"
log(integrate(f = logit_integrand, lower = -Inf, upper = Inf, X=X1.2, A = B1.2, parameters = theta_mod)$value)

# Using inferference's internal log_likelihood function:
log_likelihood(
  parameters = theta_mod,
  integrand  = logit_integrand,
  X=X1.2, A = B1.2
)

Note that your example with method = 'simple' does not give an error, but all the standard errors are NaN:

example1.2 <- interference(
  formula = Y | A | B ~ X1 + X2 + (1|group) | group, 
  allocations = c(.3, .45,  .6), 
  data = vaccinesim2,  ##large group
  randomization = 2/3,
  method = 'simple', ## simple derivative does not throw error
  runSilent = FALSE)
example1.2

I don't know what to do in this situation. One thought, but it's rather hacky: if loglihood returns -Inf replace with `log( .Machine$double.eps). I don't particularly like that idea.

[BarkleyBG]

I dont think this was an issue in version 0.4.62. (Although I should double-check that). I wonder why not?

[bsaul]

When I run the 0.4.62 score functions, I get NAs for the results:

library(inferference)
vaccinesim2 <- vaccinesim
vaccinesim2$group[1:2500] <- 1
mod <- lme4::glmer(B ~ X1 + X2 + (1|group), data = vaccinesim2, family = binomial(link = 'logit'))
theta_mod <- unlist(lme4::getME(mod, c('beta', 'theta')))
X1.2 <- model.matrix(B ~ X1 + X2, data = subset(vaccinesim2, group == 1))
B1.2 <- subset(vaccinesim2, group == 1)$B

# "by hand"
log(integrate(f = logit_integrand, lower = -Inf, upper = Inf, X=X1.2, A = B1.2, parameters = theta_mod)$value)

# Using inferference's internal log_likelihood function:
score_calc(
  parameters = theta_mod,
  integrand  = logit_integrand,
  X=X1.2, A = B1.2
)

#################
## Compare to 0.4.62
##################
logit_integrand_0.4.62 <- function(b, X, A, 
                            fixed.effects,
                            random.effects = NULL,
                            x = NULL, 
                            pos = NULL, 
                            allocation = NULL, 
                            randomization = 1, 
                            integrate.allocation = FALSE)
{
  p  <- length(fixed.effects)
  re <- random.effects[1]

  ## In the case of an intercept-only model, X needs to be converted to matrix
  # for the warning to work
  if(!is.matrix(X)){
    X <- as.matrix(X)
  }

  ## Warnings ##
  if(p != ncol(X)){
    stop('The number of fixed effect parameters is not equal to the number \n
         of columns in the covariate matrix')
  }

  if(length(A) != nrow(X)){
    stop('Length of treatment vector is not equal to number of observations')
  }

  ## For taking derivative w.r.t. a parameter ##
  params <- c(fixed.effects, re)
  if(!is.null(pos)){
    params[pos] <- x
  }

  ## Calculations ## 
  if(is.null(re) || re <= 0){
    pr.b <- randomization * (plogis(X %*% params[1:p]))
  } else {
    pr.b <- randomization * (plogis(drop(outer(X %*% params[1:p], b, '+'))))
  }
  if(integrate.allocation == FALSE){
    hh <- dbinom(A, 1, pr.b)
  } else {
    hh <- (pr.b/allocation)^A * ((1-pr.b)/(1 - allocation))^(1-A)
  }
  if(is.null(re) || re <= 0){
    # in this way dnorm integrates to one when integrating from -Inf to Inf
    out <- prod(hh) * dnorm(b, mean=0, sd = 1) 
  } else {
    hha <- apply(hh, 2, prod)
    out <- hha * dnorm(b, mean=0, sd = params[p + 1])
  }

  return(out)
  }

log_likelihood_0.4.62 <- function(x, 
                           pos, 
                           integrand,
                           ...)
{
  ## Necessary pieces ##
  integrand <- match.fun(integrand)
  dots      <- list(...)
  dot.names <- names(dots)

  ## Integrate() arguments ##
  if(!'lower' %in% dot.names){
    dots$lower <- -Inf
  }

  if(!'upper' %in% dot.names){
    dots$upper <- Inf
  }

  int.args <- append(get_args(integrate, dots),
                     get_args(integrand, dots))
  args <- append(int.args, list(f = integrand, x = x, pos = pos))

  ## Calculuation ##
  attempt <- try(do.call(integrate, args = args))
  val <- ifelse(is(attempt, 'try-error'), NA, attempt$value)

  return(log(val))
}

score_calc_0.4.62 <- function(integrand,
                       hide.errors = TRUE,
                       fixed.effects,
                       random.effects,
                       ...)
{
  ## Necessary bits ##
  params <- c(fixed.effects, random.effects)
  integrand <- match.fun(integrand)
  dots <- list(...)

  ## Function arguments ##
  int.args <- append(get_args(integrand, dots),
                     get_args(integrate, dots))
  fargs    <- append(int.args, get_args(numDeriv::grad, dots))

  ## Compute the derivative of the log likelihood for each parameter ##
  scores <- sapply(1:length(params), function(i){
    args <- append(fargs,
                   list(func = log_likelihood_0.4.62,
                        integrand = integrand, 
                        fixed.effects = fixed.effects,
                        random.effects = random.effects,
                        x = params[i], 
                        pos = i))

    attempt <- try(do.call(numDeriv::grad, args = args), silent = hide.errors)
    return(ifelse(is(attempt, 'try-error'), NA, attempt))
  })

  return(scores)
}

score_calc_0.4.62(
  integrand  = logit_integrand_0.4.62,
  fixed.effects = lme4::getME(mod, 'beta'),
  random.effects = lme4::getME(mod, 'theta'),
  X=X1.2, A = B1.2
)

[bsaul]

Closing this issue. If there's still a problem, please reopen.