Merge pull request #36 from inbo/bugfix

Bugfix

.Rbuildignore

@@ -14,4 +14,4 @@
 ^_steps$
 ^NEWS\.md$
 ^CODE_OF_CONDUCT\.md$
-$CONTRIBUTING\.md$
+CONTRIBUTING.md

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: multimput
 Type: Package
 Title: Using Multiple Imputation to Address Missing Data
-Version: 0.2.8
-Date: 2017-12-10
+Version: 0.2.8.1
+Date: 2018-01-15
 Authors@R: c(
   person(
     "Thierry", "Onkelinx", role = c("aut", "cre"), 

R/aggregate_impute.R

@@ -99,6 +99,8 @@ setMethod(
     imputation <- imputation[
       data[[id_column]] %>%
         na.omit(),
+      ,
+      drop = FALSE
     ]
 
     missing.obs <- which(is.na(data[, response]))

R/impute_inla.R

@@ -8,50 +8,67 @@ setMethod(
   definition = function(model, ..., n.imp){
     assert_that(is.count(n.imp))
 
-    if (!model$.args$control.predictor$compute) {
+
+    if (!model$.args$control.compute$config) {
       stop(
-"model must be fit with the 'compute = TRUE' argument of control.predictor"
+"model must be fit with the 'config = TRUE' argument of control.compute"
       )
     }
 
+    dots <- list(...)
+    if (is.null(dots$minimum)) {
+      dots$minimum <- ""
+    }
+
     response <- as.character(model$.args$formula)[2]
     missing.obs <- which(is.na(model$.args$data[, response]))
+    if (length(missing.obs) == 0) {
+      return(
+        new(
+          "rawImputed",
+          Data = model$.args$data,
+          Response = response,
+          Imputation = matrix(integer(0), ncol = n.imp),
+          Minimum = dots$minimum
+        )
+      )
+    }
+
+    magnitude <- ceiling(log10(nrow(model$.args$data)))
+    missing.obs <- sprintf(paste0("Predictor:%0", magnitude, "i"), missing.obs)
 
     assert_that(requireNamespace("INLA", quietly = TRUE))
     imputation <- switch(
       model$.args$family,
       poisson = {
-        linpred <- sapply(
-          model$marginals.linear.predictor[missing.obs],
-          INLA::inla.rmarginal,
-          n = n.imp
+        samples <- INLA::inla.posterior.sample(
+          n = n.imp,
+          model
         )
-        matrix(
-          rpois(length(linpred), lambda = exp(linpred)),
-          ncol = n.imp,
-          byrow = TRUE
+        sapply(
+          samples,
+          function(x) {
+            rpois(
+              n = length(missing.obs),
+              lambda = exp(x$latent[missing.obs, 1])
+            )
+          }
         )
       },
       nbinomial = {
-        linpred <- sapply(
-          model$marginals.linear.predictor[missing.obs],
-          INLA::inla.rmarginal,
-          n = n.imp
-        )
-        h <- model$marginals.hyperpar
-        h <- h[grepl("size for the nbinomial", names(h))]
-        size <- INLA::inla.rmarginal(
+        samples <- INLA::inla.posterior.sample(
           n = n.imp,
-          marginal = h[[1]]
+          model
         )
-        matrix(
-          rnbinom(
-            n = length(linpred),
-            size = rep(size, ncol(linpred)),
-            mu = exp(linpred)
-          ),
-          ncol = n.imp,
-          byrow = TRUE
+        sapply(
+          samples,
+          function(x) {
+            h <- x$hyperpar
+            rnbinom(
+              n = length(missing.obs),
+              size = h[grepl("size for the nbinomial", names(h))],
+              mu = exp(x$latent[missing.obs, 1]))
+          }
         )
       },
       stop(
@@ -61,10 +78,6 @@ a reproducible example at https://github.com/ThierryO/multimput/issues"
       )
     )
 
-    dots <- list(...)
-    if (is.null(dots$minimum)) {
-      dots$minimum <- ""
-    }
     new(
       "rawImputed",
       Data = model$.args$data,

R/model_impute.R

@@ -112,20 +112,31 @@ setMethod(
     object@Imputation <- object@Imputation[object@Covariate[[id_column]], ]
 
     form <- as.formula(paste("Imputed", rhs, sep = "~"))
-    lapply(
+    m <- lapply(
       seq_len(ncol(object@Imputation)),
       function(i){
         data <- cbind(
           Imputed = object@Imputation[, i],
           object@Covariate
         )
-        model.args <- c(list(data = data), model.args)
-        model <- do.call(model.fun, c(form, model.args))
-        do.call(extractor, c(list(model), extractor.args)) %>%
-          as.data.frame() %>%
-          rownames_to_column("Variable")
+        model <- try(
+          do.call(model.fun, c(form, list(data = data), model.args)),
+          silent = TRUE
+        )
+        if (inherits(model, "try-error")) {
+          NULL
+        } else {
+          do.call(extractor, c(list(model), extractor.args)) %>%
+            as.data.frame() %>%
+            rownames_to_column("Variable")
+        }
       }
-    ) %>%
+    )
+    failed <- sapply(m, is.null)
+    if (all(failed)) {
+      stop("model failed on all imputations")
+    }
+    m %>%
       bind_rows() %>%
       select_(Parameter = 1, Estimate = 2, SE = 3) %>%
       mutate_(

tests/testthat/test_ab_impute.R

@@ -68,6 +68,7 @@ describe("impute", {
       Count ~ factor(Year) + factor(Period) + f(Site, model = "iid"),
       data = dataset,
       family = "poisson",
+      control.compute = list(config = TRUE),
       control.predictor = list(compute = TRUE, link = 1)
     )
     expect_is(
@@ -127,6 +128,7 @@ describe("impute", {
       Count ~ factor(Year) + factor(Period) + f(Site, model = "iid"),
       data = dataset,
       family = "nbinomial",
+      control.compute = list(config = TRUE),
       control.predictor = list(compute = TRUE, link = 1)
     )
     expect_is(
@@ -217,6 +219,21 @@ describe("impute", {
       impute(model, dataset, minimum = "Junk"),
       "object@Data does not have name Junk"
     )
+
+    if (!require(INLA)) {
+      skip("INLA package not available")
+    }
+    model <- inla(
+      Count ~ Year + factor(Period) + factor(Site),
+      data = dataset,
+      family = "nbinomial",
+      control.compute = list(config = TRUE)
+    )
+    expect_is(
+      imputed <- impute(model, dataset, minimum = "Bottom"),
+      "rawImputed"
+    )
+
   })
 
 
@@ -264,21 +281,21 @@ describe("impute", {
     model <- INLA::inla(
       Mu ~ factor(Year) + factor(Period) + f(Site, model = "iid"),
       data = dataset,
-      family = "gamma"
+      family = "nbinomial"
     )
     expect_error(
       impute(model),
-"model must be fit with the 'compute = TRUE' argument of control.predictor"
+"model must be fit with the 'config = TRUE' argument of control.compute"
     )
     model <- INLA::inla(
-      Mu ~ factor(Year) + factor(Period) + f(Site, model = "iid"),
+      Count ~ factor(Year) + factor(Period) + f(Site, model = "iid"),
       data = dataset,
-      family = "gamma",
-      control.predictor = list(compute = TRUE, link = 1)
+      family = "gaussian",
+      control.compute = list(config = TRUE)
     )
     expect_error(
       impute(model),
-      "Imputations from the 'gamma' family not yet defined"
+      "Imputations from the 'gaussian' family not yet defined"
     )
 
     model <- lme4::glmer(

tests/testthat/test_ccc_model_impute.R

@@ -119,7 +119,7 @@ describe("model_impute", {
         rhs = "junk",
         extractor = extractor
       ),
-      "object 'junk' not found"
+      "model failed on all imputations"
     )
   })
 })

vignettes/Impute.Rmd

@@ -1,7 +1,7 @@
 ---
 title: "Model data with missing observations using multiple imputation"
 author: "Thierry Onkelinx"
-date: "`r Sys.Date()`"
+date: "2018-01-12"
 bibliography: multimput.bib
 output: rmarkdown::html_vignette
 vignette: >
@@ -128,13 +128,15 @@ imp.inla.p <- inla(
   Observed ~ fYear + fPeriod + f(Site, model = "iid"), 
   data = dataset, 
   family = "poisson", 
+  control.compute = list(config = TRUE),
   control.predictor = list(compute = TRUE, link = 1)
 )
 # the same model as imp.inla.p but with negative binomial distribution
 imp.inla.nb <- inla(
   Observed ~ fYear + fPeriod + f(fSite, model = "iid"), 
   data = dataset, 
   family = "nbinomial", 
+  control.compute = list(config = TRUE),
   control.predictor = list(compute = TRUE, link = 1)
 )
 # a mixed model with negative binomial distribution
@@ -153,6 +155,7 @@ imp.better <- inla(
     fPeriod, 
   data = dataset, 
   family = "nbinomial", 
+  control.compute = list(config = TRUE),
   control.predictor = list(compute = TRUE, link = 1)
 )
 ```

wercker.yml

@@ -1,4 +1,4 @@
-box: inbobmk/rstable:3.3.2
+box: inbobmk/rstable
 build:
   steps:
     - inbobmk/r-check
@@ -8,10 +8,6 @@ builddedicated:
   box: inbobmk/r-multimput
   steps:
     - inbobmk/r-check
-buildold:
-  box: inbobmk/rstable:3.3.0
-  steps:
-    - inbobmk/r-check
 
 master:
   deploy: