Annals of Human Biology

This page gives the "R" code for analyses reported  in Konigsberg (2015).  Because of some code "tweaking," a portion of the results from stochastic methods differ slightly from those reported in the AHB article.

If you do not already have "R" you can get it from:

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

You will also need to download some specific "R" libraries.  If you are unfamiliar with this, see the link below.

 
Getting Some Data

First things first, anything shown in fixed-width green text is intended to be copied and pasted into your R "Gui" (graphical user interface).

Click on the link below to get ectocranial suture closure data for 363 individuals from McKern and Stewart's (1957) Korean War dead study.

KoreaSut.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 save the data to an object, say (copy and paste the line below to your "R" Gui):

And now get comparable data from the Terry Collection:

TerrySut.csv

Sutures were scored as: 1 = "open," 2 = "minimal," 3 = "significant," and 4 = "closed."  The sutures scored here (in the order in the files) are:

  1. Mid-lambdoid
  2. Lambda
  3. Anterior sagittal
  4. Bregma
  5. Sphenofrontal
 
Getting the Data into One Big, Happy R Object

If you take a look-see at the Korean War data, you may notice that the day, month, and year for birth and death are entered, but the age-at-death is not calculated.  After getting the "chron" and "date" libraries (see the transition analysis tutorial to see how to get libraries) you can use the following script to calculate ages-at-death:

 
get.BD.DD.age = function () 
{
library(chron)
library(date)

Korea.sut = read.csv('KoreaSut.csv')
decimal.date = function (MoDayYear) 
{
sto=date.mdy(MoDayYear)
dec.year=(mdy.date(month=sto$month,day=sto$day,year=1970)-as.date(1/1/1970)+1)/365
dec.year=dec.year+sto$year
return(as.numeric(noquote(format(dec.year,digits=8))))
}

N=NROW(Korea.sut)
BD=decimal.date(julian(x=Korea.sut[,2],d=Korea.sut[,3],y=1900+Korea.sut[,4]))
DD=decimal.date(julian(x=Korea.sut[,5],d=Korea.sut[,6],y=1900+Korea.sut[,7]))
age=round(DD-BD,3)
Korea.sut2=cbind(rep('K',N),Korea.sut[,1],age,Korea.sut[,8:12])
colnames(Korea.sut2)=c('Coll','ID',colnames(Korea.sut2)[-(1:2)])
return(rbind(Korea.sut2,read.csv('TerrySut.csv')))
}
 
And then:

sutures = get.BD.DD.age()
"Assuming a Gompertz mortality model (Finch and Pike, 1996, Wilson, 1993, Pletcher, 1999) starting at age 15 years, the a3 and b3 parameters for the combined Korean War dead and Terry samples are 0.025 and 0.018, respectively."

So we will need a script to estimate these parameters (be sure to download the library "alabama"):

fit.Gomp = function(){
library(alabama)
t=sutures[,3]
N = length(t)

Gomp = function (x,age)
{
         a3 = x[1]
         b3 = x[2]
         shift = 15
         h.t = a3*exp(b3*(age-shift))
         S.t = exp(a3/b3*(1-exp(b3*(age-shift))))
         return(log(S.t*h.t))
}

lnlk = function(x){
  zlnlk = 0
  for(i in 1:N) zlnlk = zlnlk + Gomp(x,t[i])
  return(-zlnlk)
}

hin=function(par){
  h=NA
  h[1] = par[1]
  h[2] = par[2]
  return(h)
}

start.theta = c(0.001,0.1)
constrOptim.nl(par=start.theta,fn=lnlk,hin=hin,control.outer=list(trace=F))
}


fit.Gomp()

$par
[1] 0.02480333 0.01791355

$value
[1] 4864.281

$countsa3 and b3
function gradient
     164       23

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 3

$barrier.value
[1] 2.064618e-06

Table 1

This Table looks at a "composite" score (sum) from the five suture sites.  For reasons apparent later in the article, the first four suture sites were collapsed to a three-stage scoring where 1 = "open," 2 = "minimal," and 3 = "significant" or "closed."  With that scoring the minimum composite score is five (all five sutures open) and the maximum is sixteen (first four sutures "significant" or "closed" and last suture "closed."  The first part of the Table gives observed statistics for the actual data, where "count" is the number of cases within scores, and "mean," "s.d.," "median," and "range" are the statistics for age within score.  The second part of the Table gives the same statistics, but on simulated data where the sutures close on the same "schedule" as in the actual data but the age-at-death structure is radically different (a3 = 0.01 and b3 = 0.06).  The whole point of this Table is to show how "age within stage" statistics depend in part on the age structure of the sample.  This is a point often missed in the literature (for example Roberts et al. 2008, Mitchell et al. 2009, and Peiris et al. 2009).  Below are the results for Table 1 out of the "R" code.

Actual Data
 Composite Count  Mean  s.d. Median         Range
       5-6   313  26.5  12.1   22.0  14.0 -  95.0
       7-8   150  33.7  15.7   28.0  17.7 -  86.0
      9-10   116  41.1  17.3   36.0  19.0 - 102.0
     11-12   162  47.0  19.7   45.0  17.0 - 101.0
     13-14   169  49.0  18.2   48.0  17.6 -  95.0
     15-16   242  56.8  17.6   55.5  21.7 -  95.0


Simulated Data
 Composite Count Mean s.d. Median       Range
       5-6   277 33.3 12.5   31.1 15.1 - 73.6
       7-8   151 37.9 11.8   37.8 16.3 - 68.7
      9-10   125 41.4 12.7   41.2 15.6 - 73.7
     11-12   160 43.4 12.3   43.3 15.4 - 75.5
     13-14   166 46.3 12.6   45.8 16.3 - 82.5
     15-16   273 50.0 11.1   50.7 19.3 - 77.0
recode = function ()
{
dat = sutures
rec = c(1:3,3)
sutures2 = dat[,1:3]
for(i in 4:7) sutures2 = cbind(sutures2,rec[dat[,i]])
sutures2 = cbind(sutures2,dat[,8])
colnames(sutures2) = colnames(dat)
return(sutures2)
}

sutures2 = recode()


to.symm = function(x)
{
#  A utility function to take the vector of correlation coefficients and place into a
#  symmetric matrix
k=length(x)
n.dim=(sqrt(8*k+1)+1)/2
sym.mat=diag(n.dim)
ip=0
for(i in 1:(n.dim-1)){
for(j in (i+1):n.dim){
  ip=ip+1
  sym.mat[i,j]=sym.mat[j,i]=x[ip]
}
}
return(sym.mat)
}


sim.Gomp = function (N=1000,a3=0.01,b3=0.06) 
{
# Simulate mortality from Gompertz function, starting at age 15 years
U=runif(N,0,1)
return(log(1-b3*log(U)/a3)/b3+15)
}


Table.01 = function (N=1152)
{
sim.sutures=function (N=1152)
{
library(mvtnorm)
set.seed(123456)
theta=c(5.04467057878564, 5.77527005349954, 1.45692404535338, 4.72064443260498,
5.32468651156025, 1.41421360354464, 3.55510517393914, 4.08555277396893,
1.14203006862722, 5.44197491169678, 6.11694986309359, 1.55625108460126,
6.36874788009504, 6.7132844292257, 7.19071270808259, 1.73393538627762,
0.837295888826435, 0.698777441913304, 0.656195525416774, 0.466261783209169,
0.785390250744974, 0.707835949910551, 0.440554938604524, 0.825120067898404,
0.506050357944747, 0.563940211861995)
thresh=matrix(NA,nr=5,nc=5)
thresh[1,]=c(-100,theta[1:2],100,NA)
thresh[2,]=c(-100,theta[4:5],100,NA)
thresh[3,]=c(-100,theta[7:8],100,NA)
thresh[4,]=c(-100,theta[10:11],100,NA)
thresh[5,]=c(-100,theta[13:15],100)
top=c(rep(4,4),5)
slope=theta[c(3,6,9,12,16)]
r=to.symm(theta[17:26])
ages=sim.Gomp(N)

it=matrix(NA,nr=N,nc=5)
for(j in 1:N){
mu=log(ages[j])*slope
sto=rmvnorm(1,mu,r)
for(i in 1:5){
 it[j,i]=which(hist(sto[i],br=thresh[i,1:top[i]],plot=F)$counts==1)
}
}
return(cbind(ages,it))
}

stats = function ()
{
age=sutures2[,3]
N=length(age)
sen=(apply(sutures2[,4:8],1,sum)-5)%/%2 + 1
stats = matrix(NA,nr=6,nc=6)

for(i in 1:6){
  sto = age[sen==i]
  stats[i,1] = length(sto)
  stats[i,2] = mean(sto)
  stats[i,3] = sd(sto)
  stats[i,4] = median(sto)
  stats[i,5:6] = range(sto)
}
stats[,-1]=format(stats[,-1],digits=3)
lo.hi = noquote(paste(stats[,5],rep('-',6),stats[,6],sep=' '))
Composite = c('5-6','7-8','9-10','11-12','13-14','15-16')
stats = data.frame(cbind(Composite,stats[,-(5:6)],lo.hi))
colnames(stats) = c('Composite','Count','Mean','s.d.','Median','Range')
print(stats,row.names=F)
}

cat('\nActual Data\n')
stats()

sutures2=cbind(rep(NA,N),rep(NA,N),sim.sutures())

cat('\n\nSimulated Data\n')
stats()
cat('\n')

}

Table.01()
Table 2

This Table gives results from  Lagrange multiplier tests for goodness of fit to log-normal transitions in the original 4-stage coding.  As the tests for the first four sutures are significant, they are also tested for a 3-stage "collapse" (described above).  The pseudo-R2 is Cragg and Uhler's (1970) measure for the 3-stage scoring for the first four sutures and the 4-stage scoring for the last suture.  The "R" version of Table 2 is as follows:


        Suture 4-stages 3-stages pseudo R2
  Mid-lambdoid  <0.0001   0.5618     0.290
        Lambda  <0.0001   0.2190     0.294
 Ant. sagittal  <0.0001   0.1088     0.218
        Bregma   0.0006   0.8213     0.327
 Sphenofrontal   0.4394    ----      0.346
LM.test = function (dat,collapse=1:10,log.it=T)
{
library(MASS)
age=dat[,1]
if(log.it==T) age=log(age)
stage=dat[,2]
#  Number of individuals
n=NROW(stage)
stage=collapse[stage]
#  Number of states
m=NROW(table(stage))
if(m==2) stage=stage-1
if(m>=3) {

results=polr(factor(stage)~age,method='probit')
beta=as.vector(results$coef)
alpha=as.vector(results$zeta)

E1sum=0
E2sum=0
BOX=matrix(rep(0,(m+2)*(m+2)),nc=(m+2))
for(i in 1:n){
 
  vec=rep(0,m+2)

  if(stage[i]==1){
     hj=alpha[1]-beta*age[i]
     norm=dnorm(hj)/pnorm(hj)
     E1=-(1-hj^2)*norm
     E2=-(3+hj^2)*hj*norm
     vec[1]=norm
     vec[m]=-vec[1]*age[i]
  }

  if(stage[i]==m){
    hjm1=alpha[m-1]-beta*age[i]
    norm=1/(1-pnorm(hjm1))
    E1=((1-hjm1^2)*dnorm(hjm1))*norm
    E2=((3+hjm1^2)*hjm1*dnorm(hjm1))*norm
    vec[m-1]=-dnorm(hjm1)*norm
    vec[m]=-vec[m-1]*age[i]
   
  }

  if(stage[i]>1 & stage[i]<m){
     hj=alpha[stage[i]]-beta*age[i]
     hjm1=alpha[stage[i]-1]-beta*age[i]
     norm=1/(pnorm(hj)-pnorm(hjm1))
     E1=((1-hjm1^2)*dnorm(hjm1)-(1-hj^2)*dnorm(hj))*norm
     E2=((3+hjm1^2)*hjm1*dnorm(hjm1)-(3+hj^2)*hj*dnorm(hj))*norm
     vec[stage[i]-1]=-dnorm(hjm1)*norm
     vec[stage[i]]=dnorm(hj)*norm
     vec[m]=(dnorm(hjm1)-dnorm(hj))*norm*age[i]
    }

    vec[m+1]=E1
    vec[m+2]=E2
    BOX=BOX+(vec%o%vec)

    E1sum=E1sum+E1
    E2sum=E2sum+E2
}
   J=diag(m+2)[,(m+1):(m+2)]
   E=c(E1sum,E2sum)
   x2=as.vector(t(E)%*%t(J)%*%solve(BOX)%*%J%*%E)
   return(pchisq(x2,2,lower.tail=F))
}

#  For binary

cons=as.vector(glm(stage~age,family=binomial(link="probit"))$coeff)
beta=cons[2]
cons=cons[1]
vec=c(0,0)
mat=matrix(rep(0,4*4),nc=4)

for(i in 1:n){
   h=cons+beta*age[i]
   sto=dnorm(h)
   wi=sto
   p=pnorm(h)
   wi=wi/(p*(1-p))
   if(is.finite(wi)){
   vec=vec+wi*(stage[i]-p)*h*c(h,3+h^2)
   Xi=c(age[i],h,h^2-1,h*(3+h^2))
   mat=mat+sto*wi*(Xi%o%Xi)
   }
}

J=diag(4)[,3:4]
x2= as.vector(vec%*%t(J)%*%solve(mat)%*%J%*%vec)
return(pchisq(x2,2,lower.tail=F))
}

pseudoR2 = function (iwhich=1)
{
library(MASS)
library(pscl)
dat=sutures2
age=log(dat[,3])

ncats = length(table(dat[,iwhich+3]))

if(ncats>=3){
stage=dat[,iwhich+3]
op1 <- polr(factor(stage) ~ age,Hess=TRUE,method="probit")}

else{
stage=dat[,iwhich+3]-1
op1 <- glm(factor(stage) ~ age, family=binomial(link = "probit"))}


pR2(op1)
}

Table.02=function ()
{
p4=0
for (i in 1:5) p4[i] = format(round(LM.test(sutures[c(3,i+3)]),4),nsmall=4,sci=F)
for(i in 1:5) if(as.numeric(p4[i])<0.0001) p4[i]='<0.0001'
p3=rep(0,5)
for (i in 1:4) p3[i] = format(round(LM.test(sutures2[c(3,i+3)]),4),nsmall=4,sci=F)
p3[5]=' ---- '
R2= 0
for(i in 1:5) R2[i] = format(round(pseudoR2(i)[6],3),nsmall=3)
suts = c('Mid-lambdoid','Lambda','Ant. sagittal','Bregma','Sphenofrontal')
suts = data.frame(cbind(suts,p4,p3,R2))
colnames(suts) = c('Suture','4-stages','3-stages','pseudo R2')
cat('\n\n')
print(suts,row.names=F)
cat('\n\n')
}


Table.02()
Figure 1
"Fig. 1 - Cumulative probit curves for closure of the mid-lambdoid suture.  The braces show the probabilities that a 50 year-old would be in each of the three states for this suture."

Fig.01 = function (sut=1,top=50)
{
library(VGAM)
library(pBrackets)
theta = coef(vglm(sutures2[,sut+3]~sutures2[,3],cumulative(link='probit',parallel=T,reverse = F)))
print(round(theta,3))
age = seq(-50,150,0.01)
plot(age,pnorm(theta[1]+age*theta[3],0,1),type='l',xlab='Age',ylab='Probability',ylim=c(0,1),lty=2,lwd=2)
abline(h=0)
abline(h=1)
lines(age,pnorm(theta[2]+age*theta[3],0,1),lwd=2)
thet=round(theta,3)
thet[3]=-thet[3]

legend(-60,.5,legend=substitute(paste(Phi(b2 - b3%*%'Age')),list(b2=thet[2],b3=thet[3])),
    bty='n',lwd=2)
legend(-60,.4,legend=substitute(paste(Phi(b1 - b3%*%'Age')),list(b1=thet[1],b3=thet[3])),
    bty='n',lwd=2,lty=2)
top1=pnorm(theta[1]+50*theta[3])
top2=pnorm(theta[2]+50*theta[3])
print(as.vector(c(theta,top1,top2)))
brackets(50,0,50,top1,h=5,lwd=2)
brackets(50,top1,50,top2,h=-5,lwd=2)
brackets(50,top2,50,1,h=-5,lwd=2)
lines(rep(50,2),c(0,1),lty=2)
text(57,top2+(1-top2)/2,'1.0 - 0.581 = 0.419',adj=0)
lines(c(55,80),rep(top1+(top2-top1)/2,2))
text(85,top1+(top2-top1)/2,'0.581 - 0.304 = 0.277',adj=0)
text(20,top1/2,'0.304',adj=0)
}

Fig.01()
 

Figure 2
"Fig. 2 – Re-drawing of Fig. 1 using transition distributions with a common standard deviation."

Fig.02 = function (sut=1)
{
library(VGAM)
theta = coef(vglm(sutures2[,sut+3]~sutures2[,3],cumulative(link='probit',parallel=T,reverse = F)))
sdev=1/(-theta[3])
mu1=theta[1]*sdev
age = seq(-50,150,0.01)
plot(age,dnorm(age,mu1,sdev),type='l',xlab='Age',ylab='Density',lty=2,
  main='Normal transition distributions',lwd=2)
mu2=theta[2]*sdev
lines(age,dnorm(age,mu2,sdev),lwd=2)
legend(105,.012,lty=c(2,1),lwd=c(2,2),legend=c(expression(1%->%2),expression(2%->%3)),bty='n')
}

Fig.02()
  
Figure 3
"Fig. 3 - A redrawing of the model depicted in Fig. 2 but using log-normal transitions rather than normal transitions.  The modal age for the first transition (19.2) is indicated as are the left and right boundaries for the 50% highest posterior density region."

Fig.03 = function (sut=1)
{
library(VGAM)
library(TeachingDemos)
theta = coef(vglm(sutures2[,sut+3]~log(sutures2[,3]),cumulative(link='probit',parallel=T,reverse = F)))
print(round(theta,3))
sdev=1/(-theta[3])
mu1=theta[1]*sdev
print(round(sdev,3))
print(round(mu1,3))

L.R=hpd(qlnorm,meanlog=mu1,sdlog=sdev,conf=.5)

mode = exp(mu1-sdev^2)
median = exp(mu1)
mean = exp(mu1 + 0.5*sdev^2)
v = exp(2*(mu1+sdev^2))-exp(2*mu1+sdev^2)

age = seq(0.01,100,.01)
plot(age,dlnorm(age,mu1,sdev),type='l',xlab='Age',ylab='Density',lty=2,
  main='Log normal transition distributions',lwd=2)

age.50 = seq(L.R[1],L.R[2],0.001)
polygon(c(L.R[1],age.50,L.R[2],L.R[1]),c(0,dlnorm(age.50,mu1,sdev),0,0),dens=5,lty=2)

lines(rep(mode,2),c(0.003,dlnorm(mode,mu1,sdev)),lty=2)
lines(rep(mode,2),c(0,.002),lty=2)
text(mode,0.0025,round(mode,1),cex=1.2)

text(L.R[1],0.001,round(L.R[1],1),cex=1.2)
text(L.R[2],0.001,"35.0",cex=1.2)

abline(h=0)
mu2=theta[2]*sdev
print(round(mu2,3))
lines(age,dlnorm(age,mu2,sdev),lwd=2)
legend(60,.02,lty=c(2,1),lwd=c(2,2),legend=c(expression(1%->%2),expression(2%->%3)),bty='n')
}

Fig.03()
 
Figure 4A
"Fig. 4A – Latent trait representation of the model from Fig. 3.  The curved line represents the regression of the mean latent trait value on the logarithm of age.  This line has an intercept of 0 on the log scale (one year on the raw scale)."


Fig.04A = function (sut=1)
{
library(VGAM)
theta = coef(vglm(sutures2[,sut+3]~log(sutures2[,3]),cumulative(link='probit',parallel=T,reverse = F)))
threshes = theta[1:2]
slope = -theta[3]
ages = c(20,40,60)
plot(ages,slope*log(ages),ylim=c(0,9),xlim=c(0,60),
   type='n',xlab='Age',ylab='Latent trait',
   main='Zero Intercept Parameterization')
for(i in 1:3){
   mu = log(ages[i])*slope
   a = seq(mu-3,mu+3,0.01)
   dens = dnorm(a,mu,1)
   lines(ages[i]-20*dens,a)
   lower.a=seq(mu-3,threshes[1],.01)
   dens = dnorm(lower.a,mu,1)
   polygon(c(ages[i]-20*dens,ages[i]),c(lower.a,threshes[1]),den=20)
   mid.a=seq(threshes[1],threshes[2],.01)
   dens=dnorm(mid.a,mu,1)
   polygon(c(ages[i],ages[i]-20*dens,ages[i]),c(threshes[1],mid.a,threshes[2]),col='grey')
   upper.a=seq(threshes[2],mu+3,.01)
   dens = dnorm(upper.a,mu,1)
   polygon(c(ages[i]-20*dens,ages[i]),c(upper.a,threshes[2]),den=20,angle=-45)

  
 }
for(i in 1:2) lines(c(-10,120),rep(threshes[i],2),lty=2,lwd=2)
for(i in 1:3){
 lines(rep(ages[i],2),c(0,100))
 points(ages[i],log(ages[i])*slope,pch=19)
}
ages=1:60
lines(ages,log(ages)*slope)
legend(-1,9,legend=c(expression(''>=' significant'),'minimal','open'),bty='n',
     dens=c(30,-10,30),fill=c('black','grey','black'),angle=c(-45,0,45),cex=1.25)
}

Fig.04A()
 
Figure 4B

"Fig. 4B – A redrawing of Fig. 4A that uses zero for the first threshold and consequently has a non-zero intercept for the regression line."



Fig.04B = function (sut=1)
{
library(VGAM)
theta = coef(vglm(sutures2[,sut+3]~log(sutures2[,3]),cumulative(link='probit',parallel=T,reverse = F)))
threshes = theta[1:2]
intercept = threshes[1]
threshes = threshes - intercept
slope = -theta[3]
ages = c(20,40,60)
plot(ages,slope*log(ages),ylim=c(-5.4,4),xlim=c(0,60),
   type='n',xlab='Age',ylab='Latent trait',
   main='Zero First Threshold Parameterization')

for(i in 1:3){
   mu = log(ages[i])*slope - intercept
   a = seq(mu-3,mu+3,0.01)
   dens = dnorm(a,mu,1)
   lines(ages[i]-20*dens,a)
   lower.a=seq(mu-3,threshes[1],.01)
   dens = dnorm(lower.a,mu,1)
   polygon(c(ages[i]-20*dens,ages[i]),c(lower.a,threshes[1]),den=20)
   mid.a=seq(threshes[1],threshes[2],.01)
   dens=dnorm(mid.a,mu,1)
   polygon(c(ages[i],ages[i]-20*dens,ages[i]),c(threshes[1],mid.a,threshes[2]),col='grey')
   upper.a=seq(threshes[2],mu+3,.01)
   dens = dnorm(upper.a,mu,1)
   polygon(c(ages[i]-20*dens,ages[i]),c(upper.a,threshes[2]),den=20,angle=-45)
 }
for(i in 1:2) lines(c(-10,120),rep(threshes[i],2),lty=2,lwd=2)
for(i in 1:3){
 lines(rep(ages[i],2),c(-100,100))
 points(ages[i],log(ages[i])*slope-intercept,pch=19)
}
ages=1:60
lines(ages,log(ages)*slope-intercept)
legend(-1,4,legend=c(expression(''>=' significant'),'minimal','open'),bty='n',
     dens=c(30,-10,30),fill=c('black','grey','black'),angle=c(-45,0,45),cex=1.25)
}

Fig.04B()
 


Bivariate Cumulative Probit

The code below is to get the bivariate cumulative probit parameters used to draw Figure 5.  It is a pretty slow process, so if you are in a hurry just take my word for it and use the Figure 5 code further below.

full.biv = function ()
{
    library(mvtnorm)
    library(MASS)
    library(alabama)
#   Age is in 3rd column, 1st suture of interest in 4th, 2nd of interest in 8th
    sto = sutures2[,c(3,4,8)]
    sto = sto[complete.cases(sto), ]
    N = NROW(sto)
    age <- log(sto[, 1])
    t1 = sto[, 2]
    probit1 = polr(factor(t1) ~ age, method = "probit")
    beta1 = probit1$coeff
    alpha1 = probit1$zeta
    NS1 = NROW(alpha1) + 1
    t2 = sto[, 3]
    probit2 = polr(factor(t2) ~ age, method = "probit")
    beta2 = probit2$coeff
    alpha2 = probit2$zeta
    NS2 = NROW(alpha2) + 1

    lnlk <- function(xpar) {
        alpha1 = xpar[1:2]
        beta1 = xpar[3]
        alpha2 = xpar[4:6]
        beta2 = xpar[7]
        L <- rep(0, 2)
        R <- rep(0, 2)
        corr <- matrix(c(1, xpar[8], xpar[8], 1), nc = 2)
        LLK <- 0
        for (i in 1:N) {
            j <- sto[i, 2]
            if (j == 1) {
                L[1] <- -Inf
                R[1] <- alpha1[1] - beta1 * age[i]
            }
            if (j > 1 & j < NS1) {
                L[1] <- alpha1[j - 1] - beta1 * age[i]
                R[1] <- alpha1[j] - beta1 * age[i]
            }
            if (j >= NS1) {
                L[1] <- alpha1[NS1 - 1] - beta1 * age[i]
                R[1] <- Inf
            }
            j <- sto[i, 3]
            if (j == 1) {
                L[2] <- -Inf
                R[2] <- alpha2[1] - beta2 * age[i]
            }
            if (j > 1 & j < NS2) {
                L[2] <- alpha2[j - 1] - beta2 * age[i]
                R[2] <- alpha2[j] - beta2 * age[i]
            }
            if (j >= NS2) {
                L[2] <- alpha2[NS2 - 1] - beta2 * age[i]
                R[2] <- Inf
            }
            if (L[1] > R[1]) {
                L[1] <- R[1]
            }
            if (L[2] > R[2]) {
                L[2] <- R[2]
            }
            LLK <- LLK + log(pmvnorm(lower = L, upper = R, corr = corr)[1])
        }
        cat(LLK,'\r')
        flush.console()
        return(-LLK)
    }

    hin=function(xx){
    h=0
#   thresholds are ordered
    h[1] = xx[2]- xx[1]
    h[3] = xx[5]-xx[4]
    h[4] = xx[6] - xx[5]   
#   slopes are positive
    h[2] = xx[3]
    h[5] = xx[7]
#   constrain correlation between -1 and 1
    h[6] = 1 - abs(xx[8])  
    }


     start = c(alpha1,beta1,alpha2,beta2,0)
     names(start) = c(names(start)[-8],'r')
     cat('\n\n')
     return(constrOptim.nl(par=start,fn=lnlk,hin=hin))
}

full.biv()
 
Figure 5
"Fig. 5 – Latent trait representation for the cumulative probit analysis of the mid-lambdoid and sphenofrontal sutures.  Isodensity ellipses are shown at ages 20 and 60 years for the 10%, 20%, 30%, 40%, and 50% levels."

Fig.05 = function ()
{
library(ellipse)

     
threshes1=c(4.9178464, 5.6551242)
slope1=1.4252749 
threshes2=c(6.2051313, 6.5505043, 7.0309677)
slope2=1.6935090 
r =  0.4575244
ages = c(20,60)
prec=solve(matrix(c(1,r,r,1),nc=2))
print(prec)

plot(ellipse(r,centre=c(slope1*log(ages[1]),slope2*log(ages[1])),level=.1),type='l',
   xlim=c(3,8.5),ylim=c(3,8.5),xlab='Mid-lambdoid suture',
   ylab='Sphenofrontal suture',
   main='Isodensity ellipses for 20 and 60 year-olds')
for(i in 2:5){
lines(ellipse(r,centre=c(slope1*log(ages[1]),slope2*log(ages[1])),level=i/10))
}

for(i in 1:5){
lines(ellipse(r,centre=c(slope1*log(ages[2]),slope2*log(ages[2])),level=i/10))
}


for (i in 1:2) lines(rep(threshes1[i],2),c(-100,100),lty=2,lwd=2)
for (i in 1:3) lines(c(-100,100),rep(threshes2[i],2),lty=2,lwd=2)

ages=c(1,300)
lines(slope1*log(ages),slope2*log(ages))
points(c(slope1*log(20),slope1*log(60)),c(slope2*log(20),slope2*log(60)),pch=19)
text((3+threshes1[1])/2,3,'open')
text(sum(threshes1)/2,3,'minimal')
text((threshes1[2]+8.5)/2,3,expression(''>=' significant'))

text(2.9,(3+threshes2[1])/2,'open',srt=90)
text(2.9,(threshes2[1]+threshes2[2])/2,'M',srt=90)
text(2.9,(threshes2[2]+threshes2[3])/2,'S',srt=90)
text(2.9,(threshes2[3]+8.5)/2,'closed',srt=90)
}

Fig.05()
 
"To circumvent a potentially long 'burn-in' phase..."

Getting starting values that have already converged (or nearly so) is pretty slow.  The code further below ("And here's the code for getting...") will get you there, but if you want to skip the drudgery then just copy the "pars" below which willl give you good starting values.  "Pars" came from running the code below.

pars = structure(list(beta = structure(c(-4.9082290591926, 1.42329593336931,
-4.86567094256147, 1.45764581143168, -3.79215127221542, 1.21235660483456,
-5.46499515953558, 1.56457356178808, -6.1702067007965, 1.68450785578428
), .Dim = c(10L, 1L)), alpha = structure(c(-Inf, -Inf, -Inf,
-Inf, -Inf, 0, 0, 0, 0, 0, 0.74031453673205, 0.610737780982399,
0.547909573779901, 0.675160023777213, 0.347932876845134, Inf,
Inf, Inf, Inf, 0.830077370660206, Inf, Inf, Inf, Inf, Inf), .Dim = c(5L,
5L)), r = c(0.832447060371051, 0.698357639338269, 0.655002097559437,
0.456740411134408, 0.776395490197772, 0.702750897790749, 0.421615491459039,
0.817337119625701, 0.49378646434538, 0.556925328951005)), .Names = c("beta",
"alpha", "r"))
And here's the code for getting "pars"

corr.only = function (beta1,alpha1,beta2,alpha2, it1=1,it2=2)
{
    library(mvtnorm)
    library(MASS)
    library(alabama)
    sto = sutures2[,c(3,it1+3,it2+3)]
    N = NROW(sto)
    age <- log(sto[, 1])
    NS1 = NROW(alpha1) + 1
    t2 = sto[, 3]
    NS2 = NROW(alpha2) + 1
    lnlk <- function(xpar) {
        L <- rep(0, 2)
        R <- rep(0, 2)
        corr <- matrix(c(1, xpar[1], xpar[1], 1), nc = 2)
        LLK <- 0
        for (i in 1:N) {
            j <- sto[i, 2]
            if (j == 1) {
                L[1] <- -Inf
                R[1] <- alpha1[1] - beta1 * age[i]
            }
            if (j > 1 & j < NS1) {
                L[1] <- alpha1[j - 1] - beta1 * age[i]
                R[1] <- alpha1[j] - beta1 * age[i]
            }
            if (j >= NS1) {
                L[1] <- alpha1[NS1 - 1] - beta1 * age[i]
                R[1] <- Inf
            }
            j <- sto[i, 3]
            if (j == 1) {
                L[2] <- -Inf
                R[2] <- alpha2[1] - beta2 * age[i]
            }
            if (j > 1 & j < NS2) {
                L[2] <- alpha2[j - 1] - beta2 * age[i]
                R[2] <- alpha2[j] - beta2 * age[i]
            }
            if (j >= NS2) {
                L[2] <- alpha2[NS2 - 1] - beta2 * age[i]
                R[2] <- Inf
            }
            if (L[1] > R[1]) {
                L[1] <- R[1]
            }
            if (L[2] > R[2]) {
                L[2] <- R[2]
            }
            LLK <- LLK + log(pmvnorm(lower = L, upper = R, corr = corr)[1])
        }
        cat(LLK,'\r')
        flush.console()
        return(-LLK)
    }

    hin=function(xx){
    h=0
#   constrain correlation between -1 and 1
    h[1] = 1 - abs(xx[1])  
    }

     return(constrOptim.nl(par=0,fn=lnlk,hin=hin))
}

pre.MCMC = function ()
{
    library(MASS)
    NS=c(rep(3,4),4)
    beta=matrix(NA,nr=2,nc=5)
    alpha=matrix(Inf,nr=5,nc=5)
    alpha[,1]=-Inf
    alpha[,2]=0
    sto = sutures2[,c(3:8)]
    N = NROW(sto)
    age <- log(sto[, 1])
    for(i in 1:5){
    t1 = sto[, i+1]
    probit = polr(factor(t1) ~ age, method = "probit")
    beta[1,i]=-probit$zeta[1]
    beta[2,i] = probit$coeff
    alpha[i,3:NS[i]] = probit$zeta[-1]-probit$zeta[1]
}
beta=matrix(beta,nc=1)
r.ip=0
r=0
for(it1 in 1:4){
for(it2 in (it1+1):5){
  print(c(it1,it2))
  flush.console()
  r.ip=r.ip+1
  r[r.ip] = corr.only(beta[it1],alpha[it1,1:(NS[it1]-1)],beta[it2],
            alpha[it2,1:(NS[it2]-1)],it1,it2)$par
}
}
return(list(beta=beta,alpha=alpha,r=r))
}


pars = pre.MCMC()
 
"Fifty thousand iterations of the MCMC were run retaining every 50th iteration for a total of 1,000 retained iterations."

The below is from Scott M. Lynch's (2007) code with some minor modifications.

Danger!

Warning: This is 50,000 iterations if you take the default.  So kinda' slow.  If you want to cut to the chase, then paste:

theta = c(5.04592569817729, 5.77671763748464, 1.45727631004962, 4.7223052525106,
5.32676372802259, 1.41473209930088, 3.55687789153358, 4.08785339807076,
1.14266620375592, 5.44444696530223, 6.11998718034894, 1.55698058091359,
6.3695231086795, 6.71415070932476, 7.19165042526973, 1.73415205669435,
0.837029246339445, 0.697905438608117, 0.655895562057807, 0.466170985415461,
0.784327419736573, 0.707167600569747, 0.440697227321095, 0.824497203634007,
0.50457617985205, 0.563144348998157)

And jump to the Table.03 code below.


multiv.MCMC = function (N.saved.iters=1000,jump.sd=.018,x=log(sutures2[,3]),
     z=sutures2[,4:8],thin=50)
{

#  Code from Scott M. Lynch (2007) Introduction to Applied Bayesian Statistics and Estimation for
#  Social Scientists.  Springer: New York, pp. 297-299.  With some slight modifications from
#  Konigsberg

set.seed(123456)
#R program for multivariate probit model

N = length(x)
x = cbind(rep(1,N),x)

# d is the dimensionality of the probit (from # of columns in "z"
# k is the number of explanatory variables (plus 1 for the intercept).  Ordinarily
# k=2, for age and intercept.
d=NCOL(z)
k=NCOL(x)
n.stages = apply(z,2,max)
max.stages=max(n.stages)

#create variable
zstar=matrix(0,nrow(z),d)

# Starting values for intercepts and slopes from univariate MLEs

b=pars$beta

# Starting values for thresholds from univariate MLEs

tz = pars$alpha


# Starting correlations from pairwise marginal likelihood (holding
# slopes, intercepts, and thresholds at univariate MLEs)

s=cs=to.symm(pars$r)


ctz=tz
acc1=acc2=acctot=0
n.corrs=d*(d-1)/2
nc.results=sum(n.stages)+n.corrs
temp = rep(NA,sum(n.stages))
results=matrix(0,nc=nc.results,nr=N.saved.iters)


#begin Gibbs sampling
ip=0
acc=0
iter=0
for(i in 1:(thin*N.saved.iters)){

#draw latent data: one-iteration gibbs sampler for truncated mvn simulation

bb=matrix(b,2,d)
m=x%*%bb
for(j in 1:d)
{
mm = m[,j] + (zstar[,-j]-m[,-j])%*%solve(s[-j,-j])%*%s[j,-j]
ss = s[j,j]-t(s[j,-j])%*%(solve(s[-j,-j]))%*%s[j,-j]
zstar[,j]=qnorm(runif(nrow(z),
min=pnorm(tz[j,z[,j]],mm,sqrt(ss)),
max=pnorm(tz[j,(z[,j]+1)],mm,sqrt(ss))),mean=mm,sd=sqrt(ss))
}

# draw new thresholds from uniforms defined by max(zstar) for stage below
# and min(zstar) for stage above threshold.  Note that Cowles MH sampler would be
# needed for large reference samples (typically not the case that they are)

for(ii in 1:d){
  if(n.stages[ii]>2){
  for(jj in 3:(n.stages[ii])){
    tz[ii,jj]= runif(1,max(zstar[z[,ii]==(jj-1),ii]),min(zstar[z[,ii]==jj,ii]))
  }
  }
}

#draw b from mvn
vb=solve(solve(s)%x%(t(x)%*%x))
mn=vb%*%(as.vector(t(x)%*%zstar%*%t(solve(s))))
b=mn+t(rnorm((d*k),0,1)%*%chol(vb))

#use metropolis-hastings sampling to draw sigma
e=matrix((as.vector(zstar)-(diag(d)%x%x%*%b)),nrow(z),d)
v=t(e)%*%e

 for(j in 1:(d-1)){
 for(jj in (j+1):d){
   cs=s
   iter=iter+1
   like=-.5*(nrow(z)+n.corrs+1)*log(det(s))-.5*sum(diag(v%*%solve(s)))
   cs[j,jj]=cs[jj,j]=cs[j,jj]+rnorm(1,mean=0,sd=jump.sd)
   if(abs(cs[j,jj])<1){
       det.cs = det(cs)
       if(det.cs<=0) clike=-Inf
   else {clike=-.5*(nrow(z)+n.corrs+1)*log(det.cs)-.5*sum(diag(v%*%solve(cs)))}
       if((clike-like)>log(runif(1,0,1))){
       acc=acc+1
       s=cs
   }
}
}
}

if(i%%thin==0){
cat(i/thin,' ',round(acc/iter*100,2),'\r')
flush.console()
ip=ip+1
b.vec=as.vector(b)
ip2 = 0
for(i.trait in 1:d){
   ip2=ip2+1
   temp[ip2] = -b.vec[2*i.trait-1]
   if(n.stages[i.trait]>2){
   for(j.stage in 3:(n.stages[i.trait])){
     ip2=ip2+1
     temp[ip2] = tz[i.trait,j.stage]-b[2*i.trait-1]
}
}
 ip2=ip2+1function(){

cbinom=function (x=90,p=.5)
{
N=500
sbin=0
sto=dbinom(0:N,N,p)
for(i in 0:N){
 if(dbinom(i,N,p)<=dbinom(x,N,p)) sbin=sbin+dbinom(i,N,p)
}
return(sbin)
}

get.Table = function (left.most=F,i.which=1)
{
 dat = paste('hpds',i.which,'.txt',sep='')
 dat = dget(dat)
 age=simulated.sutures[,1]
 correct=0
 p.value=0
 for(i in 1:9){
   L.bound = (i-1)*2 + 1
   U.bound = L.bound + 1
   correct[i] = sum(dat[,L.bound]<=age & age <=dat[,U.bound])
   p.value[i] = round(cbinom(correct[i],i/10),3)
   p.value[i] = format(p.value[i],nsmall=3,sci=F)
   if(as.numeric(p.value[i])<0.001) p.value[i]='<0.001'

}
if(left.most==T)
 {Table = cbind(noquote(paste(seq(10,90,10),'%',sep='')),500*seq(.1,.9,.1),correct,p.value)}
else
 {Table = cbind(correct,p.value)
 }
 Table = data.frame(Table)
if(left.most==T) colnames(Table) = c('Coverage','Expect','Obs','p-value')
else{ colnames(Table) = c('Obs','p-value')}
 return(Table)
}
cat('\n\n                 Informative Uniform     Uniform\n')
cat('                 Prior       (15-120)    (15-85)\n')
print(cbind(get.Table(T,1),get.Table(F,2),get.Table(F,3)),row.names=F)
cat('\n')
}
 temp[ip2] = b.vec[i.trait*2]
}

results[ip,]=c(temp,s[lower.tri(s)])

}
}
cat('\n\n')

find.it = matrix(NA,nr=d,nc=3)

ip2 = 0
for(i.trait in 1:d){
   find.it[i.trait,1]=i.trait
   ip2=ip2+1
   find.it[i.trait,2]=ip2
   if(n.stages[i.trait]>2){
   for(j.stage in 3:(n.stages[i.trait])){
     ip2=ip2+1
   }
   }
 ip2=ip2+1
 find.it[i.trait,3]=ip2
}
ip2=ip2+1
colnames(find.it)=c('Trait #','start','stop')


print(find.it)
cat(paste('Corr matrix starts at:',ip2,'and ends at:',ip2+n.corrs-1,'\n\n'))

return(apply(results,2,mean))
}

theta = multiv.MCMC()


 Which will print the following to the screen.  The "1000" is showing the 1000th retained iteration, while the 48.33 is the acceptance percentage for the Metropolis-Hastings Sampler at that point.  The iterations will "over print" so that you can watch progress.  The "Trait # start stop" and "Corr matrix starts at:" is your secret decoder ring for columns to parameters.

1000   48.33


     Trait # start stop
[1,]       1     1    3
[2,]       2     4    6
[3,]       3     7    9
[4,]       4    10   12
[5,]       5    13   16
Corr matrix starts at: 17 and ends at: 26

But if you don't want to bother with the secret decoder ring:

Table.03_4 = function ()
{
  library(MASS)
  ests = round(theta,3)
  cat('\nSlope   T1    T2    T3\n')
  mat.ests = matrix(NA,nc=4,nr=5)
  for(i in 1:4) mat.ests[i,] = ests[c(i*3,i*3-c(2:1,1))]
  mat.ests[5,] = ests[c(16,13:15)]
  write.matrix(mat.ests)
  cat('\nResidual correlation matrix\n')
  write.matrix(to.symm(ests[17:26]))
}

Table.03_4()

Which gets you:

Slope   T1    T2    T3
1.457 5.046 5.777 5.777
1.415 4.722 5.327 5.327
1.143 3.557 4.088 4.088
1.557 5.444 6.120 6.120
1.734 6.370 6.714 7.192

Residual correlation matrix
1.000 0.837 0.698 0.656 0.466
0.837 1.000 0.784 0.707 0.441
0.698 0.784 1.000 0.824 0.505
0.656 0.707 0.824 1.000 0.563
0.466 0.441 0.505 0.563 1.000


 
Figure 6
"Fig. 6 – Posterior densities of age under a uniform prior from 15 to 120 years.  The first four sutures are closed while the last (sphenofrontal) is open.  Labelled dashed lines give the individual sutures, the heavy dashed line gives the posterior under the assumption of conditional independence, and the heavy solid line gives the posterior accounting for residual correlations between the sutures. "

Gompertz = function (t,x=c(0.025,0.018))
{
if(sum(x)==0) return(1)
         a3 = x[1]
         b3 = x[2]
         shift = 15
         h.t = a3*exp(b3*(t-shift))
         S.t = exp(a3/b3*(1-exp(b3*(t-shift))))
         return(S.t*h.t)
}

Fig.06 = function (stages=c(3,3,3,3,1),right=120,individ.traits=T)
{

library(mvtnorm)
labs = c('S1','S2','S3','S4','S5')
labs.hor = c(50,115,25,100,20)

threshes1=c(-Inf,theta[1:2],Inf)
threshes2=c(-Inf,theta[4:5],Inf)
threshes3=c(-Inf,theta[7:8],Inf)
threshes4=c(-Inf,theta[10:11],Inf)
threshes5=c(-Inf,theta[13:15],Inf)
slopes=theta[c(3,6,9,12,16)]
r=to.symm(theta[17:26])


bot = c(threshes1[stages[1]],threshes2[stages[2]],threshes3[stages[3]],threshes4[stages[4]],threshes5[stages[5]])
top = c(threshes1[stages[1]+1],threshes2[stages[2]+1],threshes3[stages[3]+1],threshes4[stages[4]+1],threshes5[stages[5]+1])

################### CORRELATION #####################################

get.densr = function(x){
  return(pmvnorm(mean=slopes*log(x),lower=bot,upper=top,corr=r)[1]*Gompertz(x,c(0,0)))
}
denomr = integrate(Vectorize(get.densr),15,right)$val
get.postr = function(x){
  return(pmvnorm(mean=slopes*log(x),lower=bot,upper=top,corr=r)[1]*Gompertz(x,c(0,0))/denomr)
}
age0=optimize(get.postr,int=c(15,right),maximum=T)$max
cat('corr ',round(age0,1),'\n')


################## NO CORRELATION ##################################

get.dens0 = function(x){
  return(pmvnorm(mean=slopes*log(x),lower=bot,upper=top,corr=diag(5))[1]*Gompertz(x,c(0,0)))
}
denom0 = integrate(Vectorize(get.dens0),15,right)$val
get.post0 = function(x){
  return(pmvnorm(mean=slopes*log(x),lower=bot,upper=top,corr=diag(5))[1]*Gompertz(x,c(0,0))/denom0)
}
age0=optimize(get.post0,int=c(15,right),maximum=T)$max
cat('No corr ',round(age0,1),'\n')
max.pdf = get.post0(age0)

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

p=0
for(i in 15:right){
p[i-14] = get.postr(i)
}
title=paste(labs[1],' = ',stages[1],', ',labs[2],' = ',stages[2],', ',
            labs[3],' = ',stages[3],', ',labs[4],' = ',stages[4],', ',
            labs[5],' = ',stages[5], sep='')
plot(15:right,p,type='l',xlab='Age',ylab='Density',lwd=2,ylim=c(0,.024),
   main=title)

for(i in 15:right){
p[i-14] = get.post0(i)
}

lines(15:right,p,lwd=2,lty=2)
abline(h=0)
legend(2/3*right,8/9*.024,
legend=c('independent','dependent'),
   bty='n',lwd=2,lty=c(2,1),
   cex=1.25,y.intersp=1.5)

if(individ.traits==F) return(cat('All done\n\n'))

get.dens1 = function(x,iwhich){
  return(pmvnorm(mean=slopes[iwhich]*log(x),lower=bot[iwhich],upper=top[iwhich],sig=1)[1])
}

get.post1 = function(x,iwhich){
  return(pmvnorm(mean=slopes[iwhich]*log(x),lower=bot[iwhich],upper=top[iwhich],sig=1)[1]/denom1)
}

for(iwhich in 1:5){
denom1 = integrate(Vectorize(get.dens1,'x'),0.01,150,iwhich=iwhich)$val
for(i in 1:right){
p[i] = get.post1(i,iwhich)
}
lines(1:right,p,lty=2)
left = which.max(p)
if(left<16) left = 16
#height = p[left]
height = p[labs.hor[iwhich]]
text(labs.hor[iwhich],height,labs[iwhich],adj=c(0.5,0.5))
}
}


Fig.06()
 

Simulating a new data set with 500 individuals and a3 = 0.001 and b3 = 0.1
This is needed below to look at properties of age estimation.

sim.sutures2 = function (N=500)
{
library(mvtnorm)
set.seed(123456)
thresh=matrix(NA,nr=5,nc=5)
thresh[1,]=c(-100,theta[1:2],100,NA)
thresh[2,]=c(-100,theta[4:5],100,NA)
thresh[3,]=c(-100,theta[7:8],100,NA)
thresh[4,]=c(-100,theta[10:11],100,NA)
thresh[5,]=c(-100,theta[13:15],100)
top=c(rep(4,4),5)
slope=theta[c(3,6,9,12,16)]
r=to.symm(theta[17:26])
ages=sim.Gomp(N=N,a3=0.001,b3=0.1)

it=matrix(NA,nr=N,nc=5)
for(j in 1:N){
mu=log(ages[j])*slope
sto=rmvnorm(1,mu,r)
for(i in 1:5){
 it[j,i]=which(hist(sto[i],br=thresh[i,1:top[i]],plot=F)$counts==1)
}
}
return(cbind(ages,it))
}


simulated.sutures = sim.sutures2()

Metropolis-Hastings Sampler for Age Estimation


Metrop = function (x=c(0.001,0.1),start.age=50,stage=c(1,1,1,1,1), right=120, n.iters = 1200, jump.sd= 30,thin=5,C.I=F)
{
library(mvtnorm)
library(tmvtnorm)
library(truncnorm)

thresh=matrix(NA,nr=5,nc=5)
thresh[1,]=c(-Inf,theta[1:2],Inf,NA)
thresh[2,]=c(-Inf,theta[4:5],Inf,NA)
thresh[3,]=c(-Inf,theta[7:8],Inf,NA)
thresh[4,]=c(-Inf,theta[10:11],Inf,NA)
thresh[5,]=c(-Inf,theta[13:15],Inf)
top=c(rep(4,4),5)
beta=theta[c(3,6,9,12,16)]
sigma=to.symm(theta[17:26])
if(C.I==T) sigma=to.symm(rep(0,10))


lo=-Inf
hi=Inf

for(i in 1:5){
  lo[i]=thresh[i,stage[i]]
  hi[i]=thresh[i,stage[i]+1]
}
age=start.age
sto=age
i.sto=0
acc=0

for(i in 1:(n.iters*thin)){
mu = beta*log(age)
z.star = rtmvnorm(1,mean=mu,sigma,lower=lo,
  upper = hi,algorithm='gibbs')
f.prev = dmvnorm(z.star,mu,sigma)*Gompertz(age,x)

age.candidate = rtruncnorm(1, a=15, b=right, mean = age, sd = jump.sd)
mu = beta*log(age.candidate)
f.cand = dmvnorm(z.star,mu,sigma)*Gompertz(age.candidate,x)
Ratio=(f.cand/f.prev)
Importance = dtruncnorm(age, a=15, b=right, mean = age.candidate, sd = jump.sd)/
             dtruncnorm(age.candidate, a=15, b=right, mean = age, sd = jump.sd)


Ratio = Ratio*Importance

u=runif(1,0,1)
if(Ratio>u) {
age=age.candidate
acc=acc+1
}


if(i%%thin==0){
  i.sto=i.sto+1
  sto[i.sto]=age
  cat(i.sto,round(acc/i*100,2),'\r')
  flush.console()
}
}
cat('\n\n')
return(sto)
}

Apply the Metropolis-Hastings Sampler to the simulated data under various assumptions

check.coverage = function (x=c(0.001,0.1),right = 120,C.I=F)
{
# x=c(0,0) will use a uniform age distribution from age 15 up until age specified by "right"
# C.I=F uses the estimated residual correlation matrix among sutures, while
# C.I=T (for "conditional independence = TRUE") assumes all residual correlations are zero.
set.seed(12345)
my.emp.hpd = function(x, conf=0.95){
        conf <- 1-conf
        n <- length(x)
        nn <- round( n*conf )
        x <- sort(x)
        xx <- x[ (n-nn+1):n ] - x[1:nn]
        m <- min(xx)
        nnn <- which(xx==m)[1]
        return( c( x[ nnn ], x[ n-nn+nnn ] ) )
}

hpds = matrix(NA,nr=500,nc=18)
for(i in 1:500){
print(i)
stages=simulated.sutures[i,-1]
junk = Metrop(stage=stages,x=x,right=right,C.I=C.I)[-(1:200)]
for(j in 1:9){
hpds[i,(2*j-1):(2*j)] = my.emp.hpd(junk,conf=j/10)
}
}
return(hpds)
}

WARNING - DANGER WILL ROBINSON!Warning: This i
Danger Will Robinson!

The four lines below will take a long, long time to run.  The first estimates ages under an informative prior, the second and third under uniform priors (15-120 and 15-85), and the last under an informative prior assuming conditional independence.  If you are impatient (I know I am), then download the four text files below and save them into the folder that contains your "R" workspace.

dput(check.coverage(),'hpds1.txt')
dput(check.coverage(c(0,0),'hpds2.txt')
dput(check.coverage(c(0,0),right=85),'hpds3.txt')
dput(check.coverage(C.I=T),'hpds4.txt)


hpds1.txt
hpds2.txt
hpds3.txt
hpds4.txt
Table 4
This Table gives "coverages" as follows:

                 Informative Uniform     Uniform
                 Prior       (15-120)    (15-85)
 Coverage Expect Obs p-value Obs p-value Obs p-value
      10%     50  49   0.941  36   0.037  37   0.052
      20%    100  87   0.162  75   0.004  72   0.001
      30%    150 144   0.591 131   0.064 112  <0.001
      40%    200 207   0.523 172   0.011 171   0.008
      50%    250 243   0.561 209  <0.001 217   0.004
      60%    300 299   0.927 258  <0.001 283   0.121
      70%    350 343   0.495 310  <0.001 349   0.922
      80%    400 383   0.065 385   0.094 405   0.615
      90%    450 434   0.021 464   0.037 463   0.052

Table.04 = function(){

cbinom=function (x=90,p=.5)
{
N=500
sbin=0
sto=dbinom(0:N,N,p)
for(i in 0:N){
 if(dbinom(i,N,p)<=dbinom(x,N,p)) sbin=sbin+dbinom(i,N,p)
}
return(sbin)
}

get.Table = function (left.most=F,i.which=1)
{
 dat = paste('hpds',i.which,'.txt',sep='')
 dat = dget(dat)
 age=simulated.sutures[,1]
 correct=0
 p.value=0
 for(i in 1:9){
   L.bound = (i-1)*2 + 1
   U.bound = L.bound + 1
   correct[i] = sum(dat[,L.bound]<=age & age <=dat[,U.bound])
   p.value[i] = round(cbinom(correct[i],i/10),3)
   p.value[i] = format(p.value[i],nsmall=3,sci=F)
   if(as.numeric(p.value[i])<0.001) p.value[i]='<0.001'

}
if(left.most==T)
 {Table = cbind(noquote(paste(seq(10,90,10),'%',sep='')),500*seq(.1,.9,.1),correct,p.value)}
else
 {Table = cbind(correct,p.value)
 }
 Table = data.frame(Table)
if(left.most==T) colnames(Table) = c('Coverage','Expect','Obs','p-value')
else{ colnames(Table) = c('Obs','p-value')}
 return(Table)
}
cat('\n\n                 Informative Uniform     Uniform\n')
cat('                 Prior       (15-120)    (15-85)\n')
print(cbind(get.Table(T,1),get.Table(F,2),get.Table(F,3)),row.names=F)
cat('\n')
}

Figures 7 - 9
These Figures are in the style of Figs. 2-6 from Milner and Boldsen (2012).

Fig.07to09 = function (which.Fig = 7,which.int=5)
{
# which.Fig specifies the Figure, which can be 7, 8, or 9
# which.int specifies the coverage, so = 5 (default) is 50%.
# = 1 is 10%, = 9 in 90%, etc.
labs=c(10,20,30,40,50,60,70,80,90)
ints=(1:9)/10
if(which.Fig==7){
 ages=cbind(simulated.sutures[,1],dget('hpds1.txt'))
 type='Gompertz Prior'
}
if(which.Fig==8){
 ages=cbind(simulated.sutures[,1],dget('hpds3.txt'))
 type='Uniform Prior'
}
if(which.Fig==9){
 ages=cbind(simulated.sutures[,1],dget('hpds4.txt'))
 type='Gompertz Prior'
}

ages=ages[sort(ages[,1],index.ret=T)$ix,]
bot = 2*which.int
top = bot + 1
covered = sum(ages[,1]>=ages[,bot] & ages[,1]<=ages[,top])
real = round(covered/5,1)
plot(ages[,1],type='l',ylab='Age',lwd=2,main=paste(type,' (',labs[which.int],'% HPD, ',real,
  '% realized)',sep=''),ylim=c(15,85))
for( i in seq(2,500,2)){
 if(ages[i,1]>=ages[i,bot] & ages[i,1]<=ages[i,top])
    {lines(rep(i,2),c(ages[i,bot],ages[i,top]))}
 else(lines(rep(i,2),c(ages[i,bot],ages[i,top]),lty=3))
}

too.low = sum(ages[,1]<ages[,bot])
too.hi = sum(ages[,1]>ages[,top])
cat(covered,' of 500 at nominal coverage of p = ',ints[which.int],'\n\n')
cat('age < HPD ',too.low,' age> HPD', too.hi,'\n')
if(which.Fig==9) text(0,80,'Conditional Independence Assumed',adj=0,cex=1.25)
}

Fig.07to09()
 

References