REML is used for estimating the proportions of phenotypic variance explained by kinship matrices. When a kinship matrix reflects allelic correlations across genome-wide SNPs (and individuals are unrelated), then the proportion of phenotypic variance it explains is an estimate of the SNP heritability of the phenotype (the total variance explained by all SNPs). Although kinship matrices are usually calculated from many tens of thousands of predictors, it is equally valid to compute kinships across small numbers of SNPs, for example, if wishing to estimate the variance explained by a gene or small genomic region. However, in this case, instead of first explicitly calculating the kinship matrix for this region, then performing REML, it is usually much more efficient to do this on-the-fly as part of the REML algorithm.
The generalized REML solver contained in LDAK is optimized for the case when some kinship matrices correspond to small subsets of predictors (referred to as regions). Also, when the phenotype is binary, it can convert estimates of variance explained from the observed scale to the liability scale.
If your aim is prediction, then having estimated the variance components, you can then use BLUP to estimate predictor effect sizes. These effect size estimates can subsequently be used for predicting the phenotypic values of new individuals.
LDAK can also perform Haseman-Elston Regression, which essentially is a faster, but slightly less precise, alternative to REML, and PCGC Regression, which we recommend when wishing to estimate variance components for binary phenotypes.