This page gives the "R" code for analyses reported in a manuscript submitted to the Journal of Forensic Sciences.
If you do not already have "R" you can get it from:
http://cran.r-project.org/
You will also need to download two specific "R" library, these being maxLik and numDeriv. If you want to use the shiny app at the end of this file you will also need the shiny library.
Test Data from Roche, Wainer, and Thissen (1975)
Roche AF, Wainer H, and Thissen D. 1975. Skeletal maturity: the knee joint as a biological indicator. New York: Plenum Medical Book Co., 374 p. (Springer is the current copyright holder)
First things first, anything shown in
fixed-width green text
is intended to be copied and pasted into your R "Gui" (graphical user
interface). Below is the data Roche, Wainer, and Thissen (1975)
gave for test cases. Because this book is under copyright to
Springer, you will need to obtain the data from page 354 of the book.
scores=structure(list(ID=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
Sex=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
Age=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEMMW=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEMEW=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEMEH=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEMWLC=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEMHLC=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIBMW=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIBEH=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIBEW=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIBEW=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIBMW=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.D=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.E=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.F=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.G=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.H=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.J=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.K=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.L=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FEM.M=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.C=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.D=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.E=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.F=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.G=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.H=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.J=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.K=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.L=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.M=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.N=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.P=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.Q=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
TIB.R=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIB.B=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIB.C=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIB.D=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIB.E=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
FIB.F=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)),
class="data.frame",row.names=c(NA,-27L))
Parameters
Below
are the parameters from Roche, Wainer, and Thissen for females and for
males. Again, because the RWT book is under copyright by
Springer, you will need to obtain these parameters from Tables 23 and
24 (pages 199 and 200). Alternatively, you can obtain the
parameters from RWT's BASIC code which will get you an extra digit.
males=structure(list(
tau1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
tau2=c(NA,NA,NA,NA,Inf,Inf,Inf,NA,Inf,NA,NA,NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,NA,Inf,NA,Inf,Inf,NA),
tau3=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,NA),
tau4=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,Inf),
d=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)),
row.names=c(NA,34L),class="data.frame")
females=structure(list(
tau1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
tau2=c(NA,NA,NA,NA,Inf,Inf,Inf,NA,Inf,NA,NA,NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,NA,Inf,NA,Inf,Inf,NA),
tau3=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,NA),
tau4=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,Inf),
d=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)),
row.names=c(NA,34L),class="data.frame")
Figure 1A and Figure 1B
Fig.01A = function ()
{
library(shape)
library(pBrackets)
options(scipen = 999)
p.vals = c(.01,.05,.1,.2,.5,1,2,5,seq(10,90,10),95,98,99,99.8,99.9,99.99)/100
plot(c(0,8),qlogis(c(0.0001,0.9999)),axes=F,xlab='Age',
ylab='Probability of individuals in stage',type='n',xaxs='i',yaxs='i',
main='Female FIB-C')
box()
par(cex=0.75)
axis(2,at=qlogis(p.vals),labels=p.vals,las=1)
par(cex=1)
axis(1)
axis(4,las=1,seq(-8,8,2))
tau.1 = females$tau1[31]
tau.2 = females$tau2[31]
slope = females$d[31]
inter.y1 = -tau.1*slope
print(round(inter.y1,4))
inter.y2 = -tau.2*slope
print(round(inter.y2,4))
print(round(slope,4))
abline(inter.y1,slope,lwd=2)
abline(inter.y2,slope,lty=2,lwd=2)
brackets(3.422,qlogis(.9999),3.422,inter.y1+slope*3.422,type=5,h=0.3)
brackets(3.422,inter.y1+slope*3.422,3.422,inter.y2+slope*3.422,type=5,h=0.3)
brackets(3.422,inter.y2+slope*3.422,3.422,qlogis(0.0001),type=5,h=0.3)
abline(h=inter.y1+slope*3.422,lty=3)
abline(h=inter.y2+slope*3.422,lty=3)
abline(v=3.422)
text(3.9,(qlogis(0.9999)+inter.y1+3.422*slope)/2,'p = 0.2954',adj=c(0,.5))
text(4.1,(inter.y2+3.422*slope+inter.y1+3.422*slope)/2,'p = 0.4095',adj=c(0,.5))
text(3.9,(qlogis(0.0001)+inter.y2+3.422*slope)/2,'p = 0.2951',adj=c(0,.5))
text(5,-7.46,expression(paste(y[1],' = ', -5.2916 + 1.8004%*%'Age')),adj=0)
text(5,-8.39,expression(paste(y[2],' = ', -7.0318 + 1.8004%*%'Age')),adj=0)
lines(c(4.25,4.75),rep(-7.46,2),lwd=2)
lines(c(4.25,4.75),rep(-8.39,2),lwd=2,lt=2)
rect(4.1,qlogis(0.0001),8,-6.9)
}
Fig.01A()
Fig.01B=function ()
{
library(maxLik)
# Female FIB-C
tau.1 = females$tau1[31]
tau.2 = females$tau2[31]
slope = females$d[31]
inter.y1 = -tau.1*slope
inter.y2 = -tau.2*slope
prob.stage2 = function(CA){
log(plogis(inter.y1+slope*CA)-plogis(inter.y2+slope*CA))
}
sto=maxNR(prob.stage2,start=3)
se=as.numeric(sqrt(1/(-sto$hess)))
z=qnorm(0.975)
left = sto$est-z*se
right = sto$est+z*se
age=seq(0,8,.01)
plot(age,dnorm(age,sto$est,se),type='l',yaxs='i',ylim=c(0,0.5),
xlim=c(0,7),xaxs='i',lwd=2,xlab='Age',ylab='Density',
main='Female FIB-C Stage 2')
lines(rep(left,2),c(0,dnorm(left,sto$est,se)),lwd=2)
lines(rep(right,2),c(0,dnorm(right,sto$est,se)),lwd=2)
age=c(left,seq(left,right,0.01),right)
n=length(age)
polygon(age,c(0,dnorm(age[-c(1,n)],sto$est,se),0),density=15,border=NA)
cat('\n Normal\n')
print(round(c(left,sto$est,right),3))
###############################
prob.stage2 = function(CA,denom=1){
(plogis(inter.y1+slope*CA)-plogis(inter.y2+slope*CA))/denom
}
denom=integrate(prob.stage2,lower=0,upper=30)$val
age=seq(0,8,.001)
lines(c(0,age),c(0,prob.stage2(age,denom=denom)),lwd=2,lty=2)
mle = optimize(prob.stage2,interval=c(0,20),maximum=T)$max
print(round(mle,3))
same.height = function(CA,h){
prob.stage2(CA,denom=denom)-h
}
find.95 = function(h){
left = uniroot(same.height,lower=0,upper=mle,h=h)$root
right = uniroot(same.height,lower=mle,upper=20,h=h)$root
integrate(prob.stage2,lower=left,upper=right,denom=denom)$val-.95
}
h = uniroot(find.95,lower=0.01,upper=0.2)$root
left = round(uniroot(same.height,lower=0,upper=mle,h=h)$root,3)
right = round(uniroot(same.height,lower=mle,upper=20,h=h)$root,3)
lines(rep(left,2),c(0,prob.stage2(left,denom=denom)),lwd=2,lty=2)
lines(rep(right,2),c(0,prob.stage2(right,denom=denom)),lwd=2,lty=2)
age=c(left,seq(left,right,0.01),right)
n=length(age)
polygon(age,c(0,prob.stage2(age[-c(1,n)],denom=denom),0),density=15,angle=-45,border=NA)
legend(x='topleft',legend=c('Normal','Posterior'),lty=c(1,2),lwd=c(2,2),bty='n')
legend(x='topright',legend=c('95% CI','95% HPD'),density=c(25,25),
angle=c(45,-45),bty='n')
cat('\n Posterior\n')
print(round(c(left,mle,right),3))
abline(v=mle,lty=2)
}
Fig.01B()
SAGE
Roche,
Wainer, and Thissen (1975) referred to their program as "SAGE"
(presumably for "skeletal age"). Here is an "R" version:
SAGE = function (id=72,sex=2,guess=0.3,scores=c(25,12,0,0,0,20.5,0,12,0,0,
2,2,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),HPD=95)
{
library(maxLik)
max.stages=c(5,5,5,3,2,2,2,3,2,3,3,3,5,5,2,2,2,2,2,2,2,2,2,2,2,2,3,
3,5,2,3,2,2,4)
if(sex==2){
tau=females[,1:4]
d=females[,5]
}
else{
tau=males[,1:4]
d=males[,5]
}
bounds=c(0,.5,.6,.7,.8,1000,
0,2,2.25,2.5,2.75,1000,
0,1,1.1,1.2,1.3,1000,
0,.6,.7,.8,.9,1000,
0,2.5,2.7,2.9,3.1,1000,
0,.3,.5,.7,.9,1000)
put.metrics=c(1,2,3,13,14,29)
bounds=matrix(bounds,nc=6,byrow=T)
grades=rep(0,34)
ratios=rep(0,6)
if(scores[1]>0) ratios[1]=scores[2]/scores[1]
if(scores[3]>0) ratios[2]=scores[2]/scores[3]
if(scores[5]>0) ratios[3]=scores[4]/scores[5]
if(scores[6]>0) ratios[4]=scores[8]/scores[6]
if(scores[7]>0) ratios[5]=scores[8]/scores[7]
if(scores[10]>0) ratios[6]=scores[9]/scores[10]
ratios=as.numeric(ratios)
for(i in 1:6){
if(ratios[i]>0) grades[put.metrics[i]]=which.max(hist(ratios[i],
br=bounds[i,],plot=F)$counts)
}
grades[c(4:12,15:28,30:34)]=scores[11:38]
grades=as.numeric(grades)
LogLike=function(Age,denom=1,dens=F){
LL=0
for(i in 1:34){
if(grades[i]>0){
if(grades[i]==1) LL=LL+log(1-1/(1+exp(-d[i]*(Age-tau[i,1]))))
else if(grades[i]==max.stages[i]) LL=LL+log(1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1]))))
else {
p1 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1])))
p2 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]])))
LL = LL + log(p1-p2)
}
}
}
if(dens==F) return(LL/denom)
return(exp(LL)/denom)
}
model=maxNR(LogLike,start=guess,dens=F)
mu=model$estimate
se=sqrt(1/as.numeric(-model$hessian))
lower=mu-8*mu
if(lower<(-10)) lower=0
upper=mu+8*se
if(upper>20) upper=30
full.height=dnorm(mu,mu,se)
find.right=function(right,denom,dens=T){
return(integrate(LogLike,lower=lower,upper=right,denom=denom,dens=T)$val-0.99)
}
denom=integrate(LogLike,lower=lower,upper=upper,denom=1,dens=T)$val
height.Like = LogLike(mu,denom=denom,dens=T)
height.Norm = dnorm(mu,mu,se)
height=1.1*max(c(height.Like,height.Norm))
top=uniroot(find.right,lower=lower,upper=upper,denom=denom,dens=T)$root
top=top+2
Age=seq(-.2,top,.001)
N=length(Age)
ll=vector()
for(i in 1:N) ll[i]=LogLike(Age[i],denom=denom,dens=T)
plot(Age,ll,type='l',ylim=c(0,height),xlim=c(mu-4*se,mu+4*se),xaxs='i',yaxs='i',ylab='Density')
legend(x='topright',bty='n',legend=c('Normal','Posterior'),lty=c(2,1),lwd=c(2,1))
find.height=function(age,try){
LogLike(age,denom=denom,dens=T)-try
}
find.hpd=function(try){
left=uniroot(find.height,interval=c(-10,mu),try)$root
right=uniroot(find.height,interval=c(mu,30),try)$root
area=integrate(LogLike,lower=left,upper=right,denom=denom,dens=T)$val-HPD/100
return(area)
}
hpd.height=uniroot(find.hpd,interval=c(0.000001*height.Like,height.Like))$root
left=uniroot(find.height,interval=c(lower,mu),try=hpd.height)$root
right=uniroot(find.height,interval=c(mu,upper),try=hpd.height)$root
L.area=integrate(LogLike,lower=lower,upper=left,denom=denom,dens=T)$val
R.area=integrate(LogLike,lower=right,upper=upper,denom=denom,dens=T)$val
Ages=seq(left,right,.001)
density=LogLike(Ages,denom=denom,dens=T)
Ages=c(Ages[1],Ages,Ages[length(Ages)])
density=c(0,density,0)
polygon(Ages,density,col='gray')
lines(Age,ll)
lines(Age,dnorm(Age,mu,se),lty=2,lwd=2)
abline(v=guess,lty=2)
lines(rep(mu-1.96*se,2),c(0,dnorm(mu-1.96*se,mu,se)),lwd=2)
lines(rep(mu+1.96*se,2),c(0,dnorm(mu+1.96*se,mu,se)),lwd=2)
return(list(est=round(mu,2),se=round(se,2),left=round(left,2),right=round(right,2)))
}
Results from Roche, Wainer, and Thissen (1975)
Paste the text below to your R Gui. Once again, you will
need to get these numbers from the RWT book. Specifically, these
are "AGE," "SA" (RWT.mu), and "SE" (RWT.se) from pages 355-356.
out=structure(list(
Age=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
RWT.mu=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
RWT.se=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)),
row.names=c(NA,27L),class="data.frame")
Run the Roche, Wainer, and Thissen (1975) examples
run.SAGE = function ()
{
out2=matrix(NA,nr=27,nc=2)
cat('\n')
for(i in 1:27){
sex=scores[i,2]
guess=scores[i,3]
curr.scores=scores[i,-(1:3)]
temp=SAGE(sex=sex,guess=guess,scores=curr.scores)
out2[i,1]=temp$est
out2[i,2]=temp$se
cat('\r',i,' of 27')
flush.console()
}
cat('\n')
out=cbind(out,out2)
colnames(out)=c(colnames(out)[1:3],'est','se')
out<<-out
}
run.SAGE()
Figure 2
Fig.02 =function ()
{
plot(out[,2],out[,4],xlim=c(0,17.5),ylim=c(0,17.5),
xaxs='i',yaxs='i',xlab='RWT Estimated skeletal age',
ylab='Estimated skeletal age')
points(out[c(18,4,15,17,8,9,13,24),2],out[c(18,4,15,17,8,9,13,24),4],pch=19)
abline(0,1,lty=2)
max(abs(out[,2]-out[,4]))
}
Fig.02()
Figure 3
Fig.03 = function ()
{
plot(out[,3],out[,5],xlim=c(0.1,1),ylim=c(0.1,1),
xaxs='i',yaxs='i',xlab='RWT Estimated standard error',
ylab='Estimated standard error')
points(out[c(4,6,8,9,13),3],out[c(4,6,8,9,13),5],pch=19)
abline(0,1,lty=2)
max(abs(out[,3]-out[,5]))
}
Fig.03()
Figures 4 & 5
Fig.04.to.05 = function (Fig=4)
{
i=23
if(Fig==5)i=14
sex=scores[i,2]
guess=scores[i,3]
curr.scores=scores[i,-(1:3)]
SAGE(sex=sex,guess=guess,scores=curr.scores)
}
Fig.04.to.05()
And take a look, then:
Fig.04.to.05(5)
Figures 6 & 7
Fig.06.to.07 = function (Fig=6)
{
max.stages=c(5,5,5,3,2,2,2,3,2,3,3,3,5,5,2,2,2,2,2,2,
2,2,2,2,2,2,3,3,5,2,3,2,2,4)
tau=males[,1:4]
d=males[,5]
if(Fig==7){
tau=females[,1:4]
d=females[,5]
}
bounds=c(0,.5,.6,.7,.8,1000,
0,2,2.25,2.5,2.75,1000,
0,1,1.1,1.2,1.3,1000,
0,.6,.7,.8,.9,1000,
0,2.5,2.7,2.9,3.1,1000,
0,.3,.5,.7,.9,1000)
put.metrics=c(1,2,3,13,14,29)
bounds=matrix(bounds,nc=6,byrow=T)
grades=max.stages
LogLike=function(Age,denom=1,dens=F){
LL=0
for(i in 1:34){
if(grades[i]>0){
if(grades[i]==1) LL=LL+log(1-1/(1+exp(-d[i]*(Age-tau[i,1]))))
else
if(grades[i]==max.stages[i])
LL=LL+log(1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1]))))
else {
p1 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1])))
p2 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]])))
LL = LL + log(p1-p2)
}
}
}
if(dens==F) return(-LL/denom)
return(exp(LL)/denom)
}
height=LogLike(30,denom=1,dens=T)
find.Age=function(Age){
LogLike(Age,denom=height,dens=T)-.5
}
med.Age=uniroot(find.Age,lower=0,upper=30)$root
Age=seq(0,30,.001)
N=length(Age)
ll=vector()
for(i in 1:N) ll[i]=LogLike(Age[i],denom=height,dens=T)
if(Fig==6) my.title=paste('Males (median age = ',round(med.Age,2),')',sep='')
if(Fig==7) my.title=paste('Females (median age = ',round(med.Age,2),')',sep='')
plot(Age,ll,type='l',xlim=c(10,30),ylim=c(0,1),
xaxs='i',yaxs='i',lwd=2,main=my.title,ylab='Probability')
lines(rep(med.Age,2),c(0,0.5))
lines(c(0,med.Age),rep(0.5,2),lty=2)
}
Fig.06.to.07()
And take a look, then:
Fig.06.to.07(7)
Figures 8 & 9
Fig.08.to.09 = function (Fig=8)
{
library(numDeriv)
max.stages=c(5,5,5,3,2,2,2,3,2,3,3,3,5,5,2,2,2,2,2,2,
2,2,2,2,2,2,3,3,5,2,3,2,2,4)
tau=males[,1:4]
d=males[,5]
if(Fig==9){
tau=females[,1:4]
d=females[,5]
}
bounds=c(0,.5,.6,.7,.8,1000,
0,2,2.25,2.5,2.75,1000,
0,1,1.1,1.2,1.3,1000,
0,.6,.7,.8,.9,1000,
0,2.5,2.7,2.9,3.1,1000,
0,.3,.5,.7,.9,1000)
put.metrics=c(1,2,3,13,14,29)
bounds=matrix(bounds,nc=6,byrow=T)
grades=max.stages
LogLike=function(Age,denom=1,dens=F){
LL=0
for(i in 1:34){
if(grades[i]>0){
if(grades[i]==1) LL=LL+log(1-1/(1+exp(-d[i]*(Age-tau[i,1]))))
else
if(grades[i]==max.stages[i])
LL=LL+log(1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1]))))
else {
p1 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1])))
p2 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]])))
LL = LL + log(p1-p2)
}
}
}
if(dens==F) return(-LL/denom)
return(exp(LL)/denom)
}
height=LogLike(30,denom=1,dens=T)
find.mode=function(Age){
grad(LogLike,Age,denom=height,dens=T)
}
find.height=function(Age,try){
grad(LogLike,Age,denom=height,dens=T)-try
}
mode=optimize(find.mode,lower=0,upper=30,maximum=T)$max
pdf=function(Age){
grad(LogLike,Age,denom=height,dens=T)
}
find.hpd=function(try){
left=uniroot(find.height,interval=c(0,mode),try)$root
right=uniroot(find.height,interval=c(mode,30),try)$root
area=integrate(pdf,lower=left,upper=right)$val-0.95
return(area)
}
hpd.height=uniroot(find.hpd,interval=c(0.005,0.1))$root
left=uniroot(find.height,interval=c(0,mode),try=hpd.height)$root
right=uniroot(find.height,interval=c(mode,30),try=hpd.height)$root
L.area=integrate(pdf,lower=0,upper=left)$val
R.area=integrate(pdf,lower=right,upper=30)$val
print(c(L.area,R.area,L.area+R.area))
Age=seq(15,25,.001)
pdf=grad(LogLike,Age,denom=height,dens=T)
plot(Age,pdf,type='l',xlim=c(15,25),xaxs='i',ylim=c(0,.5),yaxs='i',ylab='Density')
Age=seq(left,right,.001)
pdf=grad(LogLike,Age,denom=height,dens=T)
Age=c(left,Age,right)
pdf=c(0,pdf,0)
polygon(Age,pdf,density=30)
if(Fig==9){
mtext(paste('Females (',round(left,2),', ',round(mode,2),', ',round(right,2),')',sep=''),cex=1.25)
return('All done')
}
mtext(paste('Males (',round(left,2),', ',round(mode,2),', ',round(right,2),')',sep=''),cex=1.25)
}
Fig.08.to.09()
And take a look, then:
Fig.08.to.09(9)
Figures 10 & 11
Fig.10.to.11 = function (Fig=10)
{
max.stages=c(5,5,5,3,2,2,2,3,2,3,3,3,5,5,2,2,2,2,2,2,
2,2,2,2,2,2,3,3,5,2,3,2,2,4)
tau=males[,1:4]
d=males[,5]
if(Fig==11){
tau=females[,1:4]
d=females[,5]
}
bounds=c(0,.5,.6,.7,.8,1000,
0,2,2.25,2.5,2.75,1000,
0,1,1.1,1.2,1.3,1000,
0,.6,.7,.8,.9,1000,
0,2.5,2.7,2.9,3.1,1000,
0,.3,.5,.7,.9,1000)
put.metrics=c(1,2,3,13,14,29)
bounds=matrix(bounds,nc=6,byrow=T)
grades=max.stages
LogLike=function(Age,denom=1,dens=F){
LL=0
for(i in 1:34){
if(grades[i]>0){
if(grades[i]==1) LL=LL+log(1-1/(1+exp(-d[i]*(Age-tau[i,1]))))
else
if(grades[i]==max.stages[i])
LL=LL+log(1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1]))))
else {
p1 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1])))
p2 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]])))
LL = LL + log(p1-p2)
}
}
}
if(dens==F) return(-LL/denom)
return(exp(LL)/denom)
}
height=LogLike(30,denom=1,dens=T)
Age=seq(15,21,.001)
norm=dnorm(Age,18,.8067)
if(Fig==10) my.title='Males'
if(Fig==11) my.title='Females'
plot(Age,norm,type='l',xlim=c(15,21),ylim=c(0,1),
xaxs='i',yaxs='i',ylab='Probability',main=my.title)
Age=c(Age[1],Age,Age[length(Age)])
norm=c(0,norm,0)
polygon(Age,norm,col='gray')
Age=seq(15,21,.001)
N=length(Age)
ll=vector()
for(i in 1:N) ll[i]=LogLike(Age[i],denom=height,dens=T)
lines(Age,ll,lwd=2)
average=function(x){
return(LogLike(x,denom=height,dens=T)*dnorm(x,18,.8067))
}
denom=integrate(average,lower=10,upper=18)$val
num=integrate(average,lower=18,upper=30)$val
lines(c(0,18),rep(denom,2),lty=2,lwd=2)
lines(c(18,30),rep(num,2),lty=3,lwd=2)
abline(v=18,lty=3)
num/denom
}
Fig.10.to.11()
And take a look, then:
Fig.10.to.11(11)
" For the
95% HPD, 24 out of the 27 (88.9%) individuals fall within the interval. On an exact binomial test this gives a
probability-value of 0.1505, suggesting that the coverage is correct. For the 50% HPD, 12 out of 27 (44.4%)
individuals fall within the interval. On
an exact binomial test this gives a probability-value of 0.7011, again
suggesting that the coverage is correct."
HPD.SAGE = function (HPD=95)
{
cat('\n')
in.HPD=0
for(i in 1:27){
sex=scores[i,2]
guess=scores[i,3]
curr.scores=scores[i,-(1:3)]
temp=SAGE(sex=sex,guess=guess,scores=curr.scores,HPD=HPD)
if(out[i,1]>=temp$left & out[i,1]<=temp$right) in.HPD=in.HPD+1
cat('\r',i,' of 27: within HPD = ',in.HPD)
flush.console()
}
cat('\n\n')
cat(in.HPD,' of 27 in HPD')
cat('\n prob value = ',round(binom.test(in.HPD,27,HPD/100)$p.value,4))
cat('\n')
}
For the
95% HPD run with the default as follows:
HPD.SAGE()
For the 50% HPD run as follows:
HPD.SAGE(50)
shiny App
To
use this app you will need the R library "shiny." Then paste the
following into your R Gui. You will need to press the “stop”
button when you are done using the app. Alternatively, if you are
using R
studio open a new project as a shiny app and replace the text in the
upper left-hand box with the text below. Then press the button
for "run app" in the upper right-hand corner of the upper left-hand box.
library(shiny)
ui <- fluidPage(
# Application title
titlePanel("Roche, Wainer, and Thissen (1975) Skeletal Maturation: The Knee Joint as a Biological Indicator"),
fluidRow(
column(1,
radioButtons("Sex", h4("Sex"),
choices = list("Male" = 1, "Female" = 2),selected = 2)),
column(1,
radioButtons("find.HPD", h4("Find HPD?"),
choices = list("Yes" = 1, "No" = 2),selected = 2)),
column(3,
sliderInput("HPD",
h4("Percentage HPD"),min=1,max=99,value=95)),
column(2,numericInput("prior.mu",
h4("Prior mean"),value=18.0)),
column(2,numericInput("prior.sd",
h4("Prior s.d."),value=0.8067)),
column(2,numericInput("Guess",
h4("Age guesstimate"),value=4.5))
),
tags$h4("Zeroes represent missing data throughout the input"),
fluidRow(
column(1,
radioButtons("Maxed", h4("All ratios at max."),
choices = list("Yes" = 1, "No" = 2),selected = 2)),
column(2,
numericInput("FEMMW",
h4("Fem metaph W"),value=53)),
column(2,
numericInput("FEMEW",
h4("Fem epiph W"),value=53.5)),
column(2,
numericInput("FEMEH",
h4("Fem epiph H"),value=16.5)),
column(2,
numericInput("FEMWLC",
h4("Fem lat cond W"),value=0)),
column(2,
numericInput("FEMHLC",
h4("Fem lat cond H"),value=0))
),
fluidRow(
column(2,
numericInput("TIBMW",
h4("Tib metaph W"),value=45.5)),
column(2,
numericInput("TIBEH",
h4("Tib epiph H"),value=15.5)),
column(2,
numericInput("TIBEW",
h4("Tib epiph W"),value=44)),
column(2,
numericInput("FIBEW",
h4("Fib epiph W"),value=9)),
column(2,
numericInput("FIBMW",
h4("Fib metaph W"),value=14))
),
fluidRow(
column(1,
radioButtons("FEM.D", h4("FEM-D"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3" = 3),selected = 3)),
column(1,
radioButtons("FEM.E", h4("FEM-E"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("FEM.F", h4("FEM-F"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("FEM.G", h4("FEM-G"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("FEM.H", h4("FEM-H"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3" = 3),selected = 3)),
column(1,
radioButtons("FEM.J", h4("FEM-J"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("FEM.K", h4("FEM-K"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3" = 3),selected = 1)),
column(1,
radioButtons("FEM.L", h4("FEM-L"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3" = 3),selected = 1)),
column(1,
radioButtons("FEM.M", h4("FEM-M"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3" = 3),selected = 1))
),
fluidRow(
column(1,
radioButtons("TIB.C", h4("TIB-C"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("TIB.D", h4("TIB-D"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("TIB.E", h4("TIB-E"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("TIB.F", h4("TIB-F"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("TIB.G", h4("TIB-G"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("TIB.H", h4("TIB-H"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("TIB.J", h4("TIB-J"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("TIB.K", h4("TIB-K"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("TIB.L", h4("TIB-L"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("TIB.M", h4("TIB-M"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("TIB.N", h4("TIB-N"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("TIB.P", h4("TIB-P"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1))
),
fluidRow(
column(1,
radioButtons("TIB.Q", h4("TIB-Q"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3"=3),selected = 1)),
column(1,
radioButtons("TIB.R", h4("TIB-R"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3"=3),selected = 1)),
column(1,
radioButtons("FIB.B", h4("FIB-B"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 2)),
column(1,
radioButtons("FIB.C", h4("FIB-C"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3"=3),selected = 2)),
column(1,
radioButtons("FIB.D", h4("FIB-D"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("FIB.E", h4("FIB-E"),
choices = list("0" = 0, "1" = 1,
"2" = 2),selected = 1)),
column(1,
radioButtons("FIB.F", h4("FIB-F"),
choices = list("0" = 0, "1" = 1,
"2" = 2,"3"=3,"4"=4),selected = 1))
),
submitButton("Submit"),
mainPanel(
textOutput("scores"),
textOutput("scores3"),
textOutput("scores4"),
plotOutput("score")
)
)
# Define server logic required to plot posterior density
server <- shinyServer(function(input, output) {
output$score<-renderPlot({
input$Submit
scores=vector()
scores[1]=input$FEMMW
scores[2]=input$FEMEW
scores[3]=input$FEMEH
scores[4]=input$FEMWLC
scores[5]=input$FEMHLC
scores[6]=input$TIBMW
scores[7]=input$TIBEH
scores[8]=input$TIBEW
scores[9]=input$FIBEW
scores[10]=input$FIBMW
grades=rep(0,34)
ratios=rep(0,6)
bounds=c(0,.5,.6,.7,.8,1000,
0,2,2.25,2.5,2.75,1000,
0,1,1.1,1.2,1.3,1000,
0,.6,.7,.8,.9,1000,
0,2.5,2.7,2.9,3.1,1000,
0,.3,.5,.7,.9,1000)
put.metrics=c(1,2,3,13,14,29)
bounds=matrix(bounds,nc=6,byrow=T)
if(scores[1]>0) ratios[1]=scores[2]/scores[1]
if(scores[3]>0) ratios[2]=scores[2]/scores[3]
if(scores[5]>0) ratios[3]=scores[4]/scores[5]
if(scores[6]>0) ratios[4]=scores[8]/scores[6]
if(scores[7]>0) ratios[5]=scores[8]/scores[7]
if(scores[10]>0) ratios[6]=scores[9]/scores[10]
ratios=as.numeric(ratios)
for(i in 1:6){
if(ratios[i]>0) grades[put.metrics[i]]=which.max(hist(ratios[i],
br=bounds[i,],plot=F)$counts)
}
if(input$Maxed==1) grades[put.metrics]=5
grades[4]=as.numeric(input$FEM.D)
grades[5]=as.numeric(input$FEM.E)
grades[6]=as.numeric(input$FEM.F)
grades[7]=as.numeric(input$FEM.G)
grades[8]=as.numeric(input$FEM.H)
grades[9]=as.numeric(input$FEM.J)
grades[10]=as.numeric(input$FEM.K)
grades[11]=as.numeric(input$FEM.L)
grades[12]=as.numeric(input$FEM.M)
grades[15]=as.numeric(input$TIB.C)
grades[16]=as.numeric(input$TIB.D)
grades[17]=as.numeric(input$TIB.E)
grades[18]=as.numeric(input$TIB.F)
grades[19]=as.numeric(input$TIB.G)
grades[20]=as.numeric(input$TIB.H)
grades[21]=as.numeric(input$TIB.J)
grades[22]=as.numeric(input$TIB.K)
grades[23]=as.numeric(input$TIB.L)
grades[24]=as.numeric(input$TIB.M)
grades[25]=as.numeric(input$TIB.N)
grades[26]=as.numeric(input$TIB.P)
grades[27]=as.numeric(input$TIB.Q)
grades[28]=as.numeric(input$TIB.R)
grades[30]=as.numeric(input$FIB.B)
grades[31]=as.numeric(input$FIB.C)
grades[32]=as.numeric(input$FIB.D)
grades[33]=as.numeric(input$FIB.E)
grades[34]=as.numeric(input$FIB.F)
library(maxLik)
# Replace with values from book, per above
females=structure(list(
tau1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
tau2=c(NA,NA,NA,NA,Inf,Inf,Inf,NA,Inf,NA,NA,NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,NA,Inf,NA,Inf,Inf,NA),
tau3=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,NA),
tau4=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,Inf),
d=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)),
row.names=c(NA,34L),class="data.frame")
# Replace with values from book, per above
males=structure(list(
tau1=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
tau2=c(NA,NA,NA,NA,Inf,Inf,Inf,NA,Inf,NA,NA,NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,NA,Inf,NA,Inf,Inf,NA),
tau3=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,NA),
tau4=c(NA,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,NA,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,NA,Inf,Inf,Inf,Inf,Inf),
d=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)),
row.names=c(NA,34L),class="data.frame")
parms=females
if(input$Sex==1) parms=males
max.stages=parms[,1]
tau=parms[,2:5]
d=parms[,6]
LogLike=function(Age,denom=1,dens=F){
LL=0
for(i in 1:34){
if(grades[i]>0){
if(grades[i]==1) LL=LL+log(1-1/(1+exp(-d[i]*(Age-tau[i,1]))))
else
if(grades[i]==max.stages[i])
LL=LL+log(1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1]))))
else {
p1 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]-1])))
p2 = 1/(1+exp(-d[i]*(Age-tau[i,grades[i]])))
LL = LL + log(p1-p2)
}
}
}
if(dens==F) return(LL/denom)
return(exp(LL)/denom)
}
model=maxNR(LogLike,start=input$Guess,dens=F)
mu=model$estimate
if(mu<0) mu=0
print(round(mu,2))
hess=as.numeric(-model$hessian)
se=0
if(mu==0) se=0
if(hess>0) se=sqrt(1/hess)
lower=mu-8*se
if(lower>30) lower=0
if(lower<0) lower=0
upper=mu+8*se
if(upper==0) upper=30
if(upper>20) upper=30
full.height=dnorm(mu,mu,se)
print('lower<upper?')
print(c(lower,upper))
find.right=function(right,denom,dens=T){
return(integrate(LogLike,lower=lower,upper=right,denom=denom,dens=T)$val-0.99)
}
denom=integrate(LogLike,lower=lower,upper=upper,denom=1,dens=T)$val
height.Like = LogLike(mu,denom=denom,dens=T)
height.Norm = height.Like
if(se>0) height.Norm = dnorm(mu,mu,se)
height=1.1*max(c(height.Like,height.Norm))
top=uniroot(find.right,lower=lower,upper=upper,denom=denom,dens=T)$root
top=top+2
Age=seq(-.2,top,.001)
N=length(Age)
ll=vector()
for(i in 1:N) ll[i]=LogLike(Age[i],denom=denom,dens=T)
if(se>0){
plot(Age,ll,type='l',ylim=c(0,height),xlim=c(mu-4*se,mu+4*se),xaxs='i',yaxs='i',ylab='Density')}
if(se==0){
plot(Age,ll,type='l',ylim=c(0,height),xlim=c(0,30),xaxs='i',yaxs='i',ylab='Density')
}
if(mu==0){
plot(Age,ll,type='l',ylim=c(0,height),xlim=c(0,3),xaxs='i',yaxs='i',ylab='Density')
}
lines(Age,dnorm(Age,mu,se),lty=2,lwd=2)
legend(x='topright',bty='n',legend=c('Normal','Posterior'),lty=c(2,1),lwd=c(2,1))
find.height=function(Age,try){
LogLike(Age,denom=denom,dens=T)-try
}
find.hpd=function(try){
left=uniroot(find.height,interval=c(-10,mu),try)$root
right=uniroot(find.height,interval=c(mu,30),try)$root
area=integrate(LogLike,lower=left,upper=right,denom=denom,dens=T)$val-input$HPD/100
return(area)
}
if(input$find.HPD==1) {
print(find.hpd(0.00001*height.Like))
print(find.hpd(height.Like))
hpd.height=uniroot(find.hpd,interval=c(0.00001*height.Like,height.Like))$root
left=uniroot(find.height,interval=c(lower,mu),try=hpd.height)$root
right=uniroot(find.height,interval=c(mu,upper),try=hpd.height)$root
L.area=integrate(LogLike,lower=lower,upper=left,denom=denom,dens=T)$val
R.area=integrate(LogLike,lower=right,upper=upper,denom=denom,dens=T)$val
left.ages=seq(0,left,0.001)
y=LogLike(left.ages,denom=denom,dens=T)
polygon(c(0,left.ages,left),c(0,y,0),col='gray')
right.ages=seq(right,30,0.001)
y=LogLike(right.ages,denom=denom,dens=T)
polygon(c(right,right.ages,30),c(0,y,0),col='gray')
print(round(c(left,right),2))
print(c(L.area,R.area,L.area+R.area))
}
abline(v=input$Guess)
prior=function(x) dnorm(x,input$prior.mu,input$prior.sd)
BF=function(x) prior(x)*LogLike(x,denom=denom,dens=T)
BF.under.thresh=integrate(BF,lower=0,upper=input$prior.mu)$val
BF.over.thresh=integrate(BF,lower=input$prior.mu,upper=30)$val
BF.val=BF.over.thresh/BF.under.thresh
output$scores<-renderText({
if(mu>30) mu=30
print(paste('Age estimate = ',round(mu,2),' se = ',round(se,2)))
})
if(input$find.HPD==1){
output$scores3<-renderText({
print(paste(input$HPD,'% HPD = ',round(left,2),' to ',round(right,2)))
})
}
output$scores4<-renderText({
print(paste("Bayes' Factor (>threshold/<threshold) = ",round(BF.val,2)))
})
}
)
}
)
# Run the application
shinyApp(ui = ui, server = server)