Konigsberg et al. 2008, in press JFS

Calculations and Figures from Konigsberg, Herrmann, Wescott, and Kimmerle (2008)

First, “R” is your friend, so:

A Viewer for Age Estimation

All along in this project there has been a demand for “simple” methods for age estimation. The following provides same, but it is extremely important to understand some of the basics behind the calculations, as well as how to control them. First, let’s look at the defaults in the call.

ip = 1 , indicates the phase, which can be 1, 2, 3, 4, 5, or 6

top = 100, indicates the highest possible age.

area = 0.5, indicates the 50% highest posterior density region. If you specify area = 0.95 you will get the 95% region. Don’t go much above 0.95 as you’ll hit a numerical issue.

There are two sections that are “hard-wired,” but you can always edit these to fit your particular problem. The prior age-at-death here is a Gompertz starting at age 17 (the “bot=17”), and is taken from a reference sample of 212 Balkan individuals. You can change the starting age using “bot” and change the Gompertz values as well. You can also use a non-zero value for “a2” to get a Makeham model. The other “hard-wired” part is for the transition analysis parameters. These are the mean and common standard deviation on a log scale.

So if you are ready copy and paste what is below into your “Rgui.”

lrage.viewer=function (ip=1,top=100,area=.5)

#  Prior age-at-death from Gompertz fit to 212 (shifted 17 years)

bot=17

a2=0.0

a3=0.01346757

b3=0.04201534

#  "Transition analysis" parameters

mu=c(2.9177,3.0436,3.3147,3.7612,4.1916)

sdev=0.3118

##################################################################

probit=function (t)

imax = 6

if(ip==1) return(1-pnorm(log(t),mu[1],sdev))

if(ip==imax) return(pnorm(log(t),mu[imax-1],sdev))

return(pnorm(log(t),mu[ip-1],sdev)-pnorm(log(t),mu[ip],sdev))

##################################################################

#  pia()                                                         #

#  This function finds the unormalized posterior density of age  #

#  conditional on the observed phase and the Gompertz-Makeham    #

#                                                                #

##################################################################

pia<-function(t)

  p<-probit(t)

  h.t<-a2+a3*exp(b3*(t-bot))

  S.t<-exp(-a2*(t-bot)+a3/b3*(1-exp(b3*(t-bot))))

  return(S.t*h.t*p)

##################################################################

#  pia2()                                                        #

#  This function finds the normalized posterior density of age   #

#  conditional on the observed phase and the Gompertz-Makeham    #

#                                                                #

##################################################################

pia2<-function(t)

  p<-probit(t)

  shift<-bot

  h.t<-a2+a3*exp(b3*(t-bot))

  S.t<-exp(-a2*(t-bot)+a3/b3*(1-exp(b3*(t-bot))))

  return(S.t*h.t*p/denom)

##################################################################

#  pia5()                                                        #

#  Searches left and right to find cut points, such that the     #

#  integral between cut points less the desired area is near     #

#  zero                                                          #

#                                                                #

##################################################################

pia5<-function(divid)

   pia4<-function(tt)

     cd<-pia2(tt)-prob*divid

     return(cd)

  low<-uniroot(pia4,c(bot,hpd))$root

  hi<-uniroot(pia4,c(hpd,top))$root

  return(integrate(pia2,low,hi)$value-area)

##################################################################

#  pia6()                                                        #

#  Returns the cut points after they have been found by pia5()   #

#                                                                #

##################################################################

pia6<-function(divid)

  pia4<-function(tt)

    cd<-pia2(tt)-prob*divid

    return(cd)

  low<-uniroot(pia4,c(bot,hpd))$root

  hi<-uniroot(pia4,c(hpd,top))$root

  return(c(low,hi))

#################################################################

pia7<-function(t)

{return(integrate(pia2,bot,t)$value-area)}

#################END OF FUNCTIONS################################

denom<-integrate(pia,bot,top)$value

hpd<-optimize(pia,c(bot,top),maximum=T)$maximum

prob<-pia2(hpd)

prob2<-pia2(bot)

tunebot<-1.00001*prob2/prob

tunebot=max(.1,tunebot)

tune<-c(tunebot,.99999)

plot(c(bot,bot:top),c(0,pia2(bot:top)),

    type='l',xlab='Age',ylab='Density',xlim=c(bot,top))

cat('\n\n Maximum density occurs at ',round(hpd,digits=2),'\n')

lines(c(0,100),c(0,0))

hi<-uniroot(pia7,c(bot,top))$root

if(pia2(hi)<prob2)

  x<-c((bot:top)[(bot:top)<hi],hi)

  polygon(c(bot,x,x[NROW(x)]),c(0,pia2(x),0),density=10)

  cat(area*100,'% HPD from bottom [',bot,']',round(hi,digits=2),'\n\n')

   return(c(round(bot,digits=2),round(hpd,digits=2),round(hi,digits=2)))

if(pia2(hi)>prob2)

 divid<-uniroot(pia5,tune)$root

 x<-pia6(divid)

 x<-c(x[1],(bot:top)[(bot:top)>x[1] & (bot:top)<x[2]],x[2])

 polygon(c(x[1],x,x[NROW(x)]),c(0,pia2(x),0),density=10)

 cat(area*100,'% HPD \n')

 return(c(round(x[1],digits=2),round(hpd,digits=2),round(x[NROW(x)],digits=2)))

lrage.viewer()

To try other problems:

lrage.view(area=.95,top=50) # Still phase I, but now 95% and truncate graph at age 50

lrage.view(ip=6,area=.95) # Phase VI, 95%

Symphysis Data (N = 1,766)

The remainder of this file deals with analysis of a large sample, which you can get into “R.” Mark, copy and paste the yellow highlighted block into your “RGui.” This will create the data object “allmale2.”

# Start marking here

age=c(18.56810404,21.03764545,20.26830938,18.82272416,17.62080767,21.20465435,24.02190281,21.26214921,20.67624914,21.24298426,23.05270363,18.96783025,

19.39219713,17.5797399,21.07323751,18.09445585,19.40862423,18.91854894,19.72073922,18.99520876,18.05065024,20.3613963,21.76865161,19.51813826,17.59342916,

18.79534565,20.08761123,18.7871321,20.53661875,19.816564,19.77823409,19.63586585,18.31622177,21.9192334,19.90691307,18.48323066,20.03832991,21.9192334,19.12388775,20.10403833,20.95824778,18.99794661,21.2073922,18.90212183,18.31895962,19.34565366,20.07939767,

18.275154,17.67556468,20.09582478,20.66255989,17.85900068,18.35455168,23.41136208,18.68583162,22.66392882,19.86036961,20.49281314,21.31416838,20.38056126,18.89938398,18.92402464,22.00410678,18.83367556,19.57563313,18.9596167,19.60027379,19.42778919,20.41889117,

20.13689254,18.93223819,20.0054757,21.06776181,21.08418891,20.20260096,18.20944559,19.34565366,19.88501027,19.79466119,20.5119781,19.08829569,18.66392882,20.31485284,20.42162902,19.78918549,18.67761807,19.70978782,18.97604381,21.59890486,18.85010267,

18.90212183,18.39014374,18.7816564,21.99315537,18.15468857,17.96303901,19.78370979,20.99657769,19.70704997,17.54688569,21.29500342,19.41409993,17.77138946,18.68309377,22.02600958,22.40383299,20.12320329,24.42710472,30.11362081,23.32922656,20.74195756,21.15811088,

19.4770705,19.96988364,20.33675565,19.95345654,18.2532512,22.41478439,20.50102669,19.23613963,18.56810404,21.3798768,18.71868583,19.29637235,18.92676249,19.26078029,18.91854894,19.6495551,19.8329911,19.55099247,19.1430527,23.56741958,20.95003422,18.89664613,

19.72621492,20.06297057,19.29089665,19.41683778,18.32991102,18.4421629,21.42094456,19.6495551,19.69609856,20.49555099,20.81040383,19.23887748,21.95482546,23.47433265,20.0054757,18.97878166,19.71800137,21.70020534,19.07460643,20.38056126,20.17522245,19.44421629,

21.29226557,20.65708419,19.93155373,20.31759069,18.0752909,21.908282,18.20396988,17.9192334,28.47912389,17.56057495,19.03901437,20.82409309,21.27857632,23.57289528,22.29705681,23.26078029,21.56605065,22.91854894,24.92813142,20.5284052,22.25051335,19.54825462,

23.59206023,26.07802875,22.16837782,22.71868583,27.05817933,21.53593429,28.35865845,18.82272416,22.23682409,20.90075291,23.26899384,21.47296372,19.80013689,20.03832991,20.22997947,23.4661191,20.65434634,35.59479808,23.63039014,21.88637919,23.27173169,22.72963723,

23.19233402,27.93429158,24.36960986,25.24845996,20.64339493,25.85626283,24.15879535,19.63312799,27.4934976,16.93908282,25.77960301,24.26283368,24.27652293,24.5119781,24.24640657,23.08555784,21.44010951,20.20807666,23.88227242,22.76796715,28.62696783,27.69883641,

21.39904175,23.26078029,21.86721424,21.0513347,20.28199863,23.24161533,18.98151951,21.724846,22.98699521,36.67077344,20.60780287,21.14715948,26.89390828,23.04175222,25.72758385,34.17659138,27.75633128,27.00342231,26.15468857,23.24435318,22.66940452,25.48665298,

25.5523614,26.0752909,34.23408624,27.09924709,21.17453799,30.40109514,22.66392882,28.16427105,28.40793977,22.0971937,27.3045859,28.45995893,26.5927447,23.66872005,40.06844627,23.70431211,22.0971937,35.16495551,39.39493498,24.91991786,26.38740589,28.12594114,

26.26146475,25.52498289,25.01574264,22.61738535,26.55989049,34.26694045,32.73374401,34.63655031,30.07802875,23.68514716,20.65982204,25.65092402,26.38466804,23.75085558,32.07939767,25.5523614,35.96440794,25.29774127,37.94661191,33.37166324,30.11088296,30.51882272,

26.78439425,24.28199863,21.70841889,29.42915811,25.9247091,30.93497604,27.09650924,42.15742642,34.52429843,24.79123888,28.22997947,39.49897331,29.21560575,28.88980151,35.16769336,36.62422998,25.98494182,22.92128679,32.82956879,41.82888433,29.19096509,32.07665982,

28.59411362,30.2532512,28.42162902,33.0349076,35.07734428,34.55441478,33.43737166,25.93566051,31.49075975,33.69746749,31.89596167,28.41067762,31.29089665,26.75427789,30.42299795,26.98973306,30.57905544,37.71663244,29.81519507,44.24640657,41.09787817,27.52087611,

38.99520876,32.69267625,30.71047228,26.66392882,29.94113621,31.05544148,23.98631075,32.53114305,31.28542094,29.10609172,30.96235455,31.92060233,32.85420945,32.85147159,30.97878166,38.72689938,39.77275838,36.50102669,23.137577,27.95071869,39.09377139,42.00958248,

37.35797399,32.3504449,39.05544148,50.4476386,18,18,18,17,19,20,

16,18,16,17,16,18,21,18,19,17,18,18,19,18,16,16,20,20,15,19,15,22,20,21,20,17,18,15,18,17,16,16,19,19,17,19,15,19,17,18,18,17,16,22,20,17,19,19,21,23,17,19,20,19,17,17,18,18,20,17,16,18,18,19,17,17,16,18,17,19,16,17,18,23,21,23,17,18,21,20,18,22,21,21,20,19,22,18,21,21,23,20,23,18,19,20,24,24,21,19,23,20,21,22,22,19,20,22,19,21,24,21,21,19,22,23,22,23,21,24,24,19,22,27,22,20,22,22,23,19,23,24,22,24,27,22,29,21,22,21,22,27,24,34,27,24,21,21,25,23,26,19,23,35,24,25,26,45,24,23,31,20,29,26,25,25,27,20,25,21,30,26,25,33,24,21,27,28,30,30,24,23,23,24,28,26,25,20,36,22,24,25,25,23,28,22,25,25,24,29,31,22,36,25,26,30,34,27,34,29,27,25,43,22,24,24,37,27,27,31,23,51,28,39,23,33,25,27,26,29,27,26,31,27,27,27,36,50,71,46,51,36,30,43,32,32,51,33,35,41,42,27,46,46,34,29,43,47,39,29,28,27,38,41,28,31,22,41,36,32,33,41,52,37,31,59,26,33,30,40,34,48,59,52,35,35,34,42,31,24,23,47,48,26,30,20,30,24,27,58,37,25,49,30,45,27,27,31,40,44,29,49,34,31,33,29,32,40,31,39,33,46,47,22,32,31,34,43,32,35,26,37,26,25,37,32,

34,41,33,49,45,30,41,36,62,36,42,23,28,34,48,43,64,43,40,27,30,53,42,33,50,27,31,48,47,30,54,28,29,28,39,30,32,28,46,35,38,43,50,33,20,40,51,33,36,23,33,29,25,25,45,46,43,71,54,61,79,60,65,66,75,41,49,59,54,67,62,73,74,87,59,55,63,64,54,67,56,49,59,54,64,69,64,54,70,65,60,65,54,66,64,50,74,59,52,54,70,51,43,80,56,69,60,67,45,45,86,59,62,75,43,34,72,47,54,48,37,47,51,59,52,44,48,60,66,81,62,72,59,40,52,46,42,42,29,68,38,61,26,60,34,78,53,74,51,48,69,46,38,52,49,27,47,45,53,58,42,34,61,35,27,71,75,47,51,53,32,74,40,51,46,46,39,38,62,54,38,44,80,44,74,56,32,55,32,42,47,62,36,73,54,60,56,36,34,76,54,59,45,29,36,47,50,74,57,54,52,75,27,58,28,27,48,52,63,55,29,49,44,45,40,34,54,53,40,51,51,57,33,49,31,32,49,41,40,57,59,43,36,55,48,21,46,50,47,28,31,30,41,32,65,38,54,46,37,32,31,43,34,48,44,49,64,56,32,47,29,37,52,50,29,38,59,46,31,51,33,45,40,60,60,40,31,49,54,41,32,32,59,46,58,62,59,55,56,53,86,92,77,55,79,73,73,66,63,82,50,79,80,69,79,63,74,65,64,78,55,61,70,64,87,66,74,52,73,69,63,72,55,49,73,52,

59,68,67,68,91,71,51,76,49,72,69,66,47,56,75,67,68,61,77,49,65,48,55,71,31,72,55,72,66,63,50,63,53,55,59,60,69,63,48,48,53,67,68,45,60,46,42,50,55,75,64,49,53,36,78,63,62,42,26,20,23,23,23,25,31,19,22,19,20,20,20,30,18,22,24,22,30,17,28,24,27,24,34,33,39,34,36,28,31,22,26,23,34,33,21,24,55,45,50,50,30,32,32,23,59,40,41,26,59,28,25,22,25,28,29,26,27,46,32,56,27,30,22,22,32,27,24,39,29,34,29,29,60,26,37,38,70,38,26,26,47,42,26,40,35,48,59,26,61,42,63,31,23,23,62,27,42,51,61,39,46,61,29,37,26,26,53,36,20,73,53,53,36,37,53,45,32,31,54,29,35,34,58,36,33,32,29,42,62,41,30,55,60,52,33,64,48,72,32,42,45,52,54,74,45,29,57,35,69,38,27,50,48,45,32,58,44,37,60,46,30,52,57,57,56,58,57,60,43,50,58,49,62,38,77,30,84,41,74,50,49,42,23,31,44,38,57,43,44,46,35,54,82,57,28,58,46,49,49,48,40,40,54,33,46,48,32,40,31,59,47,41,32,30,44,61,43,31,42,43,51,63,43,28,52,58,53,56,35,28,88,46,29,22,25,66,79,45,85,41,36,58,63,65,56,36,42,51,36,51,54,61,35,62,37,56,35,26,31,36,31,63,26,37,51,31,51,32,81,81,83,80,70,34,53,54,

58,36,39,86,37,49,98,42,83,82,54,45,39,68,64,34,43,73,49,27,55,45,31,66,55,47,64,47,57,60,37,60,70,60,65,68,52,50,67,62,59,76,40,79,59,72,44,43,65,48,57,84,45,40,43,48,66,69,76,62,56,71,47,75,75,32,44,52,73,78,59,40,46,31,61,41,69,42,41,76,36,31,78,25,48,63,62,47,45,44,41,63,37,61,61,45,74,54,54,82,47,71,54,53,80,81,55,55,55,63,53,85,89,85,91,39,37,63,67,55,38,52,72,49,50,54,67,60,59,68,79,75,86,51,77,65,51,42,62,39,36,87,84,56,71,17,18,18.5,19,19,20,20.5,20.6,20.8,21,22,22,25.9,20,20.6,22.7,23,26,33,22,23.4,23.7,24,24.8,

25,25,25,26.3,27,27.9,28.4,28.9,31.9,35,35.6,42,43,45,45,24,25.2,25.9,27.4,28,28,28.9,29,30,30,30,31,32,32,33,33,33.1,33.1,33.9,34,35,35,35,36,36.8,37.2,38.7,38.8,39,39.4,39.4,40,40,40,41,42,42,43,45,46,46,47,47,47,49,49,49.1,50,51,52,52.1,53,53,54,55,55,56,57,62,62,63,64,66,70,74,23.7,27,27.9,30,31,31,33,35,35.6,36,36,36.3,37,37,38.8,39,40,40,40.9,41,41,41.5,42,43,43.7,44,44,45,45,46,46,47,47,47,47,49,49,49.9,50,50.7,51,51.1,52,52,52.9,53,56,57,57,57,57,57,58,58,58,59,59,60,60,60,60,60,61,62,63,

65,66,66,69,72.2,74,34,43,45,45,46,46,48,51,52.7,53,54,55,56,57,57,59,60,60,61,64,64,66,66,67,70,72,73,73,74,77,78,79,81,81,83,84,85,20,22,30,39,58,49,51,59,68,75,41,42,46,46,47,49,51,57,75,79,82,42,43,45,51,53,56,58,59,59,66,76,34,60,62,71,82);

Collection=rep(c('Korean War','LA Coroner','Terry','Balkan','Thai'),c(358,737,422,212,37));

suchey=rep(rep(1:6,5),c(167,83,22,50,28,8,119,81,43,153,241,100,22,46,27,173,122,32,13,6,20,65,71,37,0,5,5,11,11,5));allmale2=data.frame(Collection,age,suchey);rm(Collection);rm(age);rm(suchey)

# End marking here

JFS Figure 1.

fig.01=function ()

sto<-glm((c(0,0,1,1,1,1)[allmale2[,3]])~log(allmale2[,2]),family=binomial(link='probit'))

sto<-as.vector(sto$coef)

mu<-(-sto[1]/sto[2])

sdev<-1/sto[2]

print(exp(mu))

x<-seq(0,60,by=.1)

plot(x,dlnorm(x,mu,sdev),type='l',main='Age-at-Transition',

      xlab='Age',ylab='Density')

sto<-glm((c(0,0,1,1,1,1)[allmale2[,3]])~(allmale2[,2]),family=binomial(link='probit'))

sto<-as.vector(sto$coef)

mu<-(-sto[1]/sto[2])

sdev<-1/sto[2]

print(c(mu,sdev))

lines(x,dnorm(x,mu,sdev),lty=2)

legend(40,.06,lty=c(1,0,2),legend=c('Log Normal','','Normal'),bty='n')

fig.01()

JFS Figure 2. – NOTE: you will need to download the package GenKern

fig.02=function (xb=4)

library(GenKern)

x<-allmale2$age

y<-allmale2$suchey

rec<-c(0,0,1,1,1,1)

y<-rec[y]

all<-KernSec(x,range.x=15:100,xb=xb)$yden

p<-KernSec(x[y==0],range.x=15:100,xb=xb)$yden/all

#par(mfrow=c(1,2))

sto<- as.vector(glm(y~log(x),family=binomial(link='probit'))$coefficients)

sd<-1/sto[2]

mu<--sto[1]*sd

plot(15:100,p,type='l',ylim=c(0,1),xlab='Age',ylab='Probability',xlim=c(15,80))

lines(15:100,1-pnorm(log(15:100),mu,sd),lty=2)

lines(15:100,1-p)

lines(15:100,pnorm(log(15:100),mu,sd),lty=2)

text(60,.92,cex=1.25,'Suchey-Brooks III - VI')

text(60,.08,cex=1.25,'Suchey-Brooks I - II')

fig.02()

JFS Figure 3.

fig.03=function (xb=4)

library(GenKern)

x<-allmale2$age

y<-allmale2$suchey

par(mfrow=c(2,3))

library(MASS)

sto<-polr(factor(y)~log(x),method='probit')

sd<-1/as.vector(sto$coeff)

mu<-as.vector(sto$zeta)*sd

x<-seq(0,100,by=0.1)

plot(x,dlnorm(x,mu[1],sd),type='l',xlab='Age',ylab='Density',main='I / II')

plot(x,dlnorm(x,mu[2],sd),ylim=c(0,.09),

   type='l',xlab='Age',ylab='Density',main='II / III')

#=============

sto<-glm((c(0,0,1,1,1,1)[allmale2[,3]])~log(allmale2[,2]),family=binomial(link='probit'))

sto<-as.vector(sto$coef)

mu2<-(-sto[1]/sto[2])

sdev2<-1/sto[2]

lines(x,dlnorm(x,mu2,sdev2),lty=2)

#=============

plot(x,dlnorm(x,mu[3],sd),type='l',xlab='Age',ylab='Density',main='III / IV')

plot(x,dlnorm(x,mu[4],sd),type='l',xlab='Age',ylab='Density',main='IV / V')

plot(x,dlnorm(x,mu[5],sd),type='l',xlab='Age',ylab='Density',main='V / VI')

fig.03()

JFS Figure 4.

fig.04= function (xb=4)

library(GenKern)

x<-allmale2$age

y<-allmale2$suchey

p<-matrix(rep(0,86*6),nc=6)

all<-KernSec(x,range.x=15:100,xb=xb)$yden

for(i in 1:6){

  p[,i]<-KernSec(x[y==i],range.x=15:100,xb=xb)$yden/all

par(mfrow=c(2,3))

library(MASS)

sto<-polr(factor(y)~log(x),method='probit')

sd<-1/as.vector(sto$coeff)

mu<-as.vector(sto$zeta)*sd

plot(15:100,p[,1],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage I')

lines(15:100,1-pnorm(log(15:100),mu[1],sd),lty=2)

plot(15:100,p[,2],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage II')

lines(15:100,pnorm(log(15:100),mu[1],sd)-pnorm(log(15:100),mu[2],sd),lty=2)

plot(15:100,p[,3],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage III')

lines(15:100,pnorm(log(15:100),mu[2],sd)-pnorm(log(15:100),mu[3],sd),lty=2)

plot(15:100,p[,4],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage IV')

lines(15:100,pnorm(log(15:100),mu[3],sd)-pnorm(log(15:100),mu[4],sd),lty=2)

plot(15:100,p[,5],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage V')

lines(15:100,pnorm(log(15:100),mu[4],sd)-pnorm(log(15:100),mu[5],sd),lty=2)

plot(15:100,p[,6],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage VI')

lines(15:100,pnorm(log(15:100),mu[5],sd),lty=2)

fig.04()

JFS Figure 5.

fig.05= function (xb=4)

library(GenKern)

x<-allmale2$age

y<-c(1,2,3,4,5,5)[allmale2$suchey]

p<-matrix(rep(0,86*5),nc=5)

all<-KernSec(x,range.x=15:100,xb=xb)$yden

for(i in 1:5){

  p[,i]<-KernSec(x[y==i],range.x=15:100,xb=xb)$yden/all

par(mfrow=c(2,3))

library(MASS)

sto<-polr(factor(y)~log(x),method='probit')

sd<-1/as.vector(sto$coeff)

mu<-as.vector(sto$zeta)*sd

plot(15:100,p[,1],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage I')

lines(15:100,1-pnorm(log(15:100),mu[1],sd),lty=2)

plot(15:100,p[,2],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage II')

lines(15:100,pnorm(log(15:100),mu[1],sd)-pnorm(log(15:100),mu[2],sd),lty=2)

plot(15:100,p[,3],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage III')

lines(15:100,pnorm(log(15:100),mu[2],sd)-pnorm(log(15:100),mu[3],sd),lty=2)

plot(15:100,p[,4],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage IV')

lines(15:100,pnorm(log(15:100),mu[3],sd)-pnorm(log(15:100),mu[4],sd),lty=2)

plot(15:100,p[,5],type='l',ylim=c(0,1),xlab='Age',ylab='Probability',main='Stage V & VI')

lines(15:100,pnorm(log(15:100),mu[4],sd),lty=2)

fig.05()

JFS Figures 6 to 9 – The default call is for “LA Coroner.” To get other possibilities:

table(allmale2[,1])

and use the other samples, such as: fig.06(‘Korean War’)

fig.06=function (coll='LA Coroner')

library(survival)

b1<-allmale2[allmale2$Collection==coll & allmale2$suchey==1,2]

b2<-allmale2[allmale2$Collection==coll & allmale2$suchey==2,2]

b3<-allmale2[allmale2$Collection==coll & allmale2$suchey==3,2]

b4<-allmale2[allmale2$Collection==coll & allmale2$suchey==4,2]

b5<-allmale2[allmale2$Collection==coll & allmale2$suchey==5,2]

b6<-allmale2[allmale2$Collection==coll & allmale2$suchey==6,2]

par(mfrow=c(3,2))

plot(survfit(Surv(b1)),lty=c(1,1),lwd=c(1,1),xlim=c(15,90),

  xlab='Age',ylab='Survivorship',main=paste('Suchey-Brooks I (n = ',NROW(b1),')',sep=''))

abline(.5,0)

plot(survfit(Surv(b2)),lty=c(1,1),lwd=c(1,1),xlim=c(15,90),

  xlab='Age',ylab='Survivorship',main=paste('Suchey-Brooks II (n = ',NROW(b2),')',sep=''))

abline(.5,0)

plot(survfit(Surv(b3)),lty=c(1,1),lwd=c(1,1),xlim=c(15,90),

  xlab='Age',ylab='Survivorship',main=paste('Suchey-Brooks III (n = ',NROW(b3),')',sep=''))

abline(.5,0)

plot(survfit(Surv(b4)),lty=c(1,1),lwd=c(1,1),xlim=c(15,90),

  xlab='Age',ylab='Survivorship',main=paste('Suchey-Brooks IV (n = ',NROW(b4),')',sep=''))

abline(.5,0)

plot(survfit(Surv(b5)),lty=c(1,1),lwd=c(1,1),xlim=c(15,90),

  xlab='Age',ylab='Survivorship',main=paste('Suchey-Brooks V (n = ',NROW(b5),')',sep=''))

abline(.5,0)

plot(survfit(Surv(b6)),lty=c(1,1),lwd=c(1,1),xlim=c(15,90),

  xlab='Age',ylab='Survivorship',main=paste('Suchey-Brooks VI (n = ',NROW(b6),')',sep=''))

abline(.5,0)

q1<-c(NROW(b1),round(quantile(b1,p=c(.025,.25,.5,.75,.975))))

q2<-c(NROW(b2),round(quantile(b2,p=c(.025,.25,.5,.75,.975))))

q3<-c(NROW(b3),round(quantile(b3,p=c(.025,.25,.5,.75,.975))))

q4<-c(NROW(b4),round(quantile(b4,p=c(.025,.25,.5,.75,.975))))

q5<-c(NROW(b5),round(quantile(b5,p=c(.025,.25,.5,.75,.975))))

q6<-c(NROW(b6),round(quantile(b6,p=c(.025,.25,.5,.75,.975))))

print(rbind(q1,q2,q3,q4,q5,q6))

fig.06()

JFS Figure 10

fig.10=function ()

library(MASS)

#========================

coll='Balkan'

x2<-allmale2[allmale2$Collection!=coll,]

x<-x2$age

y<-x2$suchey

sto<-polr(factor(y)~log(x),method='probit')

sd<-1/as.vector(sto$coeff)

mu<-as.vector(sto$zeta)*sd

x2<-allmale2[allmale2$Collection==coll,]

x<-x2$age

y<-x2$suchey

sto<-polr(factor(y)~log(x),method='probit')

sd2<-1/as.vector(sto$coeff)

mu2<-as.vector(sto$zeta)*sd

#========================

par(mfrow=c(3,2))

labs=c('I / II','II / III','III / IV','IV / V','V / VI')

x<-seq(0,100,by=.1)

for(i in 1:3){

plot(x,dlnorm(x,mu2[i],sd2),type='l',lty=2,xlab='Age',ylab='Density',main=labs[i])

lines(x,dlnorm(x,mu[i],sd))

for(i in 4:5){

plot(x,dlnorm(x,mu[i],sd),type='l',lty=1,xlab='Age',ylab='Density',main=labs[i])

lines(x,dlnorm(x,mu2[i],sd2),lty=2)

fig.10()

JFS Figure 11

fig.11=function ()

{

library(MASS)

#========================

coll='Balkan'

x2<-allmale2[allmale2$Collection!=coll,]

x<-x2$age

y<-x2$suchey

print(NROW(x))

sto<-polr(factor(y)~log(x),method='probit')

sd<-1/as.vector(sto$coeff)

mu<-as.vector(sto$zeta)*sd

x2<-allmale2[allmale2$Collection==coll,]

x<-x2$age

y<-x2$suchey

print(NROW(x))

sto<-polr(factor(y)~log(x),method='probit')

sd2<-1/as.vector(sto$coeff)

mu2<-as.vector(sto$zeta)*sd2

#========================

ptran<-function(t,ii){

if(ii==1){p<-1-pnorm(log(t),mu[1],sd)

return(p)}

if(ii==6){p<-pnorm(log(t),mu[5],sd)

return(p)}

p<-pnorm(log(t),mu[ii-1],sd)-pnorm(log(t),mu[ii],sd)

return(p)

}

ptran2<-function(t,ii){

if(ii==1){p<-1-pnorm(log(t),mu2[1],sd2)

return(p)}

if(ii==6){p<-pnorm(log(t),mu2[5],sd2)

return(p)}

p<-pnorm(log(t),mu2[ii-1],sd2)-pnorm(log(t),mu2[ii],sd2)

return(p)

}

#========================

par(mfrow=c(3,2))

labs=c('Stage I','Stage II','Stage III','Stage IV','Stage V','Stage VI')

x<-seq(.1,100,by=.1)

for(i in 1:6){

y<-ptran(x,i)

norm<-optimize(ptran,c(.1,100),maximum=T,ii=i)$obj

plot(x,y/norm,type='l',xlab='Age',ylab='Norm Likelihood',main=labs[i])

abline(0.7965477,0,lty=3)

abline(0.1465001,0,lty=3)

y<-ptran2(x,i)

norm<-optimize(ptran2,c(.1,100),maximum=T,ii=i)$obj

lines(x,y/norm,lty=2)

}

fig.11()

JFS Figure 12 – This includes a Makeham pdf for Komar’s data. And to run this you will need her data set, so copy and paste:

komar.dat=structure(list(age=c(20,25,30,35,40,45,50,55,60,65,70),

No.Missing=c(1050,999,896,877,822,748,546,536,517,300,328)),.Names=c("age","No.Missing"),

row.names=c("1","2","3","4","5","6","7","8","9","10","11"), class = "data.frame")

fig.12=function()

{

par(mfrow=c(2,1))

makeham <- function(t) {

x <- c(0.023907712, 0.005382351, 0.063442251)

a2 <- x[1]

a3 <- x[2]

b3 <- x[3]

shift <- 20

h.t <- a2 + a3 * exp(b3 * (t - shift))

S.t <- exp(-a2 * (t - shift) + a3/b3 * (1 - exp(b3 *

(t - shift))))

return(S.t * h.t)

}

lftable <- function() {

n.ints <- NROW(komar.dat)

n.deaths <- sum(komar.dat$No.Missing)

dx <- komar.dat$No.Missing/n.deaths

dx <- dx/5

x <- komar.dat$age

x <- c(rep(x, each = 2), 75, 75, 20)

y <- c(0, rep(dx, each = 2), 0, 0)

plot(x, y, type = "l", xlab = "Age", ylab = "Frequency",main='Komar (2003)')

polygon(x, y, density = 20)

}

lftable()

lines(20:80, makeham(20:80), lwd = 2)

age=allmale2$age[allmale2$Coll=='Balkan']

age=age[age>=20 & age<=75]

sto=hist(age,breaks=c(komar.dat$age,75),plot=F)

lftable2 <- function() {

n.ints <- NROW(komar.dat)

n.deaths <- sum(sto$counts)

dx <- sto$counts/n.deaths

dx <- dx/5

x <- komar.dat$age

x <- c(rep(x, each = 2), 75, 75, 20)

y <- c(0, rep(dx, each = 2), 0, 0)

plot(x, y, type = "l", xlab = "Age", ylab = "Frequency",main='Current Study')

polygon(x, y, density = 20)

}

lftable2()

}

fig.12()

JFS Figure 13

fig.13=function ()

{

library(survival)

makeham<-function(t,a3,b3)

{

a2<-0

shift<-17

h.t<-a2+a3*exp(b3*(t-shift))

S.t<-exp(-a2*(t-shift)+a3/b3*(1-exp(b3*(t-shift))))

return(S.t)

}

ages=allmale2$age[allmale2$Coll=='Balkan']

plot(survfit(Surv(ages)~1),lty=c(0,1),lwd=c(1,1),xlim=c(1,100),

main='Gompertz (212 Individuals)',xlab='Age',ylab='Survivorship')

##################################

jfs.optim<-function (x)

{

makeham2<-function(t,x)

{

a2<-0

a3<-x[1]

b3<-x[2]

shift<-17

h.t<-a2+a3*exp(b3*(t-shift))

S.t<-exp(-a2*(t-shift)+a3/b3*(1-exp(b3*(t-shift))))

return(S.t*h.t)

}

t<-ages

lnlk<-0

for(i in 1:212){

lnlk<-lnlk+log(makeham2(t[i],x))

}

return(lnlk)

}

##################################

sto<-optim(c(.01,.05),jfs.optim,control=list(fnscale=-1),hessian=T)

library(MASS)

myenv <- new.env()

assign("a3",sto$par[1], env = myenv)

assign("b3", sto$par[2], env = myenv)

assign("t", 17:100, env = myenv)

dt<-numericDeriv(quote(makeham(17:100, a3, b3)), c("a3", "b3"), myenv)

est<-dt[1:84]

gr<-attr(dt,'gradient')

V<-ginv(-sto$hess)

lo<-0

hi<-0

for(i in 1:84){

se<-sqrt(as.vector(gr[i,]%*%V%*%gr[i,]))

lo[i]<-est[i]-1.96*se

hi[i]<-est[i]+1.96*se

if(lo[i]<0) lo[i]<-0

if(hi[i]>1) hi[i]<-1

}

lines(17:100,hi,lwd=2)

lines(17:100,lo,lwd=2)

}

par(mfrow=c(1,1))

fig.13()

JFS Figure 14

fig.14=function ()

{

library(survival)

age=allmale2$age[allmale2$Coll=='Thai']

makeham<-function(t,a3,b3)

{

a2<-0

shift<-20

h.t<-a2+a3*exp(b3*(t-shift))

S.t<-exp(-a2*(t-shift)+a3/b3*(1-exp(b3*(t-shift))))

return(S.t)

}

plot(survfit(Surv(age)),lty=c(0,1),lwd=c(1,1),xlim=c(1,100),

main='Thai Males',xlab='Age',ylab='Survivorship')

##################################

jfs.optim<-function (x)

{

makeham2<-function(t,x)

{

a2<-0

a3<-x[1]

b3<-x[2]

shift<-20

h.t<-a2+a3*exp(b3*(t-shift))

S.t<-exp(-a2*(t-shift)+a3/b3*(1-exp(b3*(t-shift))))

return(S.t*h.t)

}

t<-age

lnlk<-0

for(i in 1:37){

lnlk<-lnlk+log(makeham2(t[i],x))

}

return(lnlk)

}

##################################

sto<-optim(c(.01,.05),jfs.optim,control=list(fnscale=-1),hessian=T)

library(MASS)

myenv <- new.env()

assign("a3",sto$par[1], env = myenv)

assign("b3", sto$par[2], env = myenv)

assign("t", 20:100, env = myenv)

dt<-numericDeriv(quote(makeham(20:100, a3, b3)), c("a3", "b3"), myenv)

est<-dt[1:86]

gr<-attr(dt,'gradient')

V<-ginv(-sto$hess)

lo<-0

hi<-0

for(i in 1:81){

se<-sqrt(as.vector(gr[i,]%*%V%*%gr[i,]))

lo[i]<-est[i]-1.96*se

hi[i]<-est[i]+1.96*se

if(lo[i]<0) lo[i]<-0

if(hi[i]>1) hi[i]<-1

}

lines(20:100,hi,lwd=2)

lines(20:100,lo,lwd=2)

return(sto)

}

fig.14()

JFS Figure 15

fig.15=function ()

{

plot(allmale2$age[allmale2$Coll=='Balkan'],

jitter(allmale2$suchey[allmale2$Coll=='Balkan'],.4),

xlab='Age',ylab='Stage',main='50% Coverage (Balkans)')

rect(17,1.1,22.42,1.2,density=20)

rect(20.47,2.1,31.36,2.2,density=20)

rect(24.73,3.1,37.59,3.2,density=20)

rect(31.27,4.1,46.81,4.2,density=20)

rect(43.45,5.1,60.20,5.2,density=20)

rect(55.52,6.1,71.46,6.2,density=20)

m<-c(18.5,23.4,28.7,35.2,45.6,61.2)

s<-c(2.1,3.6,6.5,9.4,10.4,12.2)

rect(m[1]-.674*s[1],.8,m[1]+.674*s[1],.9,density=20,ang=-45)

rect(m[2]-.674*s[2],1.8,m[2]+.674*s[2],1.9,density=20,ang=-45)

rect(m[3]-.674*s[3],2.8,m[3]+.674*s[3],2.9,density=20,ang=-45)

rect(m[4]-.674*s[4],3.8,m[4]+.674*s[4],3.9,density=20,ang=-45)

rect(m[5]-.674*s[5],4.8,m[5]+.674*s[5],4.9,density=20,ang=-45)

rect(m[6]-.674*s[6],5.8,m[6]+.674*s[6],5.9,density=20,ang=-45)

legend(60,1.5,c('Trans. Analysis','Suchey-Brooks'),bty='n',

density=c(20,20),angle=c(45,-45))

}

fig.15()

JFS Figure 16

fig.16=function ()

{

plot(allmale2$age[allmale2$Coll=='Thai'],

jitter(allmale2$suchey[allmale2$Coll=='Thai'],.5),

xlab='Age',ylab='Stage',main='50% Coverage (Thailand)',

pch=19,xlim=c(0,100),ylim=c(1,6))

rect(20,1,25.44,1.1,density=20)

rect(21.89,2,32.73,2.1,density=20)

rect(26.41,3,40.02,3.1,density=20)

rect(34.92,4,51.76,4.1,density=20)

rect(49.06,5,65.8,5.1,density=20)

rect(60.58,6,75.31,6.1,density=20)

m<-c(18.5,23.4,28.7,35.2,45.6,61.2)

s<-c(2.1,3.6,6.5,9.4,10.4,12.2)

rect(m[1]-.674*s[1],.9,m[1]+.674*s[1],1,density=20,ang=-45)

rect(m[2]-.674*s[2],1.9,m[2]+.674*s[2],2,density=20,ang=-45)

rect(m[3]-.674*s[3],2.9,m[3]+.674*s[3],3,density=20,ang=-45)

rect(m[4]-.674*s[4],3.9,m[4]+.674*s[4],4,density=20,ang=-45)

rect(m[5]-.674*s[5],4.9,m[5]+.674*s[5],5,density=20,ang=-45)

rect(m[6]-.674*s[6],5.9,m[6]+.674*s[6],6,density=20,ang=-45)

legend(60,1.5,c('Trans. Analysis','Suchey-Brooks'),bty='n',

density=c(20,20),angle=c(45,-45))

}

fig.16()

JFS Figures 17 to 20 – To speed things along this is only set for 10 replicates, while the paper used 1,000. Also note that “coll” is by default “LA Coroner.” Use other sample names to get Figures 18, 19, and 20.

fig.17=function(reps=10,coll='LA Coroner')

{

x2<-allmale2[allmale2$Collection!=coll,]

x<-x2$age

y<-x2$suchey

library(MASS)

sto<-polr(factor(y)~log(x),method='probit')

sd<-1/as.vector(sto$coeff)

mu<-as.vector(sto$zeta)*sd

ptran<-function(t,i){

if(i==1){p<-1-pnorm(log(t),mu[1],sd)

return(p)}

if(i==6){p<-pnorm(log(t),mu[5],sd)

return(p)}

p<-pnorm(log(t),mu[i-1],sd)-pnorm(log(t),mu[i],sd)

return(p)

}

x2<-allmale2[allmale2$Collection==coll,]

x<-x2$age

y<-x2$suchey

N<-NROW(x2)

age<-x2$age

stage<-x2$suchey

#pi<-as.vector(table(stage))/N

pi<-hist(stage,0:6,plot=F)$counts/N

llr<-rep(0,N)

for(i in 1:N){

llr[i]<-ptran(age[i],stage[i])/pi[stage[i]]

}

scale<-max(log10(llr))

plot((1:N)/N*100,(sort(log10(llr),dec=T)),

type='l',lwd=2,ylab='Log Likelihood Ratio',xlab='Percentile',

main=coll)

prop<-sum(llr>1)/N*100

print(prop)

lines(c(prop,prop),c(-1000,0),lwd=2)

llr2<-rep(0,N)

for(j in 1:reps){

age<-sample(x2$age)

llr<-rep(0,N)

for(i in 1:N){

llr[i]<-ptran(age[i],stage[i])/pi[stage[i]]

}

llr2<-llr2+sort(log10(llr),dec=T)/reps

}

lines((1:N)/N*100,(llr2),lty=2)

abline(0,0)

legend(60,scale,bty='n',lty=c(1,2),lwd=c(2,1),legend=c('Actual','Permuted'))

prop<-sum(llr2>0)/N*100

print(prop)

lines(c(prop,prop),c(-1000,0),lty=2)

}

fig.17()