Calculations for Steadman, Adams, and Konigsberg (2006?)
First, “R” is your friend, so:
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