source("MAMMA.R") ################################################################## # Tumor Size Data # ################################################################## # Note: These data are not available on the book website but they # can be obtained from Broemeling (2009, Chapter 6) # # reading data tumor <- read.table(file = "tumor_size.csv", header = TRUE, sep = ",") # > str(tumor) # 'data.frame': 400 obs. of 4 variables: # $ lesion : int 1 2 3 4 5 6 7 8 9 10 ... # $ reader : int 1 1 1 1 1 1 1 1 1 1 ... # $ replication: int 1 1 1 1 1 1 1 1 1 1 ... # $ tumor.size : num 3.5 3.8 2.2 1.5 3.8 3.5 4.2 5.4 7.6 2.8 ... # > tumor$lesion <- as.factor(tumor$lesion) tumor$reader <- as.factor(tumor$reader) colnames(tumor) <- c("subject", "method", "rep", "meas") # > str(tumor) # 'data.frame': 400 obs. of 4 variables: # $ subject: Factor w/ 40 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ... # $ method : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ... # $ rep : int 1 1 1 1 1 1 1 1 1 1 ... # $ meas : num 3.5 3.8 2.2 1.5 3.8 3.5 4.2 5.4 7.6 2.8 ... # > # n = 40 lesions, 5 readers, balanced design with 2 replications # total number of observations = 40*5*2 = 40 # tumor size is measured in cm # > summary(tumor$meas) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 1.000 3.000 4.000 4.067 4.800 9.000 # > ################################################ # Displaying data # ################################################ # trellis plot plot(groupedData(meas ~ method | subject, data = tumor, order.groups = TRUE), ylim = c(0, 41), xlab = "tumor size (cm)", ylab = "sorted lesion ID") # boxplots bwplot(meas ~ method, data = tumor, xlab = "reader", pch = "|", ylab = "tumor size (cm)") # matrix of scatterplots and BA plots # randomly formed pairs set.seed(1234) ymat <- unlinked2paired(tumor) plot.sba.matrix(ymat) # averages of replicates ymat.avg <- unlinked2avg(tumor) plot.sba.matrix(ymat.avg) # interaction plots with(tumor, { interaction.plot(method, subject, meas, legend = FALSE, type = "o", lty = 1, pch = 19, cex = 0.5, lwd = 0.8, xlab = "reader", ylab = "average measurement", las = 1) }) ################################################ # Model fitting and evaluation # ################################################ fit.results <- lme.unlinked.repmeas.fit.multi(tumor, "all.pairwise") # > fit.results # $fit # Linear mixed-effects model fit by maximum likelihood # Data: dat # Log-likelihood: -344.1083 # Fixed: meas ~ method - 1 # method1 method2 method3 method4 method5 # 4.02625 3.59500 4.33625 4.27000 4.10875 # # Random effects: # Composite Structure: Blocked # # Block 1: (Intercept) # Formula: ~1 | subject # (Intercept) # StdDev: 1.516779 # # Block 2: method1, method2, method3, method4, method5 # Formula: ~method - 1 | subject # Structure: Diagonal # method1 method2 method3 method4 method5 Residual # StdDev: 0.3737545 0.7397127 0.5446508 0.2432421 0.2740662 0.3059836 # # Variance function: # Structure: Different standard deviations per stratum # Formula: ~1 | method # Parameter estimates: # 1 2 3 4 5 # 1.0000000 1.0752676 0.9715518 1.1783471 1.0507841 # Number of Observations: 400 # Number of Groups: 40 # # $method.names # [1] "1" "2" "3" "4" "5" # # $param.hat # alpha2 alpha3 alpha4 alpha5 mu.b # -0.4312500 0.3100000 0.2437500 0.0825000 4.0262500 # lsigmasq.b lsigmasq.bI1 lsigmasq.bI2 lsigmasq.bI3 lsigmasq.bI4 # 0.8331778 -1.9683124 -0.6029867 -1.2152209 -2.8273962 # lsigmasq.bI5 lsigmasq.e1 lsigmasq.e2 lsigmasq.e3 lsigmasq.e4 # -2.5887710 -2.3684475 -2.2233084 -2.4261688 -2.0402221 # lsigmasq.e5 # -2.2693742 # # $varpar.hat # $varpar.hat$sigmasq.b # (Intercept) # 2.300618 # # $varpar.hat$sigmasq.bI # 1 2 3 4 5 # 0.13969241 0.54717493 0.29664447 0.05916671 0.07511229 # # $varpar.hat$sigmasq.e # 1 2 3 4 5 # 0.09362597 0.10825038 0.08837477 0.12999983 0.10337686 # # # $ydist.hat # $ydist.hat$ymean # 1 2 3 4 5 # 4.02625 3.59500 4.33625 4.27000 4.10875 # # $ydist.hat$ycov # 1 2 3 4 5 # 1 2.533936 2.300618 2.300618 2.300618 2.300618 # 2 2.300618 2.956043 2.300618 2.300618 2.300618 # 3 2.300618 2.300618 2.685637 2.300618 2.300618 # 4 2.300618 2.300618 2.300618 2.489785 2.300618 # 5 2.300618 2.300618 2.300618 2.300618 2.479107 # # $ydist.hat$ycor # 1 2 3 4 5 # 1 1.0000000 0.8406036 0.8819074 0.9159374 0.9179077 # 2 0.8406036 1.0000000 0.8165174 0.8480242 0.8498484 # 3 0.8819074 0.8165174 1.0000000 0.8896926 0.8916064 # 4 0.9159374 0.8480242 0.8896926 1.0000000 0.9260107 # 5 0.9179077 0.8498484 0.8916064 0.9260107 1.0000000 # # $ydist.hat$yvar # 1 2 3 4 5 # 2.533936 2.956043 2.685637 2.489785 2.479107 # # # $contrast.matrix # [,1] [,2] [,3] [,4] [,5] # 2-1 -1 1 0 0 0 # 3-1 -1 0 1 0 0 # 4-1 -1 0 0 1 0 # 5-1 -1 0 0 0 1 # 3-2 0 -1 1 0 0 # 4-2 0 -1 0 1 0 # 5-2 0 -1 0 0 1 # 4-3 0 0 -1 1 0 # 5-3 0 0 -1 0 1 # 5-4 0 0 0 -1 1 # # $ddist.hat # $ddist.hat$dmean # 2-1 3-1 4-1 5-1 3-2 4-2 5-2 4-3 # -0.43125 0.31000 0.24375 0.08250 0.74125 0.67500 0.51375 -0.06625 # 5-3 5-4 # -0.22750 -0.16125 # # $ddist.hat$dcov # 2-1 3-1 4-1 5-1 3-2 4-2 # 2-1 0.8887437 0.2333184 0.2333184 0.2333184 -0.6554253 -0.6554253 # 3-1 0.2333184 0.6183376 0.2333184 0.2333184 0.3850192 0.0000000 # 4-1 0.2333184 0.2333184 0.4224849 0.2333184 0.0000000 0.1891665 # 5-1 0.2333184 0.2333184 0.2333184 0.4118075 0.0000000 0.0000000 # 3-2 -0.6554253 0.3850192 0.0000000 0.0000000 1.0404445 0.6554253 # 4-2 -0.6554253 0.0000000 0.1891665 0.0000000 0.6554253 0.8445919 # 5-2 -0.6554253 0.0000000 0.0000000 0.1784892 0.6554253 0.6554253 # 4-3 0.0000000 -0.3850192 0.1891665 0.0000000 -0.3850192 0.1891665 # 5-3 0.0000000 -0.3850192 0.0000000 0.1784892 -0.3850192 0.0000000 # 5-4 0.0000000 0.0000000 -0.1891665 0.1784892 0.0000000 -0.1891665 # 5-2 4-3 5-3 5-4 # 2-1 -0.6554253 0.0000000 0.0000000 0.0000000 # 3-1 0.0000000 -0.3850192 -0.3850192 0.0000000 # 4-1 0.0000000 0.1891665 0.0000000 -0.1891665 # 5-1 0.1784892 0.0000000 0.1784892 0.1784892 # 3-2 0.6554253 -0.3850192 -0.3850192 0.0000000 # 4-2 0.6554253 0.1891665 0.0000000 -0.1891665 # 5-2 0.8339145 0.0000000 0.1784892 0.1784892 # 4-3 0.0000000 0.5741858 0.3850192 -0.1891665 # 5-3 0.1784892 0.3850192 0.5635084 0.1784892 # 5-4 0.1784892 -0.1891665 0.1784892 0.3676557 # # $ddist.hat$dvar # 2-1 3-1 4-1 5-1 3-2 4-2 5-2 # 0.8887437 0.6183376 0.4224849 0.4118075 1.0404445 0.8445919 0.8339145 # 4-3 5-3 5-4 # 0.5741858 0.5635084 0.3676557 # # # $reliability.hat # 1 2 3 4 5 # 0.9630512 0.9633800 0.9670936 0.9477867 0.9583008 # # > # model checking plot(fit.results$fit) plot(fit.results$fit, abs(resid(., type = "pear")) ~ fitted(.) | method, type = c("p", "smooth", "r")) qqnorm(fit.results$fit, ~resid(., type = "pear") | method) qqnorm(fit.results$fit, ~ranef(.) | method) ################################################################################## ############################################################### # Evaluation of repeatability, similarity and agreement # ############################################################### conf.results <- conf.measures.unlinked.multi(tumor, fit.results$param.hat, fit.results$method.names, fit.results$contrast.matrix, conf.level = 0.95, prob = 0.9, seed = 12345) # > round(conf.results$param, 2) # estimate se lcl ucl # alpha2 -0.43 0.14 -0.71 -0.16 # alpha3 0.31 0.11 0.08 0.54 # alpha4 0.24 0.09 0.07 0.42 # alpha5 0.08 0.09 -0.09 0.26 # mu.b 4.03 0.25 3.54 4.51 # lsigmasq.b 0.83 0.23 0.39 1.28 # lsigmasq.bI1 -1.97 0.39 -2.73 -1.21 # lsigmasq.bI2 -0.60 0.27 -1.13 -0.08 # lsigmasq.bI3 -1.22 0.30 -1.79 -0.64 # lsigmasq.bI4 -2.83 0.73 -4.26 -1.40 # lsigmasq.bI5 -2.59 0.57 -3.71 -1.47 # lsigmasq.e1 -2.37 0.22 -2.81 -1.93 # lsigmasq.e2 -2.22 0.22 -2.66 -1.79 # lsigmasq.e3 -2.43 0.22 -2.86 -1.99 # lsigmasq.e4 -2.04 0.22 -2.48 -1.60 # lsigmasq.e5 -2.27 0.22 -2.71 -1.83 # > # > lapply(conf.results[-1], function(x){round(x$result,2)}) # $mean.diff # estimate se lcl ucl # 2-1 -0.43 0.14 -0.81 -0.05 # 3-1 0.31 0.11 0.00 0.62 # 4-1 0.24 0.09 0.01 0.48 # 5-1 0.08 0.09 -0.16 0.32 # 3-2 0.74 0.15 0.33 1.16 # 4-2 0.67 0.13 0.31 1.04 # 5-2 0.51 0.13 0.15 0.88 # 4-3 -0.07 0.11 -0.36 0.23 # 5-3 -0.23 0.11 -0.52 0.07 # 5-4 -0.16 0.08 -0.38 0.05 # $precision.ratio # estimate lcl ucl # 2/1 0.86 0.37 2.05 # 3/1 1.06 0.45 2.51 # 4/1 0.72 0.30 1.71 # 5/1 0.91 0.38 2.15 # 3/2 1.22 0.52 2.90 # 4/2 0.83 0.35 1.97 # 5/2 1.05 0.44 2.48 # 4/3 0.68 0.29 1.61 # 5/3 0.85 0.36 2.03 # 5/4 1.26 0.53 2.98 # $new.precision.ratio # estimate lcl ucl # 2/1 0.36 0.15 0.84 # 3/1 0.61 0.24 1.54 # 4/1 1.23 0.50 3.06 # 5/1 1.31 0.50 3.39 # 3/2 1.70 0.72 4.04 # 4/2 3.46 1.45 8.30 # 5/2 3.67 1.41 9.58 # 4/3 2.04 0.81 5.12 # 5/3 2.16 0.89 5.21 # 5/4 1.06 0.41 2.74 # $tdi # estimate lcl ucl # 2,1 1.71 0 2.14 # 3,1 1.39 0 1.71 # 4,1 1.14 0 1.40 # 5,1 1.06 0 1.30 # 3,2 2.07 0 2.56 # 4,2 1.87 0 2.33 # 5,2 1.72 0 2.16 # 4,3 1.25 0 1.53 # 5,3 1.29 0 1.61 # 5,4 1.03 0 1.26 # $ccc # estimate lcl ucl # 2,1 0.81 0.68 1 # 3,1 0.87 0.77 1 # 4,1 0.91 0.83 1 # 5,1 0.92 0.85 1 # 3,2 0.74 0.59 1 # 4,2 0.78 0.64 1 # 5,2 0.81 0.68 1 # 4,3 0.89 0.81 1 # 5,3 0.88 0.79 1 # 5,4 0.92 0.86 1 # $rep.tdi # estimate lcl ucl # 1 0.71 0 0.86 # 2 0.77 0 0.92 # 3 0.69 0 0.83 # 4 0.84 0 1.01 # 5 0.75 0 0.90 # $rep.ccc # estimate lcl ucl # 1 0.96 0.94 1 # 2 0.96 0.94 1 # 3 0.97 0.95 1 # 4 0.95 0.91 1 # 5 0.96 0.93 1 # > # > lapply(conf.results["ccc"], function(x){round(x$result,3)}) # $ccc # estimate lcl ucl # 2,1 0.811 0.683 1 # 3,1 0.866 0.770 1 # 4,1 0.905 0.832 1 # 5,1 0.917 0.852 1 # 3,2 0.743 0.593 1 # 4,2 0.780 0.640 1 # 5,2 0.807 0.679 1 # 4,3 0.888 0.807 1 # 5,3 0.882 0.794 1 # 5,4 0.921 0.859 1 # > # > lapply(conf.results["rep.ccc"], function(x){round(x$result,3)}) # $rep.ccc # estimate lcl ucl # 1 0.963 0.939 1 # 2 0.963 0.941 1 # 3 0.967 0.947 1 # 4 0.948 0.914 1 # 5 0.958 0.931 1 # > #####################################################################################