The function empirically estimate the variance of the score functions.
The variance-covariance matrix consists of two parts: the additive
part and the part for the individual-specific environmental effect.
Other reasonable decompositions are possible. This program has the following improvement over "score.r":
1. It works with selected nuclear families
2. Trait data on parents (one parent or two parents), if available,
are utilized.
3. Besides a statistic assuming no locus-specific dominance effect,
it also computes a statistic that allows for such effect.
It computes two statistics instead of one.
Function "merge" is used to merge the IBD data for a pair with the
transformed trait data (i.e., $w_kw_l$).