# Update history: # Posted online on Feb 18, 2020 ################################################################################# # R functions to implement the methodology proposed in the article titled # "Modeling and Analysis of Functional Method Comparison Data" # by G. S. R. de Silva, L. N. Rathnayake and P. K. Choudhary # For illustration, see the files Illustration_1.R and Illustration_2.R. # Requires: require(mvtnorm) require(mgcv) require(MASS) ################################################################################# ################################################################################# # Function for estimating model parameters for MPACE method # Arguments: # dataset = bivariate data # x0 = argument values for gam model # ncol = number of points in the common grid # nrow = number of subjects # nbasis = number of basis # pve.biv = proportion of total variability observed # prob = probability value # makePD = compute the nearest positive definite matrix # # Values: # mean = mean function of varaibles y1 and y2 # mean.diff = mean diffence function between two variables # evalues = eigen values # efunctions = eigen functions # npc = number of principal components # ltau.sq.1 = log scale of measurement error of y1, log(tau.sq.1) # ltau.sq.2 = log scale of measurement error of y2, log(tau.sq.2) # ltau.sq.ratio = log proportion of variance between two variables # fit = estimated fitted curves # scores = estimated scores # zcor.coef = zcor.coefficient between two methods # zccc = Fisher's z transformation of Corcordance Correlation Coefficient (CCC). # ltdi = log Total Deviation Index (TDI) value MPACE = function(dataset, x0, ncol, nrow, nbasis = 10, pve.biv = 0.99, prob = 0.90, makePD = FALSE) { # Mean estimate Y.1 = dataset[ ,1:ncol] Y.2 = dataset[ ,(ncol+1):(2*ncol)] d.vec.Y.1 = rep(x0, each = nrow) d.vec.Y.2 = rep(x0, each = nrow) gam_model.1 = gam(as.vector(Y.1) ~ s(d.vec.Y.1, k=nbasis)) mu.1 = predict(gam_model.1, newdata = data.frame(d.vec.Y.1 = x0)) gam_model.2 = gam(as.vector(Y.2) ~ s(d.vec.Y.2, k=nbasis)) mu.2 = predict(gam_model.2, newdata = data.frame(d.vec.Y.2 = x0)) mu = c(mu.1, mu.2) # Covariance estimate Y.tilde = dataset - matrix(mu, nrow, 2*ncol, byrow = T) cov.sum = cov.count = cov.mean = matrix(0,2*ncol,2*ncol) for(i in 1:nrow) { obs.points = which(!is.na(dataset[i,])) cov.count[obs.points, obs.points] = cov.count[obs.points, obs.points] + 1 cov.sum[obs.points, obs.points] = cov.sum[obs.points, obs.points] + tcrossprod(Y.tilde[i, obs.points]) } G0 = ifelse(cov.count == 0, NA, cov.sum/(cov.count)) # Unsmoothed covariance matrix diag.G0 = diag(G0) diag(G0) = NA G.11 = G0[1:ncol, 1:ncol] G.12 = G0[1:ncol, (ncol+1):(2*ncol)] G.22 = G0[(ncol+1):(2*ncol), (ncol+1):(2*ncol)] # Smooth the covariance matrix row.vec = rep(x0, each = ncol) col.vec = rep(x0, ncol) # We removed diagonal elements from covariance matrix gam_model.11 = gam(as.vector(G.11) ~ te(row.vec, col.vec, k=nbasis), weights = as.vector(cov.count[1:ncol, 1:ncol])) G.11 = matrix(predict(gam_model.11, newdata = data.frame(row.vec = row.vec, col.vec = col.vec)), ncol, ncol) gam_model.12 = gam(as.vector(G.12) ~ te(row.vec, col.vec, k=nbasis), weights = as.vector(cov.count[1:ncol, 1:ncol])) G.12 = matrix(predict(gam_model.12, newdata = data.frame(row.vec = row.vec, col.vec = col.vec)), ncol, ncol) gam_model.22 = gam(as.vector(G.22) ~ te(row.vec, col.vec, k=nbasis), weights = as.vector(cov.count[1:ncol, 1:ncol])) G.22 = matrix(predict(gam_model.22, newdata = data.frame(row.vec = row.vec, col.vec = col.vec)), ncol, ncol) # Make a symmetric matrix G.11 = (G.11 + t(G.11))/2 G.22 = (G.22 + t(G.22))/2 G.1 = cbind(G.11, G.12) G.2 = cbind(t(G.12), G.22) G = rbind(G.1, G.2) # To create positive definite covariance matrix if (makePD) { G = { tmp = Matrix::nearPD(G, corr = FALSE, keepDiag = FALSE, do2eigen = TRUE, trace = options()$verbose) as.matrix(tmp$mat) } } # exact number of principal components w = funData:::.intWeights(x0, method = "trapezoidal") w = rep(w, 2) Wsqrt = diag(sqrt(w), ncol = (2*ncol), nrow = (2*ncol)) Winvsqrt = diag(1/sqrt(w), ncol = (2*ncol), nrow = (2*ncol)) V = Wsqrt %*% G %*% Wsqrt eigen.V = eigen(V, symmetric = TRUE) evalues = eigen.V$values evalues = replace(evalues, which(evalues <= 0), 0) no.pc = min(which(cumsum(evalues)/sum(evalues) >= pve.biv)) evalues = evalues[1:no.pc] e.functions = Winvsqrt %*% eigen.V$vectors[, seq(len = no.pc)] cov.hat = e.functions %*% tcrossprod(diag(evalues, nrow = no.pc, ncol = no.pc), e.functions) # Error variance estimate T.len = x0[ncol] - x0[1] T.min = min(which(x0 >= x0[1] + 0.25 * T.len)) T.max = max(which(x0 <= x0[ncol] - 0.25 * T.len)) Diag.1 = (diag.G0-diag(cov.hat))[T.min:T.max] Diag.2 = (diag.G0-diag(cov.hat))[(ncol + T.min):(ncol + T.max)] w = funData:::.intWeights(x0[T.min:T.max], method = "trapezoidal") tau.sq.1 = max(1/(x0[T.max] - x0[T.min]) * sum(Diag.1*w,na.rm = T), 0) tau.sq.2 = max(1/(x0[T.max] - x0[T.min]) * sum(Diag.2*w,na.rm = T), 0) tau.sq = c(tau.sq.1, tau.sq.2) # Principal component scores Z = e.functions err.var = c(rep(tau.sq, each=ncol)) H = t(evalues*t(Z)) sigma.mat = cov.hat + diag(err.var, ncol = (2*ncol), nrow = (2*ncol)) scores = matrix(NA, nrow = nrow, ncol = no.pc) fit = matrix(NA, nrow = nrow, ncol = 2*ncol) for(i in 1:nrow) { obs.points = which(!is.na(Y.tilde[i, ])) inv.sig = chol2inv(chol(sigma.mat[obs.points, obs.points])) scores[i, ] = Y.tilde[i, obs.points] %*% inv.sig %*% H[obs.points, ] fit[i, ] = t(as.matrix(mu)) + tcrossprod(scores[i, ], Z) } # Mean difference mean.diff = as.vector(mu.1 - mu.2) # Variance ratio if(tau.sq.2 > 0) { ltau.sq.ratio = log(tau.sq.1/tau.sq.2) } # Concordance correlation coefficient cov.12 = diag(cov.hat[1:ncol,(ncol+1):(2*ncol)]) var.1 = diag(cov.hat)[1:ncol] + tau.sq.1 var.2 = diag(cov.hat)[(ncol+1):(2*ncol)] + tau.sq.2 zcor.coef = atanh(cov.12/(sqrt(var.1)*sqrt(var.2))) zccc = ZCCC(mean.diff, var.1, var.2, cov.12) # Total deviation index diff.var = var.1 + var.2 - 2*cov.12 diff.sd = sqrt(diff.var) l.tdi = mapply(log_tdi, mean.diff = mean.diff, diff.sd = diff.sd, MoreArgs = list(prob)) results = list(mean = mu, mean.diff = mean.diff, evalues = evalues, efunctions = e.functions, npc = no.pc, ltau.sq = log(tau.sq), ltau.sq.ratio = ltau.sq.ratio, fit = fit, scores = scores, zcor.coef = zcor.coef, ZCCC = zccc, log.tdi = l.tdi, sd.biv.1 = sqrt(var.1), sd.biv.2 = sqrt(var.2)) return(results) } ############################################ # Function for estimating model parameters for univariate data # Arguments: # dataset = univariate data # x0 = argument values for gam model # ncol = number of points in the common grid # nrow = number of subjects # nbasis = number of basis # pve.univ = proportion of variability # prob = probability value # makePD = compute the nearest positive definite matrix # # Values: # mean = mean function of data # evalues = eigen values # efunctions = eigen functions # npc = number of principal components # sigma.sq = measurement error of data # scores = estimated scores univariate = function(dataset, x0, ncol, nrow, nbasis = 10, pve.univ = 0.99, makePD = FALSE) { Y = dataset[ ,1:ncol] d.vec.Y = rep(x0, each = nrow) gam_model = gam(as.vector(Y) ~ s(d.vec.Y, k = nbasis)) mu = predict(gam_model, newdata = data.frame(d.vec.Y = x0)) Y.tilde = dataset - matrix(mu, nrow, ncol, byrow = T) cov.sum = cov.count = cov.mean = matrix(0, ncol, ncol) for(i in 1:nrow) { obs.points = which(!is.na(dataset[i,])) cov.count[obs.points, obs.points] = cov.count[obs.points, obs.points] + 1 cov.sum[obs.points, obs.points] = cov.sum[obs.points, obs.points] + tcrossprod(Y.tilde[i, obs.points]) } G0 = ifelse(cov.count == 0, NA, cov.sum/(cov.count)) # Unsmoothed covariance matrix diag.G0 = diag(G0) diag(G0) = NA # Smooth the covariance matrix row.vec = rep(x0, each = ncol) col.vec = rep(x0, ncol) # We removed diagonal elements from covariance matrix gam_model = gam(as.vector(G0) ~ te(col.vec, row.vec, k=nbasis), weights = as.vector(cov.count[1:ncol, 1:ncol])) G = matrix(predict(gam_model, newdata = data.frame(row.vec = row.vec, col.vec = col.vec)), ncol, ncol) # To make the symmetric covariance matrix G = (G + t(G))/2 # To create positive definite covariance matrix if (makePD) { G = { tmp = Matrix::nearPD(G, corr = FALSE, keepDiag = FALSE, do2eigen = TRUE, trace = options()$verbose) as.matrix(tmp$mat) } } w = funData:::.intWeights(x0, method = "trapezoidal") Wsqrt = diag(sqrt(w), ncol = ncol, nrow = ncol) Winvsqrt = diag(1/(sqrt(w)), ncol = ncol, nrow = ncol) V = Wsqrt %*% G %*% Wsqrt eigen.V = eigen(V, symmetric = TRUE) evalues = eigen.V$values evalues = replace(evalues, which(evalues <= 0), 0) no.pc = min(which(cumsum(evalues)/sum(evalues) >= pve.univ)) evalues = evalues[1:no.pc] e.functions = Winvsqrt %*% eigen.V$vectors[, seq(len = no.pc)] cov.hat = e.functions %*% tcrossprod(diag(evalues, nrow = no.pc, ncol = no.pc), e.functions) T.len = x0[ncol] - x0[1] T.min = min(which(x0 >= x0[1] + 0.25 * T.len)) T.max = max(which(x0 <= x0[ncol] - 0.25 * T.len)) Diag = (diag.G0-diag(cov.hat))[T.min:T.max] w <- funData:::.intWeights(x0[T.min:T.max], method = "trapezoidal") sigma.sq = max(1/(x0[T.max] - x0[T.min]) * sum(Diag*w,na.rm = T), 0) # Principal components scores Z = e.functions H = t(evalues*t(Z)) err.var = c(rep(sigma.sq, each=ncol)) sigma.mat = cov.hat + diag(err.var, ncol = ncol, nrow = ncol) scores = matrix(NA, nrow = nrow, ncol = no.pc) for(i in 1:nrow) { obs.points = which(!is.na(Y.tilde[i, ])) inv.sig = chol2inv(chol(sigma.mat[obs.points, obs.points])) scores[i, ] = Y.tilde[i, obs.points] %*% inv.sig %*% H[obs.points, ] } results = list(mean = mu, evalues = evalues,efunctions = e.functions,sigma.sq = sigma.sq, npc = no.pc, scores = scores) return(results) } ############################################ # Function for estimating model parameters for UPACE method # Arguments: # dataset.1, dataset.2 = bivariate data # x0 = argument values for gam model # ncol = number of points in the common grid # nrow = number of subjects # nbasis = number of basis # pve.univ = proportion of variability # pve.biv = proportion of total variability # prob = probability value # makePD = compute the nearest positive definite matrix # # Values: # mean = mean function of varaibles y1 and y2 # mean.diff = mean diffence function between two variables # evalues = eigen values # efunctions = eigen functions # npc = number of principal components # ltau.sq.1 = log scale of measurement error of y1, log(tau.sq.1) # ltau.sq.2 = log scale of measurement error of y2, log(tau.sq.2) # ltau.sq.ratio = log proportion of variance between two variables # fit = estimated fitted curves # scores = estimated scores # zcor.coef = zcor.coefficient between two methods # zccc = Fisher's z transformation of Corcordance Correlation Coefficient (CCC). # ltdi = log Total Deviation Index (TDI) value UPACE = function(dataset.1, dataset.2, x0, ncol, nrow, nbasis = 10, pve.univ = 0.99, pve.biv = 0.99, prob = 0.90, makePD = FALSE) { s.1 = univariate(dataset.1, x0, ncol, nrow, nbasis, pve.univ, makePD) s.2 = univariate(dataset.2, x0, ncol, nrow, nbasis, pve.univ, makePD) scores = cbind(s.1$scores, s.2$scores) Z = (t(scores) %*% scores)/(nrow-1) # Eigen values and eigen functions with dimention reduction eigen.Z = eigen(Z, symmetric = TRUE) evalues = eigen.Z$values evalues = replace(evalues, which(evalues <= 0), 0) M = min(which(cumsum(evalues)/sum(evalues) >= pve.biv)) evalues = evalues[1:M] evectors = eigen.Z$vectors[ , 1:M] v.1 = s.1$efunctions %*% evectors[(1:s.1$npc), 1:M] v.2 = s.2$efunctions %*% evectors[(s.1$npc + 1):(s.1$npc + s.2$npc), 1:M] e.functions = list(v.1, v.2) # Number of principal components npc = M # Error variances tau.sq = c(s.1$sigma.sq, s.2$sigma.sq) # Mean mu = c(s.1$mean, s.2$mean) # Mean difference mean.diff = as.vector(s.1$mean - s.2$mean) # Variance ratio if(s.2$sigma.sq > 0) { ltau.sq.ratio = log(s.1$sigma.sq/s.2$sigma.sq) } # Estimated bivariate covariance efun = rbind(v.1, v.2) cov.hat = efun %*% tcrossprod(diag(evalues, nrow = M, ncol = M), efun) # PC scores scores.biv = scores %*% evectors[ ,1:M] # Fitted values mu.mat = t(matrix(rep(mu, nrow), nrow = 2*ncol, ncol = nrow)) fit.biv = mu.mat + tcrossprod(scores.biv, efun) # Concordance correlation coefficient with dimention reduction cov.12 = diag(cov.hat[1:ncol,(ncol+1):(2*ncol)]) var.1 = diag(cov.hat)[1:ncol] + s.1$sigma.sq var.2 = diag(cov.hat)[(ncol+1):(2*ncol)] + s.2$sigma.sq zcor.coef = atanh(cov.12/(sqrt(var.1)*sqrt(var.2))) zccc = ZCCC(mean.diff, var.1, var.2, cov.12) # Total deviation index with dimention reduction diff.var = var.1 + var.2 - 2*cov.12 diff.sd = sqrt(diff.var) l.tdi = mapply(log_tdi, mean.diff = mean.diff, diff.sd = diff.sd, MoreArgs = list(prob)) results = list(mean = mu, mean.diff = mean.diff, evalues = evalues, efunctions = e.functions, npc = npc, ltau.sq = log(tau.sq), ltau.sq.ratio = ltau.sq.ratio, fit = fit.biv, scores = scores.biv, zcor.coef = zcor.coef, ZCCC = zccc, log.tdi = l.tdi, sd.univ.1 = sqrt(var.1), sd.univ.2 = sqrt(var.2)) return(results) } ################################################ # Function for computing Fisher's z-transformation of CCC # Arguments: # mean.iff = mean difference # var.1, var.2, cov.12 = variances and covariance of measurements # Values: # Fisher's z-transformation of ccc ZCCC = function(mean.diff, var.1, var.2, cov.12) { ccc = 2*cov.12/(mean.diff^2 + var.1 + var.2) return(atanh(ccc)) } ################################################ # Function for computing log(TDI) # Arguments: # mean.diff = mean difference # diff.sd = sd difference # prob = probability # l.lim, u.lim = range in which the true value lies # # Values: # log(tdi) log_tdi = function(mean.diff, diff.sd, prob, l.lim = 0, u.lim = 200) { cdf.abs.diff = function(q, mean.diff, diff.sd, prob) { result = pnorm((q - mean.diff)/diff.sd) - pnorm((-q - mean.diff)/diff.sd) return(result - prob) } result = uniroot(cdf.abs.diff, c(l.lim, u.lim), extendInt = "upX", mean.diff = mean.diff, diff.sd = diff.sd, prob = prob)$root return(log(result)) } ############################################ # Function for computing simultaneous confidence bounds or intervals for a # vector of model parameters or functions assuming given bias and variance matrix # of estimates # Arguments: # param.est = estimated model parameter vector # bias.est = estimated bias of param.hat # var.hat = estimated covariance matrix of param.hat # tail = "two.sided" or "less" (upper bound) or "greater" (lower bound) # alpha = 1-level of confidence # # Value: # A data frame containing corresponding simultaneous confidence bounds or intervals simul.bootstrap.ci = function(param.est, bias.est, var.hat, tail = c("lower", "two.sided", "upper"), alpha) { L = dim(var.hat)[1] param.sd = sqrt(diag(var.hat)) Tmu = rep(0, L) Tvar = cov2cor(var.hat) nsim = 50000 # two sided confidence interval if(tail == "two.sided") { Tobs = replicate(nsim, max(abs(mvrnorm(n = 1, mu = Tmu, Sigma = Tvar)))) calpha = quantile(Tobs, 1 - alpha, type = 1) upper.ci = param.est - bias.est + calpha * param.sd lower.ci = param.est - bias.est - calpha * param.sd } # upper confidence interval if(tail == "upper") { Tobs = replicate(nsim, max(mvrnorm(n = 1, mu = Tmu, Sigma = Tvar))) calpha = quantile(Tobs, 1 - alpha, type = 1) lower.ci = param.est - bias.est - calpha * param.sd upper.ci = rep(Inf, L) } # lower confidence interval if(tail == "lower") { Tobs = replicate(nsim, max(mvrnorm(n = 1, mu = Tmu, Sigma = Tvar))) calpha = quantile(Tobs, 1 - alpha, type = 1) lower.ci = rep(-Inf, L) upper.ci = param.est + calpha * param.sd } result = cbind(param.est = param.est, lcl = lower.ci, ucl = upper.ci, calpha = calpha) } ############################################ # Function for bootstrap calculations for parameter estimates and getting # simultaneous confidence bounds/intervals of them # # Arguments: # nboot = number of bootstraps # data.1, data.2 = bivaraite data # x0 = argument values for gam model # ncol = number of points in the union grid # nrow = number of subjects # nbasis = number of basis # pve.univ = proportion of variability # pve.biv = proportion of total variability # alpha = significance level # prob = probability value # makePD = compute the nearest positive definite matrix # # Value: # A data frame containing corresponding simultaneous confidence bounds or intervals for MPACE and UPACE method bootstrap = function(nboot, data.1, data.2, x0, ncol, nrow, nbasis = 10, pve.univ = 0.99, pve.biv = 0.99, alpha = 0.05, prob = 0.90, makePD = FALSE) { dataset = cbind(data.1, data.2) mu1.star.univ = matrix(data = NA, nrow = ncol, ncol = nboot) mu2.star.univ = matrix(data = NA, nrow = ncol, ncol = nboot) mean.diff.star.univ = matrix(data = NA, nrow = ncol, ncol = nboot) ltau.sq1.star.univ = c() ltau.sq2.star.univ = c() ltau.sq.ratio.star.univ = c() zcor.coef.star.univ = matrix(data = NA, nrow = ncol, ncol = nboot) zccc.star.univ = matrix(data = NA, nrow = ncol, ncol = nboot) ltdi.star.univ = matrix(data = NA, nrow = ncol, ncol = nboot) mu1.star.biv = matrix(data = NA, nrow = ncol, ncol = nboot) mu2.star.biv = matrix(data = NA, nrow = ncol, ncol = nboot) mean.diff.star.biv = matrix(data = NA, nrow = ncol, ncol = nboot) ltau.sq1.star.biv = c() ltau.sq2.star.biv = c() ltau.sq.ratio.star.biv = c() zcor.coef.star.biv = matrix(data = NA, nrow = ncol, ncol = nboot) zccc.star.biv = matrix(data = NA, nrow = ncol, ncol = nboot) ltdi.star.biv = matrix(data = NA, nrow = ncol, ncol = nboot) for(i in 1:nboot) { boot.sample = dataset[sample(1:nrow, nrow, replace = T), ] sample.1 = boot.sample[ ,1:ncol] sample.2 = boot.sample[ ,(ncol + 1):(2*ncol)] # Bootstrap estimates for univariate PACE method est.univ = UPACE(sample.1, sample.2, x0, ncol, nrow, nbasis, pve.univ, pve.biv, prob, makePD) mu1.star.univ[ , i] = est.univ$mean[1:ncol] mu2.star.univ[ , i] = est.univ$mean[(ncol + 1):(2*ncol)] mean.diff.star.univ[ , i] = est.univ$mean.diff ltau.sq1.star.univ[i] = est.univ$ltau.sq[1] ltau.sq2.star.univ[i] = est.univ$ltau.sq[2] ltau.sq.ratio.star.univ[i] = est.univ$ltau.sq.ratio zcor.coef.star.univ[ , i] = est.univ$zcor.coef zccc.star.univ[ , i] = est.univ$ZCCC ltdi.star.univ[ , i] = est.univ$log.tdi # Bivariate bootstrap estimates est.bi = MPACE(boot.sample, x0, ncol, nrow, nbasis, pve.biv, prob, makePD) mu1.star.biv[ , i] = est.bi$mean[1:ncol] mu2.star.biv[ , i] = est.bi$mean[(ncol + 1):(2*ncol)] mean.diff.star.biv[ , i] = est.bi$mean.diff ltau.sq1.star.biv[i] = est.bi$ltau.sq[1] ltau.sq2.star.biv[i] = est.bi$ltau.sq[2] ltau.sq.ratio.star.biv[i] = est.bi$ltau.sq.ratio zcor.coef.star.biv[ , i] = est.bi$zcor.coef zccc.star.biv[ , i] = est.bi$ZCCC ltdi.star.biv[ , i] = est.bi$log.tdi } # Bias estimates for univariate method Uni.est = UPACE(data.1, data.2, x0, ncol, nrow, nbasis, pve.univ, pve.biv, prob, makePD) mu1.bias.univ = rowMeans(mu1.star.univ) - Uni.est$mean[1:ncol] mu1.var.univ = var(t(mu1.star.univ)) mu2.bias.univ = rowMeans(mu2.star.univ) - Uni.est$mean[(ncol + 1):(2*ncol)] mu2.var.univ = var(t(mu2.star.univ)) mean.diff.bias.univ = rowMeans(mean.diff.star.univ) - Uni.est$mean.diff mean.diff.var.univ = var(t(mean.diff.star.univ)) ltau.sq1.bias.univ = mean(ltau.sq1.star.univ) - Uni.est$ltau.sq[1] ltau.sq1.var.univ = var(ltau.sq1.star.univ) ltau.sq2.bias.univ = mean(ltau.sq2.star.univ) - Uni.est$ltau.sq[2] ltau.sq2.var.univ = var(ltau.sq2.star.univ) ltau.sq.ratio.bias.univ = mean(ltau.sq.ratio.star.univ) - Uni.est$ltau.sq.ratio ltau.sq.ratio.var.univ = var(ltau.sq.ratio.star.univ) zcor.coef.bias.univ = rowMeans(zcor.coef.star.univ) - Uni.est$zcor.coef zcor.coef.var.univ = var(t(zcor.coef.star.univ)) zccc.bias.univ = rowMeans(zccc.star.univ) - Uni.est$ZCCC zccc.var.univ = var(t(zccc.star.univ)) ltdi.bias.univ = rowMeans(ltdi.star.univ) - Uni.est$log.tdi ltdi.var.univ = var(t(ltdi.star.univ)) # Bias estimates for Bivariate method Bi.est = MPACE(dataset, x0, ncol, nrow, nbasis, pve.biv, prob, makePD) mu1.bias.biv = rowMeans(mu1.star.biv) - Bi.est$mean[1:ncol] mu1.var.biv = var(t(mu1.star.biv)) mu2.bias.biv = rowMeans(mu2.star.biv) - Bi.est$mean[(ncol + 1):(2*ncol)] mu2.var.biv = var(t(mu2.star.biv)) mean.diff.bias.biv = rowMeans(mean.diff.star.biv) - Bi.est$mean.diff mean.diff.var.biv = var(t(mean.diff.star.biv)) ltau.sq1.bias.biv = mean(ltau.sq1.star.biv) - Bi.est$ltau.sq[1] ltau.sq1.var.biv = var(ltau.sq1.star.biv) ltau.sq2.bias.biv = mean(ltau.sq2.star.biv) - Bi.est$ltau.sq[2] ltau.sq2.var.biv = var(ltau.sq2.star.biv) ltau.sq.ratio.bias.biv = mean(ltau.sq.ratio.star.biv) - Bi.est$ltau.sq.ratio ltau.sq.ratio.var.biv = var(ltau.sq.ratio.star.biv) zcor.coef.bias.biv = rowMeans(zcor.coef.star.biv) - Bi.est$zcor.coef zcor.coef.var.biv = var(t(zcor.coef.star.biv)) zccc.bias.biv = rowMeans(zccc.star.biv) - Bi.est$ZCCC zccc.var.biv = var(t(zccc.star.biv)) ltdi.bias.biv = rowMeans(ltdi.star.biv) - Bi.est$log.tdi ltdi.var.biv = var(t(ltdi.star.biv)) # Bootstrap confidence intervals ltau.sq1.univ.ci = (Uni.est$ltau.sq[1] - ltau.sq1.bias.univ) + c(-1, 1) * qnorm(1 - alpha/2) * sqrt(ltau.sq1.var.univ) ltau.sq2.univ.ci = (Uni.est$ltau.sq[2] - ltau.sq2.bias.univ) + c(-1, 1) * qnorm(1 - alpha/2) * sqrt(ltau.sq2.var.univ) ltau.sq.ratio.univ.ci = (Uni.est$ltau.sq.ratio - ltau.sq.ratio.bias.univ) + c(-1, 1) * qnorm(1 - alpha/2) * sqrt(ltau.sq.ratio.var.univ) ltau.sq1.biv.ci = (Bi.est$ltau.sq[1] - ltau.sq1.bias.biv) + c(-1, 1) * qnorm(1 - alpha/2) * sqrt(ltau.sq1.var.biv) ltau.sq2.biv.ci = (Bi.est$ltau.sq[2] - ltau.sq2.bias.biv) + c(-1, 1) * qnorm(1 - alpha/2) * sqrt(ltau.sq2.var.biv) ltau.sq.ratio.biv.ci = (Bi.est$ltau.sq.ratio - ltau.sq.ratio.bias.biv) + c(-1, 1) * qnorm(1 - alpha/2) * sqrt(ltau.sq.ratio.var.biv) mu1.univ.pci = data.frame(lcl = (Uni.est$mean[1:ncol] - mu1.bias.univ) - qnorm(1 - alpha/2) * sqrt(diag(mu1.var.univ)), ucl = (Uni.est$mean[1:ncol] - mu1.bias.univ) + qnorm(1 - alpha/2) * sqrt(diag(mu1.var.univ))) mu2.univ.pci = data.frame(lcl = (Uni.est$mean[(ncol+1):(2*ncol)] - mu2.bias.univ) - qnorm(1 - alpha/2) * sqrt(diag(mu2.var.univ)), ucl = (Uni.est$mean[(ncol+1):(2*ncol)] - mu2.bias.univ) + qnorm(1 - alpha/2) * sqrt(diag(mu2.var.univ))) mean.diff.univ.pci = data.frame(lcl = (Uni.est$mean.diff - mean.diff.bias.univ) - qnorm(1 - alpha/2) * sqrt(diag(mean.diff.var.univ)), ucl = (Uni.est$mean.diff - mean.diff.bias.univ) + qnorm(1 - alpha/2) * sqrt(diag(mean.diff.var.univ))) zcor.coef.univ.pci = data.frame(lcl = (Uni.est$zcor.coef - zcor.coef.bias.univ) - qnorm(1 - alpha/2) * sqrt(diag(zcor.coef.var.univ)), ucl = (Uni.est$zcor.coef - zcor.coef.bias.univ) + qnorm(1 - alpha/2) * sqrt(diag(zcor.coef.var.univ))) zccc.univ.pci = data.frame(lcl = (Uni.est$ZCCC - zccc.bias.univ) - qnorm(1 - alpha) * sqrt(diag(zccc.var.univ))) ltdi.univ.pci = data.frame(ucl = (Uni.est$log.tdi) + qnorm(1 - alpha) * sqrt(diag(ltdi.var.univ))) mu1.biv.pci = data.frame(lcl = (Bi.est$mean[1:ncol] - mu1.bias.biv) - qnorm(1 - alpha/2) * sqrt(diag(mu1.var.biv)), ucl = (Bi.est$mean[1:ncol] - mu1.bias.biv) + qnorm(1 - alpha/2) * sqrt(diag(mu1.var.biv))) mu2.biv.pci = data.frame(lcl = (Bi.est$mean[(ncol+1):(2*ncol)] - mu2.bias.biv) - qnorm(1 - alpha/2) * sqrt(diag(mu2.var.biv)), ucl = (Bi.est$mean[(ncol+1):(2*ncol)] - mu2.bias.biv) + qnorm(1 - alpha/2) * sqrt(diag(mu2.var.biv))) mean.diff.biv.pci = data.frame(lcl = (Bi.est$mean.diff - mean.diff.bias.biv) - qnorm(1 - alpha/2) * sqrt(diag(mean.diff.var.biv)), ucl = (Bi.est$mean.diff - mean.diff.bias.biv) + qnorm(1 - alpha/2) * sqrt(diag(mean.diff.var.biv))) zcor.coef.biv.pci = data.frame(lcl = (Bi.est$zcor.coef - zcor.coef.bias.biv) - qnorm(1 - alpha/2) * sqrt(diag(zcor.coef.var.biv)), ucl = (Bi.est$zcor.coef - zcor.coef.bias.biv) + qnorm(1 - alpha/2) * sqrt(diag(zcor.coef.var.biv))) zccc.biv.pci = data.frame(lcl = (Bi.est$ZCCC - zccc.bias.biv) - qnorm(1 - alpha) * sqrt(diag(zccc.var.biv))) ltdi.biv.pci = data.frame(ucl = (Bi.est$log.tdi) + qnorm(1 - alpha) * sqrt(diag(ltdi.var.biv))) # Simultaneous Bootstrap CI mu1.univ.ci = simul.bootstrap.ci(Uni.est$mean[1:ncol], mu1.bias.univ, mu1.var.univ, "two.sided", alpha) mu2.univ.ci = simul.bootstrap.ci(Uni.est$mean[(ncol+1):(2*ncol)], mu2.bias.univ, mu2.var.univ, "two.sided", alpha) mean.diff.univ.ci = simul.bootstrap.ci(Uni.est$mean.diff, mean.diff.bias.univ, mean.diff.var.univ, "two.sided", alpha) zcor.coef.univ.ci = simul.bootstrap.ci(Uni.est$zcor.coef, zcor.coef.bias.univ, zcor.coef.var.univ, "two.sided", alpha) zccc.univ.ci = simul.bootstrap.ci(Uni.est$ZCCC, zccc.bias.univ, zccc.var.univ, "upper", alpha) ltdi.univ.ci = simul.bootstrap.ci(Uni.est$log.tdi, ltdi.bias.univ, ltdi.var.univ, "lower", alpha) mu1.biv.ci = simul.bootstrap.ci(Bi.est$mean[1:ncol], mu1.bias.biv, mu1.var.biv, "two.sided", alpha) mu2.biv.ci = simul.bootstrap.ci(Bi.est$mean[(ncol+1):(2*ncol)], mu2.bias.biv, mu2.var.biv, "two.sided", alpha) mean.diff.biv.ci = simul.bootstrap.ci(Bi.est$mean.diff, mean.diff.bias.biv, mean.diff.var.biv, "two.sided", alpha) zcor.coef.biv.ci = simul.bootstrap.ci(Bi.est$zcor.coef, zcor.coef.bias.biv, zcor.coef.var.biv, "two.sided", alpha) zccc.biv.ci = simul.bootstrap.ci(Bi.est$ZCCC, zccc.bias.biv, zccc.var.biv, "upper", alpha) ltdi.biv.ci = simul.bootstrap.ci(Bi.est$log.tdi, ltdi.bias.biv, ltdi.var.biv, "lower", alpha) results = list(mu1.univ.pci = mu1.univ.pci, mu2.univ.pci = mu2.univ.pci, mu1.univ.ci = mu1.univ.ci, mu2.univ.ci = mu2.univ.ci, mean.diff.univ.pci = mean.diff.univ.pci, mean.diff.univ.ci = mean.diff.univ.ci, ltau.sq1.univ.ci = ltau.sq1.univ.ci, ltau.sq2.univ.ci = ltau.sq2.univ.ci, ltau.sq.ratio.univ.ci = ltau.sq.ratio.univ.ci, zcor.coef.univ.pci = zcor.coef.univ.pci, zccc.univ.pci = zccc.univ.pci, ltdi.univ.pci = ltdi.univ.pci, zcor.coef.univ.ci = zcor.coef.univ.ci, zccc.univ.ci = zccc.univ.ci, ltdi.univ.ci = ltdi.univ.ci, mu1.biv.pci = mu1.biv.pci, mu2.biv.pci = mu2.biv.pci, mu1.biv.ci = mu1.biv.ci, mu2.biv.ci = mu2.biv.ci, mean.diff.biv.pci = mean.diff.biv.pci, mean.diff.biv.ci = mean.diff.biv.ci, ltau.sq1.biv.ci = ltau.sq1.biv.ci, ltau.sq2.biv.ci = ltau.sq2.biv.ci, ltau.sq.ratio.biv.ci = ltau.sq.ratio.biv.ci, zcor.coef.biv.pci = zcor.coef.biv.pci, zccc.biv.pci = zccc.biv.pci, ltdi.biv.pci = ltdi.biv.pci, zcor.coef.biv.ci = zcor.coef.biv.ci, zccc.biv.ci = zccc.biv.ci, ltdi.biv.ci = ltdi.biv.ci) return(results) } # End of file #################################################################################