List of animal names

List of animal names (Wikipedia) is a very usefull resource.

Relationships between distributions

John D Cook has created a nice graph of relationships between distributions.

Conjugate priors

John D Cook has posted a nice diagram of relationship between sampling distribution of the data i.e. the likelihood and the conjugate prior. He also posted a link to a compendium of conjugate priors by written by Fink.

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

Plot of normal distribution with shaded area

Theory of quantitative genetics much relies on the normal (Gaussian) distribution. Therefore, one would often like to plot it for presentations, class notes, etc. This is not that hard with excellent tools we have today. I played with plotting in R and came up with the following plot of normal distribution showing:

• the density of standardized breeding values with a mean of 100 units and a standard deviation of 12 unit and
  
• shaded area that roughly corresponds to upper 5 % of the distribution.

Here is the plot, followed by the R code.

Plot of a normal (Gaussian) distribution with shaded area - Upload a Document to Scribd

## Define the mean and the standard deviation
mu <- 100
sigma <- 12

## Define the % of upper are to shade
k <- 0.05

## Define the grid
x <- (seq(-5, 5, 0.01) * sigma) + mu

## Compute the density of normal distribution over the grid
y <- dnorm(x=x, mean=mu, sd=sigma)

## Compute the "threshold" for upper K % of the distribution
t <- qnorm(p=1-k, mean=mu, sd=sigma)

## Open PDF
## pdf(file="normal.pdf", width=5, height=3, pointsize=12)

## Open Windows metafile --> good for inclusion into MS Office documents
## win.metafile(filename="normal.wmf", width=5, height=3, pointsize=12)

par(mar=c(5, 4, 1, 1) + 0.1, # c(bottom, left, top, right)
    bty="l",
    pty="m")
plot(y ~ x, type="l", axes=FALSE, 
     xlab="Standardized breeding value", 
     ylab="Distribution", lwd=3)
## Mark the mean and the "threshold"
abline(v=c(mu, t), lwd=2, lty=2)
abline(h=0, lwd=3)
axis(1)
axis(2, labels=FALSE, tick=FALSE)
box()

## Add shaded polygon
testK <- x >= t
xK <- x[testK]
yK <- y[testK]
polygon(x=c(xK, rev(xK)), y=c(yK, rep(0, times=length(xK))), 
        col="black", border=NA)

## Add arrow --> this will need some manual work, i.e. modify the values
x2 <- t * 1.1
y[which(round(x) == round(x2))[1]]
max(y)
arrows(x0=(mu + sigma * 3.3), y0=y[which(round(x) == t)[1]], 
       x1=x2, y1=y[which(round(x) == x2)[1]], lwd=2)
text(x=142, y=y[which(round(x) == t)[1]] + 0.0015, labels=paste(t, "%"))
text(x=(mu + sigma * 3.3), y=0.02, labels="Good\n(+)")
text(x=(mu - sigma * 3.3), y=0.02, labels="Bad\n(-)")

## Close the device
## dev.off()

Runnning R2WinBUGS on Mac

New Mac computers/laptops (I do not know the details) are shipped with Intel processor. This brings the possibility to install a lot of UNIX/Linux software since Mac OS X is essentially a desktop Linux - I am again stating that I do not know the details about this since I can not afford a Mac.

Why is this important? WinBUGS is a great software for running Bayesian analyses. The name WinBUGS comes from two "parts": BUGS means Bayesian analysis Using Gibbs Sampling (though also other samplers are implemented), while Win corresponds to MS Windows. Therefore, one can natively run WinBUGS only on MS Windows. I personaly like the R2WinBUGS package, which allows one to run WinBUGS (and OpenBUGS, a new incarnation of BUGS software) directly from R. This approach is to a large extent promoted by Andrew Gelman in his "ARM book". Another possibilities are BRugs (using BUGS as directly linked with R - I can not find the link to BRugs package, since it is now in CRAN Extra repository) and rbugs (support for Linux systems with OpenBugs is emphasized) and the following two packages for interfacing JAGS (Just Another Gibbs Sampler): rjags and runjags.

It is quite some time since we managed to "tweak" R2WinBUGS so that it works also on Linux like systems that have installed wine emulator. One just needs to install R, wine and WinBUGS to "C:/Program Files/WinBUGS14" and the bugs() function will work out of the box on Linux!

This week as well as before we had a long discussion about running R2WinBUGS on new Macs with Intel processors. Rodney Sparapani found out that the following is needed:

schools.sim <-  bugs(data, inits, parameters, model.file,      n.chains=3, n.iter=5000,      bugs.directory="/Applications/WinBUGS14",      working.directory=".",      WINE="/Applications/Darwine/Wine.bundle/Contents/bin/wine",      WINEPATH="/Applications/Darwine/Wine.bundle/Contents/bin/winepath") 

Some notes are needed:

• bugs.directory="/Applications/WinBUGS14" - this seems to be a reasonable default for installation of WinBUGS on Mac, however, we will not "hardwire" this into the R2WinBUGS code since this place can vary a lot between different operation systems and locales


• working.directory="." - this is needed for the moment, since it seems that Mac launches wine emulator as another user, who does not have a permission to write to a temporary folder created by tempdir() by R; update on 2009-03-25, setting working.directory to "." will not be needed anymore in R2WinBUGS version 2.1-13 and later
  


• WINE="/Applications/Darwine/Wine.bundle/Contents/bin/wine" and WINEPATH="/Applications/Darwine/Wine.bundle/Contents/bin/winepath" are paths to wine and winepath binaries that we use in R2WinBUGS; unfortunatelly we can not find these programatically from R as it is the case in UNIX/Linux - if one cretaes links to /usr/bin or /usr/local/bin, there is no need to set this two arguments in bugs()

MathGL

MathGL (library for scientific graphics) has some very cool examples. I like this one:

Animal Breeding & Genomics 9

Animal Breeding & Genomics 9

Financiranje rejskega programa za lisasto pasmo govedi v Italiji

Te dni se mudim v Italiji na tečaju, ki ga vodi Ignacy Misztal. Tečaj sicer organizira ANAPRI - rejsko organizacija za lisasto pasmo govedi v Italiji. Po pogovoru z gostitelji sem izvedel, da se rejski program financira v grobem na sledeči način: 60 % iz proračuna, 15 % s strani rejcev in ostalo z "aktivnostjo" same zveze. Rejci morajo na kravo na leto plačati 11 EUR. S tem prispevajo del stroškov za izvajanje kontrole in ostale dejavnosti povezane s selekcijo.