Traditional multiple imputation method
Multiple imputation is carried out by using the conditional distribution of the missing values given the observed values. Let Y_{i} = (Y^{T}_{mis}, Y^{T}_{obs})^{T} be the quantitative phenotype values for the i^{th} family with Y_{mis} representing the vector of missing values and Y_{obs} representing the vector of observed values. Similarly, we can partition the mean and covariance matrix. Thus, for the polygenic model (including an overall mean and a gender effect), the distribution of Y_{mis} given Y_{obs} is
Y_{mis}  Y_{obs} ~ N(μ_{1}, Σ_{1}),
with μ_{1} = (μ + X_{1}β) + (σ^{2} D_{12})(σ^{2} D_{22} + τ^{2} I_{22})^{1}(Y_{obs}  (μ + X_{2} β))
and Σ_{1} = (σ^{2} D_{11} + τ^{2} I_{11})  (σ^{2} D_{12})(σ^{2}D_{22} + τ^{2} I_{22})^{1} (σ^{2} D_{21}).
The imputation is carried out by generating missing values from the conditional multivariate normal distribution, taking into account the family structure. This imputation is completed m times to produce m complete data sets to analyze. From these m analyses, the final point estimate would be the mean of the m estimates. The computation of the standard error can be done by first computing the betweenimputation variation, B_{
m
}, and the withinimputation variation, W_{
m
}. Then, the total variability is T_{
m
}= W_{
m
}+ B_{
m
}(m + 1)/m. Confidence intervals for parameters of interest use the tdistribution with degrees of freedom (m  1) (1 + 1/ (m + 1)*W_{
m
}/B_{
m
})^{2} [6].
The one problem with this imputation scheme is that values for μ, β, σ^{2}, and τ^{2} are required for the imputation. To address this issue, we ran the analysis on the observed data and used the resulting estimates to create k complete data sets, from which k sets of point estimates are found. Then, the average of the k sets of estimates is used for the m multiple imputations. Point estimates were found using the EM algorithm program PolyEM with k = 5 and m = 25 [5, 7]. Variance estimates were found using the Fisher Information matrix and the large sample properties of maximum likelihood estimates (MLEs) [4, 7, 8].
Multiple imputation via Gibbs sampler
Another approach for the imputation of missing data is through a Bayesian analysis via a Gibbs sampler. The Gibbs sampler is a particular Markov chain algorithm that is useful when working with high dimensional problems. In addition to the traditional use of the Gibbs sampler, an imputation step can also be added to impute missing values. A multiple imputation scheme can be implemented by having an imputation step at the beginning of the Gibbs sampler. For each iteration, a new imputation is done, giving multiple imputations for each missing value [9–15]. For a detailed proof of the data augmentation methodology, see Tanner and Wong [12].
A Bayesian polygenic model with noninformative prior distributions for the i^{th} family is Y_{i} = X_{i}β + a_{i} + ε_{i}, where Y_{i} is a vector containing the individual responses in family i, X_{i} is a design matrix containing covariate information, β is a vector of covariate effects, a_{i} is a vector of random family effects where a_{i} ~ MVN(0, σ^{2}D_{i}) and D_{i} is a known coefficient of relationship matrix, and ε_{i} ~ MVN(0, τ^{2}I). Noninformative priors were then placed on all other parameters in the model, i.e., p(β) proportional to 1, p(σ^{2}) proportional to 1/σ^{2}, and p(τ^{2}) proportional to 1/τ^{2} [6].
The steps for implementing the Gibbs sampler with an imputation step for the Bayesian polygenic model follow.

1.
Set starting values for β^{(0)}, σ^{2(0)}, τ^{2(0)}, a_{i}^{(0)} for all i = 1,..., k, and set m = 1 (iteration).

2.
If y_{ij} is missing, impute y_{ij} by simulating an observation from N(X_{ij} β ^{(m1)} + a_{ij}^{(m1)}, τ^{2(m1)}).

3.
Generate β^{(m)} from MVN((X^{T}X)^{1}X^{T}(y  a^{(m1)}), τ^{2(m1)}(X^{T}X)^{1}).

4.
Generate τ^{2(m)} from INGAM(N/2, 1/2 Σ(y_{i}X_{i}β^{(m)}  a_{i}^{(m1)})^{T}(y_{i}X_{i}β^{(m)}  a_{i}^{(m1)})).

5.
Generate σ^{2(m)} from INGAM(N/2, 1/2 Σ a_{i}^{(m1)T} D_{i}^{1}a_{i}^{(m1)}).

6.
Generate a_{i}^{(m)} from MVN(μ_{a}^{(m)}, V_{a}^{(m)}), where μ_{a}^{(m)} = (1/σ^{2(m)} * D_{i}^{1} + 1/τ^{2(m)}* I_{i})^{1} * (1/ τ^{2(m)} * I_{i} (y_{i}  X_{i}β^{(m)})) and V_{a}^{(m)} = (1/σ^{2(m)} * D_{i}^{1} + 1/τ^{2(m)} * I_{i})^{1}.

7.
After one iteration of the algorithm, you have (β^{(m)}, σ^{2(m)}, τ^{2(m)}, a^{(m)}). Set m = m + 1 and repeat steps 1 through 6.
Approximate (1  α)% posterior confidence intervals are found by taking the α/2 and the 1α/2 percentiles of the simulated marginal posterior distributions for parameters of interest. The simulated posterior distributions will reflect the uncertainty in the estimation of the parameter along with the uncertainty due to the imputation of the missing data.