Calculations for Steadman, Adams, and Konigsberg (2006?)

       First, “R” is your friend, so:

http://cran.r-project.org/

or, link for Rweb:  http://spider.stat.umn.edu/Rweb/Rweb.general.html

Most of the code listed below can be submitted to “Rweb,” which is a quick way to “try before you buy” (but “R” is a free download anyhow!).  The first block below (“Terry Data Analysis”) cannot be submitted to “Rweb” because the code requires file creation.


 

Terry Data Analysis

Below are two lines to get the data and save into an “R” object called “terry.”  Just copy these and paste into the R GUI (be sure to hit “enter” after the second line).

 

  download.file('http://konig.la.utk.edu/terry.txt','terry.txt')
  terry<-dget('terry.txt')
 

Now copy the bold text below and paste into the R GUI to get the MLEs for pubic symphysis transition parameters.  Note: this model already starts close to the MLEs.

 
  pube<-function() 
  {
  attach(terry)
  age<-log(Age)
  indic<-Pube3
 
  lnlk<-function(xpars) 
  {
  sd<-xpars[1]
  mu1<-xpars[2]
  mu2<-xpars[3]
  LLK<-0
  for(i in 1:746){
  if(indic[i]==1){
  f<-1-pnorm(age[i],mu1,sd)
  }
  if(indic[i]==2){
  f<-pnorm(age[i],mu1,sd)-pnorm(age[i],mu2,sd)
  }
  if(indic[i]==3){
  f<-pnorm(age[i],mu2,sd)
  }
  if(f<=0.0) f<-1.e-20
  LLK<-LLK+log(f)
  }
  return(LLK)
  }
  out<-optim(c(.45, 3.45, 3.92),lnlk,control=list(fnscale=-1))
  return(round(out$par,4))
  }
  pube()
 

Now copy below and paste into R GUI to get the MLEs for the auricular surface transition parameters.

 
  auric<- function() 
  {
  attach(terry)
  age<-log(Age)
  indic<-Aur3
  lnlk<-function(xpars) 
  {
  sd<-xpars[1]
  mu1<-xpars[2]
  mu2<-xpars[3]
  LLK<-0
  for(i in 1:746){
  if(indic[i]==1){
  f<-1-pnorm(age[i],mu1,sd)
  }
  if(indic[i]==2){
  f<-pnorm(age[i],mu1,sd)-pnorm(age[i],mu2,sd)
  }
  if(indic[i]==3){
  f<-pnorm(age[i],mu2,sd)
  }
  if(f<=0.0) f<-1.e-20
  LLK<-LLK+log(f)
  }
  return(LLK)
  }
  out<-optim(c(0.45, 3.45, 3.92),lnlk,control=list(fnscale=-1))
  return(round(out$par,4))
  }
  auric()
 
 
 
 
 


 
 
 
 

AAPA Figure 1. – This has no data so you can copy and submit to “Rweb” or paste into your own RGUI.

  aapa.fig1<-function()
  {
  sdev<-0.4536 
  mu1<- 3.4552
  mu2<- 3.9201
  plot(1:120,dnorm(log(1:120),mu1,sdev),type='l',xlab='Age at Transition',ylab='Density')
  lines(1:120,dnorm(log(1:120),mu2,sdev),lty=2)
  legend(75,.8,legend=c('Stage 1 to Stage 2','Stage 2 to Stage 3'),cex=1,lty=c(1,2),bty='n')
  }
  aapa.fig1()
 
 
 
 
 


 
 
 
 

AAPA Figure 2. – Again no data, so you can copy and submit to “Rweb” or paste into your own RGUI.

  aapa.fig2<-function()
  {
  sdev<-0.4536 
  mu1<- 3.4552
  mu2<- 3.9201
  plot(1:100,1-pnorm(log(1:100),mu1,sdev),type='l',xlab='Age',ylab='"Survival"',
     ylim=c(0,1))
  lines(1:100,1-pnorm(log(1:100),mu2,sdev))
  xs<-seq(44,46,.1)
  y1<-1-pnorm(log(xs),mu1,sdev)
  x<-c(44,xs,46,44)
  y<-c(0,y1,0,0)
  polygon(x,y)
  y2<-1-pnorm(log(xs),mu2,sdev)
  x<-c(xs,rev(xs),44)
  y<-c(y1,rev(y2),y1[1])
  polygon(x,y,density=40)
  x<-c(xs,46,44,44)
  y<-c(y2,1,1,y2[1])
  polygon(x,y)
  text(20,.4,'Stage 1')
  text(36,.6,'Stage 2')
  text(60,.75,'Stage 3')
  }
  aapa.fig2()
  
 
 
 
 


 
 
 
 

AAPA Figure 3. – And once more no data, so you can copy and submit to “Rweb” or paste into your own RGUI.

  aapa.fig3<-function()
  {
  sdev<-0.4536 
  mu1<- 3.4552
  mu2<- 3.9201
  x<-1:120
  y1<-1-pnorm(log(1:120),mu1,sdev)
  y2<-1-pnorm(log(1:120),mu2,sdev)
  plot(x,y1,type='l',xlab='Age',ylab='Density',
     ylim=c(0,1.25),axes=F)
  axis(2,at=c(0,.2,.4,.6,.8,1))
  axis(1)
  box(bty="o")
  legend(80,1.3,legend=c('f(Stage 1)','f(Stage 2)','f(Stage 3)'),cex=1,lty=c(1,2,3),bty='n')
  lines(x,y2-y1,lty=2)
  lines(x,1-y2,lty=3)
  xs<-seq(44,46,.1)
  y1<-1-pnorm(log(xs),mu1,sdev)
  y2<-1-pnorm(log(xs),mu2,sdev)
  y1<-y2-y1
  x<-c(44,xs,46,44)
  y<-c(0,y1,0,0)
  polygon(x,y,density=40)
  rect(44,0,46,1)
  }
  aapa.fig3()
 
 
 
 
 


 
 
 
 

“Using numerical integration in “R” we find the shaded area to be equal to 0.759…”

Copy, paste, and “submit.”  Note that you can specify the indicator state (1, 2, or 3) and the bone (1 for auricular surface and 2 for pubic symphysis).  If you take the default you’ll get a value of 0.759.

 

  integ1<-function (indic=2,bone=2) 
  {
# bone = 1 is auricular
# bone = 2 is pubic symphysis
   trans<-function (age) 
   {
   lage=log(age)
   sto1<-c(1/.4286,3.8129/.4286,3.9171/.4286)
   sto2<-c(1/.4536,3.4552/.4536,3.9201/.4536)
   parms<-rbind(sto1,sto2)
# Stage 1
   if(indic==1){
   L<-0
   R<-pnorm(parms[bone,2]-parms[bone,1]*lage)}
# Stage 2
   if(indic==2){
   L<-pnorm(parms[bone,2]-parms[bone,1]*lage)
   R<-pnorm(parms[bone,3]-parms[bone,1]*lage)}
# Stage 3
   if(indic==3){
   L<-pnorm(parms[bone,3]-parms[bone,1]*lage)
   R<-1}
   return(R-L)
   }
   tot<-function(tt){
   return(trans(tt))
   }
   p<-integrate(tot,44,46)$value
   print(round(p,3))
   }
   integ1()
 
 
 
 
 


 
 
 
 

NCIC Data Analysis

Copy, paste, and “submit” everything below in order to estimate the Gompertz-Makhem parameters

 

  ncic<-structure(list(open = c(20, 21, 22, 30, 40, 50, 60, 70, 80, 90), 
  close = c(21, 22, 30, 40, 50, 60, 70, 80, 90, 100), F = c(1221, 
  962, 3459, 3384, 2598, 1131, 492, 383, 272, 109), M = c(610, 
  478, 3035, 4296, 4332, 2529, 1373, 1044, 663, 257), tot = c(1831, 
  1440, 6494, 7680, 6930, 3660, 1865, 1427, 935, 366)), .Names = c("open", 
  "close", "F", "M", "tot"), row.names = c("21", "22", "23", "24", 
  "25", "26", "27", "28", "29", "30"), class = "data.frame")
 
  fit.ncic<-function(x){
  open<-ncic$open
  close<-ncic$close
  Dx<-ncic$tot
 
  makeham<-function (x,tt)
  {
  a2=x[1]
  a3=x[2] 
  b3=x[3]
  t<-tt-20
  h.t<-a2+a3*exp(b3*t)
  S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t)))
  return(S.t)
  }
 
  lnlk<-as.vector(log(makeham(x,open)-makeham(x,close))%*%Dx)
  return(lnlk)
  }
 
  optim(c(.01,.01,.01),fit.ncic,control=list(fnscale=-1)) 

AAPA Figure 4. – Copy, paste, and “submit” everything below in order to plot results from the previous analysis

 
  par(mfrow=c(1,2))
 
  ncic<-structure(list(open = c(20, 21, 22, 30, 40, 50, 60, 70, 80, 90),
  close = c(21, 22, 30, 40, 50, 60, 70, 80, 90, 100), F = c(1221, 
  962, 3459, 3384, 2598, 1131, 492, 383, 272, 109), M = c(610, 
  478, 3035, 4296, 4332, 2529, 1373, 1044, 663, 257), tot = c(1831, 
  1440, 6494, 7680, 6930, 3660, 1865, 1427, 935, 366)), .Names = c("open", 
  "close", "F", "M", "tot"), row.names = c("21", "22", "23", "24", 
  "25", "26", "27", "28", "29", "30"), class = "data.frame")
 
  makeham2<-function(tt)
  {
  x=c(0.03044,0.0066,0.0392)
  a2=x[1]
  a3=x[2] 
  b3=x[3]
  t<-tt-20
  h.t<-a2+a3*exp(b3*t)
  S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t)))
  return(S.t)
  }
 
  plot.ncic<-function () 
  {
  plot(20:100,makeham2(tt=20:100),type='l', xlab='Age',ylab='Survivorship')
  Dx<-ncic$tot
  open<-ncic$open
  n.ints<-NROW(ncic)
  print(n.ints)
  tot<-ncic$tot
  n.deaths<-sum(tot)
  lx<-c(n.deaths,n.deaths-cumsum(tot))/n.deaths
  print(lx)
  points(open,lx[1:n.ints])
  }
 
  plot.ncic()
 
  plot.ncic2<-function () 
  {
  make3<-function(tt)
  {
  a2=0.03044
  a3=0.0066
  b3=0.0392
  t<-tt-20
  h.t<-a2+a3*exp(b3*t)
  S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t)))
  return(S.t*h.t)
  }
 
  Dx<-ncic$tot
  open<-ncic$open
  close<-ncic$close
  dx<-(Dx/sum(Dx))/(close-open)
  plot(c(rep(open,each=2),100,100),c(0,rep(dx,each=2),0),
  ylim=c(0,.06),type='l',lty=2,xlim=c(20,100),
  xlab='Age',ylab='Density')
  lines(20:100,make3(20:100))
  }
 
  plot.ncic2()

“From equation 6 we find the probability that an individual would be in stage 1 is 0.384…”

Copy, paste, and “submit.”  Note that you can specify the indicator state (1, 2, or 3) and the bone (1 for auricular surface and 2 for pubic symphysis).  If you take the default you’ll get a value of 0.384.

 

integ2<-function (indic=1,bone=2) 
{
#
#  bone = 1 is auricular
#  bone = 2 is pubic symphysis
#
trans<-function (age) 
{
lage=log(age)
sto1<-c(1/.4286,3.8129/.4286,3.9171/.4286)
sto2<-c(1/.4536,3.4552/.4536,3.9201/.4536)
parms<-rbind(sto1,sto2)
# Stage 1
if(indic==1){
L<-0
R<-pnorm(parms[bone,2]-parms[bone,1]*lage)}
# Stage 2
if(indic==2){
L<-pnorm(parms[bone,2]-parms[bone,1]*lage)
R<-pnorm(parms[bone,3]-parms[bone,1]*lage)}
# Stage 3
if(indic==3){
L<-pnorm(parms[bone,3]-parms[bone,1]*lage)
R<-1}
return(R-L)
}
 
make3<-function(tt)
{
a2=0.03044
a3=0.0066
b3=0.0392
t<-tt-20
h.t<-a2+a3*exp(b3*t)
S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t)))
return(S.t*h.t)
}
 
tot<-function(tt){
return(trans(tt)*make3(tt))
}
p<-integrate(tot,20,200)$value
print(round(p,3))
}
 
integ2()
 
  
 
 
 
 


 
 
 
 

AAPA Figure 5. – Copy and submit to “Rweb” or paste into your own RGUI.

  plot.pden<-function (indic=3) 
  {
   par(mfrow=c(1,1))
#============================
  make3<-function(tt)
  {
  a2=0.03044
  a3=0.0066
  b3=0.0392
  t<-tt-20
  h.t<-a2+a3*exp(b3*t)
  S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t)))
  return(S.t*h.t)
  }
#===========================
  trans<-function (age) 
  {
  lage=log(age)
  sto<-c(1/.4538243,3.4553671/.4538243,3.9202782/.4538243)
# Stage 1
  if(indic==1){
  L<-0
  R<-pnorm(sto[2]-sto[1]*lage)}
# Stage 2
  if(indic==2){
  L<-pnorm(sto[2]-sto[1]*lage)
  R<-pnorm(sto[3]-sto[1]*lage)}
# Stage 3
  if(indic==3){
  L<-pnorm(sto[3]-sto[1]*lage)
  R<-1}
  return(R-L)
  }
#==========================
 
  tot<-function(tt){
  return(make3(tt)*trans(tt))
  }
  return(tot(20:100))
  }
 
  plot.pden2<-function () 
  {
# Integral calculated using "integrate" from 20 to 200 years
  y1<-plot.pden(1)/0.3840
  y2<-plot.pden(2)/0.2916
  y3<-plot.pden(3)/0.3244
  plot(20:100,y1,type='l',xlab='Age',ylab='Density')
  lines(20:100,y2,lty=2)
  lines(20:100,y3,lty=3)
  legend(75,.06,legend=c('Stage 1','Stage 2','Stage 3'),cex=1,lty=1:3,bty='n')
  }
 
  plot.pden2()
 
 
 


 
 

Bivariate (pubic symphysis and auricular surface) Analyses

While it would be nice to do this in “R,” there currently seems to be a certain amount of “hinkeyness” in the code.  Consequently, we do the analyses in FORTRAN (ANSI 77 version).  We compiled the source code and ran under Windows XP using GNU g77.  GNU gcc should also be able to compile, but we haven’t tried this – yet.  The source code can be obtained from:

 

http://konig.la.utk.edu/4446.txt

http://konig.la.utk.edu/atlarge.txt