# Multiple imputation of glm binomial size parameter

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  • You could treat those few missing$n_i$'s as parameters to be estimated, see stats.stackexchange.com/questions/123367/… solutions in answers to that could be adapted. – kjetil b halvorsen May 24 '17 at 10:48 • I may be wrong here but each nest has its separate$n_i$(unlike the problem in the link) so if this is estimated as a free parameter, wouldn't this in effect be equivalent to removing the observation all together, that is, if$n_i$is set equal to its MLE, then$y_i$would equal its expected value and would carry no information about the effects of the covariates in the linear predictor? With multiple imputation the idea would be to simulate a plausible values of$n_i$within the range of values seen in cases where it is non-missing while at the same time satisfying$n_i\ge y_i$. – Jarle Tufto May 24 '17 at 11:07 • Maybe. Would you say could be implemented with bayes methods, with a prior on$n_i$implementing what you know about the variability in the$n_i$? Simulation-based inference with such a (informative) prior would be a kind of multiple imputation. – kjetil b halvorsen May 24 '17 at 11:09 • What's the purpose/point of imputation in this case? Why not simply discard those data points? – Jacob Socolar May 27 '17 at 15:26 • Discarding those data points would be to throw away information. While$y_i$contains less information about the effects of the covariates when$n_i\$ is missing, these observations still contain some information. – Jarle Tufto May 27 '17 at 16:11