source("MAMMA.R") ############################################### # Unreplicated Blood Pressure Data # ############################################### # Read the data bp <- read.table(file = "unreplicated_blood_pressure.txt", sep = "\t", header = T) # > head(bp) # armsys fingsys # 1 140 154 # 2 110 112 # 3 138 156 # 4 124 152 # 5 142 142 # 6 112 104 # > # > str(bp) # 'data.frame': 200 obs. of 2 variables: # $ armsys : int 140 110 138 124 142 112 122 128 132 160 ... # $ fingsys: int 154 112 156 152 142 104 126 132 144 180 ... # > # These data are from Bland and Altman (1995). They have measurements of systolic # blood pressure from two methods --- a finger pressure method (test method) # and an arm blood pressure method (standard method) # armsys = method 1, fingsys = method 2 # data frame with two columns --- first = method1, second = method2 bp.2col <- data.frame(bp$armsys, bp$fingsys) colnames(bp.2col) <- c("arm pressure", "finger pressure") # data frame with three columns --- subject (factor), # method (factor with levels method1 and method2), and meas bp.3col <- make.3col(bp.2col) # trellis plot plot(groupedData(meas ~ method | subject, data = bp.3col, order.groups = TRUE), ylim = c(-1, 202), cex = 0.5, xlab = "systolic blood pressure (mm~Hg)", ylab = "sorted subject ID") # boxplot bwplot(meas ~ method, data = bp.3col, xlab = "method", pch = "|", ylab = "systolic blood pressure (mm~Hg)") # Scatterplot and Bland-Altman plot y1 <- bp$armsys y2 <- bp$fingsys diff <- y2 - y1 avg <- (y1 + y2)/2 par(mfrow = c(1, 2), pty = "m", oma = c(0.5, 0, 0, 0), mar = c(4.5, 4.5, 2.5, 0.5) + 0.1, cex.main = 0.8, cex.axis = 0.8, cex.lab = 0.8) plot(y1, y2, xlab = "arm pressure (mm~Hg)", ylab = "finger pressure (mm~Hg)", panel.first = grid(lty = "solid")) abline(a = 0, b = 1) title(main = "(a)") plot(avg, diff, xlab = "average", ylab = "finger $-$ arm", panel.first = grid(lty = "solid")) abline(h = 0) title(main = "(b)") par(mfrow = c(1, 1)) ################################################################################ ################################################ # Analysis based on bivariate normal model # ################################################ # model fitting fit.results.bvn <- bvn.fit(bp.2col) # > fit.results.bvn # $param.hat # mu.1 mu.2 lsigmasq.1 lsigmasq.2 zrho # 128.520000 132.815000 6.290827 6.483934 1.177170 # $bvn.hat # mean1 mean2 var1 var2 corr # 128.5200000 132.8150000 539.5996000 654.5407750 0.8265572 # $loglik # [1] -1730.104 # $aic # [1] 3470.208 # $bic # [1] 3490.165 # $D.hat # mean var # 4.295 211.698 # $loa # [1] -24.80469 33.39469 # $fit.flag # [1] "yes" # $param.hat.lme # alpha mu.b lsigmasq.b lsigmasq.e1 lsigmasq.e2 # 4.295000 128.520000 6.196895 3.879053 5.095709 # > # confidence intervals conf.results.bvn <- conf.measures(bp.2col, fit.results.bvn$param.hat, fit.method = "bvn", bvn.fit.flag = "yes", prob = 0.9) # > conf.results.bvn # $param # estimate se lcl ucl # mu.1 128.520000 1.64255837 125.300645 131.739355 # mu.2 132.815000 1.80906160 129.269304 136.360696 # lsigmasq.1 6.290827 0.10000000 6.094831 6.486824 # lsigmasq.2 6.483934 0.10000000 6.287937 6.679930 # zrho 1.177170 0.07071068 1.038580 1.315761 # $similarity # $similarity$meandiff # est se lcl ucl # 4.295000 1.028829 2.278531 6.311469 # $similarity$varpar # lest lest.se lcl.lpar ucl.lpar est lcl.par # lambda -1.2166554 0.63373960 -2.45876223 0.02545137 0.2962192 0.08554077 # varratio 0.1931065 0.07959939 0.03709457 0.34911843 1.2130120 1.03779116 # ucl.par # lambda 1.025778 # varratio 1.417817 # $tdi # ltdi.hat se.ltdi.hat ucb.ltdi tdi.hat ucb.tdi # 3.21710285 0.04984997 3.29909875 24.95571494 27.08821460 # $ccc # zccc.hat se.zccc.hat lcb.zccc ccc.hat lcb.ccc # 1.12762068 0.06883373 1.01439927 0.81020337 0.76757547 # > # model checking par(mfrow = c(2, 2)) qqnorm(y1) qqnorm(y2) qqnorm(y1 - y2) qqnorm(y1 + y2) par(mfrow = c(1, 1)) ################################################ # Nonparametric analysis # ################################################ npm.results <- npm.inference(bp.3col) # > npm.results # $moments # meanarm pressure meanfinger pressure # 128.5200000 132.8150000 # sdarm pressure sdfinger pressure # 23.2292832 25.5839945 # corrarm pressurefinger pressure # 0.8265572 # # $mean.diff # $mean.diff$ci.fun # Warning: namespace ‘multcomp’ is not available and has been replaced # by .GlobalEnv when processing object ‘’ # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = 1.96 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 == 0 4.2950 2.2785 6.3115 # # # $mean.diff$result # estimate se lcl ucl # 1 4.295 1.028829 2.278531 6.311469 # # # $varratio # $varratio$ci.fun # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = 1.96 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 == 0 0.19311 0.03731 0.34890 # # # $varratio$lresult # estimate se lcl ucl # 1 0.1931065 0.07948843 0.03731205 0.348901 # # $varratio$result # estimate lcl ucl # 1 1.213012 1.038017 1.417509 # # # $cp # $cp$ci.fun # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = 1.6449 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 >= 0 0.9000 -Inf 0.9349 # # # $cp$result # estimate se lcl ucl # 1 0.9 0.0212132 -Inf 0.9348926 # # # $tdi # $tdi$tdi.hat # tdi # 23 # # $tdi$ucb.tdi # [1] 28 # # # $ccc # $ccc$ci.fun # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = -1.6449 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 <= 0 1.1276 0.9808 Inf # # # $ccc$zresult # estimate se lcl ucl # 1 1.127621 0.08927095 0.980783 Inf # # $ccc$result # estimate lcl ucl # 1 0.8102034 0.7534047 1 # # # > ################ # Recalibrate finger by subtracting 4.30 from its measurements bp.3col.new <- bp.3col bp.3col.new[bp.3col$method == "finger pressure", "meas"] <- bp.3col[bp.3col$method == "finger pressure", "meas"] - 4.3 npm.results.new <- npm.inference(bp.3col.new) # > npm.results.new # $moments # meanarm pressure meanfinger pressure # 128.5200000 128.5150000 # sdarm pressure sdfinger pressure # 23.2292832 25.5839945 # corrarm pressurefinger pressure # 0.8265572 # # $mean.diff # $mean.diff$ci.fun # Warning: namespace ‘multcomp’ is not available and has been replaced # by .GlobalEnv when processing object ‘’ # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = 1.96 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 == 0 -0.0050 -2.0215 2.0115 # # # $mean.diff$result # estimate se lcl ucl # 1 -0.005 1.028829 -2.021469 2.011469 # # # $varratio # $varratio$ci.fun # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = 1.96 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 == 0 0.19311 0.03731 0.34890 # # # $varratio$lresult # estimate se lcl ucl # 1 0.1931065 0.07948843 0.03731205 0.348901 # # $varratio$result # estimate lcl ucl # 1 1.213012 1.038017 1.417509 # # # $cp # $cp$ci.fun # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = 1.6449 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 >= 0 0.9150 -Inf 0.9474 # # # $cp$result # estimate se lcl ucl # 1 0.915 0.01971991 -Inf 0.9474364 # # # $tdi # $tdi$tdi.hat # tdi # 23.7 # # $tdi$ucb.tdi # [1] 26.3 # # # $ccc # $ccc$ci.fun # # Simultaneous Confidence Intervals # # Fit: NULL # # Quantile = -1.6449 # 95% family-wise confidence level # # # Linear Hypotheses: # Estimate lwr upr # 1 <= 0 1.1652 1.0137 Inf # # # $ccc$zresult # estimate se lcl ucl # 1 1.165175 0.09207549 1.013725 Inf # # $ccc$result # estimate lcl ucl # 1 0.8227193 0.7672982 1 # # # > #################################################################################