Skip to content
Prev 13016 / 20628 Next
Original Raw

Best way to handle missing data?

Bonnie Dixon 
Thanks, Malcolm!  Yes, that extracts all the results I need from out of the
multiple imputation models.  This method also seems preferable over using
Zelig because it provides the flexibility to create the models using lme4,
or using a different package, such as nlme or MCMCglmm.  And, yes, the
results using this method are identical to the results from Zelig.

Bonnie

P.S.  By the way, that error that occurs when using summary() on the object
created by Zelig can be avoided by using the modified version of the
summary.MI function below (from
http://stackoverflow.com/questions/16571580/multi-level-regression-model-on-multiply-imputed-data-set-in-r-amelia-zelig-l
).

summary.MI <- function (object, subset = NULL, ...) {
  if (length(object) == 0) {
    stop('Invalid input for "subset"')
  } else {
    if (length(object) == 1) {
      return(summary(object[[1]]))
    }
  }

    getcoef <- function(obj) {
    # S4
    if (!isS4(obj)) {
      coef(obj)
    } else {
      if ("coef3" %in% slotNames(obj)) {
        obj at coef3
      } else {
        obj at coef
      }
    }
  }

  res <- list()

  # Get indices
  subset <- if (is.null(subset)) {
    1:length(object)
  } else {
    c(subset)
  }

  # Compute the summary of all objects
  for (k in subset) {
    res[[k]] <- summary(object[[k]])
  }


  # Answer
  ans <- list(
    zelig = object[[1]]$name,
    call = object[[1]]$result at call,
    all = res
  )

  coef1 <- se1 <- NULL

  for (k in subset) {
    tmp <- coef(res[[k]])
    coef1 <- cbind(coef1, tmp[, 1])
    se1 <- cbind(se1, tmp[, 2])
  }

  rows <- nrow(coef1)
  Q <- apply(coef1, 1, mean)
  U <- apply(se1^2, 1, mean)
  B <- apply((coef1-Q)^2, 1, sum)/(length(subset)-1)
  var <- U+(1+1/length(subset))*B
  nu <- (length(subset)-1)*(1+U/((1+1/length(subset))*B))^2

  coef.table <- matrix(NA, nrow = rows, ncol = 4)
  dimnames(coef.table) <- list(rownames(coef1),
                               c("Value", "Std. Error", "t-stat",
"p-value"))
  coef.table[,1] <- Q
  coef.table[,2] <- sqrt(var)
  coef.table[,3] <- Q/sqrt(var)
  coef.table[,4] <- pt(abs(Q/sqrt(var)), df=nu, lower.tail=F)*2
  ans$coefficients <- coef.table
  ans$cov.scaled <- ans$cov.unscaled <- NULL

  for (i in 1:length(ans)) {
    if (is.numeric(ans[[i]]) && !names(ans)[i] %in% c("coefficients")) {
      tmp <- NULL
      for (j in subset) {
        r <- res[[j]]
        tmp <- cbind(tmp, r[[pmatch(names(ans)[i], names(res[[j]]))]])
      }
      ans[[i]] <- apply(tmp, 1, mean)
    }
  }

  class(ans) <- "summaryMI"
  ans
}

On Tue, Mar 3, 2015 at 1:53 AM, Malcolm Fairbrother <
M.Fairbrother at bristol.ac.uk> wrote:

            

      
      
      
    

        

  
  


  


  


  
  

    

  

    
      
    
    
Thread (4 messages)

  

  

    
            
        Bonnie Dixon
        
        Best way to handle missing data?
      
      Mar 2
    
    
      
            
        Malcolm Fairbrother
        
        Best way to handle missing data?
      
      Mar 3
    
    
      
            
        Bonnie Dixon
        
        Best way to handle missing data?
      
      Mar 3
    
    
      
            
        Malcolm Fairbrother
        
        Best way to handle missing data?
      
      Mar 3