## 2008-10-10

### Lactation curve demonstration

Here is another demo plot. This time a lactation curve is showed. Lactation curve shows how the amount of produced milk is chaning over the lactation - from parturition to drying off. In this demo I used the so called Ali&Schaeffer function as applied to some goat data. Some special "things" added to the plot:
• a line to show the time of weaning the offspring and
• a line to show the time of drying off.
Here is the plot, followed by the R code.

Ali and Schaeffer. 1987. Accounting for covariances among test day milk yields in dairy cows. Can. J. Anim. Sci. v67. 637-644.

`## Days in milkdimMin <- 1dimMax <- 300dimBy <- 0.5## Days in milkAliSch <- function(x, c=305){  cbind(x / c,        (x / c) * (x / c),        log(c / x),        (log(c / x)) * (log(c / x)))}x <- seq(from=dimMin, to=dimMax, by=dimBy)XInt <- rep(1, times=length(x))XDim <- AliSch(x=x)## Design matrixX <- cbind(XInt, XDim)## Parameter estimates - from one analyis in goatsint <- 2800.222aliSch <-  c(-1678.68038, 219.29628, 53.04481, -55.26959)b <- as.matrix(c(int, aliSch))## EstimateyHat <- (X %*% b) / 1000## Plotmax(yHat)## 2.47yMax <- 3yMin <- 1.5## Open PDF## pdf(file="lact_curve.pdf", width=5, height=3.2, pointsize=11)par(bty="l",               # brez okvirja    mar=c(3.2, 3, .5, .5), # robovi od osi    mgp=c(2, 0.7, 0),      # robovi za stevilke na osi    pty="m")               # maksimalno izkoristi povrsinoplot(x=x, y=yHat, type="l", lwd=2, ylim=c(yMin, yMax),     axes=FALSE,     xlab="Days in milk", ylab="Daily milk yield (kg)")box()text(x=0, y=yMin - 0.15, xpd=TRUE,     labels="Parturition")xP <- c(80, 120, 145, 180, 210, 240, 265)points(x=xP, y=yHat[which(x %in% xP)], cex=2, pch=19)points(x=xP, y=yHat[which(x %in% xP)], cex=3, pch=21)text(x=200, y=yMax - 0.8, labels="Test-day records")     abline(v=75, lwd=2, lty=2, col="gray")arrows(x0=75 - 30, x1=75 + 30, y0=yMin + 0.005, y1=yMin + 0.005,       lty=1, lwd=2, code=3, length=0.15, col="gray")text(x=75, y=yMin - 0.15, xpd=TRUE,     labels="Weaning")text(x=30, y=yMax - 0.2,      labels="Suckling period")abline(v=275, lwd=2, lty=2, col="gray")arrows(x0=275 - 30, x1=275 + 30, y0=yMin + 0.005, y1=yMin + 0.005,       lty=1, lwd=2, code=3, length=0.15, col="gray")text(x=275, y=yMin - 0.15, xpd=TRUE,     labels="Drying off")text(x=175, y=yMax - 0.2,     labels="Milking period")## Close the plot## dev.off()`