source("MAMMA.R") ################################################ # Systolic Blood Pressure Data # ################################################ # Note: These data are not publicly available. In the code below, it is # assumed that there is a dataset called "sbp" with the following # characteristics: # > str(sbp) # 'data.frame': 912 obs. of 3 variables: # $ subject: Factor w/ 228 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ... # $ method : Factor w/ 4 levels "MS1","MS2","MS3",..: 1 1 1 1 1 1 1 1 1 1 ... # $ meas : num 110 116 118 112 108 82 82 86 92 90 ... # > ################################################ # Displaying data # ################################################ # trellis plot plot(groupedData(meas~method | subject, data=sbp, order.groups=TRUE), ylim=c(-2,231), xlab ="systolic blood pressure (mm~Hg)", ylab="sorted subject ID", pch=4) # boxplots bwplot(meas~method, data=sbp, xlab="observer", pch="|", ylab="systolic blood pressure (mm~Hg)") # matrix of scatterplots and BA plots ymat <- unlinked2paired(sbp) plot.sba.matrix (ymat) ################################################ # Model fitting and evaluation # ################################################ fit.results <- lme.unreplicated.fit.multi (sbp, "all.pairwise") # > fit.results # $fit # Linear mixed-effects model fit by maximum likelihood # Data: dat # Log-likelihood: -3148.002 # Fixed: meas ~ method - 1 # methodMS1 methodMS2 methodMS3 methodDS # 128.6930 128.9825 130.1842 131.4474 # Random effects: # Formula: ~1 | subject # (Intercept) Residual # StdDev: 30.47186 3.124265 # Variance function: # Structure: Different standard deviations per stratum # Formula: ~1 | method # Parameter estimates: # MS1 MS2 MS3 DS # 1.0000000 1.2017181 0.9425705 1.8257249 # Number of Observations: 912 # Number of Groups: 228 # $method.names # [1] "MS1" "MS2" "MS3" "DS" # $param.hat # alphaMS2 alphaMS3 alphaDS mu.b lsigmasq.b # 0.2894737 1.4912281 2.7543860 128.6929825 6.8336075 # lsigmasq.eMS1 lsigmasq.eMS2 lsigmasq.eMS3 lsigmasq.eDS # 2.2783981 2.6459026 2.1601090 3.4823523 # $varpar.hat # $varpar.hat$sigmasq.b # (Intercept) # 928.5345 # $varpar.hat$sigmasq.e # MS1 MS2 MS3 DS # 9.761031 14.096163 8.672083 32.536167 # $ydist.hat # $ydist.hat$ymean # MS1 MS2 MS3 DS # 128.6930 128.9825 130.1842 131.4474 # $ydist.hat$ycov # MS1 MS2 MS3 DS # MS1 938.2955 928.5345 928.5345 928.5345 # MS2 928.5345 942.6306 928.5345 928.5345 # MS3 928.5345 928.5345 937.2065 928.5345 # DS 928.5345 928.5345 928.5345 961.0706 # $ydist.hat$ycor # MS1 MS2 MS3 DS # MS1 1.0000000 0.9873189 0.9901718 0.9778012 # MS2 0.9873189 1.0000000 0.9878923 0.9755502 # MS3 0.9901718 0.9878923 1.0000000 0.9783691 # DS 0.9778012 0.9755502 0.9783691 1.0000000 # $ydist.hat$yvar # MS1 MS2 MS3 DS # 938.2955 942.6306 937.2065 961.0706 # $contrast.matrix # [,1] [,2] [,3] [,4] # MS2-MS1 -1 1 0 0 # MS3-MS1 -1 0 1 0 # DS-MS1 -1 0 0 1 # MS3-MS2 0 -1 1 0 # DS-MS2 0 -1 0 1 # DS-MS3 0 0 -1 1 # $ddist.hat # $ddist.hat$dmean # MS2-MS1 MS3-MS1 DS-MS1 MS3-MS2 DS-MS2 DS-MS3 # 0.2894737 1.4912281 2.7543860 1.2017544 2.4649123 1.2631579 # $ddist.hat$dcov # MS2-MS1 MS3-MS1 DS-MS1 MS3-MS2 DS-MS2 DS-MS3 # MS2-MS1 23.857194 9.761031 9.761031 -14.096163 -14.09616 0.000000 # MS3-MS1 9.761031 18.433114 9.761031 8.672083 0.00000 -8.672083 # DS-MS1 9.761031 9.761031 42.297198 0.000000 32.53617 32.536167 # MS3-MS2 -14.096163 8.672083 0.000000 22.768246 14.09616 -8.672083 # DS-MS2 -14.096163 0.000000 32.536167 14.096163 46.63233 32.536167 # DS-MS3 0.000000 -8.672083 32.536167 -8.672083 32.53617 41.208249 # $ddist.hat$dvar # MS2-MS1 MS3-MS1 DS-MS1 MS3-MS2 DS-MS2 DS-MS3 # 23.85719 18.43311 42.29720 22.76825 46.63233 41.20825 # $reliability.hat # MS1 MS2 MS3 DS # 0.9895971 0.9850459 0.9907469 0.9661459 # > # model checking plot(fit.results$fit) plot(fit.results$fit, resid(., type="pear") ~ fitted(.) | method, type=c("p","smooth")) plot(fit.results$fit, abs(resid(., type="pear")) ~ fitted(.) | method, type=c("p","smooth")) qqnorm(fit.results$fit, ~resid(., type="pear") | method) qqnorm(fit.results$fit, ~ranef(.) ) boxplot(ranef(fit.results$fit)) hist(ranef(fit.results$fit)[[1]]) ############################################################### # Evaluation of repeatability, similarity and agreement # ############################################################### conf.results <- conf.measures.unreplicated.multi(sbp, 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 # alphaMS2 0.29 0.32 -0.34 0.92 # alphaMS3 1.49 0.28 0.93 2.05 # alphaDS 2.75 0.43 1.91 3.60 # mu.b 128.69 2.03 124.72 132.67 # lsigmasq.b 6.83 0.09 6.65 7.02 # lsigmasq.eMS1 2.28 0.14 2.00 2.56 # lsigmasq.eMS2 2.65 0.12 2.41 2.89 # lsigmasq.eMS3 2.16 0.15 1.86 2.46 # lsigmasq.eDS 3.48 0.10 3.28 3.69 # > # > lapply(conf.results[-1], function(x){round(x$result,2)}) # $mean.diff # estimate se lcl ucl # MS2-MS1 0.29 0.32 -0.54 1.12 # MS3-MS1 1.49 0.28 0.76 2.22 # DS-MS1 2.75 0.43 1.65 3.85 # MS3-MS2 1.20 0.32 0.39 2.01 # DS-MS2 2.46 0.45 1.31 3.62 # DS-MS3 1.26 0.43 0.18 2.35 # $precision.ratio # estimate lcl ucl # MS2/MS1 0.69 0.42 1.14 # MS3/MS1 1.13 0.61 2.08 # DS/MS1 0.30 0.19 0.48 # MS3/MS2 1.63 0.95 2.78 # DS/MS2 0.43 0.28 0.66 # DS/MS3 0.27 0.17 0.43 # $tdi # estimate lcl ucl # MS2,MS1 8.05 0 8.94 # MS3,MS1 7.48 0 8.29 # DS,MS1 11.62 0 12.84 # MS3,MS2 8.09 0 8.98 # DS,MS2 11.94 0 13.13 # DS,MS3 10.76 0 11.92 # $ccc # estimate lcl ucl # MS2,MS1 0.99 0.98 1 # MS3,MS1 0.99 0.99 1 # DS,MS1 0.97 0.97 1 # MS3,MS2 0.99 0.98 1 # DS,MS2 0.97 0.96 1 # DS,MS3 0.98 0.97 1 # > # > lapply(conf.results["ccc"], function(x){round(x$result,3)}) # $ccc # estimate lcl ucl # MS2,MS1 0.987 0.983 1 # MS3,MS1 0.989 0.985 1 # DS,MS1 0.974 0.966 1 # MS3,MS2 0.987 0.983 1 # DS,MS2 0.972 0.964 1 # DS,MS3 0.977 0.970 1 # > #################################################################################