I came across a small bug in ameliabind. The returned amelia object does
not include an overvalues element. The following
code (a modified version of the example in moPrep) produces the error
library(Amelia)
Loading required package: foreign
##
## Amelia II: Multiple Imputation
## (Version 1.5-4, built: 2011-08-21)
## Copyright (C) 2005-2011 James Honaker, Gary King and Matthew Blackwell
## Refer to
http://gking.harvard.edu/amelia/ for more information
##
data(africa)
m.out <- moPrep(africa, trade ~ trade, error.proportion = 0.1)
## Without ameliabind
a.out1 <- amelia(m.out, m=2, ts = "year", cs = "country")
-- Imputation 1 --
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
-- Imputation 2 --
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
a.out1$overvalues
[1] 29.69 31.31
35.22 40.11 37.76 41.11 37.71 40.21 43.09 43.64
[11] 42.62 40.11 42.72 44.63 41.83 41.56 37.92 35.19 38.37 38.75
[21] 26.81 24.35 25.29 27.28 30.38 34.89 31.93 36.88 32.10 31.18
[31] 37.12 33.79 35.06 31.67 35.48 38.55 38.42 46.48 48.24 48.23
[41] 50.01 51.78 49.32 54.25 61.55 59.59 61.93 63.53 65.02 46.01
[51] 37.32 31.99 38.48 37.71 34.63 80.10 75.00 112.70 99.64 106.93
[61] 110.85 104.79 96.04 120.14 134.11 123.77 109.71 107.25 112.79 93.50
[71] 80.36 81.05 83.81 97.91 90.79 65.20 67.23 91.44 78.37 80.62
[81] 94.29 75.34 73.88 72.26 83.56 80.70 78.12 85.34 70.63 58.45
[91] 55.49 51.87 58.76 58.91 56.38 86.04 82.45 91.26 92.85 81.72
[101] 81.27 70.48 81.86 86.80 69.78 64.16 64.41 68.30 73.64 86.02
[111] 75.16 59.47 60.63 72.47 71.86
## With ameliabind
a.out2 <- do.call(ameliabind,
+ replicate(2, amelia(m.out,
m=1, ts = "year", cs =
"country"),
+ simplify=FALSE))
-- Imputation 1 --
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
-- Imputation 1 --
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21
print(a.out2$overvalues)
NULL
print(try(plot(a.out2)))
Error in
x[[jj]][iseq] <- vjj : replacement has length zero
[1] "Error in x[[jj]][iseq] <- vjj : replacement has length zero\n"
attr(,"class")
[1] "try-error"
This patch should fix the bug. Additionally, this patch adds a feature to
ameliabind. Rather than throwing an error if there is only one amelia
object as an argument, it simply returns that object. This way something
like do.call(ameliabind, replicate(n, amelia(...))) can work for n = 1,
e.g. debugging and testing, without altering the code.
diff -rupNb -x .git -x '*.pdf' -x '*.Rout' Amelia/man/ameliabind.Rd
/home/jeff/workspace/AmeliaTest2/Amelia/man/ameliabind.Rd
--- Amelia/man/ameliabind.Rd 2011-08-21 09:00:21.000000000 -0400
+++ /home/jeff/workspace/AmeliaTest2/Amelia/man/ameliabind.Rd 2011-11-11
15:59:48.498846781 -0500
@@ -10,7 +10,7 @@
ameliabind(...)
}
\arguments{
- \item{...}{two or more objects of class \code{amelia} with the same
+ \item{...}{one or more objects of class \code{amelia} with the same
arguments and created from the same data.}
}
\details{
diff -rupNb -x .git -x '*.pdf' -x '*.Rout' Amelia/R/emb.r
/home/jeff/workspace/AmeliaTest2/Amelia/R/emb.r
--- Amelia/R/emb.r 2011-08-21 09:00:21.000000000 -0400
+++ /home/jeff/workspace/AmeliaTest2/Amelia/R/emb.r 2011-11-09
12:11:39.763371625 -0500
@@ -865,7 +865,7 @@ am.inv <- function(a,tol=.Machine$double
##
## ameliabind - combines multiple Amelia outputs
##
-## INPUTS: >1 amelia output
+## INPUTS: >=1 amelia output
##
## OUTPUTS: a merged amelia output with all inputted lists
##
@@ -873,12 +873,11 @@ am.inv <- function(a,tol=.Machine$double
ameliabind <- function(...) {
args <- list(...)
- if (length(args) < 2)
- stop("We need at least two amelia outputs to bind")
-
- if (any(lapply(args, class)!="amelia"))
+ if (any(!sapply(args, is, "amelia")))
stop("All arguments must be amelia output.")
+ if (length(args) > 1) {
+
## test that data is the same. we'll just compare the missMatrices.
## this will allow datasets with the same size and missingness
## matrix to be combined unintentionally, but this seems unlikely.
@@ -905,6 +904,9 @@ ameliabind <- function(...) {
mu = matrix(NA, nrow = k, ncol = newm),
covMatrices = array(NA, dim = c(k,k,newm)),
code = integer(0),
+ ## overvalues for all args assumed to be the same
+ ## so only use overvalues from first object
+ overvalues = args[[1]]$overvalues,
message = character(0),
iterHist = list(),
arguments = list())
@@ -935,6 +937,9 @@ ameliabind <- function(...) {
}
class(out) <- "amelia"
class(out$imputations) <- c("mi","list")
+ } else {
+ out <- args
+ }
return(out)
}
Jeff
---
Jeffrey Arnold
Department of Political Science
University of Rochester
http://jrnold.me
jeffrey.arnold(a)gmail.com
jeffrey.arnold(a)rochester.edu