R/mice.impute.quasipoisson.R

In kkleinke/countimp: Multiple Imputation of incomplete count data

#' Multiple Imputation of Overdispersed Count Data based on a \code{quasipoisson} GLM
#' 
#' Imputes univariate missing data based  on a \code{quasipoisson} GLM following either the Bayesian regression or bootstrap regression (appendix \code{.boot}) MI approach.
#' 
#' Overdispersed count data (meaning that the variance of the count variable is larger than its the mean) are typically analyzed by a negative binomial (NB) or by a \code{quasipoisson} model. The \code{quasipoisson} model is identical to an ordinary Poisson model, except that it estimates an additional dispersion parameter. For details, see Zeileis, Kleiber, & Jackman (2008), or Hilbe (2007).
#' The Bayesian method consists of the following steps:
#' \enumerate{
#' \item Fit the model, and find bhat, the posterior mean, and V(bhat), the posterior variance of model parameters b.
#' \item Draw b.star from N(bhat,V(bhat)).
#' \item Compute fitted values using \code{exp(x[!ry, ] \%*\% b.star)}
#' \item Simulate imputations from a negative binomial distribution to ensure an adequate dispersion of imputed values.
#' }
#' \code{quasipoisson} imputation relies on the standard \code{glm.fit} function, using the \code{quasipoisson} family. 
#' The bootstrap method draws a bootstrap sample from \code{y[ry]} and \code{x[ry,]} and consists of the following steps:
#' \enumerate{
#' \item Fit the model to the bootstrap sample and get model parameters \code{b.star}
#' \item Compute fitted values using \code{exp(x[!ry, ] \%*\% b.star)}
#' \item Simulate imputations from a negative binomial distribution to ensure an adequate dispersion of imputed values.
#' }
#' @param y Numeric vector with incomplete data
#' @param ry Response pattern of \code{y} (\code{TRUE}=observed, \code{FALSE}=missing)
#' @param x matrix with \code{length(y)} rows containing complete covariates
#' @param wy Logical vector of length \code{length(y)}. A \code{TRUE} value indicates locations in \code{y} for which imputations are created. Default is \code{!ry}
#' @param EV should automatic outlier handling of imputed values be enabled?  Default is \code{TRUE}: extreme imputations will be identified. These values will be replaced by imputations obtained by predictive mean matching (function \code{mice.impute.midastouch()})
#' @param ... Other named arguments.
#' @return Numeric vector of length \code{sum(!ry)} with imputations
#' @aliases mice.impute.quasipoisson.boot mice.impute.qpois.boot mice.impute.qpois
#' @references 
#' \itemize{
#' \item Hilbe, J. M. (2007). \emph{Negative binomial regression}. Cambridge: Cambridge University Press.
#' \item Kleinke, K., & Reinecke, J. (2013). \emph{countimp 1.0 -- A multiple imputation package for incomplete count data} [Technical Report]. University of Bielefeld, Faculty of Sociology, available from \url{www.uni-bielefeld.de/soz/kds/pdf/countimp.pdf}.
#' \item Rubin, D. B. (1987). \emph{Multiple imputation for nonresponse in surveys}. New York: Wiley.
#' \item Zeileis, A., Kleiber, C., & Jackman, S. (2008). Regression models for count data in R. \emph{Journal of Statistical Software}, 27(8), 1–-25.
#' }
#' @importFrom stats coef glm.fit quasipoisson rnbinom rpois summary.glm
#' @importFrom MASS rnegbin
#' @examples 
#' ## simulate overdespersed count data
#' set.seed( 1234 )
#' b0 <- 1
#' b1 <- .75
#' b2 <- -.25
#' b3 <- .5
#' N <- 5000
#' x1 <- rnorm(N)
#' x2 <- rnorm(N)
#' x3 <- rnorm(N)
#' mu <- exp( b0 + b1 * x1 + b2 * x2 + b3 * x3 )
#' y <- MASS::rnegbin( N, theta = 2, mu )
#' NB <- data.frame( y, x1, x2, x3 )
#' 
#' ## introduce MAR missingness to simulated data
#' total <- round( .2 * N )  ##number of missing data in y
#' sm <- which( NB[,2] < mean( NB[,2] ) )  ##subset: cases with x2<mean(x2)
#' gr <- which( NB[,2] > mean( NB[,2] ) )	##subset: cases with x2>mean(x2)
#' sel.sm <- sample( sm, round( .2 * total ) )	##select cases to set as missing
#' sel.gr <- sample( gr, round( .8 * total ) )	##select cases to set as missing
#' sel <- c( sel.sm,sel.gr )
#' MNB <- NB
#' MNB[sel,1] <- NA	##delete selected data
#' 
#' ## impute missing data
#' imp <- countimp( MNB, method = c( "quasipoisson", "", "", "" )) 
#' @author Kristian Kleinke
#' @describeIn  mice.impute.quasipoisson Bayesian regression variant
#' @export
mice.impute.quasipoisson <-
function (y, ry, x, wy = NULL, EV=TRUE, ...) 
{
  if (is.null(wy)) 
    wy <- !ry
    x <- cbind(1, as.matrix(x))
    fit <- glm.fit(x[ry, ], y[ry], family = quasipoisson(link = log))
    fit.sum <- summary.glm(fit)
    beta <- coef(fit)
    rv <- t(chol(fit.sum$cov.unscaled))
    beta.star <- beta + rv %*% rnorm(ncol(rv))
    p <- exp((x[wy, , drop = FALSE] %*% beta.star))
	  im=rnbinom(n=length(p),size=(p/(fit.sum$dispersion-1)),mu=p) 
	  imputed.values<-im
	  if(EV){
	    outliers <- getOutliers(imputed.values, rho = c(0.3, 
	                                                    0.3), FLim = c(0.15, 0.85))
	    nans <- which(is.nan(imputed.values))
	    idx <- c(outliers$iLeft, outliers$iRight, nans)
	    if (length(idx) != 0) {
	      imputed.values[idx] <- NA
	      y[!ry] <- imputed.values
	      R = ry
	      ry <- !is.na(y)
	      new.values <- mice.impute.midastouch(y, ry, x, 
	                                           wy = NULL)
	      imputed.values[idx] <- new.values
	    }}
	  return(imputed.values)
}

#' @export
#' @describeIn  mice.impute.quasipoisson identical to \code{mice.impute.quasipoisson()}; included for backward compatibility
mice.impute.qpois <-
  function (y, ry, x, wy = NULL, EV=TRUE, ...) 
  {
    if (is.null(wy)) 
      wy <- !ry
    x <- cbind(1, as.matrix(x))
    fit <- glm.fit(x[ry, ], y[ry], family = quasipoisson(link = log))
    fit.sum <- summary.glm(fit)
    beta <- coef(fit)
    rv <- t(chol(fit.sum$cov.unscaled))
    beta.star <- beta + rv %*% rnorm(ncol(rv))
    p <- exp((x[wy, , drop = FALSE] %*% beta.star))
    im=rnbinom(n=length(p),size=(p/(fit.sum$dispersion-1)),mu=p) 
    imputed.values<-im
    if(EV){
      outliers <- getOutliers(imputed.values, rho = c(0.3, 
                                                      0.3), FLim = c(0.15, 0.85))
      nans <- which(is.nan(imputed.values))
      idx <- c(outliers$iLeft, outliers$iRight, nans)
      if (length(idx) != 0) {
        imputed.values[idx] <- NA
        y[!ry] <- imputed.values
        R = ry
        ry <- !is.na(y)
        new.values <- mice.impute.midastouch(y, ry, x, 
                                             wy = NULL)
        imputed.values[idx] <- new.values
      }}
    return(imputed.values)
  }

#' @export
#' @describeIn  mice.impute.quasipoisson Bootstrap regression variant
mice.impute.quasipoisson.boot <-
  function (y, ry, x, wy = NULL, EV=TRUE, ...) 
  {
    if (is.null(wy)) 
      wy <- !ry
    x <- cbind(1, as.matrix(x))
    xobs<-x[ry,]
    yobs<-y[ry]
    sel<-sample(1:length(yobs),length(yobs),replace=TRUE)
    xast<-xobs[sel,]
    yast<-yobs[sel]
    fit <- glm.fit(xast, yast, family = quasipoisson(link = log))
    fit.sum <- summary.glm(fit)
    beta.star <- coef(fit)
    p <- exp((x[wy, , drop = FALSE] %*% beta.star))
    if (fit.sum$dispersion>1)
    {ret=rnbinom(n=length(p),size=(p/(fit.sum$dispersion-1)),mu=p)}
    else
    {
      ret=rpois(length(p),p)
    } 
    imputed.values<-ret
    if(EV){
      outliers <- getOutliers(imputed.values, rho = c(0.3, 
                                                      0.3), FLim = c(0.15, 0.85))
      nans <- which(is.nan(imputed.values))
      idx <- c(outliers$iLeft, outliers$iRight, nans)
      if (length(idx) != 0) {
        imputed.values[idx] <- NA
        y[!ry] <- imputed.values
        R = ry
        ry <- !is.na(y)
        new.values <- mice.impute.midastouch(y, ry, x, 
                                             wy = NULL)
        imputed.values[idx] <- new.values
      }}
    return(imputed.values)
  }


#' @export
#' @describeIn  mice.impute.quasipoisson identical to \code{mice.impute.quasipoisson.boot()}; included for backward compatibility
mice.impute.qpois.boot <-
  function (y, ry, x, wy = NULL, EV=TRUE, ...) 
  {
    if (is.null(wy)) 
      wy <- !ry
    x <- cbind(1, as.matrix(x))
    xobs<-x[ry,]
    yobs<-y[ry]
    sel<-sample(1:length(yobs),length(yobs),replace=TRUE)
    xast<-xobs[sel,]
    yast<-yobs[sel]
    fit <- glm.fit(xast, yast, family = quasipoisson(link = log))
    fit.sum <- summary.glm(fit)
    beta.star <- coef(fit)
    p <- exp((x[wy, , drop = FALSE] %*% beta.star))
    if (fit.sum$dispersion>1)
    {ret=rnbinom(n=length(p),size=(p/(fit.sum$dispersion-1)),mu=p)}
    else
    {
      ret=rpois(length(p),p)
    } 
    imputed.values<-ret
    if(EV){
      outliers <- getOutliers(imputed.values, rho = c(0.3, 
                                                      0.3), FLim = c(0.15, 0.85))
      nans <- which(is.nan(imputed.values))
      idx <- c(outliers$iLeft, outliers$iRight, nans)
      if (length(idx) != 0) {
        imputed.values[idx] <- NA
        y[!ry] <- imputed.values
        R = ry
        ry <- !is.na(y)
        new.values <- mice.impute.midastouch(y, ry, x, 
                                             wy = NULL)
        imputed.values[idx] <- new.values
      }}
    return(imputed.values)
  }