Case study: Lactation curves for camels
Data provided by Fred Aloo (Institut 480, Uni Hohenheim, Prof. Valle-Zarate, June 2003)
Short description of objective of study and research question
Rendille pastoralists in northern Kenya group their camels (all belonging to the Rendille camel breed) into different performance and adaptation types. There is
-the Dabach type, said to have the highest milk yield during rainy season but loses body condition rapidly during dry season, with consequent drop in milk yield,
-the Godan type with has a medium but stable milk yield during both rainy and dry season and
-the Coitte type with a low milk yield but good body condition throughout.
-The Aitimaso type is said to be in-between the Dabach and the Godan type having a higher milk yield then the Godan but a better drought tolerance then the Dabach type.
The study aims at estimating the milk yield of these different Rendille camel types in order to assess whether the indigenous classification can be trusted and hence used as part of on-farm performance estimation.
Performance estimation of livestock kept under their natural production conditions is specifically difficult in remote areas under adverse conditions, such as the marginal dryland areas. Performance and adaptation of such livestock species and breeds can also hardly be determined on research stations, because of strong genotype-environment interactions. Therefore the development and use of methods for on-farm performance testing are necessary when aiming at assessing and improving livestock production systems in marginal areas. These methods seem to be more practical when the livestock keepers are involved in the performance recording. The results of the study will therefore contribute to the assessment of coherence and complementarities of indigenous and scientific knowledge and further be used for the fine-tuning of on-farm performance testing methods suitable for the area.
- Cross-sectional study, i.e., survey was conducted at one point in time (2,5 months data collection period from March to May 2003, this was towards the end of the dry season)
- Camels are classified into four different types: Dabach, Aitimaso, Godan and Coitte
- For each type, camels were surveyed which were in lactation months 1, 2, ..., 12. For each lactation month and camel type, there are replicate camels. Each camel was surveyed only once, i.e., there are no repeated measurements.
- Lactation curves were to be fitted to each type
- Questions: What model to fit to the lactation curve? Are regression curves different among types and if so, in what respect?
- The following model has been shown to fit well to lactation curves of camels:
where x is the number of lactation months. This model is known as the incomplete gamma model due to its relation to the incomplete gamma function and as Wood's function (Wood, 1967, 1972).
Suggestion for analysis
The model can be linearized by taking logarithms:
where
This can be fitted by multiple linear regression of y' = log(y) on z1 and z2.
To compare lactation curves, it is useful to fit curves simultaneously for all four types. The joint model can be written
where
= logged milk yield of the j-th camel of the i-th type at the t-th month
Since there are replicate camels per month and type, the lack-of-fit can be tested by adding a lack-of-fit effect:
(1)
where
it = lack-of-fit effect for the i-th type and t-th month
The ANOVA for model (1), using sequential SS (Type I SS in SAS), yields
Source DF Type I SS Mean Square F Value Pr > F
Type 3 5.21121050 1.73707017 8.86 <.0001
Lactmonth*Type 4 14.11557956 3.52889489 17.99 <.0001
LOG_Lactmonth*Type 4 4.19206207 1.04801552 5.34 0.0004
Type*lackfit 33 10.05184539 0.30460138 1.55 0.0314
The lack-of-fit is significant, so one may consider modifications of model (1). Specifically, we tried the following three models:
(2)
where
(a) ,
(b) and
(c) .
The ANOVAs are:
(a)
Source DF Type I SS Mean Square F Value Pr > F
Type 3 5.21121050 1.73707017 8.82 <.0001
Lactmonth*Type 4 14.11557956 3.52889489 17.93 <.0001
LOG_Lactmonth*Type 4 4.19206207 1.04801552 5.32 0.0004
Lactmon*Lactmon*Type 4 4.30056889 1.07514222 5.46 0.0003
Type*lackfit 30 5.75127650 0.19170922 0.97 0.5092
R2 = 0.31 (without lack-of-fit term)
(b)
Source DF Type I SS Mean Square F Value Pr > F
Type 3 5.21121050 1.73707017 8.86 <.0001
Lactmonth*Type 4 14.11557956 3.52889489 17.99 <.0001
LOG_Lactmonth*Type 4 4.19206207 1.04801552 5.34 0.0004
LOG_Lac*LOG_Lac*Type 4 2.91910479 0.72977620 3.72 0.0057
Type*lackfit 29 7.13274064 0.24595657 1.25 0.1788
R2 = 0.29 (without lack-of-fit term)
(c)
Source DF Type I SS Mean Square F Value Pr > F
Type 3 5.21121050 1.73707017 8.82 <.0001
Lactmonth*Type 4 14.11557956 3.52889489 17.93 <.0001
LOG_Lactmonth*Type 4 4.19206207 1.04801552 5.32 0.0004
Lactmon*LOG_Lac*Type 4 3.64047308 0.91011827 4.62 0.0012
Type*lackfit 30 6.41137231 0.21371241 1.09 0.3522
R2 = 0.30 (without lack-of-fit term)
While model (a) had the largest R2, it showed a slight upward twist for the highest lactation month (not shown), which is not plausible biologically. This problem did not occur with models (b) and (c), and among these, (c) fit slightly better. Thus, all subsequent analyses were based on model (c).
We tested interaction of type-by-covariate based on the model
(3)
where . Type-specific effects were partitioned into a main effect and an interaction effect as follows:
The ANOVA is as follows:
Source DF Type I SS Mean Square F Value Pr > F
Type 3 5.21121050 1.73707017 8.75 <.0001
Lactmonth 1 13.55550713 13.55550713 68.32 <.0001
Lactmonth*Type 3 0.56007242 0.18669081 0.94 0.4211
LOG_Lactmonth 1 3.81456270 3.81456270 19.22 <.0001
LOG_Lactmonth*Type 3 0.37749937 0.12583312 0.63 0.5935
Lactmonth*LOG_Lactmo 1 2.29964353 2.29964353 11.59 0.0007
Lactmon*LOG_Lac*Type 3 1.34082956 0.44694319 2.25 0.0822
All three interaction terms (bi, ci, di) were non-significant. The results suggest that curves are identical for the four groups, except for the intercept (ai significant), i.e., the curves are parallel. Thus, we fit the reduced model
(4)
For illustration, the fitted model for Type 1 is shown in Figs. 1 (log-scale) and Fig. 2 (original scale; obtained by back-transformation of model fitted on log-scale). A residual plot across all types is inconspicuous and not suggestive of variance heterogeneity (Fig. 3).
Fig. 1: Fitted model (4c) for Type=1 on log-scale, with 95% prediction limits.
Fig. 2: Model for Type 1 on original scale based on back-transformation of model (4c) fitted to log-transformed data (y0.29).
Fig. 3: Plot of studentized residual versus predicted value for fitted model (4c). All types.
Finally, we compute least squares means based on the reduced model (2). A letter display is obtained by the method of Piepho (2003).
1
Tab. 1: Means for types based on model (4) using log-transformed data. Means in a column followed by the same letter are not significantly different.
Scale
Typelog§original$
10.601a1.82a
20.476b1.61b
30.309c1.36c
40.486abc1.63abc
§ least squares mean
$ back-transformed least squares mean; this is a median estimate
A final improvement
Statistical inference (tests of effects) is expected to be fairly robust to nonnormality. Predictions of individual observations, however, are more sensitive. A Q-Q-plot of studentized residuals based on model (4) with log-transformed data shows some departure from normality (Fig. 4). Thus, we used the Box-Cox transformation, which yielded y0.29 as the optimal normalising transformation. Normality of studentized residuals was satisfactory on this scale (Fig. 5). Thus, we re-did all analyses on this new transformed scale. ANOVA results were the same as on the log-scale (ANOVA tables not shown), including the comparison of means (Tab. 2). Fitted curves and prediction intervals are also quite similar, though the intervals have become a bit more narrow.
Fig. 4: Normal Q-Q-plot based on model (4c) using log-transformed data.
Fig. 5: Normal Q-Q-plot based on model (4c) using power-transformed data (y0.29).
Fig. 6: Fitted model (4c) for Type=1 on power-transformed-scale (y0.29), with 95% prediction limits.
Fig. 7: Fitted model for Type=1 on original scale, with 95% prediction limits. Model on ori-ginal scale based on back-transformation of model (4c) fitted to power-transformed data (y0.29).
Tab. 2: Means for types based on model (4) using power-transformed data (y0.29). Means in a column followed by the same letter are not significantly different.
Scale
Typey0.29§original$
11.202a1.88a
21.161b1.67b
31.107c1.42c
41.163abc1.68abc
§ least squares mean
$ back-transformed least squares mean; this is a median estimate
Other types of model
Inverse polynomials are a flexible class of models (McCullagh and Nelder, 1989, Generalized linear models. Chapman and Hall, London) which may be used to fit lactation curves, e.g.,
and
These are easily linearized and so can be used in the same way as the incomplete gamma model. One may also consider using the power transformation to gain more flexibility. Yet another model (Emmans and Fisher, 1986):
Rerefences
Wood PDP 1967 Algebraic model of the lactation curve in cattle. Nature 216: 164-165.
Wood PDP 1972 A note on seasonal fluctuations in milk production. Anim. Prod. 15: 89-92.
SAS code
SAS macros for multiple comparisons and for Box-Cox-transformation are available at Details on their use can be found in the lecture notes "Gemischte Modelle in den Life Sciences" (see homepage under "Veranstaltungen").
/*lack of fit for regression model*/
procglmdata=s;
Class type lackfit;
Model LOG_yield=type type*lactmonth
type*LOG_lactmonth
type*lactmonth*LOG_lactmonth
type*lackfit;
run;
quit;
/*fit model with interaction*/
procglmdata=s;
Class type;
Model LOG_yield=type lactmonth type*lactmonth
log_lactmonth type*LOG_lactmonth
lactmonth*log_lactmonth type*lactmonth*log_lactmonth;
run;
/*interactions not significant, so reduce model -> curves run parallel on log-scale/not on original scale*/
procglmdata=s;
Class type;
Model LOG_yield=type lactmonth log_lactmonth lactmonth*log_lactmonth;
outputout=res student=student predicted=pred_log_scale LCL=LCL_log_scale UCL=UCL_log_scale;
run;
procsortdata=res out=res;
by type lactmonth;
data res;
set res;
pred_original_scale=exp(pred_log_scale);
LCL_original_scale=exp(LCL_log_scale);
UCL_original_scale=exp(UCL_log_scale);
filename fileref 'd:\hpp\beratung\aloo\fig.cgm';
goptionsreset=all dev=win target=cgmof97L gsfname=fileref keymap=winansi ftext=hwcgm001 gsflen=8092
gsfmode=replace hsize=17 cm vsize=12 cm htext=1.1;
axis1label=none offset=(0 cm, 0 cm) w=10minor=none major=(w=10);
axis2label=none offset=(0 cm, 0 cm) minor=none w=10major=(w=10);
/* plot fitted curves on log-scale*/
symbol1i=spline value=none color=black w=2;
symbol2i=none value=dot color=black h=1.5;
symbol3i=spline value=none color=black w=2line=3;
symbol4i=spline value=none color=black w=2line=3;
procgplot;
plot pred_log_scale*lactmonth log_yield*lactmonth
LCL_log_scale*lactmonth UCL_log_scale*lactmonth/overlayhaxis=axis1 vaxis=axis2 nolegend;
by type;
run;
quit;
/* plot fitted curves on original scale*/
procgplot;
plot pred_original_scale*lactmonth yield*lactmonth
LCL_original_scale*lactmonth UCL_original_scale*lactmonth/overlayhaxis=axis1 vaxis=axis2 nolegend;
by type;
run;
quit;
/*Residual plots*/
symboli=none value=dot;
procgplot;
plot student*pred_log_scale/overlayhaxis=axis1 vaxis=axis2;
run;
quit;
procunivariatedata=res;
qqplot student/normal nohlabel novlabel;
run;
/*pairwise comparisons*/
odsoutput diffs=diffs;
odsoutput lsmeans=lsmeans;
procmixeddata=s;
Class type;
Model LOG_yield=type lactmonth log_lactmonth lactmonth*log_lactmonth;
lsmeans type/pdiff;
%mult(trt=type);
run;
%boxcox(phimin=0.2,phimax=0.4, steps=20, model=type lactmonth log_lactmonth lactmonth*log_lactmonth, class=type, data=s, response=yield);
Data
Type Lactmonth Yield
220.76
243.00
233.00
331.52
222.00
340.65
143.13
222.08
145.06
244.00
343.46
153.60
353.20
281.00
181.20
241.40
144.00
342.80
243.80
143.60
310.90
243.00
441.60
142.80
111.70
143.90
342.20
342.40
141.80
241.80
252.20
142.40
242.40
351.40
351.60
142.00
141.96
452.20
352.10
251.00
3120.70
2120.76
2120.52
143.42
142.00
441.60
441.90
3111.40
1111.04
1121.00
341.50
242.50
440.50
441.70
132.00
133.80
132.40
231.44
442.80
144.20
242.40
242.80
151.30
141.30
311.60
232.00
152.84
121.90
2120.78
131.22
131.80
444.40
1120.80
131.00
322.50
252.00
231.60
2121.50
2123.60
134.20
252.20
1121.00
152.80
4121.46
352.40
4121.00
331.42
231.56
133.50
252.20
143.20
1101.92
291.50
151.40
272.80
131.84
3121.60
342.80
244.80
224.00
1111.00
471.20
231.24
442.50
352.80
342.70
352.50
452.90
261.00
351.20
152.60
151.46
452.30
341.20
2120.60
2121.82
2120.84
432.42
4120.68
1110.60
1110.90
381.00
162.68
433.20
451.60
1103.00
2121.70
151.88
232.00
173.00
172.00
441.64
462.00
152.42
1121.00
212.08
151.40
442.48
442.20
451.60
452.80
152.00
452.60
452.40
450.80
442.40
341.40
164.40
1122.60
132.60
151.84
152.40
251.80
242.60
142.20
152.40
252.00
251.00
221.00
1121.60
351.40
261.70
461.60
161.80
360.90
451.50
151.00
161.80
232.40
362.40
462.00
462.00
161.76
370.84
472.40
171.60
163.20
152.00
152.00
252.20
340.80
152.80
361.50
371.00
121.00
122.40
162.20
132.62
361.24
331.80
132.40
281.60
181.80
191.40
191.06
181.22
180.90
1112.40
191.00
281.00
370.76
2110.30
1111.00
490.80
181.40
191.80
491.00
1101.00
270.84
2100.82
390.20
380.22
390.50
171.60
251.52
271.60
252.00
280.40
262.00
171.00
250.90
1100.84
161.00
2101.30
270.90
171.22
271.40
271.00
271.00
261.40
111.40
2111.42
2111.00
181.40
281.80
281.60
181.58
282.00
192.60
2111.00
1100.60
291.40
291.60
390.80
1101.60
261.80
271.60
192.00
2100.60
371.20
351.00
291.70
160.72
191.40
2101.00
2111.44
2120.50
181.00
370.40
181.60
2101.40
2101.00
1101.40
111.60
211.00
271.60
171.80
3121.20
2111.20
1101.32
1111.38
391.30
2111.60
371.80
212.20
112.40
210.80
310.90
211.00
2120.80
2111.30
3110.76
261.26
361.22
1122.60
3122.00
262.40
273.20
361.00
152.80
174.00
251.60
131.20
2102.98
252.80
262.40
351.60
354.40
153.60
371.30
262.00
361.80
361.50
3121.60
3121.80
251.00
122.90
331.60
263.00
262.40
261.80
264.46
264.00
362.00
361.70
361.60
351.60
381.00
371.04
3100.90
371.30
3111.30
321.60
381.30
390.80
3101.10
391.60
382.00
3111.50
1111.60
311.00
152.60
133.80
444.00
211.50
452.26
1