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
# 
# 
# > 

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