# Illustration of the R program for implementing the methodology proposed in 
# "Semiparametric Modeling and Analysis of Paired Longitudinal Method Comparison Data" 
# by L. N. Rathnayake and P. K. Choudhary.

# Download the program:
# Right click on the link "R Code" and save the 
# associated file as longitudinal_data_analysis.R in the R working directory.

# Requires:
# The program requires "nlme," "mvtnorm,", "multcomp" and "Matrix" packages. 
# They can be obtained from the R website.

# Illustration using body fat data:

# Step 1. Launch R.

# Step 2. Change the working directory to the directory where 
# longitudinal_data_analysis.R file is stored.

# Step 3. Run source("longitudinal_data_analysis.R"). This will also load 
# the required packages (provided they are already installed)

# Step 4. Run the following commands in R

# show a snapshot of the body fat data (Note: the actual data are not being provided as we don’t have the requisite permission. The program can be used for any data that has the structure of the body fat data)
# > str(fat)
# 'data.frame':	1515 obs. of  5 variables:
 # $ subject: Factor w/ 112 levels "101","102","103",..: 1 1 1 1 1 1 1 1 1 1 ...
 # $ rep    : Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 2 ...
 # $ time   : num  12.3 12.9 13.3 13.8 14.3 ...
 # $ method : Factor w/ 2 levels "caliper","DEXA": 1 1 1 1 1 1 1 1 1 2 ...
 # $ meas   : num  20.8 21.7 23.2 26.3 23.9 ...
 # - attr(*, "na.action")=Class 'omit'  Named int [1:203] 10 28 46 51 64 66 82 87 90 100 ...
  # .. ..- attr(*, "names")= chr [1:203] "2" "20" "38" "48" ...
# > 

# order the data by subject, then by method, and then by rep
# (this is ordering is assumed for creating design and covariance matrices)

fat <- with(fat, fat[order(subject, method, rep), ])

# fit the proposed mixed-effects model without interaction and continuous AR(1) errors
method.names <- levels(droplevels(fat$method))
rep.names <- levels(droplevels(fat$rep))
subid <- unique(fat$subject)
nsub <- length(subid)
nobs <- nrow(fat)
nknots <- 35
knots <- quantile(unique(fat$time), seq(0, 1, l = (nknots + 2)), names = FALSE)[-c(1, nknots + 2)]

# fixed-effects design matrices
X <- model.matrix(~ -1 + method + method:time + method:I(time^2), data = fat)

# nobs x nknots matrix for method 1 spline	
W.1 <- outer(fat$time, knots, "-") * X[, paste0("method", method.names[1])]
W.1 <- W.1^2 * (W.1 > 0)     

# nobs x nknots matrix for method 2 spline
W.2 <-  outer(fat$time, knots, "-") * X[, paste0("method", method.names[2])]			
W.2 <- W.2^2 * (W.2 > 0)

# create a grouping variable to fit spline models
group.spline <- rep(1, nobs)

#  fit model
fit.noint <- lme(meas ~ - 1 + X, data = fat, method = "ML", 
	random = list(group.spline = pdBlocked(list(pdIdent(~ - 1 + W.1), pdIdent(~ - 1 + W.2))), 
	subject = pdBlocked(list(~1, pdIdent(~ - 1 + rep)))),
	weights = varIdent(form = ~1 | method),
	correlation = corCAR1(form = ~ time|group.spline/subject/method))

# set up the grid for inference
ngrid <- 30
x.g <- with(fat, seq(f = min(time), t = max(time), l = ngrid))
X.g <- cbind(rep(1, ngrid), x.g, x.g^2) # ngrid x 3 matrix
W.g <- outer(x.g, knots, "-")  
W.g <- W.g^2 * (W.g > 0)     # ngrid x nknots matrix

# get the estimates
prob <- 0.90
estimates.noint <- get.param.estimates.noint(fit.noint, X.g, W.g, method.names, rep.names, colnames(X), prob, err.cor = T) 	

# design matrix for random effects
# model without subject x method interaction
# nobs x nsub*(1 + # levels of rep) block diagonal matrix
Z.noint <- model.matrix(~ 1 + rep, data = fat, 
	contrasts = list(rep = contrasts(fat$rep, F)))
Z.noint <- lapply(split(data.frame(Z.noint), fat$subject), as.matrix)
Z.noint <- as.matrix(bdiag(Z.noint))

# get the bootstrap confidence intervals

set.seed(12345)	

nboot <- 500 #(takes about 200 minutes on a Windows computer with a 2.2 GHz processor and 4 GB RAM)

alpha <- 0.05

boot.noint <- bootstrap.draws.noint(fit.noint, dat = fat, prob = prob, Zmat = Z.noint, err.cor = T) 

ci.noint <- bootstrap.ci(boot.noint, level = 1 - alpha)

# The results stored ci.noint have been shown in Figures 1-3 of the manuscript.