source("MAMMA.R") ################################################ # Cholesterol Data # ################################################ # reading data chol <- read.table(file = "cholesterol.csv", header = TRUE, sep = ",") # > str(chol) # 'data.frame': 2000 obs. of 4 variables: # $ subject: int 1 1 1 1 1 1 1 1 1 1 ... # $ method : Factor w/ 2 levels "cobasb","echem": 2 2 2 2 2 2 2 2 2 2 ... # $ rep : int 1 2 3 4 5 6 7 8 9 10 ... # $ obs : int 212 212 209 211 212 210 206 209 210 212 ... # > chol$subject <- as.factor(chol$subject) colnames(chol) <- c("subject", "method", "rep", "meas") # n = 100, 2 methods (cobasb and echem), 10 replicates on each, total of 2000 obs. # cobasb = method1, echem = method2 ####################################################################### ################################################ # Displaying data # ################################################ # trellis plot plot(groupedData(meas ~ method | subject, data = chol, order.groups = TRUE), ylim = c(-1, 102), xlab = "cholesterol level (mg/dL)", ylab = "sorted subject ID") # scatter plot and B-A plot of randomly formed pairs and also of averages of replicates set.seed(12345) chol.paired <- unlinked2paired(chol) y1 <- chol.paired[, 1] y2 <- chol.paired[, 2] par(mfrow = c(2, 2), pty = "m", oma = c(0.5, 0, 0, 0), mar = c(4.5, 4.5, 2.5, 0.5) + 0.1, cex.lab = 0.8, cex.main = 0.8, cex.axis = 0.8) # randomly formed pairs # Panel a plot(y1, y2, xlab = "cholesterol, cobasb", ylab = "cholesterol, echem", xlim = c(45, 375), ylim = c(45, 375), panel.first = grid(lty = "solid")) abline(a = 0, b = 1) title(main = "(a)") # Panel b avg <- (y1 + y2)/2 diff <- y2 - y1 plot(avg, diff, xlab = "pair average", ylab = "echem - cobasb", ylim = c(-max(abs(diff)), max(abs(diff))), xlim = c(45, 375), panel.first = grid(lty = "solid")) abline(h = 0) title(main = "(b)") # plots of averages of replicates # Panel c chol.avg <- unlinked2avg(chol) # > head(chol.avg) # cobasb echem # 1 202.0 210.3 # 2 228.9 238.5 # 3 193.9 208.4 # 4 291.4 308.9 # 5 216.0 214.4 # 6 152.6 166.6 # > yy1 <- chol.avg[, "cobasb"] yy2 <- chol.avg[, "echem"] plot(yy1, yy2, xlim = c(45, 375), ylim = c(45, 375), xlab = "average cobasb", ylab = "average echem", panel.first = grid(lty = "solid")) abline(a = 0, b = 1) title(main = "(c)") # Panel d avg.avg <- (yy1 + yy2)/2 diff.avg <- yy2 - yy1 plot(avg.avg, diff.avg, ylim = c(-max(abs(diff.avg)), max(abs(diff.avg))), xlim = c(45, 375), xlab = "subject average", ylab = "avg echem - avg cobasb", panel.first = grid(lty = "solid")) abline(h = 0) title(main = "(d)") # interaction plots par(mfrow = c(1, 1)) with(chol, { interaction.plot(method, subject, meas, legend = FALSE, type = "o", lty = 1, pch = 19, cex = 0.5, lwd = 0.8, xlab = "method", ylab = "average cholesterol", las = 1) }) ################################################################## #################################################### # Analyzing residuals from the homoscedastic fit # ################################################### # > levels(chol$method) # [1] "cobasb" "echem" # > levels(chol$method) <- c("method1", "method2") # > levels(chol$method) # [1] "method1" "method2" # > homo.fit.results <- lme.unlinked.repmeas.fit(chol) # > homo.fit.results # $fit # Linear mixed-effects model fit by maximum likelihood # Data: dat # Log-likelihood: -5117.55 # Fixed: meas ~ method - 1 # methodmethod1 methodmethod2 # 184.380 189.978 # # Random effects: # Composite Structure: Blocked # # Block 1: (Intercept) # Formula: ~1 | subject # (Intercept) # StdDev: 65.76346 # # Block 2: methodmethod1, methodmethod2 # Formula: ~method - 1 | subject # Structure: Multiple of an Identity # methodmethod1 methodmethod2 Residual # StdDev: 5.185825 5.185825 2.498162 # # Variance function: # Structure: Different standard deviations per stratum # Formula: ~1 | method # Parameter estimates: # method2 method1 # 1.0000000 0.7852925 # Number of Observations: 2000 # Number of Groups: 100 # # $param.hat # alpha mu.b lsigmasq.b lsigmasq.bI lsigmasq.e1 lsigmasq.e2 # 5.598000 184.380000 8.372129 3.291858 1.347712 1.831110 # # $var.hat # sigmasq.b sigmasq.bI sigmasq.e1 sigmasq.e2 sigmasq # 4324.832354 26.892786 3.848611 6.240813 63.874996 # # $bvn.hat # mean1 mean2 var1 var2 corr # 184.3800000 189.9780000 4355.5737512 4357.9659527 0.9926695 # # $reliability.hat # method1 method2 # 0.9991164 0.9985680 # # > # the standardized residuals from the homoscedastic fit homo.res <- resid(homo.fit.results$fit, type = "response") homo.stdres <- resid(homo.fit.results$fit, type = "pear") # > levels(chol$method) # [1] "method1" "method2" # > orig.method <- chol$method levels(orig.method) <- c("cobasb", "echem") # > levels(orig.method) # [1] "cobasb" "echem" # > # residual plots from homoscedastic fit # standardized residuals against fitted values p1 <- xyplot(homo.stdres ~ fitted(homo.fit.results$fit) | orig.method, grid = T, data = chol, xlab = "fitted values", ylab = "standardized residuals", panel = function(...) { panel.abline(h = 0, lty = 1, col = "black") panel.xyplot(...) }) # absolute residuals against fitted values p2 <- xyplot(abs(homo.stdres) ~ fitted(homo.fit.results$fit) | orig.method, grid = T, data = chol, xlab = "fitted values", ylab = "| standardized residuals |") print(p1, split = c(1, 1, 1, 2), more = TRUE) # ?print.trellis print(p2, split = c(1, 2, 1, 2)) # variance covariate v <- get.var.covariate(chol, "avg.mean") # plot of absolute residuals against v xyplot(abs(homo.res) ~ v | orig.method, grid = T, data = chol, xlab = "variance covariate", ylab = "| standardized residuals |", type = c("p", "smooth", "r")) # additional plots of absolute residuals # power model - log(absolute residuals) against log(v) p1 <- xyplot(log(abs(homo.res)) ~ log(v) | orig.method, grid = T, data = chol, xlab = "log (variance covariate)", ylab = "log (abs residuals)", type = c("p", "r")) # exponential model - log(absolute residuals) against v p2 <- xyplot(log(abs(homo.res)) ~ v | orig.method, grid = T, data = chol, xlab = "variance covariate", ylab = "log (abs residuals)", type = c("p", "r")) print(p1, split = c(1, 1, 1, 2), more = TRUE) # ?print.trellis print(p2, split = c(1, 2, 1, 2)) #################################################### # Heteroscedastic fits # ################################################### # the variance covariate that we used before for exploratory analysis # v <- get.var.covariate(chol, "avg.mean") # set up a grid of variance covariate for confidence intervals, etc vgrid <- seq(min(chol$meas), max(chol$meas), l = 50) # the power variance function model fit.results.power <- lme.repmeas.hetero.fit(chol, v, "power", vgrid) # > names(fit.results.power) # [1] "varfunc" "fit" "param.hat" # [4] "anova.results" "vgrid" "err.var.hat" # [7] "var.hat" "bvn.hat" "var.D.hat" # [10] "reliability.hat" "loa" # > fit.power <- fit.results.power$fit # > fit.power # Linear mixed-effects model fit by maximum likelihood # Data: dat # Log-likelihood: -4868.601 # Fixed: meas ~ method - 1 # methodmethod1 methodmethod2 # 184.3839 189.9657 # # Random effects: # Composite Structure: Blocked # # Block 1: (Intercept) # Formula: ~1 | subject # (Intercept) # StdDev: 65.76154 # # Block 2: methodmethod1, methodmethod2 # Formula: ~method - 1 | subject # Structure: Multiple of an Identity # methodmethod1 methodmethod2 Residual # StdDev: 5.153692 5.153692 0.01433887 # # Combination of variance functions: # Structure: Different standard deviations per stratum # Formula: ~1 | method # Parameter estimates: # method2 method1 # 1.0000000 0.6246136 # Structure: Power of variance covariate, different strata # Formula: ~vtilde | method # Parameter estimates: # method2 method1 # 0.9757916 1.0153789 # Number of Observations: 2000 # Number of Groups: 100 # > hetero.stdres <- resid(fit.power, type = "pear") # plot of residuals against fitted values xyplot(hetero.stdres ~ fitted(fit.power) | orig.method, grid = T, data = chol, xlab = "fitted values", ylab = "standardized residuals", panel = function(...) { panel.abline(h = 0, lty = 1, col = "black") panel.xyplot(...) }) qqnorm(fit.power, ~resid(., type = "pear") | method) qqnorm(fit.power, ~ranef(.) | method) # the exponential variance function model fit.results.exp <- lme.repmeas.hetero.fit(chol, v, "exponential", vgrid) fit.exp <- fit.results.exp$fit # > fit.exp # Linear mixed-effects model fit by maximum likelihood # Data: dat # Log-likelihood: -4875.654 # Fixed: meas ~ method - 1 # methodmethod1 methodmethod2 # 184.3874 189.9572 # # Random effects: # Composite Structure: Blocked # # Block 1: (Intercept) # Formula: ~1 | subject # (Intercept) # StdDev: 65.75696 # # Block 2: methodmethod1, methodmethod2 # Formula: ~method - 1 | subject # Structure: Multiple of an Identity # methodmethod1 methodmethod2 Residual # StdDev: 5.12307 5.12307 0.7205942 # # Combination of variance functions: # Structure: Different standard deviations per stratum # Formula: ~1 | method # Parameter estimates: # method2 method1 # 1.0000000 0.7478053 # Structure: Exponential of variance covariate, different strata # Formula: ~vtilde | method # Parameter estimates: # method2 method1 # 0.006009060 0.006107093 # Number of Observations: 2000 # Number of Groups: 100 # > plot(fit.exp, resid(., type = "pear") ~ fitted(.) | method, xlab = "fitted values", ylab = "standardized residuals") qqnorm(fit.exp, ~resid(., type = "pear") | method) qqnorm(fit.exp, ~ranef(.) | method) # > anova(fit.power, fit.exp) # Model df AIC BIC logLik # fit.power 1 8 9753.202 9798.009 -4868.601 # fit.exp 2 8 9767.307 9812.114 -4875.654 # > # compare the fitted sd from the two models to the observed sd's sd.mat <- with(chol, by(meas, list(subject, method), sd)) v.unique <- with(chol, by(v, list(subject, method), unique))[, 1] # plot the fitted sd's dd.fitted <- data.frame(vgrid = rep(vgrid, 4), model = rep(c(rep("power", length(vgrid)), rep("exponential", length(vgrid))), 2), method = c(rep("cobasb", 2 * length(vgrid)), rep("echem", 2 * length(vgrid))), fitted.sd = sqrt(c(fit.results.power$err.var.hat[, "method1"], fit.results.exp$err.var.hat[, "method1"], fit.results.power$err.var.hat[, "method2"], fit.results.exp$err.var.hat[, "method2"]))) p1 <- xyplot(fitted.sd ~ vgrid | method, groups = model, data = dd.fitted, type = "l", lty = c(2, 1), grid = T, col = "black", xlab = "cholesterol", ylab = "standard deviation", key = list(text = list(c("power", "exponential")), lines = list(lty = 1:2), columns = 2)) # overlay the observed sd's dd.obs <- data.frame(vobs = rep(v.unique, 2), method = c(rep("cobasb", length(v.unique)), rep("echem", length(v.unique))), obs.sd = c(sd.mat[, "method1"], sd.mat[, "method2"])) p2 <- xyplot(obs.sd ~ vobs | method, data = dd.obs, col = "black", grid = T) p1 + as.layer(p2) ########################################################################### ################################################ # Inference based on the power model # ################################################ # the variance covariate that we used before for exploratory analysis # v <- get.var.covariate(chol, "avg.mean") # set up a grid of variance covariate for confidence intervals, etc # vgrid <- seq(min(chol$meas), max(chol$meas), l=50) conf.results <- conf.measures.hetero.repmeas(chol, fit.results.power$param.hat, v, "power", vgrid) # > names(conf.results) # [1] "param" "lambda" "tdi" "ccc" "rep.tdi" "rep.ccc" # > # explore impact of heteroscedasticity # error SD of method 1 # > sqrt(range(fit.results.power$err.var.hat[,"method1"])) # [1] 0.427330 3.649231 # > # error SD of method 2 # > sqrt(range(fit.results.power$err.var.hat[,"method2"])) # [1] 0.5884447 4.6219899 # > # between-subject SD # > sqrt(4351.141) # [1] 65.96318 # > # SD of method 1 # # > sqrt(range(fit.results.power$bvn.hat$var1)) # [1] 65.96456 66.06404 # > # SD of method 2 # > sqrt(range(fit.results.power$bvn.hat$var2)) # [1] 65.96580 66.12491 # > # correlation # > range(fit.results.power$bvn.hat$corr) # [1] 0.9899511 0.9938353 # > # SD of D # > sqrt(range(fit.results.power$var.D.hat)) # [1] 7.324614 9.370206 # > par(mfrow = c(2, 2)) plot(vgrid, fit.results.power$bvn.hat$var1, type = "l") lines(vgrid, fit.results.power$bvn.hat$var2, lty = 2) plot(vgrid, fit.results.power$bvn.hat$corr, type = "l") plot(vgrid, fit.results.power$var.D.hat, type = "l") plot(vgrid, fit.results.power$reliability.hat[, 1], type = "l") lines(vgrid, fit.results.power$reliability.hat[, 2], lty = 2) # the limits of agreement par(mfrow = c(1, 1)) plot(vgrid, fit.results.power$loa$loa.inter[, 1], ylim = c(-25, 25), type = "l", xlab = "cholesterol", ylab = "difference", panel.first = grid(lty = "solid")) lines(vgrid, fit.results.power$loa$loa.inter[, 2], lty = 1) lines(vgrid, fit.results.power$loa$loa.intra.1[, 1], lty = 2) lines(vgrid, fit.results.power$loa$loa.intra.1[, 2], lty = 2) lines(vgrid, fit.results.power$loa$loa.intra.2[, 1], lty = 3) lines(vgrid, fit.results.power$loa$loa.intra.2[, 2], lty = 3) legend("bottomleft", legend = c("inter-method", "cobasb", "echem"), lty = 1:3, bty = "n") # > round(conf.results$param, 2) # estimate se lcl ucl # alpha 5.58 0.74 4.14 7.02 # mu.b 184.38 6.60 171.45 197.31 # lsigmasq.b 8.37 0.14 8.09 8.65 # lsigmasq.bI 3.28 0.14 3.00 3.56 # lsigmasq.e1 -9.43 0.57 -10.55 -8.31 # lsigmasq.e2 -8.49 0.59 -9.64 -7.34 # delta.1 1.02 0.06 0.91 1.12 # delta.2 0.98 0.06 0.86 1.09 # > # ranges # > round(conf.results$param, 2) # estimate se lcl ucl # alpha 5.58 0.74 4.14 7.02 # mu.b 184.38 6.60 171.45 197.31 # lsigmasq.b 8.37 0.14 8.09 8.65 # lsigmasq.bI 3.28 0.14 3.00 3.56 # lsigmasq.e1 -9.43 0.57 -10.55 -8.31 # lsigmasq.e2 -8.49 0.59 -9.64 -7.34 # delta.1 1.02 0.06 0.91 1.12 # delta.2 0.98 0.06 0.86 1.09 # > # > apply(conf.results$rep.ccc, MAR=2, range) # vgrid zccc.1.hat se.zccc.1.hat lcb.zccc.1 ccc.1.hat lcb.ccc.1 # [1,] 45 3.241918 0.07431636 3.101790 0.9969488 0.9959638 # [2,] 372 5.385879 0.10529889 5.212678 0.9999580 0.9999407 # zccc.2.hat se.zccc.2.hat lcb.zccc.2 ccc.2.hat lcb.ccc.2 # [1,] 3.006071 0.07430588 2.865003 0.9951143 0.9935269 # [2,] 5.065962 0.10686722 4.890181 0.9999204 0.9998869 # > # > apply(conf.results$rep.tdi, MAR=2, range) # vgrid ltdi.1.hat se.ltdi.1.hat ucb.ltdi.1 tdi.1.hat ucb.tdi.1 # [1,] 45 -0.005973723 0.02358652 0.1226753 0.9940441 1.130517 # [2,] 372 2.138741516 0.07821305 2.2178024 8.4887479 9.187119 # ltdi.2.hat se.ltdi.2.hat ucb.ltdi.2 tdi.2.hat ucb.tdi.2 # [1,] 0.3139527 0.02359172 0.4460881 1.368825 1.562189 # [2,] 2.3750503 0.08033260 2.4559533 10.751554 11.657542 # > # > apply(conf.results$tdi, MAR=2, range) # vgrid ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # [1,] 45 2.712702 0.04790261 2.820059 15.06993 16.77784 # [2,] 372 2.886002 0.06526867 2.964795 17.92152 19.39073 # > # > apply(conf.results$ccc, MAR=2, range) # vgrid zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # [1,] 45 2.493227 0.08536508 2.352814 0.986433 0.9820737 # [2,] 372 2.661442 0.09682164 2.502184 0.990290 0.9866723 # > #################################################### # recalibrate Ektachem measurements by subtracting the alpha.hat and # redo the analysis #################################################### # > fit.results.power$param.hat["alpha"] # alpha # 5.581801 # > chol2 <- chol chol2[chol2$method == "method2", "meas"] <- chol2[chol2$method == "method2", "meas"] - fit.results.power$param.hat["alpha"] fit.results.power2 <- lme.repmeas.hetero.fit(chol2, v, "power", vgrid) conf.results2 <- conf.measures.hetero.repmeas(chol2, fit.results.power2$param.hat, v, "power", vgrid) # > apply(conf.results2$tdi, MAR=2, range) # vgrid ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # [1,] 45 2.488890 0.04587677 2.606896 12.04789 13.55690 # [2,] 372 2.735185 0.07174238 2.810646 15.41260 16.62065 # > # > apply(conf.results2$ccc, MAR=2, range) # vgrid zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # [1,] 45 2.644179 0.0841618 2.505745 0.9899507 0.9867662 # [2,] 372 2.889490 0.1007346 2.723797 0.9938353 0.9914236 # > # estimates and confidence bands for various measures # lambda par(mfrow = c(2, 2), pty = "m", oma = c(0.5, 0, 0, 0), mar = c(4.5, 4.5, 2.5, 0.5) + 0.1, cex.lab = 0.8, cex.main = 0.8, cex.axis = 0.8) par(mfrow = c(2, 2)) plot(vgrid, conf.results$lambda[, "est"], ylim = c(min(conf.results$lambda[, "lcl.par"]), max(conf.results$lambda[, "ucl.par"])), type = "l", xlab = "cholesterol", ylab = "precision ratio", panel.first = grid(lty = "solid")) lines(vgrid, conf.results$lambda[, "lcl.par"], lty = 2) lines(vgrid, conf.results$lambda[, "ucl.par"], lty = 2) title(main = "(a)") # CCC plot(vgrid, conf.results$ccc[, "lcb.ccc"], type = "l", xlab = "cholesterol", ylab = "CCC", ylim = c(0.98, 1), panel.first = grid(lty = "solid")) lines(vgrid, conf.results$rep.ccc[, "lcb.ccc.1"], lty = 2) lines(vgrid, conf.results$rep.ccc[, "lcb.ccc.2"], lty = 3) legend(50, 0.986, legend = c("inter-method", "cobasb", "echem"), lty = 1:3, bty = "n") title(main = "(b)") # TDI plot(vgrid, conf.results$tdi[, "ucb.tdi"], ylim = c(-20, 20), type = "l", xlab = "cholesterol", ylab = "difference", panel.first = grid(lty = "solid")) lines(vgrid, -conf.results$tdi[, "ucb.tdi"], lty = 1) lines(vgrid, conf.results$rep.tdi[, "ucb.tdi.1"], lty = 2) lines(vgrid, -conf.results$rep.tdi[, "ucb.tdi.1"], lty = 2) lines(vgrid, conf.results$rep.tdi[, "ucb.tdi.2"], lty = 3) lines(vgrid, -conf.results$rep.tdi[, "ucb.tdi.2"], lty = 3) legend(50, -5, legend = c("inter-method", "cobasb", "echem"), lty = 1:3, bty = "n") title(main = "(c)") # TDI after recalibration plot(vgrid, conf.results2$tdi[, "ucb.tdi"], ylim = c(-20, 20), type = "l", xlab = "cholesterol", ylab = "difference", panel.first = grid(lty = "solid")) lines(vgrid, -conf.results2$tdi[, "ucb.tdi"], lty = 1) lines(vgrid, conf.results2$rep.tdi[, "ucb.tdi.1"], lty = 2) lines(vgrid, -conf.results2$rep.tdi[, "ucb.tdi.1"], lty = 2) lines(vgrid, conf.results2$rep.tdi[, "ucb.tdi.2"], lty = 3) lines(vgrid, -conf.results2$rep.tdi[, "ucb.tdi.2"], lty = 3) legend(50, -5, legend = c("inter-method", "cobasb", "echem"), lty = 1:3, bty = "n") title(main = "(d)") #################################################### # Analysis assuming the homoscedastic model ####### #################################################### fit.results.homo <- lme.unlinked.repmeas.fit(chol) conf.results.homo <- conf.measures.unlinked(chol, fit.results.homo$param.hat, prob = 0.9) # > conf.results.homo # $param # estimate se lcl ucl # alpha 5.598000 0.74023317 4.147170 7.048830 # mu.b 184.380000 6.59705237 171.450015 197.309985 # lsigmasq.b 8.372129 0.14186996 8.094069 8.650189 # lsigmasq.bI 3.291858 0.14407561 3.009475 3.574241 # lsigmasq.e1 1.347712 0.04713977 1.255320 1.440105 # lsigmasq.e2 1.831110 0.04714155 1.738715 1.923506 # $similarity # lest lest.se lcl.lpar ucl.lpar est lcl.par # lambda -0.48339804 0.06666696 -0.6140629 -0.35273319 0.6166843 0.5411478 # newprecratio -0.07493764 0.01403047 -0.1024369 -0.04743843 0.9278013 0.9026351 # ucl.par # lambda 0.7027647 # newprecratio 0.9536692 # $tdi # ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # 2.77253949 0.05860543 2.86893684 15.99921231 17.61827716 # $ccc # zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # 2.60389973 0.09224109 2.45217665 0.98911219 0.98528066 # $rep.tdi # ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # method1 1.518081 0.02356989 1.556850 4.56346 4.743856 # method2 1.759780 0.02357078 1.798551 5.81116 6.040886 # $rep.ccc # zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # method1 3.862102 0.07429899 3.739891 0.9991164 0.9988719 # method2 3.620541 0.07428161 3.498358 0.9985680 0.9981719 # > # after recalibration fit.results.homo2 <- lme.unlinked.repmeas.fit(chol2) conf.results.homo2 <- conf.measures.unlinked(chol2, fit.results.homo2$param.hat, prob = 0.9) # > conf.results.homo2 # $param # estimate se lcl ucl # alpha 0.01619907 0.74023063 -1.434626 1.467024 # mu.b 184.38000000 6.59705229 171.450015 197.309985 # lsigmasq.b 8.37212865 0.14186996 8.094069 8.650189 # lsigmasq.bI 3.29185808 0.14407561 3.009475 3.574241 # lsigmasq.e1 1.34771241 0.04713977 1.255320 1.440105 # lsigmasq.e2 1.83111045 0.04714155 1.738715 1.923506 # $similarity # lest lest.se lcl.lpar ucl.lpar est lcl.par # lambda -0.48339804 0.06666696 -0.6140629 -0.35273319 0.6166843 0.5411478 # newprecratio -0.07493764 0.01403047 -0.1024369 -0.04743843 0.9278013 0.9026351 # ucl.par # lambda 0.7027647 # newprecratio 0.9536692 # $tdi # ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # 2.57611806 0.06070855 2.67597474 13.14600696 14.52650248 # $ccc # zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # 2.80258823 0.09316337 2.64934813 0.99266942 0.99005350 # $rep.tdi # ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # method1 1.518081 0.02356989 1.556850 4.56346 4.743856 # method2 1.759780 0.02357078 1.798551 5.81116 6.040886 # $rep.ccc # zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # method1 3.862102 0.07429899 3.739891 0.9991164 0.9988719 # method2 3.620541 0.07428161 3.498358 0.9985680 0.9981719 # > #######################################################################################