Chapter 10

---

FIGURES

Figures 10.1 and 10.2 use a function "number2words" written by John Fox (see below) and "capitalize" that you will need to copy to your workspace.

Figure 10.7 requires that you download the "car" package.  To do so, get into R and go to the top menu, then "packages", "install package(s)", pick a "mirror" close to you (e.g., "USA (IA)"), and then select "car" from the list.

---

number2words=function(x){
# Written by John Fox
# see http://socserv.socsci.mcmaster.ca/jfox/Courses/R-programming/Basic-programming.R

makeDigits=function(x) strsplit(as.character(x), "")[[1]]
makeNumber=function(x) as.numeric(paste(x, collapse=""))
ones=structure(c("zero", "one", "two", "three", "four", "five", "six",
"seven", "eight", "nine"), .Names = c("0", "1", "2", "3", "4",
"5", "6", "7", "8", "9"))
suffixes=c("thousand,", "million,", "billion,", "trillion,")
teens=structure(c("ten", "eleven", "twelve", "thirteen", "fourteen",
"fifteen", "sixteen", " seventeen", "eighteen", "nineteen"), .Names = c("0",
"1", "2", "3", "4", "5", "6", "7", "8", "9"))
tens=structure(c("twenty", "thirty", "forty", "fifty", "sixty", "seventy",
"eighty", "ninety"), .Names = c("2", "3", "4", "5", "6", "7",
"8", "9"))
trim=function(text){
gsub("(^\ *)|((\ *|-|,\ zero|-zero)$)",
"", text)
}

 negative <- x < 0
 x <- abs(x)
 digits <- makeDigits(x)
 nDigits <- length(digits)
 result <- if (nDigits == 1) as.vector(ones[digits])
 else if (nDigits == 2)
 if (x <= 19) as.vector(teens[digits[2]])
 else trim(paste(tens[digits[1]], "-", ones[digits[2]], sep=""))
 else if (nDigits == 3) {
 tail <- makeNumber(digits[2:3])
 if (tail == 0) paste(ones[digits[1]], "hundred")
 else trim(paste(ones[digits[1]], "hundred", number2words(tail)))
 }
 else {
 nSuffix <- ((nDigits + 2) %/% 3) - 1
 if (nSuffix > length(suffixes) || nDigits > 15)
 stop(paste(x, "is too large!"))
 pick <- 1:(nDigits - 3*nSuffix)
 trim(paste(number2words(makeNumber(digits[pick])),
 suffixes[nSuffix], number2words(makeNumber(digits[-pick]))))
 }
 if (negative) paste("minus", result) else result
 }

capitalize=function (string)
{
    capped <- grep("^[^A-Z]*$", string, perl = TRUE)
    substr(string[capped], 1, 1) <- toupper(substr(string[capped],
        1, 1))
    return(string)
}

---

Fig10.1=function (nloc=2,p=0.5,mean.x=72,sa=2)
{
# nloc is number of loci
# p is allele frequency at each locus
# mean.x is the mean of the measurement
# sa is the additive genetic standard deviation.  There is no environ var, so sa
#    is also the phenotypic standard deviation
nalls=2*nloc
F.cum=dbinom(0:nalls,nalls,p)
a=sa/sqrt(2*nloc*p*(1-p))
x=(0:nalls)*a
mu=as.vector(x%*%F.cum)
x=x+mean.x-mu
top=max(F.cum)
plot(x,F.cum,type='h',xlab='Stature (inches)',ylab='Frequency',
 main=noquote(paste(capitalize(number2words(nloc)),ifelse(nloc<2,'locus','loci'),'(B allele = ',round(a,2),')')),
 ylim=c(0,top),xlim=c(mean.x-3*sa,mean.x+3*sa),axes=F)
box();axis(1);axis(2)
}

---

Fig10.2=function (nloc=2,p=0.5,mean.x=72,s.dev=2)
{
nalls=2*nloc
Va=nalls*p*(1-p)
half.width=3*sqrt(Va)
x=0:nloc
F.cum=dbinom(x,nloc,p^2+2*p*(1-p))
mu=as.vector(x%*%F.cum)
x.cen=x-mu
x=x.cen*s.dev/sqrt(as.vector(x.cen^2%*%F.cum))+mean.x
a=x[2]-x[1]
print(x)
top=max(F.cum)
plot(x,F.cum,type='h',xlab='Stature (inches)',ylab='Frequency',
 main=noquote(paste(capitalize(number2words(nloc)),ifelse(nloc<2,'locus','loci'),
 '(B allele = ',round(a,2),')')),xlim=c(mean.x-3*s.dev,mean.x+3*s.dev),ylim=c(0,top))
}

---

Fig10.4=function (nids=500,ve=0,nclass=20,mu=72)
{
# nids is the number of individuals to simulate
# ve is the environmental variance
# nclass is the number of "bins" for the histogram
# mu is the population mean
nloc=10;p=0.5
nalls=2*nloc
dens=dbinom(0:nalls,nalls,p)
dens=round(dens*nids)
tot=sum(dens)
pheno=rep(62:82,dens)+rnorm(tot,0,sqrt(ve))
print(var(pheno))
junk=hist(pheno,nc=nclass,plot=F)
my.counts=junk$count
top=max(my.counts)
my.breaks=junk$breaks
plot(c(66.25,78),c(0,top),type='n',axes=F,xlab='Stature',ylab='')
axis(1)
y=c(0,rep(my.counts,each=2),0)
x=rep(junk$breaks,each=2)+.25
lines(x,y)
polygon(c(x,65.75),c(y,0),density=40)
rug(pheno)
}

---

Fig10.5=function (h2=1,Vp=4)
{
# h2 is the narrow-sense heritability
# Vp is the phenotypic variance
V=matrix(c(Vp,rep(.5*Vp*h2,2),Vp),nc=2)
x.y=matrix(rnorm(1000),nc=2)%*%chol(V)
x=x.y[,1]+72
y=x.y[,2]+72
sto=lm(y~x)
realized=round(2*sto$coef[2],2)
plot(x,y,xlim=c(66,78),ylim=c(66,78),
    main=noquote(paste('heritability = ',realized)),xlab='Mother',ylab='Daughter')
abline(sto,lwd=2)
abline(0,1,lty=2,lwd=2)
}

---

Fig10.7=function (h2=.35,Vp=4,Nids=500)
{
# h2 is the narrow-sense heritability
# Vp is the phenotypic variance
# Nids is the number of mid-parent/daughter pairs to simulate
library(car)
V=matrix(c(Vp,rep(Vp*h2,2),Vp),nc=2)
x.y=matrix(rnorm(2*Nids),nc=2)%*%chol(V)
x=x.y[,1]+72
y=x.y[,2]+72
sto=lm(y~x)
realized=round(sto$coef[2],2)
plot(x,y,xlim=c(66,78),ylim=c(66,78),xlab='Mid-Parent',ylab='Daughter',
    main=noquote(paste('Simulated under heritability = ',h2)))
abline(sto,lwd=2)
abline(0,1,lty=2,lwd=2)
dataEllipse(x,y,levels=(1:9)/10,plot.points=F,col=1,center.cex=0)
}