MICE: what does returned df mean? In MICE, the object returned by pool() has a component, df, which is included in the summary of the pooled analysis.  In my analysis I have about 55,000 cases, but the returned df is higher for most variables than 900,000.  What does this value of df refer to?
 A: For those without access to van Buuren's book, here is what the degrees of freedom mean. I am taking the steps pretty much directly from his book (which I think that you should buy, it is awesome!). We need to set up a lot of definitions first.
Our goal is to estimate a scientific estimand Q (this could be a regression coefficient, mean, etc.). In MICE, we do this by imputing $m$ datasets, running each analysis on each of the $m$ datasets, and then pooling the results.
To make this example easy to follow, let us suppose we are interested in a regression coefficient, and we impute $m=50$ times. We can find our average coefficient ($\bar{Q}$) by taking the average of each coefficient from our $l$ imputed datasets ($\hat{Q}_l$). So, we have that
$$\bar{Q}=\frac{1}{50}\sum_{l=1}^{m=50}\hat{Q}_l$$
Within each imputed dataset, we have a variance ($\bar{U}_l$) for the coefficient. $\bar{U}_l$ is really the var-covariance matrix of $\hat{Q}_l$, but since we have a scalar in this example, so too is $\bar{U}_l$. We may compute the average complete data var by doing
$$\bar{U}=\frac{1}{50}\sum_{l=1}^{m=50}\bar{U}_l$$
Between each dataset for the same coefficient, we also have variance. This is given by
$$B=\frac{1}{m-1}\sum_{l=1}^{m}(\hat{Q}_l-\bar{Q})(\hat{Q}_l-\bar{Q})'$$
Once again, we only have a scalar in this example, so it will just be the square of the estimated coefficient minus the average of that coefficient.
Lastly, we may define the total variance from imputation as 
$$T=\bar{U}+B+\frac{B}{M}$$
Which is the variation from within + variance between + simulation error
For a scalar Q, we may define the proportion of variance attributed to the missing data as
$$\lambda = \frac{B+\frac{B}{m}}{T}$$
To interpret this, if we see that a lot of our variation is due to our sampling, this, $\lambda$ will be close to 1. If there is no simulation of between variance, this will be 0.
We may define r as the relative increase in variance due to non response as
$$r=\frac{B+\frac{B}{M}}{\bar{U}}$$
This can also be written in terms of $\lambda$
$$r=\frac{\lambda}{1-\lambda}$$
Now, the old degrees of freedom, as defined by Rubin in 1987 is given by
$$\nu_{old}=(m-1)(1+\frac{1}{r^2})$$
Which can also be written as 
$$\frac{m-1}{\lambda^2}$$
However, this can be larger than the degrees of freedom, so in 1999, Barnard and Rubin proposed a new one. Letting
$$\nu_{com}=n-k$$
where n is the number of data points we have, and k is the number of parameters in a HYPOTHETICAL COMPLETE DATA.
We can then get a number called $\nu_{obs}$ by doing
$$\nu_{obs}=\frac{\nu_{com}+1}{\nu_{com}+3}\nu_{com}(1-\lambda)$$
Finally, we can get the degrees of freedom as
$$\nu=\frac{\nu_{old}\nu_{obs}}{\nu_{old}+\nu_{obs}}$$
So that is the theory behind what it ought to be.
However, you are getting a number that is above your total sample size. I believe this has to do with the mice package.
We can see function mice.df() will set the common degrees of freedom to 999999 if it can't find any. This is later called in the pool function on line 198. It tries to get the degrees of freedom residual from here, but it is returning a null. I guess the method that you are using doesn't have a way to get the $\nu_{com}$, so that is why this is happening.
So, I can offer few solutions, because I don't know the exacts of yours...


*

*calculate all of the degrees of freedom by yourself (this might
actually be easier) 

*Go into the code and change the default 999999
to something more reasonable. 

*Go into the df residuals code and
add code to allow you to get the df for your specific method.


Also, if @StefvanBuuren is on here, can you explain why the default for large sample/unknown $\nu_{com}$ is set at 999999?
Edits
I looked back at some code I did, and saw I was getting huge df's too. I am running a coxph model. But notice that we cannot get the df from here. because
class(coxph_obj$analyses[[1]])
[1] "coxph"

coxph_obj$analyses[[1]]$df.residual
NULL

df.residual(coxph_obj$analyses[[1]])
NULL

mn <- try(f <- fitted(coxph_obj$analyses[[1]]), silent = TRUE)
mn
NULL

So, it seems like even though the df is in fact available, the current code is not taking advantage of it, and thus defaults to using 9999999.
