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.

Lactation curve demonstration - Upload a Document to Scribd

## Days in milk
dimMin <- 1
dimMax <- 300
dimBy <- 0.5

## Days in milk
AliSch <- 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 matrix
X <- cbind(XInt, XDim)

## Parameter estimates - from one analyis in goats
int <- 2800.222
aliSch <-  c(-1678.68038, 219.29628, 53.04481, -55.26959)
b <- as.matrix(c(int, aliSch))

## Estimate
yHat <- (X %*% b) / 1000

## Plot

max(yHat)
## 2.47
yMax <- 3
yMin <- 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 povrsino

plot(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()