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

# > 


#####################################################################################