First, please visit http://cran.r-project.org/ and istall “R” and “RStudio”
1. Setting up your data
To import into R or R studio, we recommend setting up your individual cases and scores as a .csv file. The categories for which teeth are scored and how the sheet is set up is explained in Shackelford et al. 2012 and via http://faculty.las.illinois.edu/lylek/SHK2012/index.htm
You will need to set up your .csv file as follows, ensuring that the first test individual is a primer, containing no missing values. The value for ‘Neander’ will always be FALSE and the value for ‘Obs’ will always be 1. Label this file tooth.scores and save in your working directory as a .csv file.
| ID | dc | dm1 | dm2 | UI1 | UI2 | LI1 | LI2 | C | P3 | P4 | M1 | M2 | M3 | Neander | Obs |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Test | 10 | 10 | 10 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | FALSE | 1 |
| EX1 | NA | NA | 14 | 6 | 4 | NA | 7 | NA | NA | 7 | 8 | 3 | NA | FALSE | 1 |
2. Adding the Data Tables
Four data tables are used in the age calculation: the MFH, MFH2, MUS, and OBJ tables. Their function is described in Shackelford et al. 2012 and the updated version are also for reference in Kamnikar et al. 2018. These files are available below as text files. Click the green links to open the downloadable tables. Once the table appears in the window, right click the table and save it as a text document. The saved text files can then be imported using the following lines of code. Make sure the file names match the labels in quotations and do not modify the table names at the beginning of each function (MFH, MFH2, mus, obj).
MFH <- read.delim(“MFH_table.txt”)
MFH2 <- read.delim(“MFH2_table.txt”)
mus <- read.delim(“MUS_table.txt”)
obj <- read.delim(“obj_table.txt”)
Download the MFH table, which is already a .csv file to your R working directory. Start R and save the data to an object by typing the following:
MFH=read.csv('MFH.csv')[,-1]
Import your tooth.scores file
tooth.scores=read.csv('tooth.scores.csv')
Copy/Paste Option
Select and copy the following into your “RGui” to get the MFH 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, 0.92354436059, 1.13167626746, 1.40545953038, 1.55762869551, 1.55762869551, 1.70317412693, 1.82735985485, 1.95170759705, 2.08813489475, 2.11866225483, 2.53734095964, Inf, Inf, Inf, Inf),LI2 = c(-Inf, -Inf, -Inf, 1.02369116548, 1.28340301677, 1.55242301771, 1.66905496800, 1.66905496800, 1.80538291392, 1.93002754271, 2.03241544563, 2.13128655884, 2.22638955878, 2.58621019707, Inf, Inf, Inf, Inf), C = c(0.20351210936, 0.39632345019, 0.92602321722, 1.25440447762, 1.58071945700, 1.69908664188, 1.83126860232, 1.83126860232, 2.02135549184, 2.17322214318, 2.32958723170, 2.43774952648, 2.52003245176, 2.77353514924, Inf, Inf, Inf, Inf), P3 = c(0.96384203373, 1.13624520347, 1.35784344495, 1.56900126033, 1.70265499415, 1.84773014094, 1.98460874791, 1.98460874791, 2.13821150350, 2.28404426929, 2.37197122048, 2.44300334169, 2.51406083944, 2.76537526772, Inf, Inf, Inf, Inf), P4 = c(1.37215740356, 1.53303191723, 1.59959967530, 1.70635588995, 1.84683933550, 1.98817567261, 2.11911888133, 2.11911888133, 2.23569714745, 2.34307423212, 2.42630160184, 2.49899855748, 2.58787681082, 2.80922195614, Inf, Inf, Inf, Inf), M1 = c(-Inf, -0.07666508494, 0.27117152032, 0.89369741422, 1.07961633592, 1.30750815538, 1.43125577640, 1.56743814127, 1.69247306136, 1.82456542095, 1.97666795481, 2.10208327216, 2.22543403381, 2.59402260495, Inf, Inf, Inf, Inf), M2 = c(1.38676924834, 1.53925147178, 1.67808672474, 1.74071175872, 1.89198180703, 2.03869766911, 2.11671324136, 2.19900077674, 2.31789726110, 2.40463613839, 2.48648041078, 2.56690129734, 2.67376926713, 2.84945137891, Inf, Inf, Inf, Inf), M3 = c(2.22682113507, 2.33286208388, 2.40322196982, 2.49093013877, 2.51624380510, 2.59082326080, 2.68334758846, 2.70644891974, 2.78186741670, 2.82704731428, 2.86094294721, 2.86676189226, 2.96836107676, 3.01507891455, 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, 1.02785, 1.26875, 1.4821, 1.6305, NA, 1.76535, 1.8897, 2.0199, 2.10345, 2.3282,3.1092276, NA, NA, NA, NA, NA, NA, 1.1531, 1.41775, 1.6107, 1.7373, NA, 1.86755,1.9811, 2.08195, 2.1793, 2.40655,3.1333485, NA, NA, NA, 0.293750636861609, 0.6578, 1.09885, 1.4183, 1.64015, 1.76515,
1.9841805,NA, 2.2001235, 2.3346377, 2.4673005, 2.5858556,
2.7257493, 3.2431305, NA, NA, NA, 1.0505, 1.2481, 1.4635, 1.63575, 1.7757, 1.9750165,2.1568925, NA, 2.2843027, 2.4076769, 2.49345903,
2.57169881, 2.7063597, 3.2378149, NA, NA, NA, 1.4519,1.566, 1.653, 1.77675, 1.9713801,2.1484414, 2.2480678, NA, 2.3583136, 2.4667526,
2.54780958, 2.6389146, 2.755522, 3.2474988, NA, NA, NA, NA, 0.0911094, 0.57645, 0.9852,1.19335, 1.3693, 1.4984, 1.6292,
1.8091673, 2.0102497, 2.147042, 2.284937, 2.5487414, 3.2033894, NA, NA, NA, 1.48445, 1.6071,1.69935, 1.8671915, 2.069676, 2.1814644,
2.24409589, 2.3254254, 2.4310314, 2.51619333, 2.5908947, 2.6762534, 2.7902696, 3.2684917, NA, NA, NA, 2.3489374, 2.4371235, 2.49647252, 2.5513732, 2.61202925, 2.68115743, 2.72049858, 2.7501018, 2.79979092, 2.85643265, 2.95054171, 3.01871568, 3.3621931, 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")))3. Adding the functions
You will need to add the get.age() and plot.teeth() functions to produce age estimates and graphics. The tooth.test() function can be added if running a batch or more than a few individuals.
The get.age() function provides several outputs, which are explained in Shackelford et al. 2012 and Kamnikar et al. 2018. The outputs are used to calculate the maximum likelihood age estimates and the confidence intervals. The associated graphic provides the maximum likelihood age estimate (solid center line) and the within-tooth variance (blue line) and the within + between tooth variance (red line). The curved, dashed line represents the best-fitting log normal distribution. Copy and paste the following code in your R console.
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(0.75),upper=log(25.75),teeth=teeth,integ=1.0)$val
mu=optimize(tooths,lower=log(0.75),upper=log(25.75),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)
abline(v=exp(mu)-0.75)
text(exp(mu)-0.75,top,round(exp(mu)-0.75,4),pos=4)
abline(v=exp(mu-(2*(within+between)^0.5))-0.75,lty=2,col='red')
text(exp(mu-(2*(within+between)^0.5))-0.75,top-0.5,round(exp(mu-(2*(within+between)^0.5))-0.75,3),pos=2,cex=0.75)
abline(v=exp(mu+(2*(within+between)^0.5))-0.75,lty=2,col='red')
text(exp(mu+(2*(within+between)^0.5))-0.75,top-0.5,round(exp(mu+(2*(within+between)^0.5))-0.75,3),pos=4,cex=0.75)
abline(v=exp(mu-(2*(within)^0.5))-0.75,lty=2,col='blue')
abline(v=exp(mu+(2*(within)^0.5))-0.75,lty=2,col='blue')
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))
}
The plot.teeth() function will provide a density graphic for each tooth scored and used in the age calculation. Copy and paste the following code in your R console.
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))
}
}
Run a batch file of individuals using the tooth.test() function. Copy and paste in following code in your R console.
tooth.test=
function ()
{
tooth_ages <- c()
for (i in 1:nrow(tooth.scores))
{
# vector output
m<-tooth.scores[i,]
m[is.na(m)]<-0
if (m[,2] > 0) {h=3}
else if (m[,3] > 0) {h=3}
else if (m[,4] > 0) {h=3}
else h=25.75
if (m[,2] > 3) {h=5} # set to (m[,2] > 1) {h=5} for the recalculated MFH2 and MUS matrices
else if (m[,3] > 3) {h=5} # set to (m[,3] > 1) {h=5} for the recalculated MFH2 and MUS matrices
else if (m[,4] > 3) {h=5} # set to (m[,4] > 1) {h=5} for the recalculated MFH2 and MUS matrices
else h=h
if (m[,2] > 8) {h=15}
else if (m[,3] > 8) {h=15}
else if (m[,4] > 8) {h=15}
else h=h
#if ((m[,12] < 12) & (m[,2] < 1) & (m[,3] < 1) & (m[,4] < 1)) {h=15}
#else h=h
model <- get.age(i,hi=h,def.int=0.01)
scores_i <- cbind(model$lab,h,model$mu,model$within,model$between,model$p.seq)
# add vector to a dataframe
age_i <- data.frame(scores_i)
tooth_ages <- rbind(tooth_ages,age_i)
}
write.table(tooth_ages, file="tooth.test-results.csv",sep=",",col.names=c("lab","hi","mu","within","between","p.seq"),row.names=FALSE)
return(data.frame(tooth_ages))
}