Shackelford, Harris, and Konigsberg 2011

Calculations and Figures from Shackelford, Harris, and Konigsberg (2011)

First, “R” is your friend, so:

Converting the Moorees, Fanning, and Hunt graphs to useable parameters

Note: this section is here just to show you where all those numbers in various tables come from. You don't need to run this stuff if you are happy with just taking the tables at face value, but if not...

First, you will need the digitized ages which are contained in a *.csv file at:

MFH.csv

If you have the *.csv extension associated with a spreadsheet program, then this will open in a spreadsheet. If so, save the file (with the extension *.csv) in your R working directory. And then start R and and save the data to an object, say:

MFH=read.csv('MFH.csv')[,-1]

You will need this "junk" later to draw the first three Figures from the article.

The columns in this file are as follows:

Sex - "M" or "F"
Tooth - c, m1, m2, UI1, UI2, LI1, LI2, P3, P4, M1, M2, or M3. Mesial versus distal roots for molars are indicated with a lowercase "m" or "d," so m2m is the mesial root of the 2nd lower deciduous molar.
Stage - as in Moorees, Fanning and Hunt
L2SD - the digitized age minus two standard deviations below the mean
L1SD - the digitized age minus one standard deviation below the mean
Mean - digitized mean age
U1SD - the digitized age one standard deviation above the mean
U2SD - the digitized age two standard deviations above the mean
Mu - Mean age from Harris and Buck (2002)
SD - standard deviation from Harris and Buck (2002)
WithoutMean - the least squares estimate of the mean transition age (plus 0.75 years for conception age) on a natural log scale, where the estimate is made ignoring the "mean" point (variable #6)
WithMean - same as variable #11, but including the mean (variable #6)

An "NA" in a column means that a value is "not available." Note that columns 11 and 12 are what we are trying to estimate, although to make life simpler they are included in the file. These estimates were obtained using the following script.

get.estimates = function ()
{
#####################################
points2=function (x=c(NA,NA,NA,0.00324,0.0714))
{
Get.points=function(mu)
{
points=rep(0,5)
points[3]=mu
points[2]=points[3]-.0967
points[1]=points[2]-.0967
points[4]=points[3]+.0967
points[5]=points[4]+.0967
return(points)
}
Find.Min=function(mu){
try=Get.points(mu)
try=(try-x)[-(1:3)]
return(as.numeric(try%*%try))
}
x=log(x+0.75)
estimu=optimize(Find.Min,lower=-0.29,upper=3.3)$minimum
Points=Get.points(estimu)
Resid=x-Points
return(estimu)
}
#############################
points4=function (x=c(NA,0.0783,0.1568,0.2438,0.3513))
{
Get.points=function(mu)
{
points=rep(0,5)
points[3]=mu
points[2]=points[3]-.0967
points[1]=points[2]-.0967
points[4]=points[3]+.0967
points[5]=points[4]+.0967
return(points)
}
Find.Min=function(mu){
try=Get.points(mu)
try=(try-x)[-c(1,3)]
return(as.numeric(try%*%try))
}
x=log(x+0.75)
estimu=optimize(Find.Min,lower=-0.29,upper=3.3)$minimum
Points=Get.points(estimu)
Resid=x-Points
return(estimu)
}
####################################
points4.b=function (x=c(NA,0.0783,0.1568,0.2438,0.3513))
{
Get.points=function(mu)
{
points=rep(0,5)
points[3]=mu
points[2]=points[3]-.0967
points[1]=points[2]-.0967
points[4]=points[3]+.0967
points[5]=points[4]+.0967
return(points)
}
Find.Min=function(mu){
try=Get.points(mu)
try=(try-x)[-1]
return(as.numeric(try%*%try))
}
x=log(x+0.75)
estimu=optimize(Find.Min,lower=-0.29,upper=3.3)$minimum
return(estimu)
}
##################################
points5=function (x=c(.3958,.5232,.6507,.7912,.9546))
{
Get.points=function(mu)
{
points=rep(0,5)
points[3]=mu
points[2]=points[3]-.0967
points[1]=points[2]-.0967
points[4]=points[3]+.0967
points[5]=points[4]+.0967
return(points)
}
Find.Min=function(mu){
try=Get.points(mu)
try=(try-x)[-3]
return(as.numeric(try%*%try))
}
x=log(x+0.75)
estimu=optimize(Find.Min,lower=-0.29,upper=3.3)$minimum
Points=Get.points(estimu)
Resid=x-Points
return(estimu)
}
#################################
points5.b=function (x=c(.3958,.5232,.6507,.7912,.9546))
{
Get.points=function(mu)
{
points=rep(0,5)
points[3]=mu
points[2]=points[3]-.0967
points[1]=points[2]-.0967
points[4]=points[3]+.0967
points[5]=points[4]+.0967
return(points)
}
Find.Min=function(mu){
try=Get.points(mu)
try=(try-x)
return(as.numeric(try%*%try))
}
x=log(x+0.75)
estimu=optimize(Find.Min,lower=-0.29,upper=3.3)$minimum
Points=Get.points(estimu)
Resid=x-Points
return(estimu)
}

#########################
teeth=MFH[,4:8]
n=NROW(teeth)
y=matrix(rep(NA,2*n),nc=2)
for(i in 1:n){
temp=as.numeric(teeth[i,])
npnts=sum(!is.na(temp))
if(npnts==5){
y[i,1]=points5(temp)
y[i,2]=points5.b(temp)
}
if(npnts==2){
y[i,1]=points2(temp)
y[i,2]=y[i,1]
}
if(npnts==4){
y[i,1]=points4(temp)
y[i,2]=points4.b(temp)
}
}
colnames(y)=c('WithoutMean','WithMean')
return(y)
}

get.estimates()

Getting the data

Mark and copy the following to your "RGui" to get the data into an object called (drum roll, wait... for ... it): tooth.scores

tooth.scores=structure(list(ID = structure(c(1L, 2L, 7L, 8L, 9L, 10L, 11L,
12L, 13L, 14L, 15L, 16L, 17L, 18L, 20L, 22L, 23L, 24L, 26L, 27L,
28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 42L, 41L,
43L, 45L, 47L, 44L, 48L, 49L, 50L, 51L, 52L, 53L, 55L, 56L, 58L,
59L, 57L, 60L, 61L, 62L, 63L, 65L, 66L, 67L, 70L, 72L, 68L, 73L,
75L, 76L, 77L, 78L, 79L, 81L, 82L, 19L, 21L, 25L, 46L, 39L, 40L,
64L, 69L, 71L, 80L, 54L, 3L, 4L, 5L, 6L, 74L), .Label = c("Abri Pataud 26.230 B",
"Abri Pataud 26.236", "Anderson et al. (1976) \"A\"", "Anderson et al. (1976) \"B\"",
"Anderson et al. (1976) \"C\"", "Anderson et al. (1976) \"D\"",
"Arene Candide (young prince #?)", "Arene Candide VB (listed as IIB in folder)",
"Arene Candide VIII", "Arene Candide XI", "Atapuerca 2", "Badegoule #3 Mand #H-412 Coll en Côte.",
"Barma Grande adolescent (3 or 4?)", "Barma Grande Child", "Brillenhöhle",
"Chateauneuf-sur-Charente 2", "Combe Grenal I", "Ehringsdorf Juvenile",
"Engis 2 (Smith et al. 2010)", "Fontéchevade 1957-53", "Gibraltar 2 (Smith et al. 2010)",
"Grotte du Rousset", "Hohlenstein Infant #1", "Hortus 2", "Irhoud 3 (Smith et al. 2010)",
"Isturitz III (1937) 1950-7", "Isturitz III 1937 – 5 -1", "Isturitz III 1950-6",
"Kostenki XVIII", "Kostenski XV", "Krapina #45 and 45.1 (Max A)",
"Krapina #46 (Max B)", "Krapina #47 (Max C)", "Krapina #51 (Mand A)",
"Krapina #52 (Mand B)", "Krapina #53 (Mand C)", "Krapina #54 (Mand D)",
"Krapina #55 (Mand E)", "Krapina Max B (Smith et al. 2010)",
"Krapina Max C (Smith et al. 2010)", "La Chaud 3", "La Ferrassie 8",
"La Genière #3 1926", "La Madeleine", "La Quina H18 (75372)",
"La Quina H18 (Smith et al. 2010)", "La Quina Q-761 (child)",
"Laugerie- Basse 1", "Laugerie- Basse 2", "Laugerie- Basse 3",
"Laugerie- Basse 6", "Laugerie-Basse teen", "Le Fate II", "Le Moustier (Smith et al. 2010)",
"Le Placard 56029", "Le Placard 61397", "Le Placard 61401- 61397",
"Le Placard 61401 (D.G. #31 or 32?) ", "Le Placard 61401 (D.G. #40) (defaced label) ",
"Mal’ta child (older)", "Mal’ta child (Younger)", "Mas D’Azil",
"Montgaudier 3", "Obi Rakhmat (Smith et al. 2010)", "Paglicci adolescent (2)",
"Parpallo I", "Pech de l'Aze", "Qafzeh 10", "Qafzeh 10 (Smith et al. 2010)",
"Qafzeh 11", "Qafzeh 15 (Smith et al. 2010)", "Qafzeh 4", "Roc de Marsal",
"Roc de Marsal (Bayle et al. 2009b)", "Rochereil 1945-18", "Saint Germain La Rivière 1970-8 B3",
"Saint Germain La Rivière 1970-8 B4", "Saint Germain La Rivière 1970-8 B5",
"Saint Germain La Rivière 1970-8 B6 and B7", "Scladina (Smith et al. 2010)",
"Skhul I", "Skhul X"), class = "factor"), dc = c(NA, NA, NA,
NA, 15L, 11L, NA, NA, NA, NA, NA, 15L, 15L, NA, NA, NA, 10L,
NA, NA, NA, NA, NA, NA, NA, NA, NA, 14L, 16L, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
10L, NA, NA, NA, NA, NA, 10L, NA, NA, NA, 11L, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, 11L), dm1 = c(NA, 15L, NA, 14L, 15L, 14L, NA, NA, NA,
NA, 15L, NA, 15L, NA, NA, NA, 10L, NA, 13L, NA, NA, 17L, 15L,
16L, NA, NA, NA, NA, NA, NA, NA, 11L, NA, 14L, NA, NA, NA, 15L,
14L, 15L, 14L, NA, NA, NA, NA, NA, NA, 15L, 11L, 7L, NA, NA,
NA, NA, 11L, NA, NA, NA, 14L, 12L, NA, NA, 14L, 15L, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 15L
), dm2 = c(NA, 15L, NA, 13L, 14L, 10L, NA, 14L, NA, NA, 15L,
NA, 15L, NA, 14L, NA, 9L, NA, 11L, 14L, 14L, 16L, 14L, 15L, NA,
NA, NA, NA, NA, NA, NA, 10L, NA, 13L, NA, NA, NA, 14L, 14L, 14L,
13L, 16L, NA, 15L, NA, NA, 14L, 15L, 9L, 6L, NA, 16L, NA, NA,
10L, NA, 15L, NA, 12L, 10L, NA, 15L, 14L, NA, 14L, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 11L), UI1 = c(NA,
NA, 14L, 5L, 6L, 5L, NA, NA, 14L, 14L, NA, 9L, NA, NA, NA, NA,
NA, NA, NA, NA, NA, 14L, 9L, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 14L, NA, NA, NA, NA,
NA, NA, NA, NA, 7L, NA, 14L, NA, NA, NA, NA, NA, 5L, 5L, NA,
NA, NA, NA, NA, NA, 6L, 7L, NA, 9L, 9L, NA, NA, 7L, 11L, NA,
14L, NA, NA, NA, NA, 5L), UI2 = c(NA, 6L, 14L, NA, 6L, 4L, NA,
NA, 14L, 14L, NA, 7L, NA, NA, NA, NA, NA, NA, NA, NA, NA, 14L,
7L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, 5L, 14L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 14L,
NA, NA, NA, NA, NA, 4L, 4L, NA, NA, NA, NA, NA, NA, 6L, 7L, NA,
7L, 9L, NA, 11L, 6L, 10L, 14L, 14L, NA, NA, NA, NA, 5L), LI1 = c(NA,
NA, NA, 7L, 10L, NA, NA, NA, 14L, NA, 13L, NA, NA, 14L, NA, NA,
NA, NA, NA, NA, 9L, 14L, 10L, NA, NA, NA, 10L, NA, NA, NA, 14L,
NA, NA, 7L, NA, NA, NA, NA, NA, 10L, NA, NA, NA, 11L, NA, NA,
NA, NA, NA, NA, NA, NA, 14L, NA, NA, 14L, 11L, NA, 6L, NA, NA,
NA, 9L, NA, NA, 9L, NA, 7L, 11L, NA, NA, NA, NA, 7L, 11L, 14L,
14L, 9L, 14L, 7L, 13L, 5L), LI2 = c(NA, NA, NA, 6L, 9L, NA, NA,
6L, 14L, NA, NA, NA, NA, 14L, NA, NA, 5L, NA, 6L, 9L, 7L, 14L,
9L, NA, NA, NA, 10L, 11L, 14L, 14L, 14L, NA, 10L, 6L, NA, 10L,
NA, 9L, NA, 9L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
14L, NA, NA, 14L, 10L, NA, 5L, NA, NA, NA, 9L, NA, NA, 7L, NA,
7L, 10L, NA, NA, NA, NA, 7L, 10L, NA, 14L, 7L, 13L, 6L, 11L,
5L), C = c(NA, 7L, NA, 5L, 9L, 4L, NA, 5L, NA, 11L, 9L, 7L, 9L,
12L, NA, NA, NA, 11L, 4L, 7L, 6L, NA, 7L, 10L, NA, NA, 9L, 10L,
NA, 14L, 13L, NA, 7L, 5L, 7L, 7L, NA, 6L, NA, 6L, 4L, 11L, 10L,
9L, NA, 10L, NA, NA, NA, NA, 6L, 10L, 13L, NA, NA, 12L, 9L, NA,
4L, NA, 9L, 10L, 6L, 7L, NA, 5L, 5L, 6L, 9L, 7L, 7L, NA, 9L,
6L, 9L, 13L, 14L, 6L, 11L, 6L, 10L, 5L), P3 = c(NA, 7L, NA, NA,
7L, NA, NA, 4L, 14L, 11L, NA, 6L, 6L, 14L, 5L, NA, NA, 9L, NA,
6L, 6L, 11L, 6L, NA, 6L, NA, 7L, 9L, NA, 14L, 12L, NA, 6L, NA,
6L, 6L, NA, 6L, 4L, 6L, 3L, 11L, 9L, 7L, NA, 11L, 4L, 9L, NA,
NA, 5L, 9L, NA, NA, NA, 11L, 6L, NA, 2L, NA, 9L, 7L, 6L, 6L,
NA, NA, 5L, 5L, 9L, 7L, 7L, NA, 9L, 6L, 9L, 11L, 14L, 6L, 11L,
4L, 7L, 3L), P4 = c(3L, 6L, NA, NA, 6L, NA, NA, 3L, 14L, 11L,
6L, 5L, NA, NA, 4L, 4L, NA, 9L, NA, 5L, 5L, 10L, 5L, NA, 5L,
NA, NA, 9L, 10L, 14L, 13L, NA, 5L, NA, 5L, 5L, NA, 5L, NA, 5L,
NA, 10L, 7L, 6L, NA, 10L, 2L, 9L, NA, NA, NA, 7L, NA, NA, NA,
9L, 6L, 5L, NA, NA, 6L, 7L, NA, NA, NA, NA, NA, 4L, 7L, 6L, 6L,
9L, 9L, 5L, 7L, 10L, 14L, 4L, 11L, 2L, 7L, NA), M1 = c(NA, NA,
14L, 6L, 11L, 5L, NA, 8L, 14L, NA, 12L, 10L, 12L, NA, 9L, NA,
4L, 13L, NA, NA, 8L, 14L, 10L, 12L, 11L, NA, NA, 14L, 14L, 14L,
14L, 4L, 11L, 7L, 12L, NA, 14L, NA, 7L, NA, 6L, 14L, 13L, NA,
14L, NA, NA, NA, 5L, NA, NA, 13L, 14L, NA, 5L, 14L, 12L, 10L,
6L, 5L, 11L, 13L, 10L, 10L, 6L, NA, 8L, 9L, 12L, 11L, 11L, 14L,
NA, 10L, 12L, 14L, 14L, 9L, 14L, 9L, 12L, 6L), M2 = c(5L, NA,
14L, NA, 5L, NA, NA, NA, 13L, NA, NA, 5L, NA, 12L, 4L, NA, NA,
NA, NA, 6L, NA, 10L, NA, NA, NA, 10L, NA, NA, 11L, 12L, 13L,
NA, 5L, NA, 4L, NA, 14L, 5L, NA, NA, NA, 10L, 6L, 6L, 14L, NA,
NA, NA, NA, NA, NA, 7L, 12L, 14L, NA, 10L, 5L, 5L, NA, NA, NA,
6L, NA, NA, NA, NA, NA, 4L, 6L, 5L, 6L, 10L, 9L, 5L, 7L, 10L,
14L, 3L, 10L, 1L, 6L, NA), M3 = c(NA, NA, 11L, NA, NA, NA, 7L,
NA, 9L, NA, NA, NA, NA, 6L, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, 6L, NA, NA, NA, NA, NA, NA, NA, 10L, NA,
NA, NA, NA, 5L, NA, NA, 9L, NA, NA, NA, NA, NA, NA, NA, 8L, 7L,
NA, 5L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, 5L, 10L, NA, 1L, NA, NA, NA), Neander = c(FALSE,
FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE,
FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE,
FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE,
FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUE,
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE,
FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE,
FALSE, FALSE, FALSE, FALSE), Obs = c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L)), .Names = c("ID",
"dc", "dm1", "dm2", "UI1", "UI2", "LI1", "LI2", "C", "P3", "P4",
"M1", "M2", "M3", "Neander", "Obs"), class = "data.frame", row.names = c(NA,
-82L))

Getting the log scale conception-corrected mean transition ages (averaged across the sexes) and a lot of other housekeeping type stuff

Mark and copy the following to your "RGui" to get the Moorrees, Fanning, and Hunt parameters into an object called "MFH2", the mean log age within stage as "MUS", the "names" for stages into "scores," and "obj" which holds the maximum of the "objective" functions (likelihood values at the max likelihoods to scale the "plot.teeth" plots).

MFH2=structure(list(dc = c(-Inf, -0.385121087118164, -0.103360231065519,
0.0129993730872631, 0.200274234705511, 0.356347028223297, 0.463099858231227,
0.463099858231227, 0.569237250992917, 0.719588896525789, 0.940011978072852,
1.01336672773577, 1.18074017139880, 1.32534017350979, 1.82706248424488,
2.14852015367735, 2.30218770825353, Inf), dm1 = c(-Inf, -Inf,
-0.41055416872582, -0.0896085389559922, -0.0180138761278178,
0.122529585839533, 0.275467644616887, 0.302796188944659, 0.366450336422772,
0.503761988964755, 0.641972706579564, 0.711770606997633, 0.839163052228666,
0.962451404328723, 1.77697341356485, 2.09883576693258, 2.28758325145222,
Inf), dm2 = c(-Inf, -0.403379589886261, -0.0775551507179944,
0.0082416082845772, 0.202331283637113, 0.374367800593468, 0.51670013401181,
0.546027962318958, 0.727751949817333, 0.837900035950872, 0.969226311960765,
1.01774965010600, 1.14644839560940, 1.30678830067758, 1.95783771480066,
2.21621364098845, 2.39469630477986, Inf), UI1 = c(-Inf, -Inf,
-Inf, -Inf, -Inf, 1.76478075807267, 1.76478075807267, 1.76478075807267,
1.93066074274181, 2.01252753930847, 2.14594441751995, 2.21275468384625,
2.26818791637201, 2.26818792558235, Inf, Inf, Inf, Inf), UI2 = c(-Inf,
-Inf, -Inf, -Inf, -Inf, 1.87687798018485, 1.87687798018485, 1.87687798018485,
2.01367429576234, 2.09020211360232, 2.22445622326653, 2.31090685887088,
2.33241512201773, 2.33241508287378, Inf, Inf, Inf, Inf), LI1 = c(-Inf,
-Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, 1.66392006371608, 1.7857915914593,
1.94759753864067, 2.02193103433356, 2.11676348953656, 2.15571473737620,
Inf, Inf, Inf, Inf), LI2 = c(-Inf, -Inf, -Inf, -Inf, -Inf, -Inf,
-Inf, -Inf, 1.75496726704851, 1.92154219262198, 2.05230891973101,
2.14127921894058, 2.20535175303785, 2.26284850706191, Inf, Inf,
Inf, Inf), C = c(0.196163993672336, 0.391282611409868, 0.676501006212301,
1.02027416756059, 1.30061917094731, 1.55511945037896, 1.70271434432993,
1.70271434432993, 1.84030362416174, 2.12282361258594, 2.27590737434281,
2.33226195635993, 2.4525587691839, 2.56154365380084, Inf, Inf,
Inf, Inf), P3 = c(0.895753497247203, 1.08781845426543, 1.27107148030256,
1.45142782205652, 1.62001687637545, 1.75710414200012, 1.87584995422872,
1.87584995422872, 1.99917332634546, 2.20356758346562, 2.32953717847462,
2.38652598683234, 2.50390562006323, 2.60226560339250, Inf, Inf,
Inf, Inf), P4 = c(1.26665436224112, 1.43859336011361, 1.58798597895525,
1.69468916522858, 1.80545373858716, 1.92729476168307, 2.01738689877063,
2.01738689877063, 2.12285858885350, 2.28375125480169, 2.40843327876547,
2.47390989471367, 2.57101033406347, 2.68877133295038, Inf, Inf,
Inf, Inf), M1 = c(-Inf, -0.0800009566037727, 0.260904657373109,
0.550270526009671, 0.776433441019678, 1.05042616475525, 1.22306306235392,
1.43069030091526, 1.68202604785279, 1.76315195082442, 1.81941548862697,
1.87426062710447, 2.00481032175298, 2.19486922288907, Inf, Inf,
Inf, Inf), M2 = c(1.44608548846335, 1.5048249016104, 1.64467883374839,
1.72845177242483, 1.82136187992417, 1.94596970999109, 2.04573757673928,
2.15040361532396, 2.29050395400051, 2.36344617251816, 2.42905678280068,
2.46738287648498, 2.54515666917048, 2.69232326617605, Inf, Inf,
Inf, Inf), M3 = c(2.31645286993104, 2.36383348575149, 2.41644325389745,
2.46985058456016, 2.51407682458934, 2.55588483449181, 2.60509750674667,
2.66188345735770, 2.73988441618072, 2.78633685995019, 2.82786979152452,
2.8566610254016, 2.92021776431994, 3.01381543524097, Inf, Inf,
Inf, Inf)), .Names = c("dc", "dm1", "dm2", "UI1", "UI2", "LI1",
"LI2", "C", "P3", "P4", "M1", "M2", "M3"), row.names = c("C.i",
"C.co", "C.oc", "Cr.5", "Cr.75", "Cr.c", "R.i", "Cl.i", "R.25",
"R.5", "R.75", "R.c", "A.5", "A.c", "Res.25", "Res.5", "Res.75",
""), class = "data.frame")
MUS=structure(c(NA, -0.244210093821536, -0.0451769371865635, 0.106637284838301,
0.278310430906716, 0.409724607287536, 0.516132515849055, NA,
0.644412871432217, 0.829786015294346, 0.976701013285029, 1.09705847524341,
1.25304640260997, 1.57620319089673, 1.98774620927372, 2.22534718951923,
NA, NA, NA, -0.250073955898128, -0.0538364364400467, 0.0522684704241992,
0.198972699368432, 0.289126372078225, 0.334616596142382, 0.435105436959326,
0.572853634545099, 0.676882185465514, 0.775466604750228, 0.900794762740116,
1.36971171099702, 1.93791463860420, 2.19321607967430, NA, NA,
-0.240479407469982, -0.0346766011366329, 0.105287363407646, 0.288349499038704,
0.445531984634842, 0.5313513740979, 0.636877922910098, 0.782825191479169,
0.903545172946022, 0.99350187591832, 1.08209771996008, 1.22663583644906,
1.63229328460603, 2.08702720801156, 2.30545493343871, NA, NA,
NA, NA, NA, NA, 1.84770713224503, 1.84770713224503, NA, 1.97158684486789,
2.07923739085107, 2.17934761978367, 2.24045129088881, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, 1.94526460264653, 1.94526460264653,
NA, 2.05197823260549, 2.15734184638886, 2.26768201069163, 2.32167460564372,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1.72484696640806,
1.86669107211658, 1.98478105852691, 2.06934624235780, 2.13625120997274,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1.83827044062701,
1.98692565340228, 2.09677987103023, 2.17331678388273, 2.23409983219926,
NA, NA, NA, NA, 0.293750636861609, 0.533935974492477, 0.848389043842237,
1.16044939116755, 1.42786805409702, 1.62892263014194, 1.77153255579987,
NA, 1.98154358122653, 2.1993813642209, 2.30408740110175, 2.39242013760075,
2.50705375915771, NA, NA, NA, NA, 0.991787831104052, 1.17943490678436,
1.36125672523502, 1.5357317965779, 1.68855009730478, 1.81645460239049,
1.9375013302448, NA, 2.10139379772819, 2.26655267660225, 2.35808332879738,
2.44524978387237, 2.55308407896496, NA, NA, NA, NA, 1.35262416027652,
1.51328382709518, 1.64137968431062, 1.7500322821843, 1.86635867357081,
1.97233472716624, 2.07012120914100, NA, 2.20333679827773, 2.34609037210263,
2.44120832697799, 2.52246107363013, 2.62993698432684, NA, NA,
NA, NA, NA, 0.0904515639428928, 0.405617443979164, 0.663354260558872,
0.913424679693308, 1.13674449012028, 1.32687531105396, 1.55640252841111,
1.72258375282689, 1.79128258527421, 1.8468391383155, 1.93952495713922,
2.09986880643202, NA, NA, NA, NA, 1.47545125195769, 1.57472141905819,
1.68654712995320, 1.77490590832603, 1.88366736576709, 1.99580790320727,
2.09805212129396, 2.22043008531651, 2.32701353218885, 2.39625886007653,
2.44817725534085, 2.50627258338588, 2.61875249532713, NA, NA,
NA, NA, 2.34013892296218, 2.39015391370519, 2.44318922176223,
2.49198095429838, 2.53498027463377, 2.58052727403965, 2.63345973880698,
2.70085214960755, 2.76310684288337, 2.80710339184673, 2.84226609793515,
2.88843735520561, 2.96699886017695, NA, NA, NA, NA), .Dim = c(17L,
13L), .Dimnames = list(c("C.i", "C.co", "C.oc", "Cr.5", "Cr.75",
"Cr.c", "R.i", "Cl.i", "R.25", "R.5", "R.75", "R.c", "A.5", "A.c",
"Res.25", "Res.5", "row17"), c("dc", "dm1", "dm2", "UI1", "UI2",
"LI1", "LI2", "C", "P3", "P4", "M1", "M2", "M3")))
scores = c("C.i", "C.co", "C.oc", "Cr.5", "Cr.75", "Cr.c", "R.i", "Cl.i",
"R.25", "R.5", "R.75", "R.c", "A.5", "A.c", "Res.25", "Res.5",
"Res.75", "")
obj=structure(c(0.164645682728558, 0.899805458667659, 0.476377944956582,
0.702186459070747, 0.610830839190904, 0.441056898489228, 0.438759998183299,
0.438759998183299, 0.592678016549752, 0.784797440327601, 0.311055525598271,
0.645424255087385, 0.574002865155155, 1.04264536703987, 0.951037367800852,
0.603251011876375, 1.05263157894737, 0.106983302325428, 0.106983302325428,
0.950477424479285, 0.303930832342928, 0.560578693935637, 0.600937181772089,
0.118275522304579, 0.271496947011247, 0.549736107208176, 0.552765483523776,
0.296632903223354, 0.515656432319588, 0.501210545113548, 1.05260484251372,
0.951478237481924, 0.706173350611486, 1.05263157894737, 0.121513242605929,
0.955713841062925, 0.360691203578481, 0.720388706753091, 0.659203428691458,
0.566524120990583, 0.126864391815836, 0.686893327179715, 0.453655450510473,
0.529314877750986, 0.208514199691008, 0.520212555293192, 0.624091652828244,
1.05182888104397, 0.861474807272684, 0.677769827799991, 1.05263157894737,
1.05263157894737, 1.05263157894737, 1.05263157894737, 1.05263157894737,
1.05263157894737, 0.640949634036484, 0.640949634036484, 0.640949634036484,
0.34515621246639, 0.536497524792233, 0.284444003994213, 0.237453576573366,
3.99942135549742e-08, 1.05263157894737, 0, 0, 0, 1.05263157894737,
1.05263157894737, 1.05263157894737, 1.05263157894737, 1.05263157894737,
0.547995170044931, 0.547995170044931, 0.547995170044931, 0.323837464742566,
0.539359042893524, 0.363262574262047, 0.0932035581993022, -7.1645152099498e-167,
1.05263157894737, 0, 0, 0, 1.05263157894737, 1.05263157894737,
1.05263157894737, 1.05263157894737, 1.05263157894737, 1.05263157894737,
1.05263157894737, 1.05263157894737, 0.496177632957919, 0.628592287282343,
0.315006818240173, 0.395872077937051, 0.168002310435728, 1.052623,
0, 0, 0, 1.05263157894737, 1.05263157894737, 1.05263157894737,
1.05263157894737, 1.05263157894737, 1.05263157894737, 1.05263157894737,
1.05263157894737, 0.643035477268336, 0.527383460805062, 0.373134610440663,
0.273217331667863, 0.246040199683812, 1.052623, 0, 0, 0, 0.723082030657873,
0.904934157264072, 0.973149519047962, 0.897666339216396, 0.854486616674623,
0.583779895251521, 0.550672622840476, 0.550672622840476, 0.900943053328292,
0.601399885800024, 0.241290755594014, 0.490556538821616, 0.449353215443344,
1.052623, 0, 0, 0, 0.715046739332682, 0.691147727293633, 0.683061117454533,
0.64904440142063, 0.548978474863221, 0.484992191262669, 0.501334649869696,
0.501334649869696, 0.746711480302694, 0.510670698143865, 0.243929016894758,
0.480068145102355, 0.4093936349583, 1.052623, 0, 0, 0, 0.658918251102838,
0.589593596115461, 0.440871734546410, 0.455929609150581, 0.496069012419981,
0.377511075976835, 0.436271038558216, 0.436271038558216, 0.625791899557453,
0.506138346469598, 0.278981775538870, 0.404579764543550, 0.481444789432897,
1.052623, 0, 0, 0, 0, 0.970549505605535, 0.910910462052479, 0.797597605200458,
0.887786372704815, 0.660956320047176, 0.754672724997884, 0.848642228403703,
0.34221229526318, 0.240911816542405, 0.235001581059735, 0.526633297407221,
0.709698938256659, 1.052623, 0, 0, 0, 0.251198028106785, 0.558276260913264,
0.352708199899313, 0.388452077745746, 0.505876583055556, 0.414755559656101,
0.433252539962073, 0.559099730545655, 0.309387796192153, 0.279531173153755,
0.165341309821737, 0.32883415731014, 0.582388701799474, 1.052623,
0, 0, 0, 0.203702321573151, 0.225662301266584, 0.228997897145796,
0.190383895628813, 0.180139762050102, 0.211413462422148, 0.243085454368063,
0.329743744827291, 0.199788753008826, 0.178972687132590, 0.124560481010746,
0.271096232547244, 0.391110081204736, 1.052623, 0, 0, 0), .Dim = c(17L,
13L), .Dimnames = list(NULL, c("col1", "col2", "col3", "col4",
"col5", "col6", "col7", "col8", "col9", "col10", "col11", "col12",
"col13")))

Getting the "get.age" script (for Figures 4, 5, 7, 9, and 13)

Mark and copy the following to your "RGui" to get the script "get.age" into your workspace. The calling parameters are as follows.

1. id - Which case to do from the data file "tooth.scores." The default (=1) is to do the first case
2. lo - lower limit for integration. Integration is used to find the within-tooth variance (see eqn 5 in the manuscript). The default value of 0.75 years starts the integration at birth (0.75 years of conception time).
3. hi - upper limit for integration, with a default of 25.75 years (so less 0.75 years conception time, equals 25 years)
4. left - For clipping of the plot, giving the "left" (lowest) age plotted. The default is 0, so birth.
5. right - Highest age plotted (default is 25 years).
6. def.int - The interval for plotting. The default of 0.1 is generally adequate, but if you want a smoother plot then decrease this parameter (and wait longer for your plot!).
7. put.points - A logical value (i.e., either T or F) indicating whether the line fitted to the total (across teeth) likelihood should be marked with filled circles. The default is to not mark the fitted line with points, and to just show the line as a dashed line (with the total likelihood line as a solid line). But if all is well in the universe, then the fitted and the actual lines will be superimposed and you would need the filled circles to see the fitted line.
8. put.L - If points are plotted on the fitted line (see parameter #7), then this parameter is used to indicate at what lower (left) age the points should start being plotted.
9. put.R - as with parameter #8, but for the highest (right) age.
10. drop - To indicate any teeth that shouldn't be used in the function. For example, if you wanted to exclude the third molar then drop=13, as the third molar is the 13th tooth in the tooth.scores file. To drop all molars you could do drop=c(11,12,13) or drop11:13
11. retain - The opposite of "drop."
12. Known.age - If the age is known, enter the known age so that it will be plotted as a vertical line.
13. teeths - a different way to get the tooths scores for a case analyzed.
14. lab - a case label to use if you used the 13th argument ("teeths").
15. anote - an annotation to a case from the "tooth.scores" file.

Output:
The graph will show a solid line which is the total (across tooth) likelihood, a dashed line as the fitted line to the likelihood (note that this may inditinguishable, so see parameter #7), and another dashed line which is the final age density accounting for both within-tooth variance (the line fitted to the total likelihood line) and between-tooth variance. Numeric output is both "printed" to the screen and returned as a list containing $mu, $within, $between, and $p.seq. "p.seq" calculation is explained in the paper.

Examples:
Once you have copied "get.age" you can paste these pairs of lines (or just the second of each pair as the first is just a comment) to your Gui to reproduce plots from the manuscript.

# Figure 4. "Subject A" (note: retain=10 to only use the P4, anote='Only P4' to annotate, Known=5.0 as this is a five year old).
get.age(78,left=3,right=7,put.p=T,retain=10,Known=5.0,anote='Only P4')

# Figure 5. Example of poor fit to log-normal, dm2, A.c (note: using teeths=c(c(rep(NA,2),14) to 'hard enter' a 2nd deciduous molar w/ root complete)
get.age(id=NA,teeths=c(rep(NA,2),14),left=1,right=9,lab=expression(paste('Deciduous ',2^{nd},' molar, Apex complete')))

# Figure 7 "Subject A"
get.age(78,put.points=T,left=3.5,right=6.5,put.L=4,put.R=5.5,Known.age=5,def.int=.01)

# Figure 9 - Roc de Marsal (from current study)
get.age(59,put.points=T,left=1,right=4,put.L=1.8,put.R=2.8,def.int=.005)

# Figure 13A
get.age(76,put.points=T,left=7,right=15,put.L=9,put.R=12.25,def.int=.02,Known.age=8)

# Figure 13B
get.age(76,put.points=T,left=7,right=15,put.L=9,put.R=12.25,def.int=.02,Known.age=8,drop=13,anote='Third molar dropped')

get.age=function (id=1,lo=0.75,hi=25.75,left=0,right=25,def.int=0.1,
     put.points=F,put.L=0,put.R=20,drop=NA,retain=NA,Known.age=NA,
     teeths=rep(NA,13),lab=NA,anote='')
{
##################################################################
tooths <-function(tt,teeth,integ)
{
# j is the index for tooth (1:13)
# i is the score for the given tooth (2:17, as cusp initiation
# cannot be scored in paleoanth)

#   tt=log10(tt+0.75)
   cum=0
   for(j in 1:13){
     i=teeth[j]
     if(!is.na(i)){
   mu1=MFH2[i,j]
   mu2=MFH2[i+1,j]
   while(mu1==mu2)
   {
    i=i+1
    mu2=MFH2[i,j]
   }
   sto=max(pnorm(tt,mu1,.042*log(10))-pnorm(tt,mu2,.042*log(10)),1.e-200)
   cum=cum+log(sto)
   }
   }
   return(exp(cum)/integ)
}
##################################################################
tooths2 <-function(tt,teeth,integ,mu)
{
# j is the index for tooth (1:13)
# i is the score for the given tooth (2:17, as cusp initiation
# cannot be scored in paleoanth)

#   tt=log10(tt+0.75)
   cum=0
   for(j in 1:13){
     i=teeth[j]
     if(!is.na(i)){
   mu1=MFH2[i,j]
   mu2=MFH2[i+1,j]
   while(mu1==mu2)
   {
    i=i+1
    mu2=MFH2[i,j]
   }
   sto=max(pnorm(tt,mu1,.042*log(10))-pnorm(tt,mu2,.042*log(10)),1.e-200)
   cum=cum+log(sto)
   }
   }
   return(exp(cum)/integ*(tt-mu)^2)
}
##################################################################
tooths3 <-function(mu,i.tooth)
{
   mu.tran=c(-Inf,MFH2[,i.tooth])
   sto=pnorm(mu,mu.tran,.042*log(10),lower=F)
   toother=which.max(diff(sto))-1
   if(toother==0) toother=1
   return(toother)
}
##################################################################
if(!is.na(id)) teeth=as.numeric(tooth.scores[id,-1])
else(teeth=teeths)
if(nchar(anote)>0) anote=paste(', ',anote)
print('Tooth scores before any potential dropped scores')
print(teeth)

N.drop=length(drop)
for(i in 1:N.drop) teeth[drop[i]]=NA
if(!is.na(retain[1])){
sto=rep(NA,13)
N.retain=length(retain)
for(i in 1:N.retain){
   sto[retain[i]]=teeth[retain[i]]
}
teeth=sto
}

denom=integrate(Vectorize(tooths,"tt"),lower=log(lo),upper=log(hi),teeth=teeth,integ=1.0)$val
mu=optimize(tooths,lower=log(lo),upper=log(hi),maximum=T,teeth=teeth,integ=denom)$max
mus=rep(NA,13)

for(i in 1:13){
if(!is.na(teeth[i])){
mus[i]=MUS[teeth[i],i]
}
}
between=var(mus,na.rm=T)

within=integrate(Vectorize(tooths2,"tt"),lower=log(lo),upper=log(hi),teeth=teeth,integ=denom,mu=mu)$val
age=log(seq(0,25,def.int)+0.75)
N=NROW(age)
prob=0
for(i in 1:N) prob[i]=tooths(age[i],teeth,denom)
top=max(prob)
top=max(top,dnorm(mu,mu,sqrt(within)))
if(!is.na(id)){
plot(exp(age)-.75,prob,type='l',xlim=c(left,right),ylim=c(0,top),lwd=2,
     xlab='Age (years)',ylab='Density',main=paste(tooth.scores[id,1],anote),axes=F)}
else
{
plot(exp(age)-.75,prob,type='l',xlim=c(left,right),ylim=c(0,top),lwd=2,
     xlab='Age (years)',ylab='Density',main=lab,axes=F)}

box()
axis(1)
if(put.points==T){
     age=log(seq(put.L,put.R,.05)+0.75)
     points(exp(age)-.75,dnorm(age,mu,sqrt(within)),pch=19)
     }
else {lines(exp(age)-.75,dnorm(age,mu,sqrt(within)),lwd=2,lty=2)}
age=log(seq(0,25,def.int)+0.75)
if(!is.na(between)) lines(exp(age)-.75,dnorm(age,mu,sqrt(within+between)),lwd=2,lty=2)
if(!is.na(Known.age)) lines(rep(Known.age,2),c(0,1000),lwd=2)
print(noquote(paste('Mean natural log conception-corrected age = ',mu)))
print(noquote(paste('Within-tooth variance = ',within)))
print(noquote(paste('Between-tooth variance = ',between)))
print(noquote(paste('Lower limit of integration (straight scale) = ',lo)))
print(noquote(paste('Upper limit of integration (straight scale) = ',hi)))
if(!is.na(id)) (lab=as.vector(tooth.scores[id,1]))
new.teeth=rep(NA,13)
for(i in 1:13){
   if(!is.na(teeth[i])) new.teeth[i]=tooths3(mu,i)
}

if(!is.na(new.teeth[1]))
   {if(new.teeth[1]==8) new.teeth[1]=teeth[1]}
for(i in 4:5){
   if(!is.na(new.teeth[i]))
   {if(new.teeth[i]<=8) new.teeth[i]=teeth[i]}}
for(i in 6:7){
   if(!is.na(new.teeth[i]))
   {if(new.teeth[i]<=8) new.teeth[i]=teeth[i]}}
for(i in 8:10){
   if(!is.na(new.teeth[i]))
   {if(new.teeth[i]==8) new.teeth[i]=teeth[i]}}
obj=tooths(tt=mu,teeth=teeth,integ=1)
nu=tooths(tt=mu,teeth=new.teeth,integ=1)
print('Most likely scores given age (estimated after any potential drops)')
print(new.teeth)
print(obj)
print(nu)
for(i in 1:13){
   new.teeth[i]=tooths3(mu,i)
}
if(!is.na(new.teeth[1]))
   {if(new.teeth[1]==8) new.teeth[1]=teeth[1]}
for(i in 4:5){
   if(!is.na(new.teeth[i]))
   {if(new.teeth[i]<=8) new.teeth[i]=teeth[i]}}
for(i in 6:7){
   if(!is.na(new.teeth[i]))
   {if(new.teeth[i]<=8) new.teeth[i]=teeth[i]}}
for(i in 8:10){
   if(!is.na(new.teeth[i]))
   {if(new.teeth[i]==8) new.teeth[i]=teeth[i]}}

print(new.teeth)
print(obj/nu)

if(is.na(between)) return(list(lab=lab,mu=mu,within=within,between=between,p.seq=NA))
return(list(lab=lab,mu=mu,within=within,between=between,p.seq=obj/nu))
}

Getting the "plot.teeth" script (for Figures 6, 8, 10, and 12)

Mark and copy the following to your "RGui." There are vritually no calling parameters here, just the id (as in "get.age") and "Known.L" and "Known.R" which are if there is a known age. ".L" and ".R" are provided in case there is only a known age range, rather than a known point age. Oh, and there's also "teeths" and "sto.lab" to direct enter data as opposed to having it read from "tooth.scores".

Examples:
Once you have copied "plot.tooth" you can paste these pairs of lines (or just the second of each pair as the first is just a comment) to your Gui to reproduce plots from the manuscript.

# Figure 6. "Subject A"
plot.teeth(78,Known.L=5.0)

# Figure 8 - Roc de Marsal (from current study)
plot.teeth(59)

# Figure 10 - Roc de Marsal
plot.teeth(82)

# Figure 12
plot.teeth(76,8)

plot.teeth=function (id=1,Known.L=-10,Known.R=-10,teeths=c(NA,15,14,NA,5,NA,NA,9,6,5,NA,NA,NA),
          sto.lab='Expected Formation at Age 10')
{
##################################################################
tooths <-function(tt,j = 1, i = 2)
{
# j is the index for tooth (1:13)
# i is the score for the given tooth (2:17, as cusp initiation
# cannot be scored in paleoanth)
   mu1=MFH2[i,j]
   mu2=MFH2[i+1,j]
   while(mu1==mu2)
   {
    i=i+1
    mu2=MFH2[i,j]
   }
   return(pnorm(tt,mu1,.042*log(10))-pnorm(tt,mu2,.042*log(10)))
}
##################################################################
if(!is.na(id)){
plot(c(0,1),c(0,1),type='n',xlim=c(-1,22),axes=F,
      xlab='Age',ylab='Density',ylim=c(0,13),
      main=paste(tooth.scores[id,1]))}
else{plot(c(0,1),c(0,1),type='n',xlim=c(-1,22),axes=F,
      xlab='Age',ylab='Density',ylim=c(0,13),main=sto.lab)}
box()
axis(1)
for(i in 0:12) abline(i,0)
lines(c(0,0),c(-1,14))
lines(c(20,20),c(-1,14))
if(Known.L!=-10) lines(rep(Known.L,2),c(-1,14))
if(Known.R!=-10) lines(rep(Known.R,2),c(-1,14))
if(is.na(id)) teeth=teeths
else {teeth=tooth.scores[id,-1]}
print(teeth)

test.uni=c(1,4:10)
for(i in 1:8){
test.it=teeth[test.uni[i]]
if(!is.na(test.it)){
if(test.it==8) return(' No root clefts for uniradicular tooth')
}
}

logage=seq(log(.75),log(20.75),.01)

for(j in 1:13){
if(!is.na(teeth[j])){

mle=optimize(tooths,interval=c(log(.75),log(25)),j=j,i=as.numeric(teeth[j]),maximum=T)$objective/0.95

polygon(c(0,exp(logage)-.75,20,0),
    c(13-j,tooths(logage,j=j,i=as.numeric(teeth[j]))/obj[as.numeric(teeth[j]),j]+13-j,13-j,13-j),density=40)
}
}

teeth.lab=c('c','m1','m2','UI1','UI2','LI1','LI2','C','P3','P4','M1','M2','M3')
for(i in 1:13) text(-.9,13-i+.5,teeth.lab[i])
for(i in 1:13){
j=as.numeric(teeth[i])
if(!is.na(j)) text(20.2,13-i+.5,scores[j],adj=c(0,.5))
}

}

Figure 1

Mark and copy the following to your "RGui"

Fig.01=function ()
{
plot(exp(MFH[,12])-.75,MFH$Mean,xlab='Estimate without Digitized Mean (years)',
ylab='Digitized Mean (years)')
abline(0,1)
}

Fig.01()

Figure 2

Mark and copy the following to your "RGui." Note: you will need to get the package "sfsmisc (see: "packages\install package(s)" from the R menu).

Fig.02=function ()
{
library(sfsmisc)
plot(exp(MFH[,13])-.75,MFH$Mu.HB,xlab='Estimate including Digitized Mean',
ylab='Harris & Buck (2002) Mean',xlim=c(7,12),ylim=c(7,12))
points(exp(MFH[283,13])-.75,MFH$Mu.HB[283],pch=19)
abline(0,1)
p.arrows(exp(MFH[283,13])-.75, MFH$Mu.HB[283],
exp(MFH[283,13])-.75, 10.27614*MFH$SD.HB[283]-.75,fill=1)
}

Fig.02()

Figure 3

Mark and copy the following to your "RGui."

Fig.03=function ()
{

panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    x=x[x>=0]
    h <- hist(x, plot = FALSE,nc=20)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts
    max.y=max(y)
    y <- y/max.y
    rect(breaks[-nB], 0, breaks[-1], y, col="grey", ...)
    to.lab=seq(.2,1,.2)
    axis(4,las=1,at=to.lab,labels=round(to.lab*max.y,1))
}

panel.plot <- function(x, ...)
{
    points(x,...)
    abline(0,1)
    axis(2)
}

Digit = exp(MFH[,13])-.75
HB = MFH$Mu.HB ; HB[283]=9.2
Smith = MFH$Mu.Smith
pairs(cbind(Smith,HB,Digit),upper.panel=NULL,lower.panel=panel.plot,
    diag.panel=panel.hist,
    labels=c('Smith (1991)', 'Harris & Buck (2002) & Harris (2010)','Current Study'),yaxt='n')
}

Fig.03()

Figure 11

Mark and copy the following to your "RGui." Note: you will need to get the package "sfsmisc (see: "packages\install package(s)" from the R menu).

Smith=structure(list(ID = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L), .Label = c("Engis 2 (Smith et al. 2010)", "Gibraltar 2 (Smith et al. 2010)",
"Krapina Max B (Smith et al. 2010)", "Obi Rakhmat (Smith et al. 2010)",
"Scladina (Smith et al. 2010)", "Le Moustier (Smith et al. 2010)",
"Irhoud 3 (Smith et al. 2010)","Qafzeh 10 (Smith et al. 2010)"), class = "factor"), Mean = c(1.71109136883003,
1.73362953894211, 1.92086549950957, 2.19674519855111, 2.411891,
2.836648753268, 2.00777098267658, 1.78272716607327), WithinV = c(0.00228001129603942,
0.00126462693231161, 0.0018835062545257, 0.00223576071953356,
0.002196587, 0.00639174487740038, 0.00163423615783339, 0.00118876925881134
), BetweenV = c(0.0441363051371139, 0.0109695792402254, 0.00658109123119354,
0.018834293669224, 0.005374583, NA, 0.00318328582435489, 0.00996411866888163
), KnownL = c(3, 4.6, 5.9, 6, 8, 11.6, 7.8, 5.1), KnownR = c(3,
4.6, 5.9, 8.1, 8, 12.1, 7.8, 5.1)), .Names = c("ID", "Mean",
"WithinV", "BetweenV", "KnownL", "KnownR"), row.names = c(NA,
8L), class = "data.frame")

Fig.11=function ()
{

# Have ability to plot Qafzeh 10 and Irhoud 3 as well, but skipped here

Smith=Smith[1:6,]
mu=Smith[,2]
V = apply(Smith[,3:4],1,sum,na.rm=T)
Known = apply(Smith[,5:6],1,mean)
est=exp(mu)-0.75
lo=exp(mu-1.96*sqrt(V))-0.75
hi=exp(mu+1.96*sqrt(V))-0.75
plot(Known,est,xlim=c(0,20),ylim=c(0,20), xlab='Known Age',ylab='MFH Estimated Age')
abline(0,1)
for(i in 1:6) lines(rep(Known[i],2),c(lo[i],hi[i]))
for(i in 1:6) lines(c(Smith[i,5],Smith[i,6]),rep(est[i],2))
par(srt=90)
text(Known[1],0,'Engis 2',adj=0)
text(Known[2],0,'Gibraltar 2',adj=0)
text(Known[3],0,'Krapina Max. B',adj=0)
text(Known[4],0,'Obi Rakhmat',adj=0)
text(Known[5],0,'Scladina',adj=0)
text(Known[6],0,'Le Moustier',adj=0)
par(srt=0)
}

Fig.11()

Figure 14

Mark and copy the following to your "RGui", then to draw each plot:
Fig.14('A')
Fig.14('B')
Fig.14('C')
Fig.14('D')

Scladina.UC=function ()
{
mu=MFH2[1,8]
ymax=dlnorm(exp(mu),mu,log(10)*.042)
plot(seq(.1,13.75,.01)-.75,dlnorm(seq(.1,13.75,.01),mu,log(10)*.042),type='l',
     xlab='Age (years)',ylab='Density',xlim=c(0,14),ylim=c(-ymax/20,ymax),main='Canine',axes=F)
box()
axis(1)
text(exp(mu)-.75,-ymax/20,'C.i')
lines(c(0,100),rep(0,2))
lines(rep(102/325.25,2),c(0,10),lwd=2,lty=3)

mu=MFH2[6,8]
lines(seq(.1,6.75,.01)-.75,dlnorm(seq(.1,6.75,.01),mu,log(10)*.042))
lines(rep(1403/325.25,2),c(0,10),lwd=2,lty=2)
text(exp(mu)-.75,-ymax/20,'Cr.c')

mu=MFH2[13,8] # A.5 (observed stage)
lines(seq(6,15,.01)-.75,dlnorm(seq(6,15,.01),mu,log(10)*.042))
text(exp(mu)-.75,-ymax/20,'A.5')
mu=MFH2[14,8] # A.c (next higher stage)
lines(seq(6,15,.01)-.75,dlnorm(seq(6,15,.01),mu,log(10)*.042))
text(exp(mu)-.75,-ymax/20,'A.c')
lines(rep(8,2),c(0,10),lwd=2)

}

Scladina.UM1=function ()
{
mu=MFH2[6,11]
ymax=dlnorm(exp(mu),mu,log(10)*.042)
plot(seq(.1,13.75,.01)-.75,dlnorm(seq(.1,13.75,.01),mu,log(10)*.042),type='l',
     xlab='Age (years)',ylab='Density',xlim=c(0,14),ylim=c(-ymax/20,ymax),main='First Molar',axes=F)
box()
axis(1)
text(exp(mu)-.75,-ymax/20,'Cr.c')
lines(c(0,100),rep(0,2))
lines(rep(858/325.25,2),c(0,10),lwd=2,lty=2)

mu=MFH2[14,11]
lines(seq(.1,15,.01)-.75,dlnorm(seq(.1,15,.01),mu,log(10)*.042))
lines(rep(8,2),c(0,10),lwd=2)
text(exp(mu)-.75,-ymax/20,'A.c')

#mu=MFH2[13,8] # A.5 (observed stage)
#lines(seq(6,15,.01)-.75,dlnorm(seq(6,15,.01),mu,log(10)*.042))
#text(exp(mu)-.75,-ymax/20,'A.5')
#mu=MFH2[14,8] # A.c (next higher stage)
#lines(seq(6,15,.01)-.75,dlnorm(seq(6,15,.01),mu,log(10)*.042))
#text(exp(mu)-.75,-ymax/20,'A.c')
#lines(rep(8,2),c(0,10),lwd=2)

}

Scladina.UM2=function ()
{
mu=MFH2[1,12]
ymax=dlnorm(exp(mu),mu,log(10)*.042)
plot(seq(.1,15.75,.01)-.75,dlnorm(seq(.1,15.75,.01),mu,log(10)*.042),type='l',
     xlab='Age (years)',ylab='Density',xlim=c(2,14),ylim=c(-ymax/20,ymax),
     main='Second Molar',axes=F)
box()
axis(1)
text(exp(mu)-.75,-ymax/20,'C.i')
lines(c(0,100),rep(0,2))
lines(rep(1084/325.25,2),c(0,10),lwd=2,lty=3)

mu=MFH2[6,12]
lines(seq(2,13,.01)-.75,dlnorm(seq(2,13,.01),mu,log(10)*.042))
lines(rep(2051/325.25,2),c(0,10),lwd=2,lty=2)
text(exp(mu)-.75,-ymax/20,'Cr.c')

mu=MFH2[10,12]
lines(seq(2,15,.01)-.75,dlnorm(seq(2,15,.01),mu,log(10)*.042))
text(exp(mu)-.75,-ymax/20,scores[10],adj=c(1,.5))

mu=MFH2[11,12]
lines(seq(2,15,.01)-.75,dlnorm(seq(2,15,.01),mu,log(10)*.042))
text(exp(mu)-.75,-ymax/20,scores[11],adj=c(0,.5))
lines(rep(8,2),c(0,10),lwd=2)

}

Scladina.UM3=function ()
{
mu=MFH2[1,13]
ymax=dlnorm(exp(mu),mu,log(10)*.042)
plot(seq(.1,15.75,.01)-.75,dlnorm(seq(.1,15.75,.01),mu,log(10)*.042),type='l',
     xlab='Age (years)',ylab='Density',xlim=c(6,14),ylim=c(-ymax/20,ymax),
     main='Third Molar',axes=F)
box()
axis(1)
text(exp(mu)-.75,-ymax/20,'C.i')
lines(c(0,100),rep(0,2))
lines(rep(2155/325.25,2),c(0,10),lwd=2,lty=3)

mu=MFH2[5,13]
lines(seq(2,15,.01)-.75,dlnorm(seq(2,15,.01),mu,log(10)*.042))
text(exp(mu)-.75,-ymax/20,scores[5],adj=c(1,.5))

mu=MFH2[6,13]
lines(seq(2,15,.01)-.75,dlnorm(seq(2,15,.01),mu,log(10)*.042))
text(exp(mu)-.75,-ymax/20,scores[6],adj=c(0,.5))
lines(rep(8,2),c(0,10),lwd=2)

}

Fig.14=function (which='A')
{
if(which=='A') {Scladina.UC(): return()}
if(which=='B') {Scladina.UM1(): return()}
if(which=='C') {Scladina.UM2(): return()}
if(which=='D') {Scladina.UM3(): return()}

}

Well, you've been awfully patient, so...

To plot all of the "tooth.scores", copy and paste the following. This will plot all of 'em to a *.pdf file called "test.pdf."

plot.all.teeth=function ()
{
plot.teeth()
dev.copy(device=pdf,file='test.pdf')
for(i in 2:82) {
plot.teeth(i)
dev.copy()
}
dev.off()
}

plot.all.teeth()