Break between v0.4.62 and v0.5.1?
Closed this issue · 5 comments
BarkleyBG commented
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 commented
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 commented
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 commented
I dont think this was an issue in version 0.4.62. (Although I should double-check that). I wonder why not?
bsaul commented
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 commented
Closing this issue. If there's still a problem, please reopen.