Suppose we have a generalized linear model with a binomial response $y_i\sim \mathrm{bin}(n_i,p_i)$ where $p_i$ is determined by the linear predictor in the usual way via some link function. Is there some way to do multiple imputation when $n_i$ is sometimes missing (MCAR or MAR) (in addition to sometimes some of the ordinary covariates)? Obviously, we must have imputed values $n_i\ge y_i$ and it seems to me that standard methods of multiple imputation available in the mi and mice R-packages don't handle this case.
Some context: $n_i$ is the number of eggs laid and $y_i$ is the number of surviving fledglings in different nests $i=1,2,\dots,N$.
EDIT: Rather than working with $n_i$ and $y_i$ (number of Bernoulli trials and number of successes) one idea would be to work with the number of failures and successes $y_{0i}$ and $y_{1i}$ both of which are non-negative integers. $y_{0i}$ can then perhaps be imputed using predictive mean matching (recommended by van Buuren 2012, p. 78) without running into the problem of generating impossible data.
To test this the code below simulates some data nsim times from the glm $\mathrm{logit}p_i = \beta_0 + \beta_1 x_i$ with $n_i$ following a Poisson distribution with mean 5. A proportion $1-p$ of the $y_{0i}$ are set missing. Each simulated data set is analysed using multiple imputation using the MICE algorithm, the glm is fitted to the imputed datasets and estimates are pooled using standard methods.
createdata <- function(beta0=0.5, beta1=0.2,
n=100, mx=5, sdx=3, mn=5,
p=0.5) {
size <- rpois(n,lambda=mn) (number of bernoulli trials)
x <- round(rnorm(n,mean=mx,sd=sdx),1)
eta <- beta0 + beta1*x
y1 <- rbinom(n,size,prob=plogis(eta)) # number of successes
y0 <- size - y1 # number of failures
data <- data.frame(y0,y1,x)
missing <- sample(1:n, round(n*(1-p)))
data[missing,"y0"] <- NA
return(data)
}
simulate <- function(nsim=10, seed=2, maxit=5, ...) {
set.seed(seed)
cols <- c("est","se","lambda")
res <- array(NA, dim=c(2, length(cols), nsim, 2))
for (i in 1:nsim) {
data <- createdata(...)
# analysis using mice
imp <- mice(data,print=FALSE,maxit=maxit)
fit <- with(imp, glm(cbind(y1,y0) ~ x, binomial))
est <- pool(fit)
tab <- summary(est)[,cols]
res[,,i,1] <- tab
# complete case analysis
fit <- glm(cbind(y1, y0) ~ x, binomial, data = na.omit(data))
res[,1:2,i,2] <- summary(fit)$coef[,1:2]
}
dimnames(res) <- list(rownames(tab),
cols,
as.character(1:nsim),
c("mice","complete case analysis")
)
res
}
The results from following run is not very promising however, the mean of estimates obtained via multiple imputation taken over all simulated samples indicates bias and a mean square error much greater than from analysis of the data using the complete observation only (true parameter values are $\beta_0=0.5$ and $\beta_1=0.2$.
> res <- simulate(nsim=1000, p=.5, n=100)
> apply(res,c(1,2,4),mean)
, , mice
est se lambda
(Intercept) 0.5615169 0.26826391 0.3570187
x 0.1911394 0.05507577 0.3747948
, , complete case analysis
est se lambda
(Intercept) 0.5073962 0.29463861 NA
x 0.2024261 0.06022877 NA
Increasing the number of iterations of the MICE Gibbs sampler from maxit=5 to maxit=100 does not help
> res <- simulate(nsim=1000, p=.5, n=100,maxit=100)
> apply(res,c(1,2,4),mean)
, , mice
est se lambda
(Intercept) 0.5584763 0.2704222 0.3571965
x 0.1912448 0.0553893 0.3766410
, , complete case analysis
est se lambda
(Intercept) 0.5045417 0.29609451 NA
x 0.2023053 0.06034015 NA
Increasing the sample size from $100$ to $1000$, keeping the proportion of missing constant at $0.5$ do seem to help however
> res <- simulate(nsim=1000, p=.5, n=1000)
> apply(res,c(1,2,4),mean)
, , mice
est se lambda
(Intercept) 0.5064046 0.08436607 0.3925808
x 0.1983546 0.01715852 0.3977634
, , complete case analysis
est se lambda
(Intercept) 0.5030140 0.08994154 NA
x 0.1992938 0.01826588 NA
reducing the mean square error by about 7% relative to complete case analysis:
> apply((res[,1,,] - c(.5,.2))^2, c(1,3),mean)
mice complete case analysis
(Intercept) 0.0074264798 0.0079707049
x 0.0003164675 0.0003255542
