title: Can I pool imputed random effect model estimates using the mi package?
body: <p>It appears that the <code>mi</code> package has had a pretty big rewrite at some point within the past couple of years. </p>
  
  <p>The "old" way of doing things is well-outlined in the following tutorial: <a href="http://thomasleeper.com/Rcourse/Tutorials/mi.html" rel="nofollow">http://thomasleeper.com/Rcourse/Tutorials/mi.html</a></p>
  
  <p>The "new" way of doing things (sticking with Leeper's simulation demo) looks something like this: </p>
  
  <pre><code>#load mi
  library(mi)
  #set seed
  set.seed(10)
  #simulate some data (with some observations missing)
  x1 &lt;- runif(100, 0, 5)
  x2 &lt;- rnorm(100)
  y &lt;- 2*x1 + 20*x2 + rnorm(100)
  mydf &lt;- cbind.data.frame(x1, x2, y)
  mydf$x1[sample(1:nrow(mydf), 20, FALSE)] &lt;- NA
  mydf$x2[sample(1:nrow(mydf), 10, FALSE)] &lt;- NA
  
  # Convert to a missing_data.frame
  mydf_mdf &lt;- missing_data.frame(mydf)
  
  # impute
  mydf_imp &lt;- mi(mydf_mdf)
  </code></pre>
  
  <p>Although function names have changed, this is actually pretty similar to the "old" way of doing things. </p>
  
  <p>The biggest change (from my vantage) is the replacement of the following "old" functions </p>
  
  <p><code>lm.mi(formula, mi.object, ...)</code></p>
  
  <p><code>glm.mi(formula, mi.object, family = gaussian, ...)</code></p>
  
  <p><code>bayesglm.mi(formula, mi.object, family = gaussian, ...)</code></p>
  
  <p><code>polr.mi(formula, mi.object, ...)</code></p>
  
  <p><code>bayespolr.mi(formula, mi.object, ...)</code></p>
  
  <p><code>lmer.mi(formula, mi.object, rescale=FALSE, ...)</code></p>
  
  <p><code>glmer.mi(formula, mi.object, family = gaussian, rescale=FALSE, ...)</code>. </p>
  
  <p>Previously, a user could compute a model for each imputed dataset using one of these functions and then pool the results using <code>mi.pooled()</code> (or <code>coef.mi()</code> if we are following the Leeper example).</p>
  
  <p>In the current version of <code>mi</code> (I have v1.0 installed), these last steps appear to have been combined into a single function, <code>pool()</code>. The <code>pool()</code> function appears to read the family and link function that was assigned to a variable during the imputation process above and then estimate a model with <code>bayesglm</code> using the specified formula as shown below. </p>
  
  <pre><code># run models on imputed data and pool the results
  summary(pool(y ~ x1 + x2, mydf_imp))
  
  ## 
  ## Call:
  ## pool(formula = y ~ x1 + x2, data = mydf_imp)
  ## 
  ## Deviance Residuals: 
  ##      Min        1Q    Median        3Q       Max  
  ## -1.98754  -0.40923   0.03393   0.46734   2.13848  
  ## 
  ## Coefficients:
  ##             Estimate Std. Error t value Pr(&gt;|t|)    
  ## (Intercept) -0.34711    0.25979  -1.336    0.215    
  ## x1           2.07806    0.08738  23.783 1.46e-13 ***
  ## x2          19.90544    0.11068 179.844  &lt; 2e-16 ***
  ## ---
  ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ## 
  ## (Dispersion parameter for gaussian family taken to be 0.7896688)
  ## 
  ##     Null deviance: 38594.916  on 99  degrees of freedom
  ## Residual deviance:    76.598  on 97  degrees of freedom
  ## AIC: 264.74
  ## 
  ## Number of Fisher Scoring iterations: 7
  </code></pre>
  
  <p>This looks like we are coming close to recovering our simulated beta values (2 and 20). In other words, it is behaving as expected. </p>
  
  <p>Let's take a slightly larger set of data with a naively simulated random effect just for the sake of getting a grouping variable. </p>
  
  <pre><code>mydf2 &lt;- data.frame(x1 = rep(runif(100, 0, 5), 20)
                     ,x2 = rep(rnorm(100, 0, 2.5), 20)
                     ,group_var = rep(1:20, each = 100)
                     ,noise = rep(rnorm(100), 20))
  
  mydf2$y &lt;- 2*mydf2$x1 + 20*mydf2$x2 + mydf2$noise
  
  mydf2$x1[sample(1:nrow(mydf2), 200, FALSE)] &lt;- NA
  mydf2$x2[sample(1:nrow(mydf2), 100, FALSE)] &lt;- NA
  
  # Convert to a missing_data.frame
  mydf2_mdf &lt;- missing_data.frame(mydf2)
  
  show(mydf2_mdf)
  
  ## Object of class missing_data.frame with 2000 observations on 5 variables
  ## 
  ## There are 4 missing data patterns
  ## 
  ## Append '@patterns' to this missing_data.frame to access the corresponding pattern for every observation or perhaps use table()
  ## 
  ##                 type missing method  model
  ## x1        continuous     200    ppd linear
  ## x2        continuous     100    ppd linear
  ## group_var continuous       0   &lt;NA&gt;   &lt;NA&gt;
  ## noise     continuous       0   &lt;NA&gt;   &lt;NA&gt;
  ## y         continuous       0   &lt;NA&gt;   &lt;NA&gt;
  ## 
  ##             family     link transformation
  ## x1        gaussian identity    standardize
  ## x2        gaussian identity    standardize
  ## group_var     &lt;NA&gt;     &lt;NA&gt;    standardize
  ## noise         &lt;NA&gt;     &lt;NA&gt;    standardize
  ## y             &lt;NA&gt;     &lt;NA&gt;    standardize
  </code></pre>
  
  <p>Since <code>missing_data.frame()</code> appears to be intepreting <code>group_var</code> as continuous, I use the <code>change()</code> function from <code>mi</code> to reassign to <code>"un"</code> for "unordered categorical" and then proceed as above. </p>
  
  <pre><code>mydf2_mdf &lt;- change(mydf2_mdf, y = "group_var", what = "type", to = "un"  )
  
  # impute
  mydf2_imp &lt;- mi(mydf2_mdf)
  </code></pre>
  
  <p>Now, unless version 1.0 of <code>mi</code> has removed the functionality of previous versions (i.e. functionality available with <code>lmer.mi</code> and <code>glmer.mi</code>), I would assume that the addition of a random effect in the formula should point <code>pool()</code> to the appropriate <code>lme4</code> function. However, the initial error message suggests that this is not the case. </p>
  
  <pre><code># run models on imputed data and pool the results
  summary(pool(y ~ x1 + x2 + (1|group_var), mydf2_imp))
  ## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
  ## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
  ## Error in if (prior.scale[j] &lt; min.prior.scale) {: missing value where TRUE/FALSE needed
  </code></pre>
  
  <p>Following my warning message and extracting the integers out of my factor does get me an estimate, but the results suggest that <code>pool()</code> is still estimating a fixed-effect model with <code>bayesglm</code> and holding my attempted random-effect constant. </p>
  
  <pre><code>summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), mydf2_imp))
  
  ## 
  ## Call:
  ## pool(formula = y ~ x1 + x2 + (1 | as.numeric(as.character(group_var))), 
  ##     data = mydf2_imp)
  ## 
  ## Deviance Residuals: 
  ##      Min        1Q    Median        3Q       Max  
  ## -1.93633  -0.69923   0.01073   0.56752   2.12167  
  ## 
  ## Coefficients:
  ##                                               Estimate Std. Error  t value
  ## (Intercept)                                  1.383e-01  2.596e+02    0.001
  ## x1                                           1.995e+00  1.463e-02  136.288
  ## x2                                           2.000e+01  8.004e-03 2499.077
  ## 1 | as.numeric(as.character(group_var))TRUE -3.105e-08  2.596e+02    0.000
  ##                                             Pr(&gt;|t|)    
  ## (Intercept)                                        1    
  ## x1                                            &lt;2e-16 ***
  ## x2                                            &lt;2e-16 ***
  ## 1 | as.numeric(as.character(group_var))TRUE        1    
  ## ---
  ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ## 
  ## (Dispersion parameter for gaussian family taken to be 0.8586836)
  ## 
  ##     Null deviance: 5384205.2  on 1999  degrees of freedom
  ## Residual deviance:    1713.9  on 1996  degrees of freedom
  ## AIC: 5377
  ## 
  ## Number of Fisher Scoring iterations: 4
  </code></pre>
  
  <p>My questions are: </p>
  
  <ol>
  <li>Is it possible to easily generate pooled random effects estimates using <code>mi</code>?, and </li>
  <li>If yes, how?</li>
  </ol>