3

I'm trying to multiply impute some data, and it's failing in a very odd way. I have options(error=recover) set, and mi() doesn't throw an error, but it does stop imputing after the first attempt. It also doesn't return a value. consequently I have no idea where to even start with debugging. minimal reproducible example below.

> library(mi)
> temp <- mi(dat)
Beginning Multiple Imputation ( Wed Dec 14 10:44:44 2011 ):
Iteration 1 
 Chain 1 : HLTHA5.fac*  BMI*  INCOME*  
> temp
Error: object 'temp' not found


dat<-structure(list(treat = c(FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, 
FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, 
FALSE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TRUE), NUMADULT = c(2, 
1, 2, 1, 2, 1, 2, 1, 2, 4, 1, 1, 2, 2, 1, 2, 1, 1, 2, 1, 2, 2, 
3, 3, 1), HLTHA5.fac = structure(c(3L, NA, 3L, 2L, 4L, 5L, 5L, 
4L, 4L, 3L, 3L, 5L, 3L, 4L, 5L, 4L, 2L, 2L, 3L, 5L, 4L, 5L, 4L, 
3L, 3L), .Label = c("0", "1", "2", "3", "4"), class = "factor"), 
    SOURCEA = structure(c(1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L), .Label = c("Yes", "No", "Don't know", "Refused"), class = "factor"), 
    BMI = c(27.363941459046, 24.0265857515842, 34.3236939308346, 
    27.0907152026518, 32.6101901381975, 34.1643655360753, 21.4628674188624, 
    29.1751398412094, 22.5924920198551, 39.6719545438681, 38.5220574557939, 
    20.1156133421915, 30.6612391698034, 35.7332536282609, 26.5664872147956, 
    25.6016897082437, 19.3649931598758, 28.1868713091175, NA, 
    32.4438116170843, 32.5507197719099, 21.1090717674633, 32.2340044872853, 
    24.3699149340904, 27.3369153440247), SMOKE2 = structure(c(2L, 
    1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L), .Label = c("Yes", "No"
    ), class = "factor"), INCOME = structure(c(16L, 4L, 13L, 
    11L, 13L, 7L, 22L, 6L, NA, 1L, 13L, 18L, 12L, 20L, NA, NA, 
    2L, 13L, 17L, NA, 12L, 21L, 9L, 15L, 13L), .Label = c("Less than $10,800", 
    "$10,800-$14,600", "$14,601-$16,250", "$16,251-$18,300", 
    "$18,301-$21,800", "$21,801-$25,000", "$25,001-$27,500", 
    "$27,501-$29,300", "$29,301-$33,100", "$33,101-$36,700", 
    "$36,701-$38,700", "$38,701-$44,200", "$44,201-$50,000", 
    "$50,001-$58,000", "$58,001-$66,500", "$66,501-$73,500", 
    "$73,501-$80,000", "$80,001-88,200", "$88,201-$100,000", 
    "$100,001-$120,000", "$120,001-$130,000", "$130,001-$150,000", 
    "$150,001-$250,000", "Over $250,000", "Don't know", "Refused"
    ), class = "factor"), RESPMAR = structure(c(1L, 5L, 1L, 4L, 
    3L, 6L, 1L, 1L, 1L, 1L, 4L, 4L, 1L, 1L, 4L, 1L, 3L, 3L, 1L, 
    3L, 1L, 1L, 1L, 1L, 1L), .Label = c("Married", "Living w partner", 
    "Widowed", "Divorced", "Separated", "Single", "Other", "Don't know", 
    "Refused"), class = "factor"), RESPGRAD = structure(c(5L, 
    1L, 2L, 5L, 3L, 3L, 5L, 2L, 4L, 2L, 4L, 3L, 4L, 4L, 2L, 3L, 
    2L, 5L, 2L, 4L, 4L, 5L, 2L, 2L, 3L), .Label = c("< HS 0-11", 
    "HS graduate", "Some colge 13-15", "Collge grad 16", "Post college 16+", 
    "Don't know", "Refused"), class = "factor"), RACEA2 = structure(c(1L, 
    1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 
    3L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L), .Label = c("White (Not-Latino)", 
    "Black (Not-Latino)", "Latino (total)", "Asian", "Biracial/Multi", 
    "Native American", "Other", "Don't know", "Refused"), class = "factor"), 
    INSUREDA = structure(c(1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 
    1L), .Label = c("Insured", "Not insured", "Don't know", "Refused"
    ), class = "factor"), PAP.adj = c(TRUE, FALSE, FALSE, FALSE, 
    FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, TRUE, 
    TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, 
    TRUE, TRUE)), .Names = c("treat", "NUMADULT", "HLTHA5.fac", 
"SOURCEA", "BMI", "SMOKE2", "INCOME", "RESPMAR", "RESPGRAD", 
"RACEA2", "INSUREDA", "PAP.adj"), row.names = c(1L, 13L, 15L, 
23L, 26L, 33L, 38L, 53L, 56L, 60L, 62L, 85L, 109L, 116L, 138L, 
217L, 240L, 262L, 264L, 269L, 277L, 295L, 328L, 334L, 338L), class = "data.frame")

Any ideas as to where to begin?

Update

Thanks to the diagnostic technique below, I tracked down the error, which I am summarizing here since it appears I am not the only one to have this problem.

The error occurs when you have an unordered categorical variable with levels that have no values. mi.default calls .initializeConvCheckArray to fill in AveVar with NAs. That function uses the levels, regardless of whether those levels are used. By contrast, to fill in AveVar it calls .getmean, which drops unused levels. Therefore the dimensions don't match up.

The easy solution on the user end is of course to drop the extra levels before calling mi.info and mi. I'm going to submit a bug fix to the package authors, though, having spent way too much time tracking this down already.

Ari B. Friedman
  • 71,271
  • 35
  • 175
  • 235

1 Answers1

7

Since the error=recover option is not working, an alternative that does work is to set options(error=dump.frames). That will get you some information about the error, which you can either print out or, more usefully, examine with debugger()

ls()
# [1] "dat"
options(error=dump.frames)
mi(dat)
ls()
# [1] "dat"       "last.dump"  # Apparently there WAS an error

# INVESTIGATE WITH debugger()
debugger(dump=last.dump)

# ALTERNATIVELY, PRINT last.dump TO CONSOLE
last.dump

$`mi(dat)`
<environment: 0x05155c44>

$`mi(dat)`
<environment: 0x05158f30>

$`.local(object, ...)`
<environment: 0x05158cac>

$`mi.default(object, info, n.imp, n.iter, R.hat, max.minutes, rand.imp.method`
<environment: 0x047dc3a0>

attr(,"error.message")
[1] "Error in aveVar[s, i, ] <- c(avevar.mean, avevar.sd) :
\n  number of items to replace is not a multiple of replacement length\n"
attr(,"class")
[1] "dump.frames"
Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
  • @PaulHiemstra and gsk3 -- Yeah, I had read about years ago, but never found a need for it. Even Chambers (2008) says "Usually, there is no advantage to dump.frames(), since recover() behaves like dump.frames() if the computations are not interactive." – Josh O'Brien Dec 14 '11 at 18:51
  • Well now there's a use for it! I added it to my list of debugging suggestions here: http://stackoverflow.com/a/5156351/636656 – Ari B. Friedman Dec 14 '11 at 18:53
  • It's prob. worth noting that both the relevant help files and Chambers (2008) have mild warnings that `debugger(dump=last.dump)` is not always exactly equivalent to `recover()`. Things like promise objects, I/O on connections, and interactive graphics aren't saved well, so are harder to debug with `dump.frames`. – Josh O'Brien Dec 14 '11 at 18:55