The increasing use of genomic selection (GS) in plant and animal breeding programs has led to the development of software that fits models based on unique scenarios. Accordingly, several R packages have been developed for GS. The lmekin function from the coxme R package was one of the first functions implemented in R to fit models with random family effects using the pedigree–based relationship matrix. The function allows the user to provide the covariance structures for the random effects; thus, the GBLUP model can be fitted. This fitting process consists of replacing, in the traditional BLUP model, the additive relationship matrix derived from a pedigree by the additive relationship matrix derived from markers. Thus, the objective of this study was to employ the lmekin function in the context of genomic prediction by comparing the results of this function with those obtained using five R packages for GS: rrBLUP, BGLR, sommer, lme4qtl, and lme4GS. The comparisons were performed considering the computational times and predicted values for a wheat dataset and simulated big data. In addition, we implemented a 5-fold cross-validation scheme through considering the values predicted by the lmekin function for the wheat dataset. The results indicated that the lmekin function was effective in predicting genomic breeding values considering multiple random effects and relatively small sample sizes. The rrBLUP package processed the fastest for the scenario with only one genetic random effect, and the high temporal efficiency of the sommer package was confirmed for the scenario with more than one genetic random effect. Differences in computational times occurred because of the different algorithms implemented in the packages to estimate the variance components.