A Simple Model to Deal with Sire by Treatment Interactions When Sires Are Related J. L. FOULLEY
Station de Genetique Quantitative and Appliquee Institut National Recherche Agronomique Jouy-€n·Josas, 78350, France C. R. HENDERSON'
University of Illinois Urbana 61801
ABSTRACT
The assumption that var(hs) = la~s is not valid when sires are related. A more logical assumption is that the variance is [I 0 A] a~s' This result is derived from the use of an equivalent multiple-trait model. This procedure was extended to a more complicated model involving treatments and random herds within treatments. A small simulation trial demonstrated that a~s is underestimated if the la~s model is used when the true model should be [I 0 A] a~s'
2a;s' This is not reasonable if sires are related because one should expect related sires to have correlated differences. McLaren et al. (17) addressed this problem when treatments are types of dams. This method was applied by Cunningham and Magee (2) to analyze preweaning growth rate in beef cattle. The formulation of McLaren et al. (17) requires cov(s, st) =1= O. A simpler method not requiring this nonzero covariance is described in the present paper. This is accomplished by invoking a simple multiple-trait model as suggested by Falconer (6).
INTRODUCTION
The most widely used model when BLUP was introduced assumed that the interaction between herds (or herd-year-seasons) and sires is negligible. Studies of this interaction were usually based on a model, var(s) = la;, var(hs) = lah ' var(e) = la~, and these vectors uncorrelateJ [see, for example, (14)]. After the easy method for inversion of the numerator relationship matrix (12) was reported, estimates of 2 a hs were usually computed under the assumption that var(s) = Aat and, as before, var(hs) = la~s' Is var(hs) = la~s an appropriate model when sires are related? Hereafter we shall refer to herds (or herd-year-seasons) as treatments to show the relationship between interaction and a multiple-trait model. Consider a case with two treatments. The interactions in the la;, la;s model can be characterized as a set of sire differences, uncorrelated and with vanance
A GENERAL MULTIPLE·TRAIT MODEL
A general multiple-trait model applied to the treatment by sire interaction case follows.
where t refers to treatments, s to sires, and ~ to fixed effects other than t, for example, groups. We assume a single trait and a single record for progeny, and we order the records, progeny within sires within treatments. The sire vector can be partitioned as s' = [s~, ... , Sk] , where Sj refers to sire effects in the i th treatment. Similarly the error effects can be partitioned, e'
1989 J Dairy Sci 72:167-172
... , ek).
Then under the additive genetic model with no selection, G
Received January 7,1988. Accepted May 2, 1988. , Department of Anim aI Sciences.
= [e~,
=var(s) = Go 0
A
[2a]
where Go is one-fourth times the usual additive genetic covariance matrix and 0 denotes Kroneker product. 167
FOULLEV AND HENDERSON
168
R=var(e) =
[
Ia t 0
[2bl
Note that R is diagonal because a different set of progeny is recorded in each treatment, and is the residual variance for progeny within the ith treatment. The mixed model solution is illustrated with a simple set of data. The records are shown in Table 1.
at
Note that there is a nontrivial solution to sire values for treatments in which there were no progeny. This multiple-trait approach is practical if k, the number of treatments, is small, because then Go and can probably be estimated satisfactorily. The number of such parameters is k(k + 1)/2 + k. But if herds represent treatments, the number of parameters required would be very large, and there would be few observations per sire by treatment subclass. Consequently, a model with fewer parameters is needed. Such a model is described in the next section.
at
A SIMPLIFIED MULTIPLE-TRAIT MODEL
The assumption for a simplified model is that Go has only two distinct parameters, namely, v in all diagonals and c in all offdiagonals. Thus in matrix notation:
at = 50, a~ = 61, a~ = 70. 1.0
.5
.5
0
.5
1.0
.25
0
Go
= (v -
c)I + cJ
[3al
LetA = .5
.25
0
0
1.0
0
0
1.0
It was assumed that t is the only fixed effect. Then the mixed model solution is t' = [4.527, 5.762, 5.4891, s with rows representing treatments is:
[-M8
.316
-.161
-.060
.274
-.122
-.097
.283
-.404
-088] -.081
where J is a square matrix of 1's. That is, sire covariance is the same for all treatments, and the covariances for all pairs of treatments are identical. Also, it is assumed that all a~ equal 2 . 1 a e . That IS: R
= IU~
[3bl
Then only three parameters are required. An equivalent interaction model can be written 'h as, 2 u 2 ' an d u 2 as t h e parameters. For a Wit e ts definition of equivalent models see Henderson (13). The model is:
.113 Now s has as many elements as there are sires rather than k such subvectors. Var(s) = AU;, var(ts) = 1(9 Auis'
TABLE 1. Records.
ment
1 2 3
var(e)
Sires
Treat2
3,4,5 6,4
6, 7 8
9
4 4, 3
5,4,6,7 2, 3,4,6. 5
Journal of Dairy Science Vol. 72, No.1, 1989
2,8,9
= IU~
[4bl
and these vectors have null covariances with each other. The identity matrix in var(ts) has order equal to the number of treatments. For 2 '1ence Us2 = c an d U st ' eqUiva = v-c. Th e matrIx, I (9 A, is block diagonal with blocks of order
169
SIRE BY TREATMENT INTERACTIONS
equal to the number of sires. The mixed model equations are easy to write because
and A- 1 is simple to compute. These two simple equivalent models are illustrated with the same data as before and with:
and a~ == 40. Then the multiple-trait solution is t = [4.472, 5.944, 5.610] '.
s=
[~:::: -.199
.635
-.398
.525
-.439
.006
.592
-.699
.138
a;
-.122]
Under the equivalent interaction model
ais = 7 -
a;
==
4,
a;t
=
B oOB BOO] [:
a;t
where B == 3A. The 0 are matrices with the same order as A. Then the mixed model solution is:
i = same as before,
[023
.168
.012
-.057
.058
-.029
-.066
.124
-.289
TREATMENT AND HERDS WITHIN TREATMENTS BY SIRES
Another possible model is one with fixed treatments and random herds (or herd-yearseasons) within treatments. A model is:
S = [-.133, .467, -.410, .006]',
ts =
a;
a;.
4 == 3,
var(ts)
example, 51 + t531 = -.133 - .066 = -.199 = value for sire 1 in trait 3 of the equivalent multiple-trait solution. Caution must be used in applying the interaction model. Formally, the only condition to satisfy in applying [2a] is that Go, a variancecovariance matrix, is positive definite. This property is met when the eigenvalues of the matrix are positive. Due to the form of Go, those are v - c and v + (k - l)c [(22), p. 322]. Imposing these two quantities to be positive results in the following conditions: i) c < v for c > 0 or ii) lei < v/(k - 1) for c < O. Notice that strictly speaking the interaction model is defined only under i), because negative variances are by definition not possible. On the contrary, the parameter space for the simplified multiple-trait model is larger including also negative values of c satisfying ii). We illustrate this with the same data as before and with v = 5, c == -2. This is valid because 2 < 5/(3 - 1). Then for the interaction model = -2, = 7. Solutions are consistent with each other under both models provided the parameter space of the interaction model is extended to permit negative This is i == [4.589, 5.688, 5.491] , in both, and Sj + fS j are equal to the corresponding elements of s of the multipletrait model. There is a potential problem in < 0 in that we subtract a positive using number from the diagonals of the ts equations, and this could result in an ill-behaved coefficient matrix. If the idea of a negative variance is unacceptable, regard < 0 as merely a computing strategy, not an equivalent model. Note that c < 0 implies negative genetic correlations.
~12'] .132
Now if the prediction of the merit of the jth sire on the ith treatment is wanted, its BLUP is Sj + tSij, and these are equal to the intratreatment S of the multiple-trait model. For
y =
XI3 + Xtt + Zhh/t + Zss + Zhshs/t + e
[5a]
Var(h/ti) == la~(i)' i = 1, ... , k, that is a, different a~ within each treatment. Var(s)
=G o 0
var(hslt j ) ==
A,
10 Aa~S(i)
[5b]
Journal of Dairy Science Vol. 72, No.1, 1989
FOULLEV AND HENDERSON
170
that is, a different a~s for each treatment. Var(e) = diagonal matrix with elements Iai, Ia~, ... , Ia~. Note that the error vector has k subvectors, one for each treatment. We illustrate with data in Table 2. a~(i) = 5, 3. Thus, var(h) = diag(5, 5,3,3,3).
s ordered sires in treatments = -.275
.123
[ .016
.154
.101] -.113
h8ft ordered sires in herd in treatments = Go =
[: :l
[, :}
-.229
-.074
.193
-.082
.096
-.162
-.159
.194
.077
.154
-.164
.199
.099
-.170
0
.5
1
A=
0
2
2
2 (2) = 2. T h e a12=3 0, a2 = 2 4 ,ahs(l) = 4, a hs
first 16 diagonal elements of Rare 30 and the last 15 are 24. Then:
var(hslt) =
B
o
o
o
o
o
B
o
o
o
0
o
c
o
o
o
o
o
c
o
o
o
o
o
c
where B = 4A, C = 2A. It is assumed that X{3 = IJ.l.. Then a mixed model solution is:
fl = o.
i=
[5.495,5.2801',
hit = [-.012, .012, -.029, -.014, .0431',
TABLE 2. Records.
Treatments
If the correct variance is in fact I (8) Aa~ but ais is estimated under the model var(ts) = Ia~, are the estimators of the variances biased? The answer is "yes," the magnitude of bias depending in part upon how far A departs from I. We illustrate this by using MIVQUE in our example with true values a; = 6, a~ = 2, a~ = 20 and taking as guessed values for MIVQUE estimation 4, 3, 40 that were used in BLUP. Using the la~ model, the estimates have expected values, 6.155,1.474,19.986, but with the proper form of interaction the expected values are 6, 2, 20 as they should be. In general, a~ is underestimated with the Ia~ model. Because most estimates of a~s have been done using the I model, can we conclude that interaction may be more important than we have thought? If that is the case, perhaps we should include a herd X sire interaction term in sire evaluation. This in fact is done in the USDA program through use of a c 2 term. See Henderson (10) for a discussion of the relationship between c 2 and interaction.
Sires Herd
1
2
3
1
3,4,5,6 5,4
5.4,8 6, 5,8,9, 2
7,6
3,6,5
4 8, 5
6, 7 4,3,6 3
2 2
BIAS IN ESTIMATION OF VARIANCES USING IUts
3 4 5
8,9,2
Journal of Dairy Science Vol. 72, No.1, 1989
DISCUSSION
The general model shown in [1), [2al, and [2b 1 is obviously more flexible than the interaction model in [4al and [4b) , as already pointed out by Mallard et al. (15). The matrix Go has a general variance covariance structure that allows heterogeneity in variances and covariances among sire effects in different
171
SIRE BY TREATMENT INTERACTIONS
environments (4, 8) as well as possibly negative covariances contrarily to the form based on an intraclass correlation structure. These two models have been discussed by Dickerson (5) and Yamada (24) among others and, as stressed by Robertson (21), a~s may indicate changes in ranking as well as scaling effects. Therefore, caution is recommended by users in choosing the interaction model. A general discussion of the relative merits of the different statistical models available for analyzing sire (or genotype) by environment interactions is beyond the scope of this paper. Readers interested in that area are referred to the specialized literature [see (1, 3, 7, 9, 16, 18,19,20)]. A model for treatment (herds. etc.) x sire interaction that has variance la;s is not logical when sires are related. Very little additional computation is required to utilize a more logical assumption, that is, variance is [I (8) A] a;s' Use of the first model when sires are related yields predictors unbiased but with larger than necessary prediction error variances and to less progress through selection (11). Further, the estimates of a;s in supposed translation invariant quadratic unbiased methods are biased downward, leading possibly to concluding that a;s is unimportant when it in fact may be important. Finally the approach developed in this paper for a strict sire model can be easily employed in other models. The animal model proposed by Wiggans et al. (23) includes fixed herd-yearseason effects (h) and random herd x sire (hs), animal (a), and permanent environmental effects (p). According to formulae [4a] and [4b] , the variance covariance structure of random effects in this model may be written as: hs
o Aiaru,
0
0
0
a
0
A*ai
0
0
P
0
0
Ina p
0
e
0
0
0
INa~
var
The notation E!l denotes a set of direct sums. If all interaction~ between herds and sires are included var(hs) = I (8) Asa~s where As is the relationship among sires. The expression E!lAia~s results when rows and columns for misstng subclasses are deleted. Similar interaction effects (i.e., herd x dam) might be envisioned to take into account in this animal model preferential treatment applied within herds to daughter progeny born out of superior donors via multiple ovulation and egg transfer to recipient cows. REFERENCES
2
3
4
5
6 7
8
9
10 2
where Ai ,A * are the additive genetic relationships among sires present in herd i and among . d"d I . I 2 2 2 2 In IV) ua s respective y; a hs ' aa' a p , a e are variance components for hs, a, p, and e effects; n,N sizes of the p and e vectors, respectively.
11
12
13
Bucio Alanis, L. 1966. Environment and genotypeenvironmental components of variability. I. Inbred lines. Heredity 21 :387. Cunningham, B. E., and W. T. Magee. 1986. Evidence for a sire X dam type interaction in calves from Simmental sires mated to two crossbred dam types.]. Anim. Sci. (Suppl. 1):63:91. (Abstr;) Denis, J. B., J. P. Gouet, and J. Tranchefort. 1980. Methodes d'etude de la structure de l'interaction "genotype X milieu." Application a l'analyse d'essais de varietes de ble tendre. Pages 98- I 00 in Biometrie et genetique. Vol. 4. J. M. Legay, ]. P. Masson, and R. Tomassone, ed. Soc. Fr. Biomet., Inst. Natl. Rech. Agron., Dep_ Biomer., Paris, Fr. De Veer, ]. C., and L. D. Van Vleck. 1987. Genetic parameters for first lactation milk yields at three different levels of herd production. J. Dairy Sci. 70: 1434. Dickerson, G. E. 1962. Implications of geneticenvironmental interaction in animal breeding. Anim. Prod. 4:47. Falconer, D. S. 1952. The problem of environment and selection. Am. Nature 86:293. Freeman, G. H. 1973. Statistical methods for the analysis of genotype-environment interactions. Heredity 31: 339. Gianola, D. 1986. On selection criteria and estimation of parameters when the variance is heterogeneous. Theor. Appl. Genet. 72:671. Goffinet, B., and P. Vincourt. 1980. Prise en compte de I'interaction genotype X millieu en selection vegetale. Pages 90-97 in Biometrie et genetique. Vol. 4. J. M. Legay, ]. P. Masson, and R. Tomassone, ed. Soc. Fr. Biomet., Inst. Natl. Rech. Agron., Dep. Biomer., Paris, Fr. Henderson, C. R. 1973. Sire evaluation and genetic trends. Pages 10-41 in Proc. of the Animal Breeding Gener. Sympos. in Honor of Dr. J. L. Lush. Am. Soc. Anim. Sci.-Am. Dairy Sci. Assoc., Champaign,IL. Henderson, C. R. 1975. Best linear unbiased estimation and prediction under a selection model. Biometrics 31:423. Henderson, C. R. 1976. A simple method for computing the inverse of a numerator relationship matrix. Biometrics 32:69. Henderson, C. R. 1985. Equivalent linear models Journal of Dairy Science Vol. 72, No.1, 1989
172
FOULLEY AND HENDERSON
to reduce computations. J. Dairy Sci. 68:2267. 14 Hickman, C. G., and C. R. Henderson.-l956 .. Sire by herd interaction in production traits in dairy cattle. J. Dairy Sci. 39:1055. 15 Mallard, J., J. P. Masson, and M. Douaire. 1983. Interaction genotype X millieu et modele mixte. I. Modclisation. Genet. Sel. Evol. 15: 379. 16 Mather, K., and P.D.S. Valigari. 1974. Genotype X envitonment interactions. I. Regression of interaction on overall effect of the environment. Heredity 33 :43. 17 McLaren, D. G., D. S. Buchanan, and R. L. Hintz. 1985. Sire ranking based upon pure-bred versus crossbred progeny performance in swine. J. Anim. Sci. 60:902. 18 Meyer, K. 1987. Variances due to sire X herd interactions and common environmental covariances between half sisters for first lactation production. J. Dairy Sci. (Suppl. 1) 70: 126. (Abstr.)
J oumal of Dairy Science Vol. 72, No. I, 1989
19 Moav, R., and G. Wohlfarth. 1974. Magnification through competition oigenetic differences in yield capacity in carp. Heredity 33: 181. 20 Perkins, J. M., and J. L. Jinks. 1968. Environmental and genotype X environmental components of variability. III. Multiple lines and crosses. Heredity 23:339. 21 Robertson, A. 1959. The sampling variance of the genetic correlation coefficient. Biometrics 15 :469. 22 Searle, S. R. 1982. Matrix algebra useful for statistics. Wiley & Sons, New York, NY. 23 Wiggans, G. R., I. Misztal, and L. D. Van Vleck. 1987. Investigation on an animal model for na· tional genetic evaluation of milk yield. J. Dairy Sci. (SuppJ. 1) 70: 124. (Abstr.) 24 Yamada, Y. 1962. Genotype by environment interaction and genetic correlation of the same trait under different environments. Jpn. J. Genet. 37: 498.