Getting Some Data

It is debatable what is the one most terrifying aspect when first using "R," but importing data has to be pretty high on the list.  The good news is that "R" reads *.csv (comma-delimited files) quite easily.  So click on the following oh intrepid traveler:

refbones.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):

refbones=read.csv('refbones.csv')

Now check to see that this really exists in your workspace (where anything in green is what I actually typed) :

> refbones=read.csv('refbones.csv')
> ls()
[1] "refbones"
> refbones[1:10,]
   Sample ID Age Sex Weight  Hum  Rad  Uln  Fem  Tib  Fib
1     F&K  1  12   M     12  8.5  5.5  5.5  7.0  5.1  5.0
2     F&K  2  12   F     40  9.0  8.0  9.0 10.0  7.0  7.0
3     F&K  3  14   M     60 11.6  9.8 10.9 11.9  9.5  9.3
4     F&K  4  14   F    100 11.7  9.1 10.2 11.5  9.1  8.5
5     F&K  5  14   F     60 14.0 11.5 12.5 13.8 12.0 12.0
6     F&K  6  16   M    110 19.5 16.0 17.5 18.0 16.0 14.5
7     F&K  7  16   M    100 18.0 16.5 18.0 19.0 16.5 15.5
8     F&K  8  16   M    100 20.0 19.0 20.7 22.5 18.9 17.8
9     F&K  9  16   M    180 20.0 20.5 22.0 24.0 21.0 20.5
10    F&K 10  16   M    190 20.3 18.2 20.7 22.5 19.5 19.2

OK, all's good.  Let's try to get a handle on what we have here:

> table(refbones[,1])

  F&K M&D39
  138   150

So this is data on the 138 cases from Fazekas and Kósa (1978) and 150 girls from Maresh and Deming (1939).  OK, except that I lied about that being 150 girls.  This is semi-longitudinal data, so it is actually:

> length(unique(refbones$ID[-(1:138)]))
[1] 41

Forty-one girls.   OK, that looked kinda' nasty, so here's the explanation on what just happened.  refbones$ID is going to just pull out the "ID" variable.  The first 138 cases are from Fazekas and Kosa, so the -(1:138) drops those cases, "unique" is, well, unique, and length is length.

Now let's get another dataset to play with:

uk.csv

> UK.bones=read.csv('uk.csv')
> ls()
[1] "refbones" "UK.bones"
> UK.bones[1:10,]
   ID Site  Hum  Rad  Uln  Fem  Tib      fem
1   1   RF   NA   NA 48.6 61.2   NA 61.20000
2   2   RF   NA 46.6   NA 63.2 54.6 63.20000
3   3   RF 59.9 49.0   NA 65.9 58.7 65.90000
4   4   RF 58.2   NA   NA   NA 60.0 67.55579
5   5   RF 60.0   NA   NA 66.8   NA 66.80000
6   6   RF 58.5 47.8 53.9 68.4 58.9 68.40000
7   7   RF 60.9   NA   NA 68.7 60.6 68.70000
8   8   RF 62.1 49.6 56.1 68.9 59.9 68.90000
9   9   RF   NA   NA   NA 69.9   NA 69.90000
10 10   RF   NA   NA 55.5 70.1 60.2 70.10000

These data are from Lewis and Gowland (2007), and give long bone lengths on:

> NROW(UK.bones)
[1] 118

118 individuals.  The "NA" values are "not available."  To make life simpler for at least a bit, we will just work with femur length.  Because this is missing for many of the individuals, specifically for:

> sum(is.na(UK.bones$Fem))
[1] 36

36 of them, there's a final created column called "fem" (case matters with "R", so "Fem" and "fem" are different variables) which is the regression estimate of femur length on all available bones (or the actual value if available).  The single imputation for missing femoral lengths is just that.  That single imputation was done using the R package norm.

Now, let the wild rumpus begin!

Then paste the following into your Rgui:

FemRef.mfp = function ()
{
# uncomment line below if you actually want to run the mfp
# library(mfp)
set.seed(1234)
RF = data.frame(cbind(refbones[,3]*7,refbones[,9]))
colnames(RF)=c('Age','fem')

# Comment out the line below and uncomment the line below that to run
# the mfp
  sto = lm(fem~I(sqrt(Age)),data=RF)
# sto=mfp(formula = fem ~ fp(Age,df=2), data = RF)

betas=as.vector(sto$coeff)

mu = function(t) {
  age=t
  return(betas[1]+betas[2]*sqrt(age))
}

absResid=abs(RF$fem-mu(RF$Age))
# uncomment line below, then comment out the line below that if you want to use mfp
# betas2=as.vector(mfp(formula = absResid~fp(RF$Age,df=1))$coef)
betas2 = lm(absResid ~ RF$Age)$coeff
SD = function(t) {(betas2[1]+betas2[2]*t)*sqrt(pi/2)}

age=seq(10*7,70*7,.1)
est.mu = mu(age)
est.sd = SD(age)

plot(age,est.mu+1.96*est.sd,type='l',lty=2,xlab='Age (LMP days)',
    ylab='Femur Length (mm.)',ylim=c(0,120))
text(75,109,substitute(paste(mu[y],' = ',b1,' + ',b2,sqrt(x)),
   list(b1=round(betas[1],4),b2=signif(betas[2],4))),adj=0,cex=1.25)
text(75,99,substitute(paste(sigma[y],' = ',b1,' + ',b2,x),
   list(b1=signif(betas2[1]*sqrt(pi/2),5),b2=signif(betas2[2]*sqrt(pi/2),4))),adj=0,cex=1.25)
lines(age,est.mu)
lines(age,est.mu-1.96*est.sd,lty=2)

dat = cbind(jitter(RF$Age,factor=1),RF$fem)
points(dat,cex=.5)

est.mu = mu(RF$Age)
est.sd = SD(RF$Age)

test.it=cbind(est.mu-1.96*est.sd,RF[,2],est.mu+1.96*est.sd)
perc=round(sum(test.it[,2]>=test.it[,1]&test.it[,2]<=test.it[,3])/NROW(RF)*100,2)
cat(noquote(paste(sep='','est.mu +/- 1.96*est.sd femur length conditioned on age contains ',perc,'% of data','\n')))
}


And now the wild rumpus really begins:

> ls()
[1] "FemRef.mfp" "refbones"   "UK.bones" 
> FemRef.mfp()
Loading required package: survival
Loading required package: splines
est.mu +/- 1.96*est.sd femur length conditioned on age contains 95.83% of data

and you'll get the following:

FemReg
The regression technique used here was a  "fractional polynomial" (with cherries on top?) and regression of absolute residuals.  See Altman (1993).  Now make a function to "invert" this, so that given a femur length we can read the graph "backwards" to get an age in days LMP:

est.age=function (fem=60) return(((fem+73.8884)/8.829)^2)

And now we can find the age in days LMP for someone with a femur length of 60 mm:

> est.age()
[1] 229.9656

OK, a point estimate is not really going to cut it.  So let's think about this a bit more.   Our lovely graph above shows us how to predict the mean and standard deviation of femoral length given known age, so the distribution of femur length given known age.  But we want the reverse, or the distribution of age given a known femur length.  This screams Bayes' Theorem.  So let's start down that path.

First we'll want a function for the likelihood:

Likelihood = function (x=280,y=80)
{
 mu.y = -73.8884 + 8.829*sqrt(x)
 sd.y = 1.6216 + 0.009508*x
 return(dnorm(y,mu.y,sd.y))
}

And now we can try this in a very rough example:

plot(100:300,Likelihood(x=100:300,y=60),type='l')
lines(rep(est.age(60),2),c(0,2))

But that's not "normalized" (in other words, it doesn't integrate to 1.0).  So we need to integrate across age and then divide.  The following code does this, and a lot more.

Posterior = function (y = 60)
{

# Integrate across likelihood
denom = integrate(Likelihood,lower=0,upper=600,y=y)$val

# Posterior density is likelihood at age x divided by integral across x
post = function(x){
   return(Likelihood(x,y=y)/denom)
}

# Get age estimate from function and posterior density at that age
age = est.age(y)
max.pdf = post(age)

# Following three functions are for finding left and right tails of posterior
bound = function(x,xx) return(post(x)-xx)
lo = function(xx) uniroot(bound,lower=0,upper=age,tol=.Machine$double.eps,xx=xx)$root
hi = function(xx) uniroot(bound,lower=age,upper=600,tol=.Machine$double.eps,xx=xx)$root

# For finding specific value of posterior density that contains "want.hdr" HDR
how.hi = function(xx, want.hdr){
bot = lo(xx); top = hi(xx)
 return(integrate(post,bot,top,abs.tol=.Machine$double.eps)$val - want.hdr)
}

hdr.height = uniroot(how.hi,lower=max.pdf/10e4,upper=max.pdf,want.hdr=0.999,tol=.Machine$double.eps)$root
bot = lo(hdr.height); top = hi(hdr.height)
hdr = c(bot,top)

plot(0:600,post(0:600),type='l',xlab='Age (LMP Days)',ylab='Posterior Density',xlim=hdr)
lines(c(0,600),rep(0,2))

# 95% interval using QUANTILES
left = function(x) integrate(post,lower=0,upper=x,abs.tol=.Machine$double.eps)$val-0.025
med = function(x) integrate(post,lower=0,upper=x,abs.tol=.Machine$double.eps)$val-0.5
right = function(x) integrate(post,lower=0,upper=x,abs.tol=.Machine$double.eps)$val-0.975
bot.pi=uniroot(left,lower=0,upper=age,tol=.Machine$double.eps)$root
top.pi=uniroot(right,lower=age,upper=600,tol=.Machine$double.eps)$root
med.est = uniroot(med,lower=bot.pi,upper=top.pi,tol=.Machine$double.eps)$root
PI = c(bot.pi,top.pi)
for(i in 1:2) lines(rep(PI[i],2),c(0,post(PI[i])))
horizon = max(post(bot.pi),post(top.pi))
lines(c(0,1000),rep(horizon,2),lty=2,lwd=2)
polygon(c(PI[1],seq(PI[1],PI[2],.5),rep(PI[2],2)),
      c(0,post(c(seq(PI[1],PI[2],.5),PI[2])),0),density=10)

cat('\n\n')
cat(noquote(paste(sep=' ','2.5%, 50%, and 97.% at',round(PI[1],1),round(med.est,1),round(PI[2],1),'\n\n')))
cat('This is the quantile-based probability interval\n\n')
junk=readline('Hit anything except your neighbor: ')

plot(0:600,post(0:600),type='l',xlab='Age (LMP Days)',ylab='Posterior Density',xlim=hdr)
lines(c(0,600),rep(0,2))

# Previous was 99.9% HPD for plotting.  Now find 95% HIGHEST POSTERIOR DENSITY
hdr.height = uniroot(how.hi,lower=max.pdf/10e4,upper=max.pdf,want.hdr=0.95,tol=.Machine$double.eps)$root
bot = lo(hdr.height); top = hi(hdr.height)
hdr = c(bot,top)

lines(c(0,600),rep(hdr.height,2),lty=2,lwd=2)
lines(rep(bot,2),c(0,hdr.height))
lines(rep(top,2),c(0,hdr.height))

polygon(c(hdr[1],seq(hdr[1],hdr[2],.5),rep(hdr[2],2)),
      c(0,post(c(seq(hdr[1],hdr[2],.5),hdr[2])),0),density=10)

l.area = round(integrate(post,lower=0,upper=bot,abs.tol=.Machine$double.eps)$val*100,4)
r.area = round(integrate(post,lower=top,upper=600,abs.tol=.Machine$double.eps)$val*100,4)
cat('\n')
cat(noquote(paste(sep=' ','Left, Mode, and Right at',round(hdr[1],1),round(age,1),round(hdr[2],1),'\n\n')))
cat(noquote(paste(sep='','Tail to left of 95% HDR contains ',l.area,'% of posterior density','\n')))
cat(noquote(paste(sep='','Tail to right of 95% HDR contains ',r.area,'% of posterior density','\n')))
}



Maybe you get no joy from calculus, so let's move over to OpenBUGS instead:

1. Start OpenBUGS and you'll get the following:

Opening menu

2.  Then open a new file:

OpenBUGS2

3.  Copy and paste the boxed material below:

model
{
   y <- 60
   y.mu <- -73.8884 + 8.829*sqrt(x)
# Note: In Bayesian lingo you use "precision" rather
# than variance (or standard deviation). It is
# just the inverse of variance.
   y.prec <- 1/pow(1.6216 + 0.009508*x,2)
   y ~ dnorm(y.mu, y.prec)
   x ~ dunif(0,600)
}

OpenBUGS3

4.  Go to "model" and then "specification"

OpenBUGS4



5.  And do "check model" after you have highlighted the block:

OpenBUGS5

6.  At the very bottom of the screen in the bottom left corner you should get some good news.

OpenBUGS6

7.  Now do "compile" and then "gen inits" from the menu (hey, no data to load) and you are good to go.

8.  Go to "model" and then "update" and do the default 1000 updates.

9.  Now go to "inference" and then "samples."  In "node" type "x" (without quotes) and hit set.  Then do, say, 5,000 updates.  Then click on "stats," "density" and "auto cor" to get the goods:

OpenBUGS7

Now it would be fun to see what method was used for the updating:

OpenBUGS8

So OpenBUGS was using the slice sampler (Neal 2003).   Slice sampling can be done in "R" using "stepout.slice.sample" from the SamplerCompare package.  And we can go quickly over how slice  sampling works in the following graphic:

slice1slice2

But it won't be very instructive for us to use slice sampling in R, so let's instead use Metropolis sampling within "R."  

Metrop = function (burn.in=1000, jump.sd=25)
{
# Comment out the line below if you want a random seed
set.seed(12345)
y = 60
x.curr=100

mu.y=function (x) -73.8884+8.829*sqrt(x)
sd.y=function (x) 1.6216+0.009508*x

# Called "One.Step" because this is one step in the Metropolis Sampler
One.Step = function(x.curr){
curr.dens = dnorm(y,mu.y(x.curr),sd.y(x.curr))
x.trial = x.curr + rnorm(1,0,jump.sd)
# Can't take square root of a negative, so return current value
if(x.trial<=0) return(x.curr)
trial.dens = dnorm(y,mu.y(x.trial),sd.y(x.trial))
R = trial.dens/curr.dens
if(R>=1) return(x.trial)
else{
  U = runif(1)
  if(U<=R) return(x.trial)
}
return(x.curr)
}

# 1000 iterations for burn-in
for(i in 1:burn.in){
 x.curr=One.Step(x.curr)
}

accept=0
# And now 5000 iterates
x = x.curr
for(i in 2:5000){
 sto = x.curr
 x.curr=One.Step(x.curr)
 if(x.curr!=sto) accept=accept+1
 x[i]=x.curr
}

cat('Acceptance rate = ',accept/5000,'\n\n')


return(x)
}


There are lots of "fun" things we can do with this, such as:

 plot(Metrop(burn=0),type='l') # Shows the effect of not having any burn-in

plot(Metrop(jump=1),type='l')  # Shows the effect of having a "jump" that is too small (lots o' baby-steps)

acf(Metrop(jump=1)) # Shows the disastrous effect on the autocorrelation function

acf(Metrop()) # As versus a reasonable "jump" (the default value is =25)

quantile(Metrop(),c(0.025,0.5,0.975)) # Our 2.5%, 50%, and 97.5% quantiles

acf(Metrop()[seq(1,5000,5)]) # thinning the series to get rid of autocorrelation


But it is time to actually apply this to Lewis and Gowland's data, for which we will turn back to OpenBUGS:


model
{

# Diffuse priors for log normal age distribution
log.mu ~ dnorm(0,.001)
log.prec ~ dunif(0.001,100)
log.sd <- sqrt(1/log.prec)

for(i in 1:n) {
x[i] ~ dlnorm(log.mu,log.prec)
y.mu[i] <- -73.8884 + 8.829*sqrt(x[i])
y.prec[i] <- 1/pow(1.6216 + 0.009508*x[i],2)
fem[i] ~ dnorm(y.mu[i], y.prec[i])
}
}

#Data
list(n= 118,fem = c(61.2, 63.2, 65.9, 67.6, 66.8, 68.4, 68.7, 68.9,
69.9, 70.1, 71.1, 72.3, 73.2, 74.5, 74.6, 75.3, 76.3, 76.3, 76.8,
77, 77.6, 78, 79.1, 90, 90.2, 47.3, 49.5, 54.9, 53.6, 57.9, 58.9,
60.1, 59.3, 60.5, 63, 62.3, 64.6, 65.1, 66.3, 65.7, 65.3, 66.2, 67,
67.3, 68, 68.4, 68.7, 68.3, 70.9, 71.7, 72.7, 75, 75, 75.6, 76.4,
79.9, 84.9, 87.7, 85.9, 90.9, 90.9, 91.4, 92.2, 94.3, 94.1, 94.3,
94, 95, 96.4, 95.5, 98.8, 96, 97.2, 101.9, 102.7, 102.8, 52.6,
54.8, 61, 63, 67, 71, 68.8, 71.5, 72, 70.3, 72, 74, 74, 77, 78,
78, 78, 79, 77.2, 82, 86.8, 83, 83, 84, 86.2, 87, 89.9, 91, 91.5,
92, 93, 94, 97, 102, 102.9, 104, 104, 104, 102.5, 105, 110, 111))

In OpenBUGS we can get "log.mu" (when I ran 5,000 iterations after 1,000 iterations of burn-in I got 5.689) and "log.sd" (0.1808).  First let's see how this does in terms of recovering the distribution of femur lengths:

did.it.work = function (y=80)
{
log.mu=5.689
log.sd = 0.1808

pred.den = function(x)dlnorm(x,log.mu,log.sd)*Likelihood(x,y=y)
# integrate across age from 0 to 600 days LMP
integrate(pred.den,0,600)$val
}

plot.it1 = function (br=20)
{
dens=did.it.work(1)
for(i in 2:130) dens[i]=did.it.work(i)
max.dens = max(c(dens,hist(UK.bones$fem,br=br,plot=F)$density))

plot(1:130,dens,type='l',xlab='Femur Length',ylab='Density',
 lty=2,ylim=c(0,max.dens),main=paste('Desired Number of Bins =',br))
hist(UK.bones$fem,br=br,prob=T,add=T)
}

plot.it2 = function (bw=5)
{
dens=did.it.work(1)
for(i in 2:130) dens[i]=did.it.work(i)
plot(1:130,dens,type='l',xlab='Femur Length',ylab='Density',lty=2,
    main=paste('Kernel band width =',bw))
lines(density(UK.bones$fem,bw=bw))
}


Now we can use "plot.it1" to compare to the actual histogram of femur lengths.  Looks OK-ish, but we can change the number of desired bins to get different looking results.  So maybe we should use kernel density estimation for the actual femur lengths. Use "plot.it2" to compare the kernel density plot for the 118 femoral lengths to the predicted probability density function for femoral lengths.  Looks OK, but now monkeying around with the band-width gets us different results.

So let's use the empirical cumulative density function as versus the predicted cumulative density.  From the "archives," take a gander at: https://stat.ethz.ch/pipermail/r-help/2005-March/067052.html
(my problem at that point was that nobody realized the problem that you needed "Vectorize" within "integrate"):

plot.did.it.work = function ()
{
dens = 0
for(i in 1:80) dens[i]=integrate(Vectorize(did.it.work),0,i+39)$val
plot(40:(40+79),dens,type='l',xlab='Femur Length',
       ylab='Cumulative density',lty=2)
lines(ecdf(UK.bones$fem),do.points=F,verticals=T)
}


But now, let's go multivariate:

This time, instead of transforming the long bone lengths, let's make life simpler by transforming age.  Specifically, we will use (natural) log age.  The following "R" script does multivariate regression of five long bones (all except fibula) on log age.  Note that this is not the multiple regression of age onto long bones, but rather the regression of each long bone onto age.  It is just like simple bivariate regression, but we have to keep track of residual covariances between bones.  The following script does that and "dumps" some stuff we will use in OpenBUGS.

multivar.reg = function ()
{

lnAge = log(refbones[,3]*7)
sto=lm(as.matrix(refbones[,6:10])~lnAge)
theta.vcv = vcov(sto)
reord = c(seq(1,9,2),seq(2,10,2))
theta.mu = as.vector(round(coef(sto),2))[reord]
theta.prec = round(solve(matrix(theta.vcv[reord,reord],nc=10)),2)

dump('theta.mu','thetamu.txt')
dump('theta.prec','thetaprec.txt')

labs = c('Humerus','Radius','Ulna','Femur','Tibia')
new =  data.frame(lnAge = log(80:500))
par(mfrow=c(2,3))

for(i in 1:5){
plot(exp(lnAge),refbones[,5+i],ylab=labs[i],main=labs[i],xlab='Age')
sto2 = lm(refbones[,5+i]~lnAge)
preds=predict(sto2,newdata=new,interval='predict')
lines(80:500,preds[,1])
lines(80:500,preds[,2])
lines(80:500,preds[,3])

# lines(80:500,sto$coef[1,i]+sto$coef[2,i]*log(80:500))
}

cat('\n\n')
junk=readline('Hit anything except your neighbor: ')
par(mfrow=c(1,1))

y = predict(sto) - refbones[,6:10]

pairs(y,upper.panel=NULL)

N = length(lnAge)


resid.cov = cov(y)*(N-1)/(N-2)
tau = round(solve(resid.cov),5)
dump('tau','tau.txt')

}

And now it's off to OpenBUGS


model
{
# Posterior for Multivariate Regression

theta[1:10] ~ dmnorm(theta.mu[1:10],prec.theta[1:10,1:10])

# Prior for age

log.mu ~ dnorm(0,.001)
log.prec ~ dunif(0.0001,100)
log.sd <- sqrt(1/log.prec)

for(i in 1:N){
 x[i] ~ dnorm(log.mu, log.prec)
 age[i] <- exp(x[i])/7
 for(j in 1:5){
 Y.mu[i,j] <- theta[j] + theta[5+j] * x[i]
}
 Y[i,1:5] ~ dmnorm(Y.mu[i,], tau[,])
}
}

# Data
list(N = 118,
theta.mu=c(-214.22, -158.06, -181.92, -299.41, -228.11, 49.33, 37.05, 42.48, 66.63, 51.66),
prec.theta= structure(.Data=c(93.1, -29.62, -11.74, -20.28, -19.86, 513.67, -163.43, -64.76, -111.87, -109.59, -29.62, 372.48, -248.4,
 1.38, -32.31, -163.43, 2055.12, -1370.51, 7.64, -178.28, -11.74,
 -248.4, 255.79, 0.63, -15.57, -64.76, -1370.51, 1411.33, 3.47, -85.9,
 -20.28, 1.38, 0.63, 40.9, -34.11, -111.87, 7.64, 3.47, 225.64,
 -188.2, -19.86, -32.31, -15.57, -34.11, 105.88, -109.59, -178.28,
 -85.9, -188.2, 584.16, 513.67, -163.43, -64.76, -111.87, -109.59,  2849.48, -906.58, -359.27, -620.58, -607.95, -163.43, 2055.12,
 -1370.51, 7.64, -178.28, -906.58, 11400.43, -7602.68, 42.37, -988.98, -64.76, -1370.51, 1411.33, 3.47, -85.9, -359.27, -7602.68, 7829.11, 19.27, -476.52, -111.87, 7.64, 3.47, 225.64, -188.2, -620.58, 42.37,
19.27, 1251.68, -1043.99, -109.59, -178.28, -85.9, -188.2, 584.16,
-607.95, -988.98, -476.52, -1043.99, 3240.54),
.Dim=c(10,10)),
tau = structure(.Data=c(0.32326, -0.10285, -0.04076, -0.0704,
 -0.06897, -0.10285, 1.29333, -0.86249, 0.00481, -0.1122, -0.04076,
-0.86249, 0.88818, 0.00219, -0.05406, -0.0704, 0.00481, 0.00219,
 0.142, -0.11844, -0.06897, -0.1122, -0.05406, -0.11844, 0.36762),
.Dim=c(5,5)))

# WARNING! - load this data after you have loaded the list above.  Don't try one fell swoop.

Y[,1] Y[,2] Y[,3] Y[,4] Y[,5]
   NA   NA 48.6  61.2   NA
   NA 46.6   NA  63.2 54.6
 59.9 49.0   NA  65.9 58.7
 58.2   NA   NA    NA 60.0
 60.0   NA   NA  66.8   NA
 58.5 47.8 53.9  68.4 58.9
 60.9   NA   NA  68.7 60.6
 62.1 49.6 56.1  68.9 59.9
   NA   NA   NA  69.9   NA
   NA   NA 55.5  70.1 60.2
 63.8   NA 58.3  71.1 61.6
 64.7   NA   NA  72.3 62.0
 65.4 52.7 61.8  73.2 64.1
 63.9 51.7 58.6    NA 65.4
 64.5   NA   NA    NA   NA
   NA   NA   NA  75.3   NA
 65.7   NA   NA    NA   NA
 66.6 56.5   NA  76.3   NA
 66.0   NA 61.7  76.8 69.8
 64.5 51.2 59.2  77.0 67.1
 66.5 52.8   NA  77.6 65.3
 67.4 54.6   NA    NA   NA
 68.2 55.6 62.9  79.1   NA
 73.5 56.5 64.7  90.0   NA
 76.5 62.4 70.6    NA   NA
 44.8   NA 41.8    NA   NA
   NA   NA 42.3  49.5 44.8
 47.4   NA   NA  54.9 47.1
 50.0   NA   NA  53.6   NA
 53.1 43.0 49.2  57.9 52.3
 54.3 43.2 49.0  58.9 54.4
   NA   NA 49.9    NA   NA
 54.3   NA 49.6  59.3   NA
 54.6   NA   NA    NA   NA
 55.3 44.7 51.3  63.0 56.0
 56.1   NA 53.0    NA   NA
 57.5   NA   NA    NA   NA
 58.4   NA 53.1  65.1   NA
 57.8 46.9 52.1    NA 59.2
   NA 48.8 54.3  65.7 59.4
 58.3   NA 51.8    NA 56.2
 59.6 46.1 52.7  66.2 57.9
 59.2   NA   NA    NA   NA
   NA 48.9 55.2  67.3 60.5
 60.5 47.5 53.9  68.0 57.9
   NA   NA   NA  68.4   NA
 57.4 46.8 54.5  68.7 58.9
 60.4 49.2   NA    NA   NA
 61.5   NA 57.3  70.9 62.8
 64.2   NA   NA  71.7 61.1
 65.0 54.0 62.3  72.7   NA
 64.6 51.4 59.3    NA   NA
 64.5   NA 59.4  75.0 62.0
   NA   NA   NA  75.6 62.3
 68.5   NA   NA  76.4   NA
   NA 52.9   NA  79.9 68.1
 71.8   NA   NA    NA   NA
 74.8 61.9   NA  87.7 74.8
 74.4 59.4 66.7    NA 73.2
   NA 60.1   NA    NA   NA
 76.9   NA   NA  90.9   NA
   NA   NA   NA  91.4   NA
 78.0   NA   NA  92.2 79.8
 78.1   NA 67.2    NA   NA
 78.3   NA   NA    NA   NA
 78.6   NA 69.6    NA   NA
 80.4 65.3 74.0  94.0   NA
 79.0   NA   NA  95.0   NA
 79.3 60.8   NA    NA   NA
 79.3   NA   NA    NA   NA
 81.3 65.5 73.5  98.8   NA
 82.0 65.2 70.4  96.0 85.0
 83.2   NA   NA  97.2 83.6
 83.8   NA   NA    NA   NA
 84.7   NA   NA 102.7 87.0
 85.0   NA 76.1    NA   NA
 48.0 38.0 42.0    NA 47.0
 51.0 42.0 48.0    NA 49.0
 52.0 42.0 49.0  61.0   NA
   NA   NA   NA  63.0   NA
 61.0   NA   NA  67.0 60.0
 61.0 48.0 54.0  71.0 61.0
   NA   NA   NA    NA 60.0
 62.0 49.0   NA    NA   NA
 62.0 50.0 58.0  72.0   NA
 62.0 51.0 59.0    NA   NA
   NA 54.0 61.5  72.0 63.0
 64.0 50.0 57.0  74.0 63.0
 65.0 53.0 59.0  74.0 63.0
 67.0 52.0 60.0  77.0 66.0
 67.0 54.0 63.0  78.0   NA
 69.0 54.0 61.0  78.0 67.0
 67.0 54.0 61.0  78.0 69.0
 70.0 58.0 65.0  79.0 70.0
 68.0 55.0 61.0    NA 67.0
   NA 52.0   NA  82.0 69.0
 71.0 56.0 64.0    NA 75.0
 71.0   NA   NA  83.0 71.0
 69.0 57.0 64.0  83.0 72.5
 71.0 56.0 63.0  84.0 71.0
   NA 58.0 67.0    NA   NA
 74.0 57.0 64.0  87.0 75.0
 73.5 55.0 62.0  89.9 73.0
 77.0 61.5 67.0  91.0 80.0
   NA 60.0 67.0  91.5 76.0
 79.0 61.0 68.0  92.0 80.0
 76.0 58.0   NA  93.0 79.0
   NA 62.0 69.0  94.0 82.5
   NA 60.0 66.0  97.0 84.0
 80.5 62.0 68.0 102.0 85.0
 84.5   NA   NA    NA   NA
 83.0   NA 67.0 104.0 84.0
 82.0 64.0 71.0 104.0 90.0
 80.0 61.5 70.0 104.0 85.0
   NA   NA   NA    NA 87.0
 85.0 64.0 73.0 105.0 88.0
 87.0 64.0 71.0 110.0 89.0
 89.0 68.0 72.0 111.0 92.5
END

Now you can monitor Y, Y.mu, log.mu, log.sd, and age.  Fun times!  Note how Y will contain estimates for missing long bone data.  Y.mu  are predicted long bone lengths, log.mu and log.sd are as before, and "age" contains estimated age in weeks LMP.
References