############ ### DATA ### ############ MincerCrossdata <- matrix(c( 96, 34, 13, 8, 11, 54, 50, 45, 14, 6, 11, 86, 83, 38, 10, 9, 18, 74, 130, 58, 14, 27, 11, 87, 375), 5, 5, byrow=T) #NOTE: Reading in transposed matrix: MincerPercentiles <- matrix(c( 14.19, 15.14, 15.15, 15.87, 18.27, 14.40, 15.66, 16.13, 16.89, 19.47, 15.02, 16.74, 16.98, 17.91, 20.30, 16.52, 17.90, 18.19, 19.55, 22.00, 16.91, 21.42, 20.58, 20.79, 23.28), 5, 5) # # For UK counts, the lower (mandibular) left wisdom tooth stage # is used: LL8xm in the terminology of the database, where # x ranges from D to H: UKcounts <- c(164, 148, 145, 80, 54) UKPercentiles <- matrix(c( 9.20, 9.70, 10.54, 11.07, 12.38, 13.61, 15.15, 16.00, 16.45, 17.02, 18.71, 12.03, 12.84, 13.12, 13.54, 14.28, 15.25, 16.09, 16.68, 17.05, 17.39, 18.20, 12.54, 14.24, 14.55, 14.99, 15.69, 16.44, 17.20, 18.19, 18.67, 18.87, 19.36, 14.14, 15.32, 15.45, 15.89, 16.88, 17.74, 18.36, 19.27, 20.32, 21.10, 21.64, 15.47, 16.10, 16.66, 17.03, 18.10, 19.65, 20.91, 21.38, 21.50, 21.55, 21.63), 5, 11, byrow=T) haglundMature <- c(5, 15, 65, 215, 439, 656, 774, 823, 719, 553, 345) haglundTotal <- c(1069, 1427, 1222, 1312, 1359, 1285, 1215, 1070, 858, 634, 381) OttowCountdata <- c(84, 27, 1, 32, 182) OttowPercentiles <- matrix(c( 12.05, 13.05, 13.84, 15.02, 19.15, 13.68, 15.18, 15.94, 16.98, 17.88, 17.77, 17.77, 17.77, 17.77, 17.77, 16.13, 16.79, 17.28, 17.89, 24.84, 17.46, 19.85, 21.56, 23.15, 24.98), 5, 5, byrow=T) tableDataLucas <- c(0, 0, 0, 4, 25, 27, 39, 44, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50) blekingeData <- matrix(c(0, 0, 1, 2, 12, 19, 31, 31, 22, 26, 29, 24, 23, 24, 34, 33), 2, 8, byrow=T) getData <- function () { #returns the table of actual #observations # rows: teeth # columns: knees, # 1: mature # 2: not mature # 3: "cannot tell" matrix(c(4176, 1735, 1364, 348, 1087, 237, 187, 83, 63), 3, 3) } getDataFemales <- function () { #returns the table of actual #observations # rows: teeth # columns: knees, # 1: mature # 2: not mature # 3: "cannot tell" matrix(c(190, 79, 40, 1, 12, 1, 6, 5, 3), 3, 3) } # Global parameter alpha = 3 #alpha = 10000 ######################### ## Utility functions #### ######################### rDirichlet <- function (alpha) { # Assuming alpha is the parameter vector y <- rgamma(length(alpha), alpha, 1) y/sum(y) } makePriorAge <- function() { nAge <- 100 # THE FOLLOWING IS THE STANDARD PRIOR: mmax <- pgamma(15, 4, 1) steps <- mmax *(1:nAge - 0.5)/nAge agevals <- qgamma(steps, 4, 1) + 15 # THE FOLLOWING IS THE FIRST ALTERNATIVE PRIOR: # mmin <- pnorm(15, 22.5, 4) # mmax <- pnorm(30, 22.5, 4) # steps <- mmin + (mmax - mmin)*(1:nAge - 0.5)/nAge # agevals <- qnorm(steps, 22.5, 4) # THE FOLLOWING IS THE SECOND ALTERNATIVE PRIOR: # agevals <- seq(15, 25, length.out = 100) cbind(agevals, 1/nAge) } getPriorValue <- function(x, prior) { sum(dnorm(x[1:2], prior[1:2], prior[3:4], log=TRUE)) } getProbs <- function (param, ages) { result <- matrix(0, 3, length(ages)) result[3,] <- param[3] + (ages-20)*param[4] result[1,] <- pnorm((ages - param[1])/param[2]) result[2,] <- (1-result[1,])*(1-result[3,]) result[1,] <- result[1,]*(1-result[3,]) result } plotACF <- function(result) { N <- dim(result$teeth)[1] teeth <- result$teeth[(1:(N/100)*100),1:2] knees <- result$knees[(1:(N/100)*100),1:2] par(mfrow=c(3,3)) acf(teeth[,1]) acf(teeth[,2]) acf(knees[,1]) acf(knees[,2]) acf(result$popul[,20]) acf(result$popul[,40]) acf(result$popul[,50]) acf(result$popul[,60]) acf(result$popul[,80]) } ################################################## ## Functions for computations in paper #### ################################################## ParameterEstimationFromPapers <- function(UKc = UKcounts, UCP = UKPercentiles, MincerCd = MincerCrossdata, MincerPerc = MincerPercentiles, OttowCountd = OttowCountdata, OttowPerc = OttowPercentiles, tableDataL = tableDataLucas, bleDat = blekingeData) { resultparams <- matrix(0, 7, 2) likelihood <- function (theta, data) { # This LOGlikelihoodfunction # is for numerical optimization x <- pnorm(data[,1], theta[1], exp(theta[2])) x[data[,2]<5] <- 1 - x[data[,2]<5] -sum(log(x)) } # # LUCAS: # dataLucas <- matrix(1, 1000, 2) for (i in 1:20) { dataLucas[((i-1)*50+1):((i-1)*50+50),1] <- 16 + (i-1)*0.5 + (1:50 - 0.5)/100 ind <- sample(((i-1)*50+1):((i-1)*50+50), tableDataL[i]) dataLucas[ind, 2] <- 5 } res <- nlm(likelihood, c(19, log(1)), dataLucas, stepmax = 2) resultparams[1,1] <- round(res$estimate[1],1) resultparams[1,2] <- round(exp(res$estimate[2]),1) # # Haglund # dataHaglund <- matrix(1, 11832, 2) N <- 11832 start <- cumsum(c(1, haglundTotal)) for (i in 1:11) { dataHaglund[start[i]:(start[i+1]-1), 1] <- 14 + i + ((1:haglundTotal[i] - 0.5)/haglundTotal[i]) ind <- sample(start[i]:(start[i+1]-1), haglundMature[i]) dataHaglund[ind, 2] <- 5 } res <- nlm(likelihood, c(19, log(1)), dataHaglund, stepmax = 2) resultparams[2,1] <- round(res$estimate[1],1) resultparams[2,2] <- round(exp(res$estimate[2]),1) # x <- seq(15, 26, 0.1) # y <- pnorm(x, res$estimate[1], exp(res$estimate[2])) # plot(x,y, type="l") # points(15:25 + 0.5, haglundMature/haglundTotal) # # UK database DARLInG # N <- sum(UKc) start <- cumsum(c(1, UKc)) result <- matrix(0, N, 2) perp <- c(0, 0.005, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.995, 1) for (i in 1:5) { result[start[i]:(start[i+1]-1),2] <- i xout <- qnorm(((1:UKc[i])-0.5)/UKc[i]) pp <- quantile(xout, probs = perp) ppy <- UCP[i,] result[start[i]:(start[i+1]-1),1] <- approx(x=pp, y=ppy, xout = xout)$y } res <- nlm(likelihood, c(19, log(1)), result, stepmax = 2) resultparams[3,1] <- round(res$estimate[1],1) resultparams[3,2] <- round(exp(res$estimate[2]),1) # # Construct Mincer Rawdata # N <- 271 grsize <- apply(MincerCd, 2, sum) grsize <- grsize/sum(grsize)*0.999 gr <- round(N*grsize, 0) start <- cumsum(c(1, gr)) result <- matrix(0, N, 2) perp <- c(0.1, 0.25, 0.5, 0.75, 0.9) linExtend <- function (x0, x1, x2, y1, y2) { #A linear extrapolation (x0-x1)/(x2-x1)*(y2-y1) + y1 } for (i in 1:5) { result[start[i]:(start[i+1]-1),2] <- i xout <- qnorm(((1:gr[i])-0.5)/gr[i]) pp <- quantile(xout, probs=c(0, perp, 1)) ymin <- linExtend(0, perp[1], perp[2], MincerPercentiles[i,1], MincerPerc[i,2]) ymax <- linExtend(1, perp[4], perp[5], MincerPercentiles[i,4], MincerPerc[i,5]) ymin <- max(ymin, 14.1) ymax <- min(ymax, 24.5) ppy <- c(ymin, MincerPercentiles[i,], ymax) result[start[i]:(start[i+1]-1),1] <- approx(x=pp, y=ppy, xout = xout)$y } res <- nlm(likelihood, c(19, log(1)), result, stepmax = 2) resultparams[4,1] <- round(res$estimate[1],1) resultparams[4,2] <- round(exp(res$estimate[2]),1) # # OTTOW: # gr <- OttowCountd start <- cumsum(c(1, gr)) N <- sum(gr) result <- matrix(0, N, 2) perp <- c(0, 0.25, 0.5, 0.75, 1) for (i in 1:5) { result[start[i]:(start[i+1]-1),2] <- i xout <- qnorm(((1:gr[i])-0.5)/gr[i]) pp <- quantile(xout, probs = perp) ppy <- OttowPerc[i,] if (gr[i]>1) result[start[i]:(start[i+1]-1),1] <- approx(x=pp, y=ppy, xout = xout)$y else result[start[i]] <- ppy[3] } res <- nlm(likelihood, c(19, log(1)), result, stepmax = 2) resultparams[5,1] <- round(res$estimate[1],1) resultparams[5,2] <- round(exp(res$estimate[2]),1) # # Adjusted Ottow results: # # There are 32 persons assigned category 4. # In this run, we re-assign 16 of them to category 5. ind <- sample(start[4]:(start[5]-1), 16) result[ind,2] <- 5 res <- nlm(likelihood, c(19, log(1)), result, stepmax = 2) resultparams[6,1] <- round(res$estimate[1],1) resultparams[6,2] <- round(exp(res$estimate[2]),1) # # BLEKINGE # dataBlekinge <- matrix(1, sum(bleDat[2,]), 2) start <- cumsum(c(1, bleDat[2,])) for (i in 1:8) { dataBlekinge[start[i]:(start[i+1]-1),1] <- 13 + i + (1:bleDat[2,i] - 0.5)/(2*bleDat[2,i]) if (i>2) { ind <- sample(start[i]:(start[i+1]-1), bleDat[1,i]) dataBlekinge[ind,2] <- 5 } } res <- nlm(likelihood, c(19, log(1)), dataBlekinge, stepmax = 2) resultparams[7,1] <- round(res$estimate[1],1) resultparams[7,2] <- round(exp(res$estimate[2]),1) resultparams } MakeFigure1 <- function(prior1, prior2) { age <- makePriorAge() x <- age[,1] y <- pnorm((x - prior1[1])/prior1[2]) plot(x, y, xlab="Age", ylab="Probability", type="l", main="Prior age indicator models", lwd=2) for (i in 1:5) { theta1 <- rnorm(1, prior1[1], prior1[3]) theta2 <- -1 while (theta2<=0) theta2 <- rnorm(1, prior1[2], prior1[4]) y <- pnorm((x - theta1)/theta2) lines(x, y, lty=2) } y <- pnorm((x - prior2[1])/prior2[2]) lines(x, y, lwd = 2) for (i in 1:5) { theta1 <- rnorm(1, prior2[1], prior2[3]) theta2 <- -1 while (theta2<=0) theta2 <- rnorm(1, prior2[2], prior2[4]) y <- pnorm((x - theta1)/theta2) lines(x, y, lty=3) } } MakeFigure2 <- function() { age <- makePriorAge() nAge <- dim(age)[1] # Plot prior: plot(age[,1], qbeta(0.5, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), type="l", lty=1, xlab="Age", ylab="Probability", main="Population prior") lines(age[,1], qbeta(0.25, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=2) lines(age[,1], qbeta(0.75, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=2) lines(age[,1], qbeta(0.025, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=3) lines(age[,1], qbeta(0.975, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=3) } GibbsSampler <- function(priorT, priorK) { # alpha <- 3 # Set up the age structure; draw an initial age distribution, and # compute (as age[,2]) the probabilities at the discretized ages # using this age profile: age <- makePriorAge() Nage <- dim(age)[1] #set up starting variables for Gibbs sampling teeth <- c(19, 1.5, 0.2, 0) # was c(19, 1.5, 0.2, 0) knees <- c(18, 1.5, 0.05, 0) # was c(21, 2, 0.05, 0) counts <- array(0, c(3, 3, Nage)) data <- getData() # get value for prior density at this parameter dpriorT <- getPriorValue(teeth, priorT) # get value for prior density at this parameter dpriorK <- getPriorValue(knees, priorK) probsT <- getProbs(teeth, age[,1]) # returns 3xnAge matrix probsK <- getProbs(knees, age[,1]) # returns 3xnAge matrix #Gibbs sampling maxstep = 1000000 result <- list(teeth = matrix(0, maxstep, 4), knees = matrix(0, maxstep, 4), popul = matrix(0, as.integer(maxstep/100), Nage), cumprobs = matrix(0, 4, Nage), tableprobs = array(0, c(maxstep, 3, 3, 2))) result$cumprobs[1:2,] <- probsT[2:3,] result$cumprobs[3:4,] <- probsK[2:3,] counter <- 1 countT <- 0 countK <- 0 for (step in 1:maxstep) { # simulate counts given obs model and population model for (i in 1:3) for (j in 1:3) counts[i,j,] <- rmultinom( 1, data[i,j], prob = probsT[i,]*probsK[j,]*age[,2]) # monitoring convergence: if (step/1000 == round(step/1000)) { mytab <- matrix(0, 3, 3) for (i in 1:3) for (j in 1:3) mytab[i,j] <- sum(probsT[i,]*probsK[j,]*age[,2]) mytab <- mytab/sum(mytab)*sum(data) print(round(mytab)) print(paste("====", step)) } # Simulate a new teeth model: # 0.05, 0.05, 0.01, 0.001 newteeth <- teeth + rnorm(4, 0, c(0.05, 0.07, 0.002, 0.0002)) newProbsT <- getProbs(newteeth, age[,1]) if (all(newProbsT > 0)) { # reject directly nonsensical proposals newpriorT <- getPriorValue(newteeth, priorT) A <- apply(counts, c(1,3), sum) B <- sum((log(newProbsT) - log(probsT)) * A, na.rm=T) if (runif(1) < exp(B + newpriorT - dpriorT)) { countT <- countT + 1 teeth <- newteeth probsT <- newProbsT dpriorT <- newpriorT } } # Simulate a new knee model: # 0.05, 0.05, 0.01, 0.001 newknees <- knees + rnorm(4, 0, c(0.05, 0.05, 0.002, 0.0002)) newProbsK <- getProbs(newknees, age[,1]) if (all(newProbsK > 0)) { # reject directly nonsensical proposals newpriorK <- getPriorValue(newknees, priorK) A <- apply(counts, c(2,3), sum) B <- sum((log(newProbsK) - log(probsK)) * A, na.rm=T) if (runif(1) < exp(B + newpriorK - dpriorK)) { countK <- countK + 1 knees <- newknees probsK <- newProbsK dpriorK <- newpriorK } } # simulate population given counts ct <- apply(counts, 3, sum) age[,2] <- rDirichlet(ct + alpha/Nage) # Collect results result$teeth[step,] <- teeth result$knees[step,] <- knees result$cumprobs <- result$cumprobs + rbind(probsT[-1,], probsK[-1,]) result$tableprobs[step,,,1] <- apply(counts[,,age[,1]<=18], 1:2, sum) result$tableprobs[step,,,2] <- apply(counts[,,age[,1]> 18], 1:2, sum) A <- step/100 if (A == as.integer(A)) { result$popul[counter,] <- cumsum(age[,2]) counter <- counter + 1 } } print(paste("Accep. rate teeth", countT/maxstep)) print(paste("Accep. rate knees", countK/maxstep)) result } MakeTables3and4 <- function(result) { tp <- result$tableprobs N <- dim(tp)[1] res1 <- array(0, c(N, 3, 4)) res1[,1,1] <- tp[,1,1,2] + tp[,1,2,2] + tp[,1,3,2] + tp[,2,1,2] + tp[,3,1,2] res1[,1,2] <- tp[,2,2,2] + tp[,2,3,2] + tp[,3,2,2] res1[,1,3] <- tp[,3,3,2] res1[,1,4] <- res1[,1,1] + res1[,1,2] + res1[,1,3] res1[,2,1] <- tp[,1,1,1] + tp[,1,2,1] + tp[,1,3,1] + tp[,2,1,1] + tp[,3,1,1] res1[,2,2] <- tp[,2,2,1] + tp[,2,3,1] + tp[,3,2,1] res1[,2,3] <- tp[,3,3,1] res1[,2,4] <- res1[,2,1] + res1[,2,2] + res1[,2,3] res1[,3,1] <- res1[,1,1] + res1[,2,1] res1[,3,2] <- res1[,1,2] + res1[,2,2] res1[,3,3] <- res1[,1,3] + res1[,2,3] res1[,3,4] <- res1[,1,4] + res1[,2,4] print(round(apply(res1, 2:3, mean))) print(round(apply(res1, 2:3, quantile, probs = 0.025))) print(round(apply(res1, 2:3, quantile, probs = 0.975))) percentchil <- res1[,2,4]/res1[,3,4] percentchil <- percentchil[!is.na(percentchil)] sensitivity <- res1[,1,1]/(res1[,1,1] + res1[,1,2]) sensitivity <- sensitivity[!is.na(sensitivity)] specificity <- res1[,2,2]/(res1[,2,2] + res1[,2,1]) specificity <- specificity[!is.na(specificity)] PPV <- res1[,1,1]/(res1[,1,1] + res1[,2,1]) PPV <- PPV[!is.na(PPV)] NPV <- res1[,2,2]/(res1[,2,2] + res1[,1,2]) NPV <- NPV[!is.na(NPV)] print("Proportion children:") print(round(mean(percentchil)*100)) print(round(quantile(percentchil, probs=0.025)*100)) print(round(quantile(percentchil, probs=0.975)*100)) print("Sensitivity:") print(round(mean(sensitivity)*100)) print(round(quantile(sensitivity, probs=0.025)*100)) print(round(quantile(sensitivity, probs=0.975)*100)) print("Specificity:") print(round(mean(specificity)*100)) print(round(quantile(specificity, probs=0.025)*100)) print(round(quantile(specificity, probs=0.975)*100)) print("Positive predictive value:") print(round(mean(PPV)*100)) print(round(quantile(PPV, probs=0.025)*100)) print(round(quantile(PPV, probs=0.975)*100)) print("Negative predictive value:") print(round(mean(NPV)*100)) print(round(quantile(NPV, probs=0.025)*100)) print(round(quantile(NPV, probs=0.975)*100)) tpFULL <- array(0, c(N, 4, 4, 2)) tpFULL[,-4,-4,] <- tp tpFULL[, 4,-4,] <- apply(tp,c(1,3,4), sum) tpFULL[, , 4,] <- apply(tpFULL, c(1,2,4), sum) propc <- tpFULL[,,,1]/(tpFULL[,,,1] + tpFULL[,,,2]) print("Proportion of children:") print(round(apply(propc, 2:3, mean)*100)) print(round(apply(propc, 2:3, quantile, probs = 0.025)*100)) print(round(apply(propc, 2:3, quantile, probs = 0.975)*100)) } MakeFigure3 <- function(result) { print("Mean data for teeth parameters") print(round(apply(result$teeth, 2, mean), 2)) print("Sd data for teeth parameters") print(round(apply(result$teeth, 2, sd), 2)) print("Mean data for knees parameters") print(round(apply(result$knees, 2, mean), 2)) print("Sd data for knees parameters") print(round(apply(result$knees, 2, sd), 2)) # diff <- result$teeth[,1] - result$knees[,1] print("Mean for difference teeth maturity and knee maturity") print(mean(diff)) print(quantile(diff, probs = c(0.025, 0.25, 0.75 ,0.975))) # # age <- makePriorAge() x <- age[,1] t11 <- mean(result$teeth[,1]) t12 <- mean(result$teeth[,2]) y <- pnorm((x - t11)/t12) plot(x, y, xlab="Age", ylab="Probability", type="l", main="Posterior age indicator models", lwd=2) N <- dim(result$teeth)[1] ind <- sample(N, 5) for (i in 1:5) { t11 <- result$teeth[ind[i],1] t12 <- result$teeth[ind[i],2] y <- pnorm((x - t11)/t12) lines(x, y, lty=2) } t21 <- mean(result$knees[,1]) t22 <- mean(result$knees[,2]) y <- pnorm((x - t21)/t22) lines(x, y, lwd = 2) for (i in 1:5) { t21 <- result$knees[ind[i], 1] t22 <- result$knees[ind[i], 2] y <- pnorm((x - t21)/t22) lines(x, y, lty=3) } } MakeFigure4 <- function(result) { pop <- result$popul[-(1:200),] age <- makePriorAge() x <- age[,1] y <- apply(pop, 2, median) plot(x, y, type="l", xlab="Age", ylab="Probaility", main = "Population posterior") y <- apply(pop, 2, quantile, probs = 0.25, names=FALSE) lines(x, y, lty=2) y <- apply(pop, 2, quantile, probs = 0.75, names=FALSE) lines(x, y, lty=2) y <- apply(pop, 2, quantile, probs = 0.025, names=FALSE) lines(x, y, lty=3) y <- apply(pop, 2, quantile, probs = 0.975, names=FALSE) lines(x, y, lty=3) } MakeFigure5 <- function(diffMM, diffLL, diffLM, diffLH, diffML, diffMH, diffHL, diffHM, diffHH) { plot(density(diffMM), xlab="tooth maturation age minus knee maturation age", main="Posterior densities", lwd = 2) lines(density(diffLL)) lines(density(diffLM)) lines(density(diffLH)) lines(density(diffML)) lines(density(diffMH)) lines(density(diffHL)) lines(density(diffHM)) lines(density(diffHH)) } MakeTable1Suppl <- function(teeth0, knees0) { age <- makePriorAge() data <- getData() loglik <- function(theta) { # computes the loglikelihood of data given the remaining parameters teeth <- c(teeth0, theta[1:2]) knees <- c(knees0, theta[3:4]) probsT <- getProbs(teeth, age[,1]) probsK <- getProbs(knees, age[,1]) result <- 0 for (i in 1:3) for (j in 1:3) result <- result - data[i,j]* log(sum(probsT[i,]*probsK[j,]*age[,2])) result } res <- nlm(loglik,c(0.2, 0, 0.05, 0), stepmax = 0.01)$estimate print("Maximum Likelihood estimates for missing parameters:") print(res) teeth <- c(teeth0, res[1:2]) knees <- c(knees0, res[3:4]) probsT <- getProbs(teeth, age[,1]) probsK <- getProbs(knees, age[,1]) myprobs <- matrix(0, 3, 3) for (i in 1:3) for (j in 1:3) myprobs[i,j] <- sum(probsT[i,]*probsK[j,]*age[,2]) print(chisq.test(as.vector(data), p = as.vector(myprobs))) predtable <- myprobs*sum(data) predtable <- round(predtable) print("Data:") print(data) print("Predicted data:") print(predtable) print(apply(predtable, 1, sum)) print(apply(predtable, 2, sum)) print(sum(predtable)) pvals <- matrix(0, 3, 3) for (i in 1:3) for (j in 1:3) pvals[i,j] <- binom.test(data[i,j], sum(data), myprobs[i,j])$p.value print(pvals) } #################################################################### ### RUNNING COMPUTATIONS ### #################################################################### #ParameterEstimationFromPapers() # Set the priors based on estimated parameters TMAIN <- c(19.5, 1.6, 0.8, 0.4) THIGH <- c(20.5, 1.6, 0.8, 0.4) TLOW <- c(18.5, 1.6, 0.8, 0.4) THAGLUND <- c(20.9, 2.5, 0.8, 0.4) KMAIN <- c(18.5, 1.4, 0.8, 0.4) KHIGH <- c(19.5, 1.4, 0.8, 0.4) KLOW <- c(17.5, 1.4, 0.8, 0.4) prior <- list(TMAIN = TMAIN, KMAIN = KMAIN, THIGH = THIGH, TLOW = TLOW, KHIGH = KHIGH, KLOW = KLOW) #MakeFigure1(prior$TMAIN, prior$KMAIN) #MakeFigure2() #print(system.time(resMAIN <- GibbsSampler(prior$TMAIN, prior$KMAIN))) #print(system.time(resLL <- GibbsSampler(prior$TLOW, prior$KLOW))) #print(system.time(resLM <- GibbsSampler(prior$TLOW, prior$KMAIN))) #print(system.time(resLH <- GibbsSampler(prior$TLOW, prior$KHIGH))) #print(system.time(resML <- GibbsSampler(prior$TMAIN, prior$KLOW))) #print(system.time(resMH <- GibbsSampler(prior$TMAIN, prior$KHIGH))) #print(system.time(resHL <- GibbsSampler(prior$THIGH, prior$KLOW))) #print(system.time(resHM <- GibbsSampler(prior$THIGH, prior$KMAIN))) #print(system.time(resHH <- GibbsSampler(prior$THIGH, prior$KHIGH))) #MakeTables3and4(resMAIN) #MakeFigure3(resMAIN) #MakeFigure4(resMAIN) #MakeFigure5(resMAIN$teeth[,1] - resMAIN$knees[,1], # resLL$teeth[,1] - resLL$knees[,1], # resLM$teeth[,1] - resLM$knees[,1], # resLH$teeth[,1] - resLH$knees[,1], # resML$teeth[,1] - resML$knees[,1], # resMH$teeth[,1] - resMH$knees[,1], # resHL$teeth[,1] - resHL$knees[,1], # resHM$teeth[,1] - resHM$knees[,1], # resHH$teeth[,1] - resHH$knees[,1]) #MakeTable1Suppl(TMAIN[1:2], KMAIN[1:2]) #print(system.time(resHAG <- GibbsSampler(THAGLUND, prior$KMAIN)))