Nothing
HSROCSummary <-
function (data, burn_in = 0, iter.keep = NULL, Thin = 1, sub_rs = NULL,
point_estimate = c("median", "mean"), summary.path = getwd(),
chain = getwd(), tv = NULL, digit = 6, print_plot = FALSE,
plot.ind.studies = TRUE, cred_region = TRUE, predict_region = TRUE,
col.pooled.estimate = "red", col.predict.region = "blue",
lty.cred.region = "dotdash", lty.predict.region = "dotted",
region_level = 0.95, trunc_low = 0.025, trunc_up = 0.025)
{
setwd(summary.path)
if (missing(data))
stop("You must provide a valid 'data' argument", call. = FALSE)
N = length(data[, 1])
if (burn_in < 0) {
cat("The 'burn_in' argument must be greater or equal than zero. \n")
stop("Please respecify and call HSROCSummary() again.\n")
}
if (is.null(iter.keep) == FALSE) {
if (iter.keep < 0) {
cat("The 'iter.keep' argument must be greater or equal than zero. \n")
stop("Please respecify and call HSROCSummary() again.\n")
}
}
if (Thin < 1) {
cat("The 'Thin' argument must be greater or equal than 1. \n")
stop("Please respecify and call HSROCSummary() again.\n")
}
if (is.null(sub_rs) == TRUE) {
sub_rs = list(1, 1:N)
}
if (sub_rs[[1]] != (length(sub_rs) - 1)) {
cat(paste("The value of the first element of 'sub_rs' (sub_rs[[1]] = ",
sub_rs[[1]], " ) does not match the number of remaining elements (length(sub_rs[[2:",
length(sub_rs), "]])) = ", length(2:length(sub_rs)),
"\n", sep = ""))
stop("Please respecify and call HSROCSummary() again.\n")
}
if (is.logical(print_plot) == FALSE) {
cat("The 'print_plot' argument must be a logical object. \n")
stop("Please respecify and call HSROCSummary() again.\n")
}
file.prior = "Prior.information.txt"
prior_dist_PI = "beta"
if (is.null(chain) == FALSE) {
setwd(chain[[1]])
nb_chains = length(chain)
}
else {
nb_chains = 1
}
point_estimate = match.arg(point_estimate)
prior.exist = file.exists(file.prior)
if (prior.exist == FALSE) {
stop(paste("Make sure the \"", file.prior, "\" file created by the 'HSROC' function is in the \"",
summary.path, "\" working directory. \n", sep = ""))
}
prior = read.table(file.prior, header = TRUE)
model = read.table("model.txt", header = FALSE)
Gold_se = (read.table("S2.txt", header = FALSE) == 1)
Gold_sp = (read.table("C2.txt", header = FALSE) == 1)
prior_sig_t = read.table("Prior on sigma_theta.txt", header = FALSE)
prior_sig_a = read.table("Prior on sigma_alpha.txt", header = FALSE)
if (length(prior[, 1]) == 7) {
Gold_Std = TRUE
condInd = TRUE
}
else {
if (length(prior[, 1]) > 7 & length(prior[, 1]) <= 7 +
2 * sub_rs[[1]]) {
Gold_Std = FALSE
condInd = TRUE
}
else {
if (length(prior[, 1]) > 7 + 2 * sub_rs[[1]]) {
Gold_Std = FALSE
condInd = FALSE
}
}
}
if (is.null(tv) == FALSE) {
real_life = FALSE
if (sum(dim(tv[[1]])) != N + 5) {
cat(paste("The true value for the within-study parameters were misspecified. Make sure the ordering described in the help file is preserved. \n"))
stop("Please respecify and call HSROCSummary() again.\n")
}
if (length(tv[[2]]) != 7) {
cat(paste("The true value for the between-study parameters were misspecified. Make sure the ordering described in the help file is preserved. \n"))
stop("Please respecify and call HSROCSummary() again.\n")
}
if (Gold_Std == FALSE) {
if (sum(dim(tv[[3]])) != sub_rs[[1]] + 2) {
cat(paste("The true value for the test under evaluation were misspecified. Make sure the ordering described in the help file is preserved. \n"))
stop("Please respecify and call HSROCSummary() again.\n")
}
}
}
else {
real_life = TRUE
}
beta.a = prior[1, 1]
beta.b = prior[1, 2]
prior.THETA.lower = prior[2, 1]
prior.THETA.upper = prior[2, 2]
prior.LAMBDA.lower = prior[3, 1]
prior.LAMBDA.upper = prior[3, 2]
l.disp.alpha = prior[4, 1]
u.disp.alpha = prior[4, 2]
l.disp.theta = prior[5, 1]
u.disp.theta = prior[5, 2]
low.pi = prior[6, 1]
up.pi = prior[6, 2]
low.rj = prior[7, 1]
up.rj = prior[7, 2]
SCO = FALSE
rs.length = sub_rs[[1]]
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
low.sp = prior[8:(8 + (rs.length - 1)), 1]
up.sp = prior[8:(8 + (rs.length - 1)), 2]
prior_dist_S2 = NULL
prior_dist_C2 = "beta"
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
low.se = prior[8:(8 + (rs.length - 1)), 1]
up.se = prior[8:(8 + (rs.length - 1)), 2]
prior_dist_C2 = NULL
prior_dist_S2 = "beta"
}
else {
low.se = prior[8:(8 + (rs.length - 1)), 1]
up.se = prior[8:(8 + (rs.length - 1)), 2]
low.sp = prior[(8 + rs.length):(8 + (2 * rs.length -
1)), 1]
up.sp = prior[(8 + rs.length):(8 + (2 * rs.length -
1)), 2]
prior_dist_S2 = "beta"
prior_dist_C2 = "beta"
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model == 2) {
mean.a1 = prior[8:(8 + (rs.length - 1)), 1]
sd.a1 = prior[8:(8 + (rs.length - 1)), 2]
mean.a0 = prior[(8 + rs.length):(8 + (2 * rs.length -
1)), 1]
sd.a0 = prior[(8 + rs.length):(8 + (2 * rs.length -
1)), 2]
}
else {
if (condInd == FALSE) {
low.d1 = prior[8, 1]
up.d1 = prior[8, 2]
low.d0 = prior[9, 1]
up.d0 = prior[9, 2]
mean.a1 = prior[10:(10 + (rs.length - 1)), 1]
sd.a1 = prior[10:(10 + (rs.length - 1)), 2]
mean.a0 = prior[(10 + rs.length):(10 + (2 * rs.length -
1)), 1]
sd.a0 = prior[(10 + rs.length):(10 + (2 * rs.length -
1)), 2]
low.b1 = prior[(10 + (2 * rs.length)):(10 + (3 *
rs.length - 1)), 1]
up.b1 = prior[(10 + (2 * rs.length)):(10 + (3 *
rs.length - 1)), 2]
low.b0 = prior[(10 + (3 * rs.length)):(10 + (4 *
rs.length - 1)), 1]
up.b0 = prior[(10 + (3 * rs.length)):(10 + (4 *
rs.length - 1)), 2]
}
}
}
if (prior_dist_PI == "beta") {
alpha.PI = beta.parameter(low = low.pi, up = up.pi)[1,
]
beta.PI = beta.parameter(low = low.pi, up = up.pi)[2,
]
}
else {
if (prior_dist_PI == "uniform") {
alpha.PI = low.pi
beta.PI = up.pi
}
}
if (model == 1) {
if (Gold_se == TRUE) {
Sens2.alpha = Sens2.beta = NULL
}
else {
Sens2.alpha = beta.parameter(low = low.se, up = up.se)[1,
]
Sens2.beta = beta.parameter(low = low.se, up = up.se)[2,
]
}
}
if (model == 1) {
if (Gold_sp == TRUE) {
Spec2.alpha = Spec2.beta = NULL
}
else {
Spec2.alpha = beta.parameter(low = low.sp, up = up.sp)[1,
]
Spec2.beta = beta.parameter(low = low.sp, up = up.sp)[2,
]
}
}
file.theta = "theta.txt"
file.alpha = "alpha.txt"
file.capital.THETA = "capital_THETA.txt"
file.LAMBDA = "LAMBDA.txt"
file.beta = "beta.txt"
file.PI = "PI.txt"
file.sigma.alpha = "sigma.alpha.txt"
file.sigma.theta = "sigma.theta.txt"
file.Sens2 = "Sens2.txt"
file.Spec2 = "Spec2.txt"
file.Sens1 = "Sens1.txt"
file.Spec1 = "Spec1.txt"
file.result = "estimate.txt"
file.C_overall = "C_overall.txt"
file.S_overall = "S_overall.txt"
file.d1 = "d1.txt"
file.d0 = "d0.txt"
file.a1 = "a1.txt"
file.a0 = "a0.txt"
file.b1 = "b1.txt"
file.b0 = "b0.txt"
numb.iter = scan("iter.txt", quiet = TRUE)
if (is.null(iter.keep) == TRUE) {
iter.num = numb.iter
}
else {
iter.num = iter.keep
}
if ((iter.num - burn_in)/Thin < 100)
stop("You don't have enough iterations to estimate the MC error. After taking into account the \"burn in\" and \"thinning interval\", you need at least 100 iterations to proceed.")
if (is.null(chain) == TRUE) {
N1 = N * iter.num
t1 <- file("alpha.txt", "rb")
t2 <- file("theta.txt", "rb")
t3 <- file("PI.txt", "rb")
t4 <- file("Sens1.txt", "rb")
t5 <- file("Spec1.txt", "rb")
alpha_bin = readBin(t1, double(), n = N1, endian = "little")
theta_bin = readBin(t2, double(), n = N1, endian = "little")
pi_bin = readBin(t3, double(), n = N1, endian = "little")
Sens1_bin = readBin(t4, double(), n = N1, endian = "little")
Spec1_bin = readBin(t5, double(), n = N1, endian = "little")
alpha = matrix(0, ncol = N, nrow = iter.num)
theta = matrix(0, ncol = N, nrow = iter.num)
PI = matrix(0, ncol = N, nrow = iter.num)
S1 = matrix(0, ncol = N, nrow = iter.num)
C1 = matrix(0, ncol = N, nrow = iter.num)
for (i in 1:N) {
sequence = seq(i, iter.num * N, N)
alpha[, i] = alpha_bin[sequence]
theta[, i] = theta_bin[sequence]
PI[, i] = pi_bin[sequence]
S1[, i] = Sens1_bin[sequence]
C1[, i] = Spec1_bin[sequence]
}
close(t1)
close(t2)
close(t3)
close(t4)
close(t5)
t6 <- file("capital_THETA.txt", "rb")
THETA = readBin(t6, double(), n = iter.num, endian = "little")
close(t6)
t7 <- file("LAMBDA.txt", "rb")
LAMBDA = readBin(t7, double(), n = iter.num, endian = "little")
close(t7)
t8 <- file("beta.txt", "rb")
beta = readBin(t8, double(), n = iter.num, endian = "little")
close(t8)
t9 <- file("Sens1_new.txt", "rb")
S1_new = readBin(t9, double(), n = iter.num, endian = "little")
close(t9)
t10 <- file("Spec1_new.txt", "rb")
C1_new = readBin(t10, double(), n = iter.num, endian = "little")
close(t10)
t11 <- file("sigma.alpha.txt", "rb")
sigma.alpha = readBin(t11, double(), n = iter.num, endian = "little")
close(t11)
t12 <- file("sigma.theta.txt", "rb")
sigma.theta = readBin(t12, double(), n = iter.num, endian = "little")
close(t12)
t14 <- file("S_overall.txt", "rb")
S_overall = readBin(t14, double(), n = iter.num, endian = "little")
close(t14)
t15 <- file("C_overall.txt", "rb")
C_overall = readBin(t15, double(), n = iter.num, endian = "little")
close(t15)
total = iter.num
q = burn_in
alpha = alpha[(q + 1):total, ]
THETA = THETA[(q + 1):total]
LAMBDA = LAMBDA[(q + 1):total]
beta = beta[(q + 1):total]
PI = PI[(q + 1):total, ]
sigma.alpha = sigma.alpha[(q + 1):total]
S1 = S1[(q + 1):total, ]
C1 = C1[(q + 1):total, ]
S1_new = S1_new[(q + 1):total]
C1_new = C1_new[(q + 1):total]
C_overall = C_overall[(q + 1):total]
S_overall = S_overall[(q + 1):total]
theta = theta[(q + 1):total, ]
sigma.theta = sigma.theta[(q + 1):total]
taille = length((q + 1):total)
thin.interval = Thin
thin = seq(1, taille, by = thin.interval)
alpha = as.matrix(alpha)
alpha = alpha[thin, ]
THETA = THETA[thin]
LAMBDA = LAMBDA[thin]
beta = beta[thin]
PI = as.matrix(PI)
PI = PI[thin, ]
sigma.alpha = sigma.alpha[thin]
S1 = as.matrix(S1)
S1 = S1[thin, ]
C1 = as.matrix(C1)
C1 = C1[thin, ]
S1_new = S1_new[thin]
C1_new = C1_new[thin]
C_overall = C_overall[thin]
S_overall = S_overall[thin]
theta = as.matrix(theta)
theta = theta[thin, ]
sigma.theta = sigma.theta[thin]
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
C2 = read.table(file.Spec2)
C2 = C2[q:total, ]
C2 = as.matrix(C2)
C2 = C2[thin, ]
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
S2 = read.table(file.Sens2)
S2 = S2[(q + 1):total, ]
S2 = as.matrix(S2)
S2 = S2[thin, ]
}
else {
N2 = rs.length * iter.num
t16 <- file("Sens2.txt", "rb")
t17 <- file("Spec2.txt", "rb")
Sens2_bin = readBin(t16, double(), n = N1,
endian = "little")
Spec2_bin = readBin(t17, double(), n = N1,
endian = "little")
S2 = matrix(0, ncol = rs.length, nrow = iter.num)
C2 = matrix(0, ncol = rs.length, nrow = iter.num)
for (j in 1:rs.length) {
sequence = seq(j, iter.num * rs.length, rs.length)
S2[, j] = Sens2_bin[sequence]
C2[, j] = Spec2_bin[sequence]
}
close(t16)
close(t17)
S2 = S2[(q + 1):total, ]
C2 = C2[(q + 1):total, ]
S2 = as.matrix(S2)
C2 = as.matrix(C2)
S2 = S2[thin, ]
C2 = C2[thin, ]
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
a1 = read.table(file.a1)
a0 = read.table(file.a0)
S2 = read.table(file.Sens2)
C2 = read.table(file.Spec2)
a1 = a1[(q + 1):total, ]
a0 = a0[(q + 1):total, ]
S2 = S2[(q + 1):total, ]
C2 = C2[(q + 1):total, ]
a1 = as.matrix(a1)
a0 = as.matrix(a0)
a1 = a1[thin, ]
a0 = a0[thin, ]
S2 = as.matrix(S2)
C2 = as.matrix(C2)
S2 = S2[thin, ]
C2 = C2[thin, ]
}
else {
if (condInd == FALSE) {
d1 = read.table(file.d1)
d0 = read.table(file.d0)
a1 = read.table(file.a1)
a0 = read.table(file.a0)
b1 = read.table(file.b1)
b0 = read.table(file.b0)
S2 = read.table(file.Sens2)
C2 = read.table(file.Spec2)
d1 = d1[(q + 1):total, ]
d0 = d0[(q + 1):total, ]
a1 = a1[(q + 1):total, ]
a0 = a0[(q + 1):total, ]
b1 = b1[(q + 1):total, ]
b0 = b0[(q + 1):total, ]
S2 = S2[(q + 1):total, ]
C2 = C2[(q + 1):total, ]
d1 = as.matrix(d1)
d0 = as.matrix(d0)
a1 = as.matrix(a1)
a0 = as.matrix(a0)
b1 = as.matrix(b1)
b0 = as.matrix(b0)
S2 = as.matrix(S2)
C2 = as.matrix(C2)
d1 = d1[thin, ]
d0 = d0[thin, ]
a1 = a1[thin, ]
a0 = a0[thin, ]
b1 = b1[thin, ]
b0 = b0[thin, ]
S2 = S2[thin, ]
C2 = C2[thin, ]
}
}
}
}
else {
if (is.null(chain) == FALSE) {
K = length(chain)
theta = alpha = THETA = LAMBDA = beta = PI = sigma.alpha = sigma.theta = S1 = C1 = S1_new = C1_new = S2 = C2 = S_overall = C_overall = a1 = a0 = b1 = b0 = d1 = d0 = numeric()
for (k in 1:K) {
setwd(chain[[k]])
N1 = N * iter.num
t1 <- file("alpha.txt", "rb")
t2 <- file("theta.txt", "rb")
t3 <- file("PI.txt", "rb")
t4 <- file("Sens1.txt", "rb")
t5 <- file("Spec1.txt", "rb")
alpha_bin = readBin(t1, double(), n = N1, endian = "little")
theta_bin = readBin(t2, double(), n = N1, endian = "little")
pi_bin = readBin(t3, double(), n = N1, endian = "little")
Sens1_bin = readBin(t4, double(), n = N1, endian = "little")
Spec1_bin = readBin(t5, double(), n = N1, endian = "little")
a = matrix(0, ncol = N, nrow = iter.num)
t = matrix(0, ncol = N, nrow = iter.num)
p = matrix(0, ncol = N, nrow = iter.num)
S.1 = matrix(0, ncol = N, nrow = iter.num)
C.1 = matrix(0, ncol = N, nrow = iter.num)
for (i in 1:N) {
sequence = seq(i, iter.num * N, N)
a[, i] = alpha_bin[sequence]
t[, i] = theta_bin[sequence]
p[, i] = pi_bin[sequence]
S.1[, i] = Sens1_bin[sequence]
C.1[, i] = Spec1_bin[sequence]
}
close(t1)
close(t2)
close(t3)
close(t4)
close(t5)
t6 <- file("capital_THETA.txt", "rb")
T = readBin(t6, double(), n = iter.num, endian = "little")
close(t6)
t7 <- file("LAMBDA.txt", "rb")
L = readBin(t7, double(), n = iter.num, endian = "little")
close(t7)
t8 <- file("beta.txt", "rb")
b = readBin(t8, double(), n = iter.num, endian = "little")
close(t8)
t9 <- file("Sens1_new.txt", "rb")
S.1_new = readBin(t9, double(), n = iter.num,
endian = "little")
close(t9)
t10 <- file("Spec1_new.txt", "rb")
C.1_new = readBin(t10, double(), n = iter.num,
endian = "little")
close(t10)
t11 <- file("sigma.alpha.txt", "rb")
sig.a = readBin(t11, double(), n = iter.num,
endian = "little")
close(t11)
t12 <- file("sigma.theta.txt", "rb")
sig.t = readBin(t12, double(), n = iter.num,
endian = "little")
close(t12)
t14 <- file("S_overall.txt", "rb")
S.ov = readBin(t14, double(), n = iter.num, endian = "little")
close(t14)
t15 <- file("C_overall.txt", "rb")
C.ov = readBin(t15, double(), n = iter.num, endian = "little")
close(t15)
total = iter.num
q = burn_in
a = a[(q + 1):total, ]
T = T[(q + 1):total]
L = L[(q + 1):total]
b = b[(q + 1):total]
p = p[(q + 1):total, ]
sig.a = sig.a[(q + 1):total]
S.1 = S.1[(q + 1):total, ]
C.1 = C.1[(q + 1):total, ]
S.1_new = S.1_new[(q + 1):total]
C.1_new = C.1_new[(q + 1):total]
C.ov = C.ov[(q + 1):total]
S.ov = S.ov[(q + 1):total]
t = t[(q + 1):total, ]
sig.t = sig.t[(q + 1):total]
taille = length((q + 1):total)
thin.interval = Thin
thin = seq(1, taille, by = thin.interval)
a = as.matrix(a)
a = a[thin, ]
T = T[thin]
L = L[thin]
b = b[thin]
p = as.matrix(p)
p = p[thin, ]
sig.a = sig.a[thin]
S.1 = as.matrix(S.1)
S.1 = S.1[thin, ]
C.1 = as.matrix(C.1)
C.1 = C.1[thin, ]
S.1_new = S.1_new[thin]
C.1_new = C.1_new[thin]
C.ov = C.ov[thin]
S.ov = S.ov[thin]
t = as.matrix(t)
t = t[thin, ]
sig.t = sig.t[thin]
alpha = rbind(alpha, as.matrix(a))
THETA = rbind(THETA, as.matrix(T))
LAMBDA = rbind(LAMBDA, as.matrix(L))
beta = rbind(beta, as.matrix(b))
PI = rbind(PI, as.matrix(p))
sigma.alpha = rbind(sigma.alpha, as.matrix(sig.a))
S1 = rbind(S1, as.matrix(S.1))
C1 = rbind(C1, as.matrix(C.1))
S1_new = rbind(S1_new, as.matrix(S.1_new))
C1_new = rbind(C1_new, as.matrix(C.1_new))
C_overall = rbind(C_overall, as.matrix(C.ov))
S_overall = rbind(S_overall, as.matrix(S.ov))
theta = rbind(theta, as.matrix(t))
sigma.theta = rbind(sigma.theta, as.matrix(sig.t))
if (condInd == TRUE & Gold_Std == FALSE & model ==
1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
C.2 = read.table(file.Spec2)
C.2 = C.2[(q + 1):total, ]
C.2 = as.matrix(C.2)
C.2 = C.2[thin, ]
C2 = rbind(C2, as.matrix(C.2))
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
S.2 = read.table(file.Sens2)
S.2 = S.2[(q + 1):total, ]
S.2 = as.matrix(S.2)
S.2 = S.2[thin, ]
S2 = rbind(S2, as.matrix(S.2))
}
else {
N2 = rs.length * iter.num
t16 <- file("Sens2.txt", "rb")
t17 <- file("Spec2.txt", "rb")
Sens2_bin = readBin(t16, double(), n = N1,
endian = "little")
Spec2_bin = readBin(t17, double(), n = N1,
endian = "little")
S.2 = matrix(0, ncol = rs.length, nrow = iter.num)
C.2 = matrix(0, ncol = rs.length, nrow = iter.num)
for (i in 1:rs.length) {
sequence = seq(i, iter.num * rs.length,
rs.length)
S.2[, i] = Sens2_bin[sequence]
C.2[, i] = Spec2_bin[sequence]
}
close(t16)
close(t17)
S.2 = S.2[(q + 1):total, ]
C.2 = C.2[(q + 1):total, ]
S.2 = as.matrix(S.2)
C.2 = as.matrix(C.2)
S.2 = S.2[thin, ]
C.2 = C.2[thin, ]
S2 = rbind(S2, as.matrix(S.2))
C2 = rbind(C2, as.matrix(C.2))
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
a.1 = read.table(file.a1)
a.0 = read.table(file.a0)
S.2 = read.table(file.Sens2)
C.2 = read.table(file.Spec2)
a.1 = a.1[(q + 1):total, ]
a.0 = a.0[(q + 1):total, ]
S.2 = S.2[(q + 1):total, ]
C.2 = C.2[(q + 1):total, ]
a.1 = as.matrix(a.1)
a.0 = as.matrix(a.0)
a.1 = a.1[thin, ]
a.0 = a.0[thin, ]
S.2 = as.matrix(S.2)
C.2 = as.matrix(C.2)
S.2 = S.2[thin, ]
C.2 = C.2[thin, ]
a1 = rbind(a1, as.matrix(a.1))
a0 = rbind(a0, as.matrix(a.0))
S2 = rbind(S2, as.matrix(S.2))
C2 = rbind(C2, as.matrix(C.2))
}
else {
if (condInd == FALSE) {
d.1 = read.table(file.d1)
d.0 = read.table(file.d0)
a.1 = read.table(file.a1)
a.0 = read.table(file.a0)
b.1 = read.table(file.b1)
b.0 = read.table(file.b0)
S.2 = read.table(file.Sens2)
C.2 = read.table(file.Spec2)
d.1 = d.1[(q + 1):total, ]
d.0 = d.0[(q + 1):total, ]
a.1 = a.1[(q + 1):total, ]
a.0 = a.0[(q + 1):total, ]
b.1 = b.1[(q + 1):total, ]
b.0 = b.0[(q + 1):total, ]
S.2 = S.2[(q + 1):total, ]
C.2 = C.2[(q + 1):total, ]
d.1 = as.matrix(d.1)
d.0 = as.matrix(d.0)
a.1 = as.matrix(a.1)
a.0 = as.matrix(a.0)
b.1 = as.matrix(b.1)
b.0 = as.matrix(b.0)
d.1 = d.1[thin, ]
d.0 = d.0[thin, ]
a.1 = a.1[thin, ]
a.0 = a.0[thin, ]
b.1 = b.1[thin, ]
b.0 = b.0[thin, ]
S.2 = as.matrix(S.2)
C.2 = as.matrix(C.2)
S.2 = S.2[thin, ]
C.2 = C.2[thin, ]
a1 = rbind(a1, as.matrix(a.1))
a0 = rbind(a0, as.matrix(a.0))
b1 = rbind(b1, as.matrix(b.1))
b0 = rbind(b0, as.matrix(b.0))
d1 = rbind(d1, as.matrix(d.1))
d0 = rbind(d0, as.matrix(d.0))
S2 = rbind(S2, as.matrix(S.2))
C2 = rbind(C2, as.matrix(C.2))
}
}
}
}
}
}
setwd(summary.path)
alpha = as.mcmc(alpha)
THETA = as.mcmc(THETA)
LAMBDA = as.mcmc(LAMBDA)
beta = as.mcmc(beta)
PI = as.mcmc(PI)
sigma.alpha = as.mcmc(sigma.alpha)
S1 = as.mcmc(S1)
C1 = as.mcmc(C1)
S1_new = as.mcmc(S1_new)
C1_new = as.mcmc(C1_new)
C_overall = as.mcmc(C_overall)
S_overall = as.mcmc(S_overall)
theta = as.mcmc(theta)
sigma.theta = as.mcmc(sigma.theta)
if (point_estimate == "mean") {
mean_OR_med = 2
}
else {
if (point_estimate == "median") {
mean_OR_med = 3
}
}
iter.size = length(beta)
theta.est = apply(as.matrix(theta), 2, point_estimate)
theta.HPD = HPDinterval(theta)
theta.sd = apply(theta, 2, sd)
sigma.theta.est = apply(as.matrix(sigma.theta), 2, point_estimate)
sigma.theta.HPD = HPDinterval(sigma.theta)
sigma.theta.sd = apply(sigma.theta, 2, sd)
alpha.est = apply(as.matrix(alpha), 2, point_estimate)
alpha.HPD = HPDinterval(alpha)
alpha.sd = apply(alpha, 2, sd)
THETA.est = apply(as.matrix(THETA), 2, point_estimate)
THETA.HPD = HPDinterval(THETA)
THETA.sd = apply(THETA, 2, sd)
LAMBDA.est = apply(as.matrix(LAMBDA), 2, point_estimate)
LAMBDA.HPD = HPDinterval(LAMBDA)
LAMBDA.sd = apply(LAMBDA, 2, sd)
beta.est = apply(as.matrix(beta), 2, point_estimate)
beta.HPD = HPDinterval(beta)
beta.sd = apply(beta, 2, sd)
PI.est = apply(as.matrix(PI), 2, point_estimate)
PI.HPD = HPDinterval(PI)
PI.sd = apply(PI, 2, sd)
sigma.alpha.est = apply(as.matrix(sigma.alpha), 2, point_estimate)
sigma.alpha.HPD = HPDinterval(sigma.alpha)
sigma.alpha.sd = apply(sigma.alpha, 2, sd)
S1.est = apply(as.matrix(S1), 2, point_estimate)
S1.HPD = HPDinterval(S1)
S1.sd = apply(S1, 2, sd)
C1.est = apply(as.matrix(C1), 2, point_estimate)
C1.HPD = HPDinterval(C1)
C1.sd = apply(C1, 2, sd)
S1_new.est = apply(as.matrix(S1_new), 2, point_estimate)
S1_new.HPD = HPDinterval(S1_new)
S1_new.sd = apply(S1_new, 2, sd)
C1_new.est = apply(as.matrix(C1_new), 2, point_estimate)
C1_new.HPD = HPDinterval(C1_new)
C1_new.sd = apply(C1_new, 2, sd)
C_overall.est = apply(as.matrix(C_overall), 2, point_estimate)
C_overall.HPD = HPDinterval(C_overall)
C_overall.sd = apply(C_overall, 2, sd)
S_overall.est = apply(as.matrix(S_overall), 2, point_estimate)
S_overall.HPD = HPDinterval(S_overall)
S_overall.sd = apply(S_overall, 2, sd)
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
C2 = as.mcmc(C2)
C2.est = apply(as.matrix(C2), 2, point_estimate)
C2.HPD = HPDinterval(C2)
C2.sd = apply(C2, 2, sd)
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
S2 = as.mcmc(S2)
S2.est = apply(as.matrix(S2), 2, point_estimate)
S2.HPD = HPDinterval(S2)
S2.sd = apply(S2, 2, sd)
}
else {
S2 = as.mcmc(S2)
C2 = as.mcmc(C2)
S2.est = apply(as.matrix(S2), 2, point_estimate)
S2.HPD = HPDinterval(S2)
S2.sd = apply(S2, 2, sd)
C2.est = apply(as.matrix(C2), 2, point_estimate)
C2.HPD = HPDinterval(C2)
C2.sd = apply(C2, 2, sd)
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model == 2) {
a1 = as.mcmc(a1)
a0 = as.mcmc(a0)
S2 = as.mcmc(S2)
C2 = as.mcmc(C2)
a1.est = apply(as.matrix(a1), 2, point_estimate)
a1.HPD = HPDinterval(a1)
a1.sd = apply(a1, 2, sd)
a0.est = apply(as.matrix(a0), 2, point_estimate)
a0.HPD = HPDinterval(a0)
a0.sd = apply(a0, 2, sd)
S2.est = apply(as.matrix(S2), 2, point_estimate)
S2.HPD = HPDinterval(S2)
S2.sd = apply(S2, 2, sd)
C2.est = apply(as.matrix(C2), 2, point_estimate)
C2.HPD = HPDinterval(C2)
C2.sd = apply(C2, 2, sd)
}
else {
if (condInd == FALSE) {
a1 = as.mcmc(a1)
a0 = as.mcmc(a0)
b1 = as.mcmc(b1)
b0 = as.mcmc(b0)
d1 = as.mcmc(d1)
d0 = as.mcmc(d0)
S2 = as.mcmc(S2)
C2 = as.mcmc(C2)
a1.est = apply(as.matrix(a1), 2, point_estimate)
a1.HPD = HPDinterval(a1)
a1.sd = apply(a1, 2, sd)
a0.est = apply(as.matrix(a0), 2, point_estimate)
a0.HPD = HPDinterval(a0)
a0.sd = apply(a0, 2, sd)
b1.est = apply(as.matrix(b1), 2, point_estimate)
b1.HPD = HPDinterval(b1)
b1.sd = apply(b1, 2, sd)
b0.est = apply(as.matrix(b0), 2, point_estimate)
b0.HPD = HPDinterval(b0)
b0.sd = apply(b0, 2, sd)
d1.est = apply(as.matrix(d1), 2, point_estimate)
d1.HPD = HPDinterval(d1)
d1.sd = apply(d1, 2, sd)
d0.est = apply(as.matrix(d0), 2, point_estimate)
d0.HPD = HPDinterval(d0)
d0.sd = apply(d0, 2, sd)
S2.est = apply(as.matrix(S2), 2, point_estimate)
S2.HPD = HPDinterval(S2)
S2.sd = apply(S2, 2, sd)
C2.est = apply(as.matrix(C2), 2, point_estimate)
C2.HPD = HPDinterval(C2)
C2.sd = apply(C2, 2, sd)
}
}
}
batch = 50
ssize = iter.size/batch
moy.t = moy.a = moy.T = moy.L = moy.b = moy.p = moy.sig.a = moy.sig.t = moy.s2 = moy.c2 = moy.s1 = moy.c1 = moy.s1_new = moy.c1_new = moy.s.ov = moy.c.ov = moy.a1 = moy.a0 = moy.b1 = moy.b0 = moy.d1 = moy.d0 = numeric()
if (N == 1) {
for (i in 1:batch) {
moy.a = c(moy.a, mean(alpha[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)])/ssize)
moy.p = c(moy.p, mean(PI[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)])/ssize)
moy.s1 = c(moy.s1, mean(as.matrix(S1[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))/ssize)
moy.c1 = c(moy.c1, mean(as.matrix(C1[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))/ssize)
moy.L = c(moy.L, mean(LAMBDA[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)]))
moy.T = c(moy.T, mean(THETA[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)]))
moy.b = c(moy.b, mean(beta[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)]))
moy.sig.a = c(moy.sig.a, mean(sigma.alpha[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.s.ov = c(moy.s.ov, mean(S_overall[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.c.ov = c(moy.c.ov, mean(C_overall[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.s1_new = c(moy.s1_new, mean(S1_new[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.c1_new = c(moy.c1_new, mean(C1_new[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.t = c(moy.t, mean(theta[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)])/ssize)
moy.sig.t = c(moy.sig.t, mean(sigma.theta[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
}
}
else {
for (i in 1:batch) {
moy.a = cbind(moy.a, colSums(alpha[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0), ])/ssize)
moy.p = cbind(moy.p, colSums(PI[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0), ])/ssize)
moy.s1 = cbind(moy.s1, colSums(as.matrix(S1[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0), ]))/ssize)
moy.c1 = cbind(moy.c1, colSums(as.matrix(C1[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0), ]))/ssize)
moy.L = c(moy.L, mean(LAMBDA[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)]))
moy.T = c(moy.T, mean(THETA[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)]))
moy.b = c(moy.b, mean(beta[round((1 + ssize * (i -
1)), 0):round((ssize * i), 0)]))
moy.sig.a = c(moy.sig.a, mean(sigma.alpha[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.s.ov = c(moy.s.ov, mean(S_overall[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.c.ov = c(moy.c.ov, mean(C_overall[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.s1_new = c(moy.s1_new, mean(S1_new[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.c1_new = c(moy.c1_new, mean(C1_new[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
moy.t = cbind(moy.t, colSums(theta[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0), ])/ssize)
moy.sig.t = c(moy.sig.t, mean(sigma.theta[round((1 +
ssize * (i - 1)), 0):round((ssize * i), 0)]))
}
}
if (N == 1) {
alpha.MCerror = sqrt(sum((moy.a - mean(alpha)/iter.size)^2)/(batch -
1))/sqrt(batch)
PI.MCerror = sqrt(sum((moy.p - mean(PI)/iter.size)^2)/(batch -
1))/sqrt(batch)
S1.MCerror = sqrt(sum((moy.s1 - mean(as.matrix(S1))/iter.size)^2)/(batch -
1))/sqrt(batch)
C1.MCerror = sqrt(sum((moy.c1 - mean(as.matrix(C1))/iter.size)^2)/(batch -
1))/sqrt(batch)
theta.MCerror = sqrt(sum((moy.t - mean(theta)/iter.size)^2)/(batch -
1))/sqrt(batch)
sigma.theta.MCerror = sqrt(sum((moy.sig.t - mean(sigma.theta))^2)/(batch -
1))/sqrt(batch)
}
else {
alpha.MCerror = sqrt(rowSums((moy.a - colSums(alpha)/iter.size)^2)/(batch -
1))/sqrt(batch)
PI.MCerror = sqrt(rowSums((moy.p - colSums(PI)/iter.size)^2)/(batch -
1))/sqrt(batch)
S1.MCerror = sqrt(rowSums((moy.s1 - colSums(as.matrix(S1))/iter.size)^2)/(batch -
1))/sqrt(batch)
C1.MCerror = sqrt(rowSums((moy.c1 - colSums(as.matrix(C1))/iter.size)^2)/(batch -
1))/sqrt(batch)
theta.MCerror = sqrt(rowSums((moy.t - colSums(theta)/iter.size)^2)/(batch -
1))/sqrt(batch)
sigma.theta.MCerror = sqrt(sum((moy.sig.t - mean(sigma.theta))^2)/(batch -
1))/sqrt(batch)
}
THETA.MCerror = sqrt(sum((moy.T - mean(THETA))^2)/(batch -
1))/sqrt(batch)
LAMBDA.MCerror = sqrt(sum((moy.L - mean(LAMBDA))^2)/(batch -
1))/sqrt(batch)
beta.MCerror = sqrt(sum((moy.b - mean(beta))^2)/(batch -
1))/sqrt(batch)
sigma.alpha.MCerror = sqrt(sum((moy.sig.a - mean(sigma.alpha))^2)/(batch -
1))/sqrt(batch)
S_overall.MCerror = sqrt(sum((moy.s.ov - mean(S_overall))^2)/(batch -
1))/sqrt(batch)
C_overall.MCerror = sqrt(sum((moy.c.ov - mean(C_overall))^2)/(batch -
1))/sqrt(batch)
S1_new.MCerror = sqrt(sum((moy.s1_new - mean(S1_new))^2)/(batch -
1))/sqrt(batch)
C1_new.MCerror = sqrt(sum((moy.c1_new - mean(C1_new))^2)/(batch -
1))/sqrt(batch)
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
if (sub_rs[[1]] == 1) {
for (i in 1:batch) {
moy.c2 = c(moy.c2, mean(C2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
}
C2.MCerror = sqrt(sum((moy.c2 - mean(C2))^2)/(batch -
1))/sqrt(batch)
}
else {
if (sub_rs[[1]] > 1) {
for (i in 1:batch) {
moy.c2 = cbind(moy.c2, colSums(as.matrix(C2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
}
C2.MCerror = sqrt(rowSums((moy.c2 - colSums(as.matrix(C2))/iter.size)^2)/(batch -
1))/sqrt(batch)
}
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
if (sub_rs[[1]] == 1) {
for (i in 1:batch) {
moy.s2 = c(moy.s2, mean(S2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
}
S2.MCerror = sqrt(sum((moy.s2 - mean(S2))^2)/(batch -
1))/sqrt(batch)
}
else {
if (sub_rs[[1]] > 1) {
for (i in 1:batch) {
moy.s2 = cbind(moy.s2, colSums(as.matrix(S2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
}
S2.MCerror = sqrt(rowSums((moy.s2 - colSums(as.matrix(S2))/iter.size)^2)/(batch -
1))/sqrt(batch)
}
}
}
else {
if (sub_rs[[1]] == 1) {
for (i in 1:batch) {
moy.s2 = c(moy.s2, mean(S2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.c2 = c(moy.c2, mean(C2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
}
S2.MCerror = sqrt(sum((moy.s2 - mean(S2))^2)/(batch -
1))/sqrt(batch)
C2.MCerror = sqrt(sum((moy.c2 - mean(C2))^2)/(batch -
1))/sqrt(batch)
}
else {
if (sub_rs[[1]] > 1) {
for (i in 1:batch) {
moy.s2 = cbind(moy.s2, colSums(as.matrix(S2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.c2 = cbind(moy.c2, colSums(as.matrix(C2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
}
S2.MCerror = sqrt(rowSums((moy.s2 - colSums(as.matrix(S2))/iter.size)^2)/(batch -
1))/sqrt(batch)
C2.MCerror = sqrt(rowSums((moy.c2 - colSums(as.matrix(C2))/iter.size)^2)/(batch -
1))/sqrt(batch)
}
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model == 2) {
if (sub_rs[[1]] == 1) {
for (i in 1:batch) {
moy.a1 = c(moy.a1, mean(a1[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.a0 = c(moy.a0, mean(a0[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.s2 = c(moy.s2, mean(S2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.c2 = c(moy.c2, mean(C2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
}
a1.MCerror = sqrt(sum((moy.a1 - mean(a1))^2)/(batch -
1))/sqrt(batch)
a0.MCerror = sqrt(sum((moy.a0 - mean(a0))^2)/(batch -
1))/sqrt(batch)
S2.MCerror = sqrt(sum((moy.s2 - mean(S2))^2)/(batch -
1))/sqrt(batch)
C2.MCerror = sqrt(sum((moy.c2 - mean(C2))^2)/(batch -
1))/sqrt(batch)
}
else {
if (sub_rs[[1]] > 1) {
for (i in 1:batch) {
moy.a1 = cbind(moy.a1, colSums(as.matrix(a1[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.a0 = cbind(moy.a0, colSums(as.matrix(a0[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.s2 = cbind(moy.s2, colSums(as.matrix(S2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.c2 = cbind(moy.c2, colSums(as.matrix(C2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
}
a1.MCerror = sqrt(rowSums((moy.a1 - colSums(as.matrix(a1))/iter.size)^2)/(batch -
1))/sqrt(batch)
a0.MCerror = sqrt(rowSums((moy.a0 - colSums(as.matrix(a0))/iter.size)^2)/(batch -
1))/sqrt(batch)
S2.MCerror = sqrt(rowSums((moy.s2 - colSums(as.matrix(S2))/iter.size)^2)/(batch -
1))/sqrt(batch)
C2.MCerror = sqrt(rowSums((moy.c2 - colSums(as.matrix(C2))/iter.size)^2)/(batch -
1))/sqrt(batch)
}
}
}
else {
if (condInd == FALSE) {
if (sub_rs[[1]] == 1) {
for (i in 1:batch) {
moy.a1 = c(moy.a1, mean(a1[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.a0 = c(moy.a0, mean(a0[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.b1 = c(moy.b1, mean(b1[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.b0 = c(moy.b0, mean(b0[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.d1 = c(moy.d1, mean(d1[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.d0 = c(moy.d0, mean(d0[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.s2 = c(moy.s2, mean(S2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
moy.c2 = c(moy.c2, mean(C2[round((1 + ssize *
(i - 1)), 0):round((ssize * i), 0)]))
}
a1.MCerror = sqrt(sum((moy.a1 - mean(a1))^2)/(batch -
1))/sqrt(batch)
a0.MCerror = sqrt(sum((moy.a0 - mean(a0))^2)/(batch -
1))/sqrt(batch)
b1.MCerror = sqrt(sum((moy.b1 - mean(b1))^2)/(batch -
1))/sqrt(batch)
b0.MCerror = sqrt(sum((moy.b0 - mean(b0))^2)/(batch -
1))/sqrt(batch)
d1.MCerror = sqrt(sum((moy.d1 - mean(d1))^2)/(batch -
1))/sqrt(batch)
d0.MCerror = sqrt(sum((moy.d0 - mean(d0))^2)/(batch -
1))/sqrt(batch)
S2.MCerror = sqrt(sum((moy.s2 - mean(S2))^2)/(batch -
1))/sqrt(batch)
C2.MCerror = sqrt(sum((moy.c2 - mean(C2))^2)/(batch -
1))/sqrt(batch)
}
else {
if (sub_rs[[1]] > 1) {
for (i in 1:batch) {
moy.a1 = cbind(moy.a1, colSums(as.matrix(a1[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.a0 = cbind(moy.a0, colSums(as.matrix(a0[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.b1 = cbind(moy.b1, colSums(as.matrix(b1[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.b0 = cbind(moy.b0, colSums(as.matrix(b0[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.d1 = cbind(moy.d1, colSums(as.matrix(d1[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.d0 = cbind(moy.d0, colSums(as.matrix(d0[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.s2 = cbind(moy.s2, colSums(as.matrix(S2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
moy.c2 = cbind(moy.c2, colSums(as.matrix(C2[round((1 +
ssize * (i - 1)), 0):round((ssize * i),
0), ]))/ssize)
}
a1.MCerror = sqrt(rowSums((moy.a1 - colSums(as.matrix(a1))/iter.size)^2)/(batch -
1))/sqrt(batch)
a0.MCerror = sqrt(rowSums((moy.a0 - colSums(as.matrix(a0))/iter.size)^2)/(batch -
1))/sqrt(batch)
b1.MCerror = sqrt(rowSums((moy.b1 - colSums(as.matrix(b1))/iter.size)^2)/(batch -
1))/sqrt(batch)
b0.MCerror = sqrt(rowSums((moy.b0 - colSums(as.matrix(b0))/iter.size)^2)/(batch -
1))/sqrt(batch)
d1.MCerror = sqrt(rowSums((moy.d1 - colSums(as.matrix(d1))/iter.size)^2)/(batch -
1))/sqrt(batch)
d0.MCerror = sqrt(rowSums((moy.d0 - colSums(as.matrix(d0))/iter.size)^2)/(batch -
1))/sqrt(batch)
S2.MCerror = sqrt(rowSums((moy.s2 - colSums(as.matrix(S2))/iter.size)^2)/(batch -
1))/sqrt(batch)
C2.MCerror = sqrt(rowSums((moy.c2 - colSums(as.matrix(C2))/iter.size)^2)/(batch -
1))/sqrt(batch)
}
}
}
}
}
data = list(data)
if (real_life == FALSE) {
pp = data[[1]][, 1]
pn = data[[1]][, 2]
np = data[[1]][, 3]
nn = data[[1]][, 4]
Sample.size = pp + pn + np + nn
true.alpha = tv[[1]][, 1]
true.theta = tv[[1]][, 2]
true.S1 = tv[[1]][, 3]
true.C1 = tv[[1]][, 4]
true.PI = tv[[1]][, 5]
true.THETA = tv[[2]][1]
true.sigma.theta = tv[[2]][2]
true.LAMBDA = tv[[2]][3]
true.sigma.alpha = tv[[2]][4]
true.beta = tv[[2]][5]
if (Gold_Std == TRUE) {
true.S_overall = tv[[2]][6]
true.C_overall = tv[[2]][7]
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
1) {
true.S2 = tv[[3]][1, ]
true.C2 = tv[[3]][2, ]
true.S_overall = tv[[2]][6]
true.C_overall = tv[[2]][7]
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
true.a1 = tv[[3]][3, ]
true.a0 = tv[[3]][4, ]
true.S2 = tv[[3]][1, ]
true.C2 = tv[[3]][2, ]
true.S_overall = tv[[2]][8]
true.C_overall = tv[[2]][9]
}
else {
if (condInd == FALSE) {
true.S2 = tv[[3]][1, ]
true.C2 = tv[[3]][2, ]
true.a1 = tv[[3]][3, ]
true.a0 = tv[[3]][4, ]
true.b1 = tv[[3]][5, ]
true.b0 = tv[[3]][6, ]
true.d1 = tv[[2]][6]
true.d0 = tv[[2]][7]
true.S_overall = tv[[2]][8]
true.C_overall = tv[[2]][9]
}
}
}
}
test.file = paste("Summary for N =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0), ".txt")
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("Number of chains =", nb_chains), file = test.file,
append = TRUE)
write(paste("Number of iteration within a chain =", iter.num,
" Burn in within each chain =", burn_in), file = test.file,
append = TRUE)
write(paste("Thinning interval =", Thin), file = test.file,
append = TRUE)
write(paste("Total number of iteration kept =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0)), file = test.file,
append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("File location : ", summary.path), file = test.file,
append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Date :", Sys.time()), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (Gold_Std == TRUE) {
write("Perfect reference standard", file = test.file,
append = TRUE)
}
else {
write("Imperfect reference standard", file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tSAMPLE SIZE \t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Total ++ +- -+ --"), file = test.file,
append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", Sample.size[i], "",
pp[i], "", pn[i], "", np[i], "", nn[i]), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPRIOR INFORMATION \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of prevalence (pi) is ", prior_dist_PI,
"(", round(alpha.PI, digits = 4), ",", round(beta.PI,
digits = 4), "), <=> pi in [", low.pi, ",", up.pi,
"]"), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of beta is Uniform(", round(beta.a,
4), ",", round(beta.b, 4), ")"), file = test.file,
append = TRUE)
write(paste("Prior of THETA is Uniform(", prior.THETA.lower,
",", prior.THETA.upper, ")"), file = test.file, append = TRUE)
write(paste("Prior of LAMBDA is Uniform(", prior.LAMBDA.lower,
",", prior.LAMBDA.upper, ")"), file = test.file,
append = TRUE)
if (prior_sig_a == 1) {
write(paste("Prior of sigma_alpha is uniform(", l.disp.alpha,
",", u.disp.alpha, ")"), file = test.file, append = TRUE)
}
else {
if (prior_sig_a == 2) {
write(paste("Prior of sigma_alpha^2 is uniform(",
l.disp.alpha, ",", u.disp.alpha, ")"), file = test.file,
append = TRUE)
}
else {
if (prior_sig_a == 3) {
write(paste("Prior of precision of sigma_alpha is gamma(",
l.disp.alpha, ",", u.disp.alpha, ")"), file = test.file,
append = TRUE)
}
}
}
if (prior_sig_t == 1) {
write(paste("Prior of sigma_theta is uniform(", l.disp.theta,
",", u.disp.theta, ")"), file = test.file, append = TRUE)
}
else {
if (prior_sig_t == 2) {
write(paste("Prior of sigma_theta^2 is uniform(",
l.disp.theta, ",", u.disp.theta, ")"), file = test.file,
append = TRUE)
}
else {
if (prior_sig_t == 3) {
write(paste("Prior of precision of sigma_theta is gamma(",
l.disp.theta, ",", u.disp.theta, ")"), file = test.file,
append = TRUE)
}
}
}
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
write(paste(""), file = test.file, append = TRUE)
if (Gold_se == TRUE & Gold_sp == FALSE) {
write(paste("Prior of S2 (Sensitivity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "assumed to be perfect."), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of C2 (Specificity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_C2, "(", round(Spec2.alpha[i],
digits = 4), ",", round(Spec2.beta[i],
digits = 4), "), <=> C2 in [", low.sp[i],
",", up.sp[i], "]"), file = test.file, append = TRUE)
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
write(paste("Prior of S2 (Sensitivity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_S2, "(", round(Sens2.alpha[i],
digits = 4), ",", round(Sens2.beta[i],
digits = 4), "), <=> S2 in [", low.se[i],
",", up.se[i], "]"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of C2 (Specificity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "assumed to be perfect."), file = test.file,
append = TRUE)
}
}
else {
write(paste("Prior of S2 (Sensitivity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_S2, "(", round(Sens2.alpha[i],
digits = 4), ",", round(Sens2.beta[i],
digits = 4), "), <=> S2 in [", low.se[i],
",", up.se[i], "]"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of C2 (Specificity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_C2, "(", round(Spec2.alpha[i],
digits = 4), ",", round(Spec2.beta[i],
digits = 4), "), <=> C2 in [", low.sp[i],
",", up.sp[i], "]"), file = test.file,
append = TRUE)
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a1 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a1[i], digits = 4),
",", round(sd.a1[i], digits = 4), "), <=> S2 in []"),
file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a0 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a0[i], digits = 4),
",", round(sd.a0[i], digits = 4), "), <=> C2 in []"),
file = test.file, append = TRUE)
}
}
else {
if (condInd == FALSE) {
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a1 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a1[i],
digits = 4), ",", round(sd.a1[i], digits = 4),
"), <=> S2 in []"), file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a0 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a0[i],
digits = 4), ",", round(sd.a0[i], digits = 4),
"), <=> C2 in []"), file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of b1 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Uniform (", round(low.b1[i],
digits = 4), ",", round(up.b1[i], digits = 4),
"), <=> S2 in []"), file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of b0 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Uniform (", round(low.b0[i],
digits = 4), ",", round(up.b0[i], digits = 4),
"), <=> C2 in []"), file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of d1 is Uniform(", round(low.d1,
4), ",", round(up.d1, 4), ")"), file = test.file,
append = TRUE)
write(paste("Prior of d0 is Uniform(", round(low.d0,
4), ",", round(up.d0, 4), ")"), file = test.file,
append = TRUE)
}
}
}
write(paste(), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tBETWEEN-STUDY parameters (Point estimate =",
point_estimate, ")\t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("THETA ", round(true.THETA, digits = digit),
"", round(THETA.est, digits = digit), "", round(THETA.sd,
digits = digit), "", round(THETA.MCerror, digits = digit),
"", round(THETA.HPD[1], digits = digit), "", round(THETA.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("LAMBDA ", round(true.LAMBDA, digits = digit),
"", round(LAMBDA.est, digits = digit), "", round(LAMBDA.sd,
digits = digit), "", round(LAMBDA.MCerror, digits = digit),
"", round(LAMBDA.HPD[1], digits = digit), "", round(LAMBDA.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("beta ", round(true.beta, digits = digit),
"", round(beta.est, digits = digit), "", round(beta.sd,
digits = digit), "", round(beta.MCerror, digits = digit),
"", round(beta.HPD[1], digits = digit), "", round(beta.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("sigma.alpha ", round(true.sigma.alpha, digits = digit),
"", round(sigma.alpha.est, digits = digit), "", round(sigma.alpha.sd,
digits = digit), "", round(sigma.alpha.MCerror,
digits = digit), "", round(sigma.alpha.HPD[1],
digits = digit), "", round(sigma.alpha.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("sigma.theta ", round(true.sigma.theta, digits = digit),
"", round(sigma.theta.est, digits = digit), "", round(sigma.theta.sd,
digits = digit), "", round(sigma.theta.MCerror,
digits = digit), "", round(sigma.theta.HPD[1],
digits = digit), "", round(sigma.theta.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("S overall ", round(true.S_overall, digits = digit),
"", round(S_overall.est, digits = digit), "", round(S_overall.sd,
digits = digit), "", round(S_overall.MCerror,
digits = digit), "", round(S_overall.HPD[1],
digits = digit), "", round(S_overall.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("C overall ", round(true.C_overall, digits = digit),
"", round(C_overall.est, digits = digit), "", round(C_overall.sd,
digits = digit), "", round(C_overall.MCerror,
digits = digit), "", round(C_overall.HPD[1],
digits = digit), "", round(C_overall.HPD[2],
digits = digit)), file = test.file, append = TRUE)
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tReference standard (Point estimate =",
point_estimate, ")\t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (Gold_se == TRUE & Gold_sp == FALSE) {
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.S2[i], digits = digit),
"(It was assumed to be perfect)"), file = test.file,
append = TRUE)
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.C2[i], digits = digit),
"", round(C2.est[i], digits = digit), "",
round(C2.sd[i], digits = digit), "", round(C2.MCerror[i],
digits = digit), "", round(C2.HPD[i,
1], digits = digit), "", round(C2.HPD[i,
2], digits = digit)), file = test.file,
append = TRUE)
}
}
else {
write(paste("S2 ", round(true.S2, digits = digit),
"(It was assumed to be perfect)"), file = test.file,
append = TRUE)
write(paste("C2 ", round(true.C2, digits = digit),
"", round(C2.est, digits = digit), "", round(C2.sd,
digits = digit), "", round(C2.MCerror,
digits = digit), "", round(C2.HPD[1], digits = digit),
"", round(C2.HPD[2], digits = digit)), file = test.file,
append = TRUE)
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.S2[i], digits = digit),
"", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.C2[i], digits = digit),
"(It was assumed to be perfect)"), file = test.file,
append = TRUE)
}
}
else {
write(paste("S2 ", round(true.S2, digits = digit),
"", round(S2.est, digits = digit), "",
round(S2.sd, digits = digit), "", round(S2.MCerror,
digits = digit), "", round(S2.HPD[1],
digits = digit), "", round(S2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("C2 ", round(true.C2, digits = digit),
"(It was assumed to be perfect)"), file = test.file,
append = TRUE)
}
}
else {
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.S2[i], digits = digit),
"", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.C2[i], digits = digit),
"", round(C2.est[i], digits = digit),
"", round(C2.sd[i], digits = digit),
"", round(C2.MCerror[i], digits = digit),
"", round(C2.HPD[i, 1], digits = digit),
"", round(C2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
}
else {
write(paste("S2 ", round(true.S2, digits = digit),
"", round(S2.est, digits = digit), "",
round(S2.sd, digits = digit), "", round(S2.MCerror,
digits = digit), "", round(S2.HPD[1],
digits = digit), "", round(S2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("C2 ", round(true.C2, digits = digit),
"", round(C2.est, digits = digit), "",
round(C2.sd, digits = digit), "", round(C2.MCerror,
digits = digit), "", round(C2.HPD[1],
digits = digit), "", round(C2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tReference standard (Point estimate =",
point_estimate, ")\t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("a1 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.a1[i], digits = digit),
"", round(a1.est[i], digits = digit), "",
round(a1.sd[i], digits = digit), "", round(a1.MCerror[i],
digits = digit), "", round(a1.HPD[i,
1], digits = digit), "", round(a1.HPD[i,
2], digits = digit)), file = test.file,
append = TRUE)
}
for (i in 1:rs.length) {
write(paste("a0 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.a0[i], digits = digit),
"", round(a0.est[i], digits = digit), "",
round(a0.sd[i], digits = digit), "", round(a0.MCerror[i],
digits = digit), "", round(a0.HPD[i,
1], digits = digit), "", round(a0.HPD[i,
2], digits = digit)), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.S2[i], digits = digit),
"", round(S2.est[i], digits = digit), "",
round(S2.sd[i], digits = digit), "", round(S2.MCerror[i],
digits = digit), "", round(S2.HPD[i,
1], digits = digit), "", round(S2.HPD[i,
2], digits = digit)), file = test.file,
append = TRUE)
}
for (i in 1:rs.length) {
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.C2[i], digits = digit),
"", round(C2.est[i], digits = digit), "",
round(C2.sd[i], digits = digit), "", round(C2.MCerror[i],
digits = digit), "", round(C2.HPD[i,
1], digits = digit), "", round(C2.HPD[i,
2], digits = digit)), file = test.file,
append = TRUE)
}
}
else {
write(paste("a1 ", round(true.a1, digits = digit),
"", round(a1.est, digits = digit), "", round(a1.sd,
digits = digit), "", round(a1.MCerror,
digits = digit), "", round(a1.HPD[1], digits = digit),
"", round(a1.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("a0 ", round(true.a0, digits = digit),
"", round(a0.est, digits = digit), "", round(a0.sd,
digits = digit), "", round(a0.MCerror,
digits = digit), "", round(a0.HPD[1], digits = digit),
"", round(a0.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("S2 ", round(true.S2, digits = digit),
"", round(S2.est, digits = digit), "", round(S2.sd,
digits = digit), "", round(S2.MCerror,
digits = digit), "", round(S2.HPD[1], digits = digit),
"", round(S2.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("C2 ", round(true.C2, digits = digit),
"", round(C2.est, digits = digit), "", round(C2.sd,
digits = digit), "", round(C2.MCerror,
digits = digit), "", round(C2.HPD[1], digits = digit),
"", round(C2.HPD[2], digits = digit)), file = test.file,
append = TRUE)
}
}
else {
if (condInd == FALSE) {
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tReference standard (Point estimate =",
point_estimate, ")\t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("d1 ", round(true.d1, digits = digit),
"", round(d1.est, digits = digit), "", round(d1.sd,
digits = digit), "", round(d1.MCerror,
digits = digit), "", round(d1.HPD[1], digits = digit),
"", round(d1.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("d0 ", round(true.d0, digits = digit),
"", round(d0.est, digits = digit), "", round(d0.sd,
digits = digit), "", round(d0.MCerror,
digits = digit), "", round(d0.HPD[1], digits = digit),
"", round(d0.HPD[2], digits = digit)), file = test.file,
append = TRUE)
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("a1 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.a1[i], digits = digit),
"", round(a1.est[i], digits = digit),
"", round(a1.sd[i], digits = digit),
"", round(a1.MCerror[i], digits = digit),
"", round(a1.HPD[i, 1], digits = digit),
"", round(a1.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("a0 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.a0[i], digits = digit),
"", round(a0.est[i], digits = digit),
"", round(a0.sd[i], digits = digit),
"", round(a0.MCerror[i], digits = digit),
"", round(a0.HPD[i, 1], digits = digit),
"", round(a0.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.S2[i], digits = digit),
"", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.C2[i], digits = digit),
"", round(C2.est[i], digits = digit),
"", round(C2.sd[i], digits = digit),
"", round(C2.MCerror[i], digits = digit),
"", round(C2.HPD[i, 1], digits = digit),
"", round(C2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("b1 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.b1[i], digits = digit),
"", round(b1.est[i], digits = digit),
"", round(b1.sd[i], digits = digit),
"", round(b1.MCerror[i], digits = digit),
"", round(b1.HPD[i, 1], digits = digit),
"", round(b1.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("b0 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(true.b0[i], digits = digit),
"", round(b0.est[i], digits = digit),
"", round(b0.sd[i], digits = digit),
"", round(b0.MCerror[i], digits = digit),
"", round(b0.HPD[i, 1], digits = digit),
"", round(b0.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
}
else {
write(paste("a1 ", round(true.a1, digits = digit),
"", round(a1.est, digits = digit), "",
round(a1.sd, digits = digit), "", round(a1.MCerror,
digits = digit), "", round(a1.HPD[1],
digits = digit), "", round(a1.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("a0 ", round(true.a0, digits = digit),
"", round(a0.est, digits = digit), "",
round(a0.sd, digits = digit), "", round(a0.MCerror,
digits = digit), "", round(a0.HPD[1],
digits = digit), "", round(a0.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("b1 ", round(true.b1, digits = digit),
"", round(b1.est, digits = digit), "",
round(b1.sd, digits = digit), "", round(b1.MCerror,
digits = digit), "", round(b1.HPD[1],
digits = digit), "", round(b1.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("b0 ", round(true.b0, digits = digit),
"", round(b0.est, digits = digit), "",
round(b0.sd, digits = digit), "", round(b0.MCerror,
digits = digit), "", round(b0.HPD[1],
digits = digit), "", round(b0.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("S2 ", round(true.S2, digits = digit),
"", round(S2.est, digits = digit), "",
round(S2.sd, digits = digit), "", round(S2.MCerror,
digits = digit), "", round(S2.HPD[1],
digits = digit), "", round(S2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("C2 ", round(true.C2, digits = digit),
"", round(C2.est, digits = digit), "",
round(C2.sd, digits = digit), "", round(C2.MCerror,
digits = digit), "", round(C2.HPD[1],
digits = digit), "", round(C2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
}
}
}
}
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tWITHIN-STUDY PARAMETERS \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\ttheta \t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(true.theta[i],
digits = digit), "", round(theta.est[i], digits = digit),
"", round(theta.sd[i], digits = digit), "", round(theta.MCerror[i],
digits = digit), "", round(theta.HPD[i, 1],
digits = digit), "", round(theta.HPD[i, 2],
digits = digit)), file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\talpha \t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(true.alpha[i],
digits = digit), "", round(alpha.est[i], digits = digit),
"", round(alpha.sd[i], digits = digit), "", round(alpha.MCerror[i],
digits = digit), "", round(alpha.HPD[i, 1],
digits = digit), "", round(alpha.HPD[i, 2],
digits = digit)), file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPrevalence \t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(true.PI[i], digits = digit),
"", round(PI.est[i], digits = digit), "", round(PI.sd[i],
digits = digit), "", round(PI.MCerror[i], digits = digit),
"", round(PI.HPD[i, 1], digits = digit), "",
round(PI.HPD[i, 2], digits = digit)), file = test.file,
append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tSensitivity of test 1 (S1) \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(true.S1[i], digits = digit),
"", round(S1.est[i], digits = digit), "", round(S1.sd[i],
digits = digit), "", round(S1.MCerror[i], digits = digit),
"", round(S1.HPD[i, 1], digits = digit), "",
round(S1.HPD[i, 2], digits = digit)), file = test.file,
append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tSpecificity of test 1 (C1) \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(true.C1[i], digits = digit),
"", round(C1.est[i], digits = digit), "", round(C1.sd[i],
digits = digit), "", round(C1.MCerror[i], digits = digit),
"", round(C1.HPD[i, 1], digits = digit), "",
round(C1.HPD[i, 2], digits = digit)), file = test.file,
append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPosterior predictive value of Sensitivity of test under evaluation (S1) \t "),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Sensitivity ------ ", round(S1_new.est,
digits = digit), "", round(S1_new.sd, digits = digit),
"", round(S1_new.MCerror, digits = digit), "", round(S1_new.HPD[1],
digits = digit), "", round(S1_new.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPosterior predictive value of Specificity of test under evaluation (C1) \t "),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" True_value Estimate Standard_dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Specificity ------ ", round(C1_new.est,
digits = digit), "", round(C1_new.sd, digits = digit),
"", round(C1_new.MCerror, digits = digit), "", round(C1_new.HPD[1],
digits = digit), "", round(C1_new.HPD[2], digits = digit)),
file = test.file, append = TRUE)
Num_study = c()
for (i in 1:N) {
Num_study = c(Num_study, paste("Study", i))
}
if (condInd == TRUE & Gold_Std == FALSE & model == 1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
4, 1), dimnames = list(1:rs.length, c("True Value",
paste(point_estimate, "estimate"), "HPD lower",
"HPD upper"), "C2"))
refstd_parameters[, 1, 1] <- true.C2
refstd_parameters[, 2, 1] <- C2.est
refstd_parameters[, 3, 1] <- C2.HPD[, 1]
refstd_parameters[, 4, 1] <- C2.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
1))
refstd_Parameters[, , 1] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[longueur,
i, 1], type = "n", col = 1, ylab = paste(param,
" of reference standard ", i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[,
i, 1], col = l)
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, 1], type = "l", col = "grey", ylab = paste(param,
" of reference standard ", i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i, 1]),
lwd = 4, type = "l", col = "grey", main = paste(param,
" of reference standard ", i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 4, nrow = 1)
rownames(refstd_parameters) = c("C2")
colnames(refstd_parameters) = c("True Value",
paste(point_estimate, "estimate"), "HPD.low",
"HPD.high")
refstd_parameters[1, 1] <- true.C2
refstd_parameters[1, 2] <- C2.est
refstd_parameters[1, 3] <- C2.HPD[1]
refstd_parameters[1, 4] <- C2.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 1)
colnames(refstd_Parameters) = c("C2")
refstd_Parameters[, 1] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Specificity of reference standard"
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
plot(x = longueur, y = refstd_Parameters[longueur,
1], type = "n", col = 1, ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), 1], col = l)
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Specificity of reference standard"
par(mfcol = c(5, 2))
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
1], type = "l", col = "grey", ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
longueur = 1:long
plot(density(refstd_Parameters[, 1]), lwd = 4,
type = "l", col = "grey", main = paste(param,
" of reference standard \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin))
dev.off()
}
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
4, 1), dimnames = list(1:rs.length, c("True Value",
paste(point_estimate, "estimate"), "HPD lower",
"HPD upper"), "S2"))
refstd_parameters[, 1, 1] <- true.S2
refstd_parameters[, 2, 1] <- S2.est
refstd_parameters[, 3, 1] <- S2.HPD[, 1]
refstd_parameters[, 4, 1] <- S2.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
1))
refstd_Parameters[, , 1] <- S2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Sensitivity"
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[longueur,
i, 1], type = "n", col = 1, ylab = paste(param,
" of reference standard ", i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), i, 1],
col = l)
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Sensitivity"
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, 1], type = "l", col = "grey",
ylab = paste(param, " of reference standard ",
i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = "Sensitivity"
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i, 1]),
lwd = 4, type = "l", col = "grey",
main = paste(param, " of reference standard ",
i, " \n Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin))
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 4, nrow = 1)
rownames(refstd_parameters) = c("S2")
colnames(refstd_parameters) = c("True Value",
paste(point_estimate, "estimate"), "HPD.low",
"HPD.high")
refstd_parameters[1, 1] <- true.S2
refstd_parameters[1, 2] <- S2.est
refstd_parameters[1, 3] <- S2.HPD[1]
refstd_parameters[1, 4] <- S2.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 1)
colnames(refstd_Parameters) = c("S2")
refstd_Parameters[, 1] <- S2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Sensitivity of reference standard"
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
plot(x = longueur, y = refstd_Parameters[longueur,
1], type = "n", col = 1, ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), 1], col = l)
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Sensitivity of reference standard"
par(mfcol = c(5, 2))
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
1], type = "l", col = "grey", ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = "Sensitivity of reference standard"
par(mfcol = c(5, 2))
longueur = 1:long
plot(density(refstd_Parameters[, 1]), lwd = 4,
type = "l", col = "grey", main = paste(param,
" \n Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin))
dev.off()
}
}
}
else {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
4, 2), dimnames = list(1:rs.length, c("True Value",
paste(point_estimate, "estimate"), "HPD lower",
"HPD upper"), c("S2", "C2")))
refstd_parameters[, 1, 1] <- true.S2
refstd_parameters[, 2, 1] <- S2.est
refstd_parameters[, 3, 1] <- S2.HPD[, 1]
refstd_parameters[, 4, 1] <- S2.HPD[, 2]
refstd_parameters[, 1, 2] <- true.C2
refstd_parameters[, 2, 2] <- C2.est
refstd_parameters[, 3, 2] <- C2.HPD[, 1]
refstd_parameters[, 4, 2] <- C2.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
2))
refstd_Parameters[, , 1] <- S2
refstd_Parameters[, , 2] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity")
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (j in 1:2) {
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[longueur,
i, j], type = "n", col = 1, ylab = paste(param[j],
" of reference standard ", i),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), i, j],
col = l)
}
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity")
par(mfcol = c(5, 2))
for (j in 1:2) {
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, j], type = "l", col = "grey",
ylab = paste(param[j], " of reference standard ",
i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:2) {
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i,
j]), lwd = 4, type = "l", col = "grey",
main = paste(param[j], " of reference standard ",
i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 4, nrow = 2)
rownames(refstd_parameters) = c("S2", "C2")
colnames(refstd_parameters) = c("True Value",
paste(point_estimate, "estimate"), "HPD.low",
"HPD.high")
refstd_parameters[1, 1] <- true.S2
refstd_parameters[1, 2] <- S2.est
refstd_parameters[1, 3] <- S2.HPD[1]
refstd_parameters[1, 4] <- S2.HPD[2]
refstd_parameters[2, 1] <- true.C2
refstd_parameters[2, 2] <- C2.est
refstd_parameters[2, 3] <- C2.HPD[1]
refstd_parameters[2, 4] <- C2.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 2)
colnames(refstd_Parameters) = c("S2", "C2")
refstd_Parameters[, 1] <- S2
refstd_Parameters[, 2] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
param = c("Sensitivity of reference standard",
"Specificity of reference standard")
par(mfcol = c(5, 2))
for (j in 1:2) {
plot(x = longueur, y = refstd_Parameters[longueur,
j], type = "n", col = 1, ylab = paste(param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), j], col = l)
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
Param = c("Sensitivity of reference standard",
"Specificity of reference standard")
par(mfcol = c(5, 2))
for (j in 1:2) {
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
j], type = "l", col = "grey", ylab = paste(Param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity of reference standard",
"Specificity of reference standard")
par(mfcol = c(5, 2))
for (j in 1:2) {
longueur = 1:long
plot(density(refstd_Parameters[, j]),
lwd = 4, type = "l", col = "grey",
main = paste(param[j], " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
}
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
4, 4), dimnames = list(1:rs.length, c("True Value",
paste(point_estimate, "estimate"), "HPD lower",
"HPD upper"), c("S2", "C2", "a1", "a0")))
refstd_parameters[, 1, 1] <- true.S2
refstd_parameters[, 2, 1] <- S2.est
refstd_parameters[, 3, 1] <- S2.HPD[, 1]
refstd_parameters[, 4, 1] <- S2.HPD[, 2]
refstd_parameters[, 1, 2] <- true.C2
refstd_parameters[, 2, 2] <- C2.est
refstd_parameters[, 3, 2] <- C2.HPD[, 1]
refstd_parameters[, 4, 2] <- C2.HPD[, 2]
refstd_parameters[, 1, 3] <- true.a1
refstd_parameters[, 2, 3] <- a1.est
refstd_parameters[, 3, 3] <- a1.HPD[, 1]
refstd_parameters[, 4, 3] <- a1.HPD[, 2]
refstd_parameters[, 1, 4] <- true.a0
refstd_parameters[, 2, 4] <- a0.est
refstd_parameters[, 3, 4] <- a0.HPD[, 1]
refstd_parameters[, 4, 4] <- a0.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
4))
refstd_Parameters[, , 1] <- S2
refstd_Parameters[, , 2] <- C2
refstd_Parameters[, , 3] <- a1
refstd_Parameters[, , 4] <- a0
if (print_plot == TRUE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity", "a1",
"a0")
par(mfcol = c(5, 2))
for (j in 1:4) {
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, j], type = "l", col = "grey", ylab = paste(param[j],
" of reference standard ", i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
abline(a = refstd_parameters[i, 2, j],
b = 0, col = "black", lwd = 3)
abline(a = refstd_parameters[i, 1, j],
b = 0, col = "red", lwd = 3)
abline(a = refstd_parameters[i, 3, j],
b = 0, col = "green", lwd = 3)
abline(a = refstd_parameters[i, 4, j],
b = 0, col = "green", lwd = 3)
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity", "a1",
"a0")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:4) {
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i, j]),
lwd = 4, type = "l", col = "grey",
main = paste(param[j], " of reference standard ",
i, " \n Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin))
}
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 4, nrow = 4)
rownames(refstd_parameters) = c("S2", "C2",
"a1", "a0")
colnames(refstd_parameters) = c("True Value",
paste(point_estimate, "estimate"), "HPD.low",
"HPD.high")
refstd_parameters[1, 1] <- true.S2
refstd_parameters[1, 2] <- S2.est
refstd_parameters[1, 3] <- S2.HPD[1]
refstd_parameters[1, 4] <- S2.HPD[2]
refstd_parameters[2, 1] <- true.C2
refstd_parameters[2, 2] <- C2.est
refstd_parameters[2, 3] <- C2.HPD[1]
refstd_parameters[2, 4] <- C2.HPD[2]
refstd_parameters[3, 1] <- true.a1
refstd_parameters[3, 2] <- a1.est
refstd_parameters[3, 3] <- a1.HPD[1]
refstd_parameters[3, 4] <- a1.HPD[2]
refstd_parameters[4, 1] <- true.a0
refstd_parameters[4, 2] <- a0.est
refstd_parameters[4, 3] <- a0.HPD[1]
refstd_parameters[4, 4] <- a0.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 4)
colnames(refstd_Parameters) = c("S2", "C2",
"a1", "a0")
refstd_Parameters[, 1] <- S2
refstd_Parameters[, 2] <- C2
refstd_Parameters[, 3] <- a1
refstd_Parameters[, 4] <- a0
if (print_plot == TRUE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
Param = c("Sensitivity of reference standard",
"Specificity of reference standard", "a1",
"a0")
par(mfcol = c(5, 2))
for (j in 1:4) {
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
j], type = "l", col = "grey", ylab = paste(Param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
abline(a = refstd_parameters[j, 2], b = 0,
col = "black", lwd = 3)
abline(a = refstd_parameters[j, 1], b = 0,
col = "red", lwd = 3)
abline(a = refstd_parameters[j, 3], b = 0,
col = "green", lwd = 3)
abline(a = refstd_parameters[j, 4], b = 0,
col = "green", lwd = 3)
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity of reference standard",
"Specificity of reference standard", "a1",
"a0")
par(mfcol = c(5, 2))
for (j in 1:4) {
longueur = 1:long
plot(density(refstd_Parameters[, j]), lwd = 4,
type = "l", col = "grey", main = paste(param[j],
" \n Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin))
}
dev.off()
}
}
}
}
parameters = array(0, dim = c(N, 4, 5), dimnames = list(Num_study,
c("True value", paste(point_estimate, "estimate"),
"HPD lower", "HPD upper"), c("theta", "alpha",
"pi", "S1", "C1")))
parameters[, 1, 1] <- true.theta
parameters[, 2, 1] <- theta.est
parameters[, 3, 1] <- theta.HPD[, 1]
parameters[, 4, 1] <- theta.HPD[, 2]
parameters[, 1, 2] <- true.alpha
parameters[, 2, 2] <- alpha.est
parameters[, 3, 2] <- alpha.HPD[, 1]
parameters[, 4, 2] <- alpha.HPD[, 2]
parameters[, 1, 3] <- true.PI
parameters[, 2, 3] <- PI.est
parameters[, 3, 3] <- PI.HPD[, 1]
parameters[, 4, 3] <- PI.HPD[, 2]
parameters[, 1, 4] <- true.S1
parameters[, 2, 4] <- S1.est
parameters[, 3, 4] <- S1.HPD[, 1]
parameters[, 4, 4] <- S1.HPD[, 2]
parameters[, 1, 5] <- true.C1
parameters[, 2, 5] <- C1.est
parameters[, 3, 5] <- C1.HPD[, 1]
parameters[, 4, 5] <- C1.HPD[, 2]
long = length(alpha[, 1])
Parameters = array(0, c(long, N, 5))
Parameters[, , 1] <- theta
Parameters[, , 2] <- alpha
Parameters[, , 3] <- PI
Parameters[, , 4] <- S1
Parameters[, , 5] <- C1
parameter = matrix(0, ncol = 4, nrow = 9)
rownames(parameter) = c("THETA", "LAMBDA", "beta", "sigma.alpha",
"sigma.theta", "S Overall", "C Overall", "S1_new",
"C1_new")
colnames(parameter) = c("True.value", paste(point_estimate,
"estimate"), "HPD.low", "HPD.high")
parameter[1, 1] <- true.THETA
parameter[1, 2] <- THETA.est
parameter[1, 3] <- THETA.HPD[1]
parameter[1, 4] <- THETA.HPD[2]
parameter[2, 1] <- true.LAMBDA
parameter[2, 2] <- LAMBDA.est
parameter[2, 3] <- LAMBDA.HPD[1]
parameter[2, 4] <- LAMBDA.HPD[2]
parameter[3, 1] <- true.beta
parameter[3, 2] <- beta.est
parameter[3, 3] <- beta.HPD[1]
parameter[3, 4] <- beta.HPD[2]
parameter[4, 1] <- true.sigma.alpha
parameter[4, 2] <- sigma.alpha.est
parameter[4, 3] <- sigma.alpha.HPD[1]
parameter[4, 4] <- sigma.alpha.HPD[2]
parameter[5, 1] <- true.sigma.theta
parameter[5, 2] <- sigma.theta.est
parameter[5, 3] <- sigma.theta.HPD[1]
parameter[5, 4] <- sigma.theta.HPD[2]
parameter[6, 1] <- true.S_overall
parameter[6, 2] <- S_overall.est
parameter[6, 3] <- S_overall.HPD[1]
parameter[6, 4] <- S_overall.HPD[2]
parameter[7, 1] <- true.C_overall
parameter[7, 2] <- C_overall.est
parameter[7, 3] <- C_overall.HPD[1]
parameter[7, 4] <- C_overall.HPD[2]
parameter[8, 1] <- S1_new.est
parameter[8, 2] <- S1_new.est
parameter[8, 3] <- S1_new.HPD[1]
parameter[8, 4] <- S1_new.HPD[2]
parameter[9, 1] <- C1_new.est
parameter[9, 2] <- C1_new.est
parameter[9, 3] <- C1_new.HPD[1]
parameter[9, 4] <- C1_new.HPD[2]
long = length(THETA)
Parameter = matrix(0, nrow = long, ncol = 9)
colnames(Parameter) = c("THETA", "LAMBDA", "beta", "sigma.alpha",
"sigma.theta", "S overall", "C overall", "S1_new",
"C1_new")
Parameter[, 1] <- THETA
Parameter[, 2] <- LAMBDA
Parameter[, 3] <- beta
Parameter[, 4] <- sigma.alpha
Parameter[, 5] <- sigma.theta
Parameter[, 6] <- S_overall
Parameter[, 7] <- C_overall
Parameter[, 8] <- S1_new
Parameter[, 9] <- C1_new
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
no_chains = length(chain)
}
else {
if (is.null(chain) == TRUE) {
no_chains = 1
}
}
file.pdf5 = paste("Trace plots for N =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0),
".pdf")
pdf(file.pdf5, paper = "a4", height = 20)
param = c("theta", "alpha", "PI", "S1", "C1")
Param = c("Capital Theta", "Capital Lambda", "beta",
"~sigma[alpha]", "~sigma[theta]", "S Overall",
"C Overall", "S1_new", "C1_new")
iter_chain = round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0)/no_chains
min_param = c(min(Parameters[, , 1]), min(Parameters[,
, 2]), min(Parameters[, , 3]), min(Parameters[,
, 4]), min(Parameters[, , 5]))
max_param = c(max(Parameters[, , 1]), max(Parameters[,
, 2]), max(Parameters[, , 3]), max(Parameters[,
, 4]), max(Parameters[, , 5]))
dlag = (max_param - min_param)/100
range_param = numeric()
for (j in 1:5) {
range_param = cbind(range_param, seq(min_param[j] +
dlag[j]/2, max_param[j] - dlag[j]/2, by = dlag[j]))
}
par(mfcol = c(5, 2))
longueur = 1:iter_chain
for (j in 1:5) {
for (i in 1:N) {
plot(x = longueur, y = Parameters[longueur,
i, j], type = "n", col = 1, ylab = paste(param[j],
" of study ", i), xlab = "iteration number",
main = paste("Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin),
ylim = range(range_param[, j]))
for (l in 1:no_chains) {
lines(x = longueur, y = Parameters[longueur +
(iter_chain * (l - 1)), i, j], col = l)
}
}
}
min_Param = c(min(Parameter[, 1]), min(Parameter[,
2]), min(Parameter[, 3]), min(Parameter[, 4]),
min(Parameter[, 5]), min(Parameter[, 6]), min(Parameter[,
7]), min(Parameter[, 8]), min(Parameter[, 9]))
max_Param = c(max(Parameter[, 1]), max(Parameter[,
2]), max(Parameter[, 3]), max(Parameter[, 4]),
max(Parameter[, 5]), max(Parameter[, 6]), max(Parameter[,
7]), max(Parameter[, 8]), max(Parameter[, 9]))
dlag = (max_Param - min_Param)/100
range_Param = numeric()
for (j in 1:9) {
range_Param = cbind(range_Param, seq(min_Param[j] +
dlag[j]/2, max_Param[j] - dlag[j]/2, by = dlag[j]))
}
for (j in 1:9) {
plot(x = longueur, y = Parameter[longueur, j],
type = "n", col = 1, ylab = paste(Param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin),
ylim = range(range_Param[, j]))
for (l in 1:no_chains) {
lines(x = longueur, y = Parameter[longueur +
(iter_chain * (l - 1)), j], col = l)
}
}
dev.off()
file.pdf3 = paste("Density plots for N =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0),
".pdf")
pdf(file.pdf3, paper = "a4", height = 20)
param = c("theta", "alpha", "PI", "S1", "C1")
Param = c("Capital Theta", "Capital Lambda", "beta",
"~sigma[alpha]", "~sigma[theta]", "S Overall",
"C Overall", "S1_new", "C1_new")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:5) {
for (i in 1:N) {
plot(density(Parameters[, i, j]), lwd = 4,
type = "l", col = "grey", main = paste(param[j],
" of study ", i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin))
}
}
for (j in 1:9) {
plot(density(Parameter[, j]), lwd = 4, type = "l",
col = "grey", main = paste(Param[j], " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin))
}
dev.off()
Sensi1 = apply(as.matrix(Parameters[, , 4]), 2, median)
Speci1 = apply(as.matrix(Parameters[, , 5]), 2, median)
Ov_Se = 1 - pnorm((median(Parameter[, 1]) - median(Parameter[,
2])/2)/exp(median(Parameter[, 3])/2))
Ov_Sp = pnorm((median(Parameter[, 1]) + median(Parameter[,
2])/2)/exp(-median(Parameter[, 3])/2))
thet = qnorm((1 - as.matrix(Parameters[, , 4])) +
1e-14) * exp(Parameter[, 3]/2) + Parameter[,
2]/2
min_TH = quantile(thet, trunc_low)
max_TH = quantile(thet, (1 - trunc_up))
dTH = 5e-05
TH_range = seq(min_TH + dTH/2, max_TH - dTH/2, dTH)
S_sroc = 1 - pnorm((TH_range - median(Parameter[,
2])/2)/exp(median(Parameter[, 3])/2))
C_sroc = pnorm((TH_range + median(Parameter[, 2])/2)/exp(-median(Parameter[,
3])/2))
if (cred_region == TRUE) {
min_t = 0
max_t = 2 * pi
dTt = 5e-05
range_t = seq(min_t + dTt/2, max_t - dTt/2, dTt)
bound_cte = sqrt(qchisq(1 - region_level, 2,
lower.tail = FALSE))
hat_mu_A = exp(-median(Parameter[, 3])/2) * (-median(Parameter[,
1]) + median(Parameter[, 2])/2)
s_A = sd(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2))
cos_fun_A = cos(range_t)
hat_mu_B = exp(median(Parameter[, 3])/2) * (median(Parameter[,
1]) + median(Parameter[, 2])/2)
s_B = sd(exp(Parameter[, 3]/2) * (Parameter[,
1] + Parameter[, 2]/2))
r = cor(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2), exp(Parameter[, 3]/2) *
(Parameter[, 1] + Parameter[, 2]/2))
cos_fun_B = cos(range_t + acos(r))
probit_S_credible = hat_mu_A + s_A * bound_cte *
cos_fun_A
probit_C_credible = hat_mu_B + s_B * bound_cte *
cos_fun_B
S_credible = pnorm(probit_S_credible)
C_credible = pnorm(probit_C_credible)
}
if (predict_region == TRUE) {
min_t = 0
max_t = 2 * pi
dTt = 5e-05
range_t = seq(min_t + dTt/2, max_t - dTt/2, dTt)
bound_cte = sqrt(qchisq(1 - region_level, 2,
lower.tail = FALSE))
hat_mu_A = exp(-median(Parameter[, 3])/2) * (-median(Parameter[,
1]) + median(Parameter[, 2])/2)
s_A = sd(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2))
hat_sigma_A = sqrt((exp(-median(Parameter[, 3])) *
(median((Parameter[, 5])^2) + 0.25 * median((Parameter[,
4])^2))))
cos_fun_A = cos(range_t)
hat_mu_B = exp(median(Parameter[, 3])/2) * (median(Parameter[,
1]) + median(Parameter[, 2])/2)
s_B = sd(exp(Parameter[, 3]/2) * (Parameter[,
1] + Parameter[, 2]/2))
hat_sigma_B = sqrt((exp(median(Parameter[, 3])) *
(median((Parameter[, 5])^2) + 0.25 * median((Parameter[,
4])^2))))
r = cor(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2), exp(Parameter[, 3]/2) *
(Parameter[, 1] + Parameter[, 2]/2))
hat_sigma_AB = -median((Parameter[, 5])^2) +
0.25 * median((Parameter[, 4])^2)
cos_fun_B = cos(range_t + acos(r))
probit_S_prediction = hat_mu_A + (sqrt(s_A^2 +
hat_sigma_A^2)) * bound_cte * cos_fun_A
probit_C_prediction = hat_mu_B + (sqrt(s_B^2 +
hat_sigma_B^2)) * bound_cte * cos(range_t +
acos(r))
S_prediction = pnorm(probit_S_prediction)
C_prediction = pnorm(probit_C_prediction)
}
pdf("Summary ROC curve.pdf")
default.x = range(1, 0)
default.y = range(0, 1)
plot(x = default.x, y = default.y, type = "n", xlim = rev(range(default.x)),
xlab = "", ylab = "")
title(xlab = "Specificity", ylab = "Sensitivity",
cex.lab = 1.5, main = "Summary ROC curve")
if (plot.ind.studies == TRUE) {
if (Gold_Std == TRUE) {
Scale_factor = 10
SENSi1 = data[[1]][, 1]/(data[[1]][, 1] + data[[1]][,
3])
SPECi1 = data[[1]][, 4]/(data[[1]][, 4] + data[[1]][,
2])
symbols(SPECi1, SENSi1, circles = rowSums(as.matrix(data[[1]])),
inches = 0.1 * Scale_factor/7, add = TRUE,
fg = 1)
}
else {
Scale_factor = 10
symbols(Speci1, Sensi1, circles = rowSums(as.matrix(data[[1]])),
inches = 0.1 * Scale_factor/7, add = TRUE,
fg = 1)
}
}
if (cred_region == TRUE) {
lines(C_credible, S_credible, lwd = 3, lty = lty.cred.region,
col = col.pooled.estimate)
}
if (predict_region == TRUE) {
lines(C_prediction, S_prediction, lwd = 3, lty = lty.predict.region,
col = col.predict.region)
}
lines(C_sroc, S_sroc, lwd = 3, col = "black", lty = 1)
points(Ov_Sp, Ov_Se, pch = 19, cex = 2.5, col = col.pooled.estimate)
dev.off()
}
}
else {
if (real_life == TRUE) {
d = as.matrix(data[[1]])
Sample.size = d[, 1] + d[, 2] + d[, 3] + d[, 4]
pp = d[, 1]
pn = d[, 2]
np = d[, 3]
nn = d[, 4]
test.file = paste("Summary for N =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0),
".txt")
write(paste("Number of chains =", nb_chains), file = test.file,
append = TRUE)
write(paste("Number of iteration within a chain =",
iter.num, " Burn in within each chain =",
burn_in), file = test.file, append = TRUE)
write(paste("Thinning interval =", Thin), file = test.file,
append = TRUE)
write(paste("Total number of iteration kept =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0)),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("File location : ", summary.path), file = test.file,
append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Date :", Sys.time()), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (Gold_Std == TRUE) {
write("Perfect reference standard", file = test.file,
append = TRUE)
}
else {
write("Imperfect reference standard", file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tSAMPLE SIZE \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Total ++ +- -+ --"), file = test.file,
append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", Sample.size[i],
"", pp[i], "", pn[i], "", np[i], "", nn[i]),
file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPRIOR INFORMATION \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of prevalence (pi) is ", prior_dist_PI,
"(", round(alpha.PI, digits = 4), ",", round(beta.PI,
digits = 4), "), <=> pi in [", low.pi, ",",
up.pi, "]"), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of beta is Uniform(", round(beta.a,
4), ",", round(beta.b, 4), ")"), file = test.file,
append = TRUE)
write(paste("Prior of THETA is Uniform(", prior.THETA.lower,
",", prior.THETA.upper, ")"), file = test.file,
append = TRUE)
write(paste("Prior of LAMBDA is Uniform(", prior.LAMBDA.lower,
",", prior.LAMBDA.upper, ")"), file = test.file,
append = TRUE)
if (prior_sig_a == 1) {
write(paste("Prior of sigma_alpha is uniform(",
l.disp.alpha, ",", u.disp.alpha, ")"), file = test.file,
append = TRUE)
}
else {
if (prior_sig_a == 2) {
write(paste("Prior of sigma_alpha^2 is uniform(",
l.disp.alpha, ",", u.disp.alpha, ")"), file = test.file,
append = TRUE)
}
else {
if (prior_sig_a == 3) {
write(paste("Prior of precision of sigma_alpha is gamma(",
l.disp.alpha, ",", u.disp.alpha, ")"),
file = test.file, append = TRUE)
}
}
}
if (prior_sig_t == 1) {
write(paste("Prior of sigma_theta is uniform(",
l.disp.theta, ",", u.disp.theta, ")"), file = test.file,
append = TRUE)
}
else {
if (prior_sig_t == 2) {
write(paste("Prior of sigma_theta^2 is uniform(",
l.disp.theta, ",", u.disp.theta, ")"), file = test.file,
append = TRUE)
}
else {
if (prior_sig_t == 3) {
write(paste("Prior of precision of sigma_theta is gamma(",
l.disp.theta, ",", u.disp.theta, ")"),
file = test.file, append = TRUE)
}
}
}
if (condInd == TRUE & Gold_Std == FALSE & model ==
1) {
write(paste(""), file = test.file, append = TRUE)
if (Gold_se == TRUE & Gold_sp == FALSE) {
write(paste("Prior of S2 (Sensitivity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "assumed to be perfect."), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of C2 (Specificity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_C2, "(", round(Spec2.alpha[i],
digits = 4), ",", round(Spec2.beta[i],
digits = 4), "), <=> C2 in [", low.sp[i],
",", up.sp[i], "]"), file = test.file,
append = TRUE)
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
write(paste("Prior of S2 (Sensitivity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_S2, "(", round(Sens2.alpha[i],
digits = 4), ",", round(Sens2.beta[i],
digits = 4), "), <=> S2 in [", low.se[i],
",", up.se[i], "]"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of C2 (Specificity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "assumed to be perfect."), file = test.file,
append = TRUE)
}
}
else {
write(paste("Prior of S2 (Sensitivity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_S2, "(", round(Sens2.alpha[i],
digits = 4), ",", round(Sens2.beta[i],
digits = 4), "), <=> S2 in [", low.se[i],
",", up.se[i], "]"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of C2 (Specificity of reference test) is "),
file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", prior_dist_C2, "(", round(Spec2.alpha[i],
digits = 4), ",", round(Spec2.beta[i],
digits = 4), "), <=> C2 in [", low.sp[i],
",", up.sp[i], "]"), file = test.file,
append = TRUE)
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a1 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a1[i],
digits = 4), ",", round(sd.a1[i], digits = 4),
"), <=> S2 in []"), file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a0 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i + 1]][1],
"to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a0[i],
digits = 4), ",", round(sd.a0[i], digits = 4),
"), <=> C2 in []"), file = test.file, append = TRUE)
}
}
else {
if (condInd == FALSE) {
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a1 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a1[i],
digits = 4), ",", round(sd.a1[i], digits = 4),
"), <=> S2 in []"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of a0 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Normal (", round(mean.a0[i],
digits = 4), ",", round(sd.a0[i], digits = 4),
"), <=> C2 in []"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of b1 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Uniform (", round(low.b1[i],
digits = 4), ",", round(up.b1[i], digits = 4),
"), <=> S2 in []"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of b0 is "), file = test.file,
append = TRUE)
for (i in 1:rs.length) {
write(paste("Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "Uniform (", round(low.b0[i],
digits = 4), ",", round(up.b0[i], digits = 4),
"), <=> C2 in []"), file = test.file,
append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
write(paste("Prior of d1 is Uniform(", round(low.d1,
4), ",", round(up.d1, 4), ")"), file = test.file,
append = TRUE)
write(paste("Prior of d0 is Uniform(", round(low.d0,
4), ",", round(up.d0, 4), ")"), file = test.file,
append = TRUE)
}
}
}
write(paste(), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tBETWEEN_STUDY parameters (Point estimate =",
point_estimate, ")\t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("THETA ", round(THETA.est, digits = digit),
"", round(THETA.sd, digits = digit), "", round(THETA.MCerror,
digits = digit), "", round(THETA.HPD[1], digits = digit),
"", round(THETA.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("LAMBDA ", round(LAMBDA.est, digits = digit),
"", round(LAMBDA.sd, digits = digit), "", round(LAMBDA.MCerror,
digits = digit), "", round(LAMBDA.HPD[1], digits = digit),
"", round(LAMBDA.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("beta ", round(beta.est, digits = digit),
"", round(beta.sd, digits = digit), "", round(beta.MCerror,
digits = digit), "", round(beta.HPD[1], digits = digit),
"", round(beta.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("sigma.alpha ", round(sigma.alpha.est,
digits = digit), "", round(sigma.alpha.sd, digits = digit),
"", round(sigma.alpha.MCerror, digits = digit),
"", round(sigma.alpha.HPD[1], digits = digit),
"", round(sigma.alpha.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("sigma.theta ", round(sigma.theta.est,
digits = digit), "", round(sigma.theta.sd, digits = digit),
"", round(sigma.theta.MCerror, digits = digit),
"", round(sigma.theta.HPD[1], digits = digit),
"", round(sigma.theta.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("S overall ", round(S_overall.est,
digits = digit), "", round(S_overall.sd, digits = digit),
"", round(S_overall.MCerror, digits = digit),
"", round(S_overall.HPD[1], digits = digit),
"", round(S_overall.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("C overall ", round(C_overall.est,
digits = digit), "", round(C_overall.sd, digits = digit),
"", round(C_overall.MCerror, digits = digit),
"", round(C_overall.HPD[1], digits = digit),
"", round(C_overall.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (condInd == TRUE & Gold_Std == FALSE & model ==
1) {
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tReference standard (Point estimate =",
point_estimate, ")\t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (Gold_se == TRUE & Gold_sp == FALSE) {
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "(It was assumed to be perfect)"),
file = test.file, append = TRUE)
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(C2.est[i], digits = digit),
"", round(C2.sd[i], digits = digit),
"", round(C2.MCerror[i], digits = digit),
"", round(C2.HPD[i, 1], digits = digit),
"", round(C2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
}
else {
write(paste("S2 (It was assumed to be perfect)"),
file = test.file, append = TRUE)
write(paste("C2 ", round(C2.est, digits = digit),
"", round(C2.sd, digits = digit), "", round(C2.MCerror,
digits = digit), "", round(C2.HPD[1],
digits = digit), "", round(C2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "(It was assumed to be perfect)"),
file = test.file, append = TRUE)
}
}
else {
write(paste("S2 ", round(S2.est,
digits = digit), "", round(S2.sd, digits = digit),
"", round(S2.MCerror, digits = digit),
"", round(S2.HPD[1], digits = digit),
"", round(S2.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("C2 (It was assumed to be perfect)"),
file = test.file, append = TRUE)
}
}
else {
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(C2.est[i], digits = digit),
"", round(C2.sd[i], digits = digit),
"", round(C2.MCerror[i], digits = digit),
"", round(C2.HPD[i, 1], digits = digit),
"", round(C2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
}
else {
write(paste("S2 ", round(S2.est,
digits = digit), "", round(S2.sd, digits = digit),
"", round(S2.MCerror, digits = digit),
"", round(S2.HPD[1], digits = digit),
"", round(S2.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("C2 ", round(C2.est,
digits = digit), "", round(C2.sd, digits = digit),
"", round(C2.MCerror, digits = digit),
"", round(C2.HPD[1], digits = digit),
"", round(C2.HPD[2], digits = digit)),
file = test.file, append = TRUE)
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tReference standard (Point estimate =",
point_estimate, ")\t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("a1 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(a1.est[i], digits = digit),
"", round(a1.sd[i], digits = digit),
"", round(a1.MCerror[i], digits = digit),
"", round(a1.HPD[i, 1], digits = digit),
"", round(a1.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("a0 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(a0.est[i], digits = digit),
"", round(a0.sd[i], digits = digit),
"", round(a0.MCerror[i], digits = digit),
"", round(a0.HPD[i, 1], digits = digit),
"", round(a0.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(C2.est[i], digits = digit),
"", round(C2.sd[i], digits = digit),
"", round(C2.MCerror[i], digits = digit),
"", round(C2.HPD[i, 1], digits = digit),
"", round(C2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
}
else {
write(paste("a1 ", round(a1.est, digits = digit),
"", round(a1.sd, digits = digit), "", round(a1.MCerror,
digits = digit), "", round(a1.HPD[1],
digits = digit), "", round(a1.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("a0 ", round(a0.est, digits = digit),
"", round(a0.sd, digits = digit), "", round(a0.MCerror,
digits = digit), "", round(a0.HPD[1],
digits = digit), "", round(a0.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("S2 ", round(S2.est, digits = digit),
"", round(S2.sd, digits = digit), "", round(S2.MCerror,
digits = digit), "", round(S2.HPD[1],
digits = digit), "", round(S2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("C2 ", round(C2.est, digits = digit),
"", round(C2.sd, digits = digit), "", round(C2.MCerror,
digits = digit), "", round(C2.HPD[1],
digits = digit), "", round(C2.HPD[2],
digits = digit)), file = test.file, append = TRUE)
}
}
else {
if (condInd == FALSE) {
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tReference standard (Point estimate =",
point_estimate, ")\t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("d1 ", round(d1.est, digits = digit),
"", round(d1.sd, digits = digit), "", round(d1.MCerror,
digits = digit), "", round(d1.HPD[1],
digits = digit), "", round(d1.HPD[2],
digits = digit)), file = test.file, append = TRUE)
write(paste("d0 ", round(d0.est, digits = digit),
"", round(d0.sd, digits = digit), "", round(d0.MCerror,
digits = digit), "", round(d0.HPD[1],
digits = digit), "", round(d0.HPD[2],
digits = digit)), file = test.file, append = TRUE)
if (rs.length != 1) {
for (i in 1:rs.length) {
write(paste("a1 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(a1.est[i], digits = digit),
"", round(a1.sd[i], digits = digit),
"", round(a1.MCerror[i], digits = digit),
"", round(a1.HPD[i, 1], digits = digit),
"", round(a1.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("a0 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(a0.est[i], digits = digit),
"", round(a0.sd[i], digits = digit),
"", round(a0.MCerror[i], digits = digit),
"", round(a0.HPD[i, 1], digits = digit),
"", round(a0.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
write(paste(""), file = test.file, append = TRUE)
for (i in 1:rs.length) {
write(paste("S2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(S2.est[i], digits = digit),
"", round(S2.sd[i], digits = digit),
"", round(S2.MCerror[i], digits = digit),
"", round(S2.HPD[i, 1], digits = digit),
"", round(S2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("C2 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(C2.est[i], digits = digit),
"", round(C2.sd[i], digits = digit),
"", round(C2.MCerror[i], digits = digit),
"", round(C2.HPD[i, 1], digits = digit),
"", round(C2.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("b1 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(b1.est[i], digits = digit),
"", round(b1.sd[i], digits = digit),
"", round(b1.MCerror[i], digits = digit),
"", round(b1.HPD[i, 1], digits = digit),
"", round(b1.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
for (i in 1:rs.length) {
write(paste("b0 of Study(ies) ", sub_rs[[i +
1]][1], "to", sub_rs[[i + 1]][length(sub_rs[[i +
1]])], "", round(b0.est[i], digits = digit),
"", round(b0.sd[i], digits = digit),
"", round(b0.MCerror[i], digits = digit),
"", round(b0.HPD[i, 1], digits = digit),
"", round(b0.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
}
else {
write(paste("a1 ", round(a1.est,
digits = digit), "", round(a1.sd, digits = digit),
"", round(a1.MCerror, digits = digit),
"", round(a1.HPD[1], digits = digit),
"", round(a1.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("a0 ", round(a0.est,
digits = digit), "", round(a0.sd, digits = digit),
"", round(a0.MCerror, digits = digit),
"", round(a0.HPD[1], digits = digit),
"", round(a0.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("b1 ", round(b1.est,
digits = digit), "", round(b1.sd, digits = digit),
"", round(b1.MCerror, digits = digit),
"", round(b1.HPD[1], digits = digit),
"", round(b1.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("b0 ", round(b0.est,
digits = digit), "", round(b0.sd, digits = digit),
"", round(b0.MCerror, digits = digit),
"", round(b0.HPD[1], digits = digit),
"", round(b0.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("S2 ", round(S2.est,
digits = digit), "", round(S2.sd, digits = digit),
"", round(S2.MCerror, digits = digit),
"", round(S2.HPD[1], digits = digit),
"", round(S2.HPD[2], digits = digit)),
file = test.file, append = TRUE)
write(paste("C2 ", round(C2.est,
digits = digit), "", round(C2.sd, digits = digit),
"", round(C2.MCerror, digits = digit),
"", round(C2.HPD[1], digits = digit),
"", round(C2.HPD[2], digits = digit)),
file = test.file, append = TRUE)
}
}
}
}
write(paste(""), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tWITHIN-STUDY PARAMETERS \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\ttheta \t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(theta.est[i],
digits = digit), "", round(theta.sd[i], digits = digit),
"", round(theta.MCerror[i], digits = digit),
"", round(theta.HPD[i, 1], digits = digit),
"", round(theta.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\talpha \t "), file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(alpha.est[i],
digits = digit), "", round(alpha.sd[i], digits = digit),
"", round(alpha.MCerror[i], digits = digit),
"", round(alpha.HPD[i, 1], digits = digit),
"", round(alpha.HPD[i, 2], digits = digit)),
file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPrevalence \t "), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(PI.est[i],
digits = digit), "", round(PI.sd[i], digits = digit),
"", round(PI.MCerror[i], digits = digit), "",
round(PI.HPD[i, 1], digits = digit), "", round(PI.HPD[i,
2], digits = digit)), file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tSensitivity of test 1 (S1) \t "),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(S1.est[i],
digits = digit), "", round(S1.sd[i], digits = digit),
"", round(S1.MCerror[i], digits = digit), "",
round(S1.HPD[i, 1], digits = digit), "", round(S1.HPD[i,
2], digits = digit)), file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tSpecificity of test 1 (C1) \t "),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
for (i in 1:N) {
write(paste("Study ", i, "", round(C1.est[i],
digits = digit), "", round(C1.sd[i], digits = digit),
"", round(C1.MCerror[i], digits = digit), "",
round(C1.HPD[i, 1], digits = digit), "", round(C1.HPD[i,
2], digits = digit)), file = test.file, append = TRUE)
}
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPosterior predictive value of Sensitivity of test under evaluation (S1) \t "),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Sensitivity ", round(S1_new.est, digits = digit),
"", round(S1_new.sd, digits = digit), "", round(S1_new.MCerror,
digits = digit), "", round(S1_new.HPD[1], digits = digit),
"", round(S1_new.HPD[2], digits = digit)), file = test.file,
append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste("\tPosterior predictive value of Specificity of test under evaluation (C1) \t "),
file = test.file, append = TRUE)
write(paste("______________________________________________________"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste(" Estimate Standard_Dev MC_error C.I._lower C.I._upper"),
file = test.file, append = TRUE)
write(paste(""), file = test.file, append = TRUE)
write(paste("Specificity", round(C1_new.est, digits = digit),
"", round(C1_new.sd, digits = digit), "", round(C1_new.MCerror,
digits = digit), "", round(C1_new.HPD[1], digits = digit),
"", round(C1_new.HPD[2], digits = digit)), file = test.file,
append = TRUE)
Num_study = c()
for (i in 1:N) {
Num_study = c(Num_study, paste("Study", i))
}
if (condInd == TRUE & Gold_Std == FALSE & model ==
1) {
if (Gold_se == TRUE & Gold_sp == FALSE) {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
3, 1), dimnames = list(1:rs.length, c(paste(point_estimate,
"estimate"), "HPD lower", "HPD upper"),
"C2"))
refstd_parameters[, 1, 1] <- C2.est
refstd_parameters[, 2, 1] <- C2.HPD[, 1]
refstd_parameters[, 3, 1] <- C2.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
1))
refstd_Parameters[, , 1] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[longueur,
i, 1], type = "n", col = 1, ylab = paste(param,
" of reference standard ", i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), i, 1],
col = l)
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, 1], type = "l", col = "grey",
ylab = paste(param, " of reference standard ",
i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = "Specificity"
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i, 1]),
lwd = 4, type = "l", col = "grey",
main = paste(param, " of reference standard ",
i, " \n Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin))
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 3, nrow = 1)
rownames(refstd_parameters) = c("C2")
colnames(refstd_parameters) = c(paste(point_estimate,
"estimate"), "HPD.low", "HPD.high")
refstd_parameters[1, 1] <- C2.est
refstd_parameters[1, 2] <- C2.HPD[1]
refstd_parameters[1, 3] <- C2.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 1)
colnames(refstd_Parameters) = c("C2")
refstd_Parameters[, 1] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Specificity of reference standard")
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
plot(x = longueur, y = refstd_Parameters[longueur,
1], type = "n", col = 1, ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), 1], col = l)
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Specificity of reference standard")
par(mfcol = c(5, 2))
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
1], type = "l", col = "grey", ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Specificity")
par(mfcol = c(5, 2))
longueur = 1:long
plot(density(refstd_Parameters[, 1]), lwd = 4,
type = "l", col = "grey", main = paste(param,
" of reference standard \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
dev.off()
}
}
}
else {
if (Gold_sp == TRUE & Gold_se == FALSE) {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
3, 1), dimnames = list(1:rs.length, c(paste(point_estimate,
"estimate"), "HPD lower", "HPD upper"),
"S2"))
refstd_parameters[, 1, 1] <- S2.est
refstd_parameters[, 2, 1] <- S2.HPD[, 1]
refstd_parameters[, 3, 1] <- S2.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
1))
refstd_Parameters[, , 1] <- S2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity")
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[longueur,
i, 1], type = "n", col = 1, ylab = paste(param,
" of reference standard ", i),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), i, 1],
col = l)
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity")
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, 1], type = "l", col = "grey",
ylab = paste(param, " of reference standard ",
i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity")
par(mfcol = c(5, 2))
longueur = 1:long
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i,
1]), lwd = 4, type = "l", col = "grey",
main = paste(param, " of reference standard ",
i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 3,
nrow = 1)
rownames(refstd_parameters) = c("S2")
colnames(refstd_parameters) = c(paste(point_estimate,
"estimate"), "HPD.low", "HPD.high")
refstd_parameters[1, 1] <- S2.est
refstd_parameters[1, 2] <- S2.HPD[1]
refstd_parameters[1, 3] <- S2.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 1)
colnames(refstd_Parameters) = c("S2")
refstd_Parameters[, 1] <- S2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity of reference standard")
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
plot(x = longueur, y = refstd_Parameters[longueur,
1], type = "n", col = 1, ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), 1], col = l)
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity of reference standard")
par(mfcol = c(5, 2))
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
1], type = "l", col = "grey", ylab = paste(param),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity of reference standard")
par(mfcol = c(5, 2))
longueur = 1:long
plot(density(refstd_Parameters[, 1]),
lwd = 4, type = "l", col = "grey",
main = paste(param, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
dev.off()
}
}
}
else {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
3, 2), dimnames = list(1:rs.length, c(paste(point_estimate,
"estimate"), "HPD lower", "HPD upper"),
c("S2", "C2")))
refstd_parameters[, 1, 1] <- S2.est
refstd_parameters[, 2, 1] <- S2.HPD[, 1]
refstd_parameters[, 3, 1] <- S2.HPD[, 2]
refstd_parameters[, 1, 2] <- C2.est
refstd_parameters[, 2, 2] <- C2.HPD[, 1]
refstd_parameters[, 3, 2] <- C2.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, c(long, rs.length,
2))
refstd_Parameters[, , 1] <- S2
refstd_Parameters[, , 2] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity")
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (j in 1:2) {
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[longueur,
i, j], type = "n", col = 1, ylab = paste(param[j],
" of reference standard ", i),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), i,
j], col = l)
}
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:2) {
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, j], type = "l", col = "grey",
ylab = paste(param[j], " of reference standard ",
i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:2) {
for (i in 1:rs.length) {
plot(density(refstd_Parameters[,
i, j]), lwd = 4, type = "l", col = "grey",
main = paste(param[j], " of reference standard ",
i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 3,
nrow = 2)
rownames(refstd_parameters) = c("S2", "C2")
colnames(refstd_parameters) = c(paste(point_estimate,
"estimate"), "HPD.low", "HPD.high")
refstd_parameters[1, 1] <- S2.est
refstd_parameters[1, 2] <- S2.HPD[1]
refstd_parameters[1, 3] <- S2.HPD[2]
refstd_parameters[2, 1] <- C2.est
refstd_parameters[2, 2] <- C2.HPD[1]
refstd_parameters[2, 3] <- C2.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 2)
colnames(refstd_Parameters) = c("S2", "C2")
refstd_Parameters[, 1] <- S2
refstd_Parameters[, 2] <- C2
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity of reference standard",
"Specificity of reference standard")
par(mfcol = c(5, 2))
no_chains = length(chain)
iter_chain = round((iter.num * nb_chains -
(burn_in) * nb_chains)/Thin, 0)/no_chains
longueur = 1:iter_chain
for (j in 1:2) {
plot(x = longueur, y = refstd_Parameters[longueur,
j], type = "n", col = 1, ylab = paste(param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
for (l in 1:length(chain)) {
lines(x = longueur, y = refstd_Parameters[longueur +
(iter_chain * (l - 1)), j], col = l)
}
}
}
else {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity of reference standard",
"Specificity of reference standard")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:2) {
plot(x = longueur, y = refstd_Parameters[,
j], type = "l", col = "grey", ylab = paste(param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity of reference standard",
"Specificity of reference standard")
par(mfcol = c(5, 2))
for (j in 1:2) {
longueur = 1:long
plot(density(refstd_Parameters[, j]),
lwd = 4, type = "l", col = "grey",
main = paste(param[j], " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
dev.off()
}
}
}
}
}
else {
if (condInd == TRUE & Gold_Std == FALSE & model ==
2) {
if (rs.length != 1) {
refstd_parameters = array(0, dim = c(rs.length,
3, 4), dimnames = list(1:rs.length, c(paste(point_estimate,
"estimate"), "HPD lower", "HPD upper"),
c("S2", "C2", "a1", "a0")))
refstd_parameters[, 1, 1] <- S2.est
refstd_parameters[, 2, 1] <- S2.HPD[, 1]
refstd_parameters[, 3, 1] <- S2.HPD[, 2]
refstd_parameters[, 1, 2] <- C2.est
refstd_parameters[, 2, 2] <- C2.HPD[, 1]
refstd_parameters[, 3, 2] <- C2.HPD[, 2]
refstd_parameters[, 1, 3] <- a1.est
refstd_parameters[, 2, 3] <- a1.HPD[, 1]
refstd_parameters[, 3, 3] <- a1.HPD[, 2]
refstd_parameters[, 1, 4] <- a0.est
refstd_parameters[, 2, 4] <- a0.HPD[, 1]
refstd_parameters[, 3, 4] <- a0.HPD[, 2]
long = length(alpha[, 1])
refstd_Parameters = array(0, dim = c(long,
rs.length, 4))
refstd_Parameters[, , 1] <- S2
refstd_Parameters[, , 2] <- C2
refstd_Parameters[, , 3] <- a1
refstd_Parameters[, , 4] <- a0
if (print_plot == TRUE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity",
"a1", "a0")
par(mfcol = c(5, 2))
for (j in 1:4) {
longueur = 1:long
for (i in 1:rs.length) {
plot(x = longueur, y = refstd_Parameters[,
i, j], type = "l", col = "grey",
ylab = paste(param[j], " of reference standard ",
i), xlab = "iteration number",
main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
abline(a = refstd_parameters[i, 2,
j], b = 0, col = "green", lwd = 3)
abline(a = refstd_parameters[i, 3,
j], b = 0, col = "green", lwd = 3)
}
}
dev.off()
file.pdf_RS2 = paste("RefStd density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity",
"a1", "a0")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:4) {
for (i in 1:rs.length) {
plot(density(refstd_Parameters[, i,
j]), lwd = 4, type = "l", col = "grey",
main = paste(param[j], " of reference standard ",
i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
}
}
dev.off()
}
}
else {
refstd_parameters = matrix(0, ncol = 3, nrow = 4)
rownames(refstd_parameters) = c("S2", "C2",
"a1", "a0")
colnames(refstd_parameters) = c(paste(point_estimate,
"estimate"), "HPD.low", "HPD.high")
refstd_parameters[1, 1] <- S2.est
refstd_parameters[1, 2] <- S2.HPD[1]
refstd_parameters[1, 3] <- S2.HPD[2]
refstd_parameters[2, 1] <- C2.est
refstd_parameters[2, 2] <- C2.HPD[1]
refstd_parameters[2, 3] <- C2.HPD[2]
refstd_parameters[3, 1] <- a1.est
refstd_parameters[3, 2] <- a1.HPD[1]
refstd_parameters[3, 3] <- a1.HPD[2]
refstd_parameters[4, 1] <- a0.est
refstd_parameters[4, 2] <- a0.HPD[1]
refstd_parameters[4, 3] <- a0.HPD[2]
long = length(THETA)
refstd_Parameters = matrix(0, nrow = long,
ncol = 4)
colnames(refstd_Parameters) = c("S2", "C2",
"a1", "a0")
refstd_Parameters[, 1] <- S2
refstd_Parameters[, 2] <- C2
refstd_Parameters[, 3] <- a1
refstd_Parameters[, 4] <- a0
if (print_plot == TRUE) {
file.pdf_RS = paste("RefStd trace plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS, paper = "a4", height = 20)
Param = c("Sensitivity of reference standard",
"Specificity of reference standard",
"a1", "a0")
par(mfcol = c(5, 2))
for (j in 1:4) {
longueur = 1:long
plot(x = longueur, y = refstd_Parameters[,
j], type = "l", col = "grey", ylab = paste(Param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) *
nb_chains)/Thin))
abline(a = refstd_parameters[j, 2], b = 0,
col = "green", lwd = 3)
abline(a = refstd_parameters[j, 3], b = 0,
col = "green", lwd = 3)
}
dev.off()
file.pdf_RS2 = paste("Density plots for N =",
round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0), ".pdf")
pdf(file.pdf_RS2, paper = "a4", height = 20)
param = c("Sensitivity", "Specificity",
"a1", "a0")
par(mfcol = c(5, 2))
for (j in 1:4) {
longueur = 1:long
plot(density(refstd_Parameters[, j]),
lwd = 4, type = "l", col = "grey",
main = paste(param[j], " of reference standard ",
j, " \n Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin))
}
dev.off()
}
}
}
}
parameters = array(0, dim = c(N, 3, 5), dimnames = list(Num_study,
c(paste(point_estimate, "estimate"), "HPD lower",
"HPD upper"), c("theta", "alpha", "pi", "S1",
"C1")))
parameters[, 1, 1] <- theta.est
parameters[, 2, 1] <- theta.HPD[, 1]
parameters[, 3, 1] <- theta.HPD[, 2]
parameters[, 1, 2] <- alpha.est
parameters[, 2, 2] <- alpha.HPD[, 1]
parameters[, 3, 2] <- alpha.HPD[, 2]
parameters[, 1, 3] <- PI.est
parameters[, 2, 3] <- PI.HPD[, 1]
parameters[, 3, 3] <- PI.HPD[, 2]
parameters[, 1, 4] <- S1.est
parameters[, 2, 4] <- S1.HPD[, 1]
parameters[, 3, 4] <- S1.HPD[, 2]
parameters[, 1, 5] <- C1.est
parameters[, 2, 5] <- C1.HPD[, 1]
parameters[, 3, 5] <- C1.HPD[, 2]
long = length(alpha[, 1])
Parameters = array(0, c(long, N, 5))
Parameters[, , 1] <- theta
Parameters[, , 2] <- alpha
Parameters[, , 3] <- PI
Parameters[, , 4] <- S1
Parameters[, , 5] <- C1
parameter = matrix(0, ncol = 3, nrow = 9)
rownames(parameter) = c("THETA", "LAMBDA", "beta",
"sigma.alpha", "sigma.theta", "S Overall", "C Overall",
"S1_new", "C1_new")
colnames(parameter) = c(paste(point_estimate, "estimate"),
"HPD.low", "HPD.high")
parameter[1, 1] <- THETA.est
parameter[1, 2] <- THETA.HPD[1]
parameter[1, 3] <- THETA.HPD[2]
parameter[2, 1] <- LAMBDA.est
parameter[2, 2] <- LAMBDA.HPD[1]
parameter[2, 3] <- LAMBDA.HPD[2]
parameter[3, 1] <- beta.est
parameter[3, 2] <- beta.HPD[1]
parameter[3, 3] <- beta.HPD[2]
parameter[4, 1] <- sigma.alpha.est
parameter[4, 2] <- sigma.alpha.HPD[1]
parameter[4, 3] <- sigma.alpha.HPD[2]
parameter[5, 1] <- sigma.theta.est
parameter[5, 2] <- sigma.theta.HPD[1]
parameter[5, 3] <- sigma.theta.HPD[2]
parameter[6, 1] <- S_overall.est
parameter[6, 2] <- S_overall.HPD[1]
parameter[6, 3] <- S_overall.HPD[2]
parameter[7, 1] <- C_overall.est
parameter[7, 2] <- C_overall.HPD[1]
parameter[7, 3] <- C_overall.HPD[2]
parameter[8, 1] <- S1_new.est
parameter[8, 2] <- S1_new.HPD[1]
parameter[8, 3] <- S1_new.HPD[2]
parameter[9, 1] <- C1_new.est
parameter[9, 2] <- C1_new.HPD[1]
parameter[8, 3] <- C1_new.HPD[2]
long = length(THETA)
Parameter = matrix(0, nrow = long, ncol = 9)
colnames(Parameter) = c("THETA", "LAMBDA", "beta",
"sigma.alpha", "sigma.theta", "S Overall", "C Overall",
"S1_new", "C1_new")
Parameter[, 1] <- THETA
Parameter[, 2] <- LAMBDA
Parameter[, 3] <- beta
Parameter[, 4] <- sigma.alpha
Parameter[, 5] <- sigma.theta
Parameter[, 6] <- S_overall
Parameter[, 7] <- C_overall
Parameter[, 8] <- S1_new
Parameter[, 9] <- C1_new
if (print_plot == TRUE) {
if (is.null(chain) == FALSE) {
no_chains = length(chain)
}
else {
if (is.null(chain) == TRUE) {
no_chains = 1
}
}
file.pdf5 = paste("Trace plots for N =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0),
".pdf")
pdf(file.pdf5, paper = "a4", height = 20)
param = c("theta", "alpha", "PI", "S1", "C1")
Param = c("Capital Theta", "Capital Lambda",
"beta", "~sigma[alpha]", "~sigma[theta]", "S Overall",
"C Overall", "S1_new", "C1_new")
iter_chain = round((iter.num * nb_chains - (burn_in) *
nb_chains)/Thin, 0)/no_chains
min_param = c(min(Parameters[, , 1]), min(Parameters[,
, 2]), min(Parameters[, , 3]), min(Parameters[,
, 4]), min(Parameters[, , 5]))
max_param = c(max(Parameters[, , 1]), max(Parameters[,
, 2]), max(Parameters[, , 3]), max(Parameters[,
, 4]), max(Parameters[, , 5]))
dlag = (max_param - min_param)/100
range_param = numeric()
for (j in 1:5) {
range_param = cbind(range_param, seq(min_param[j] +
dlag[j]/2, max_param[j] - dlag[j]/2, by = dlag[j]))
}
par(mfcol = c(5, 2))
longueur = 1:iter_chain
for (j in 1:5) {
for (i in 1:N) {
plot(x = longueur, y = Parameters[longueur,
i, j], type = "n", col = 1, ylab = paste(param[j],
" of study ", i), xlab = "iteration number",
main = paste("Thinning interval = ", thin.interval,
"\n Total samplesize kept = ", (iter.num *
nb_chains - (burn_in) * nb_chains)/Thin),
ylim = range(range_param[, j]))
for (l in 1:no_chains) {
lines(x = longueur, y = Parameters[longueur +
(iter_chain * (l - 1)), i, j], col = l)
}
}
}
min_Param = c(min(Parameter[, 1]), min(Parameter[,
2]), min(Parameter[, 3]), min(Parameter[, 4]),
min(Parameter[, 5]), min(Parameter[, 6]), min(Parameter[,
7]), min(Parameter[, 8]), min(Parameter[,
9]))
max_Param = c(max(Parameter[, 1]), max(Parameter[,
2]), max(Parameter[, 3]), max(Parameter[, 4]),
max(Parameter[, 5]), max(Parameter[, 6]), max(Parameter[,
7]), max(Parameter[, 8]), max(Parameter[,
9]))
dlag = (max_Param - min_Param)/100
range_Param = numeric()
for (j in 1:9) {
range_Param = cbind(range_Param, seq(min_Param[j] +
dlag[j]/2, max_Param[j] - dlag[j]/2, by = dlag[j]))
}
for (j in 1:9) {
plot(x = longueur, y = Parameter[longueur,
j], type = "n", col = 1, ylab = paste(Param[j]),
xlab = "iteration number", main = paste("Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin),
ylim = range(range_Param[, j]))
for (l in 1:no_chains) {
lines(x = longueur, y = Parameter[longueur +
(iter_chain * (l - 1)), j], col = l)
}
}
dev.off()
file.pdf3 = paste("Density plots for N =", round((iter.num *
nb_chains - (burn_in) * nb_chains)/Thin, 0),
".pdf")
pdf(file.pdf3, paper = "a4", height = 20)
param = c("theta", "alpha", "PI", "S1", "C1")
Param = c("Capital Theta", "Capital Lambda",
"beta", "~sigma[alpha]", "~sigma[theta]", "S Overall",
"C Overall", "S1_new", "C1_new")
par(mfcol = c(5, 2))
longueur = 1:long
for (j in 1:5) {
for (i in 1:N) {
plot(density(Parameters[, i, j]), lwd = 4,
type = "l", col = "grey", main = paste(param[j],
" of study ", i, " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin))
}
}
for (j in 1:9) {
plot(density(Parameter[, j]), lwd = 4, type = "l",
col = "grey", main = paste(Param[j], " \n Thinning interval = ",
thin.interval, "\n Total samplesize kept = ",
(iter.num * nb_chains - (burn_in) * nb_chains)/Thin))
}
dev.off()
Sensi1 = apply(as.matrix(Parameters[, , 4]),
2, median)
Speci1 = apply(as.matrix(Parameters[, , 5]),
2, median)
Ov_Se = 1 - pnorm((median(Parameter[, 1]) - median(Parameter[,
2])/2)/exp(median(Parameter[, 3])/2))
Ov_Sp = pnorm((median(Parameter[, 1]) + median(Parameter[,
2])/2)/exp(-median(Parameter[, 3])/2))
thet = qnorm((1 - as.matrix(Parameters[, , 4])) +
1e-14) * exp(Parameter[, 3]/2) + Parameter[,
2]/2
min_TH = quantile(thet, trunc_low)
max_TH = quantile(thet, (1 - trunc_up))
dTH = 5e-05
TH_range = seq(min_TH + dTH/2, max_TH - dTH/2,
dTH)
S_sroc = 1 - pnorm((TH_range - median(Parameter[,
2])/2)/exp(median(Parameter[, 3])/2))
C_sroc = pnorm((TH_range + median(Parameter[,
2])/2)/exp(-median(Parameter[, 3])/2))
if (cred_region == TRUE) {
min_t = 0
max_t = 2 * pi
dTt = 5e-05
range_t = seq(min_t + dTt/2, max_t - dTt/2,
dTt)
bound_cte = sqrt(qchisq(1 - region_level, 2,
lower.tail = FALSE))
hat_mu_A = exp(-median(Parameter[, 3])/2) *
(-median(Parameter[, 1]) + median(Parameter[,
2])/2)
s_A = sd(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2))
cos_fun_A = cos(range_t)
hat_mu_B = exp(median(Parameter[, 3])/2) *
(median(Parameter[, 1]) + median(Parameter[,
2])/2)
s_B = sd(exp(Parameter[, 3]/2) * (Parameter[,
1] + Parameter[, 2]/2))
r = cor(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2), exp(Parameter[, 3]/2) *
(Parameter[, 1] + Parameter[, 2]/2))
cos_fun_B = cos(range_t + acos(r))
probit_S_credible = hat_mu_A + s_A * bound_cte *
cos_fun_A
probit_C_credible = hat_mu_B + s_B * bound_cte *
cos_fun_B
S_credible = pnorm(probit_S_credible)
C_credible = pnorm(probit_C_credible)
}
if (predict_region == TRUE) {
min_t = 0
max_t = 2 * pi
dTt = 5e-05
range_t = seq(min_t + dTt/2, max_t - dTt/2,
dTt)
bound_cte = sqrt(qchisq(1 - region_level, 2,
lower.tail = FALSE))
hat_mu_A = exp(-median(Parameter[, 3])/2) *
(-median(Parameter[, 1]) + median(Parameter[,
2])/2)
s_A = sd(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2))
hat_sigma_A = sqrt((exp(-median(Parameter[,
3])) * (median((Parameter[, 5])^2) + 0.25 *
median((Parameter[, 4])^2))))
cos_fun_A = cos(range_t)
hat_mu_B = exp(median(Parameter[, 3])/2) *
(median(Parameter[, 1]) + median(Parameter[,
2])/2)
s_B = sd(exp(Parameter[, 3]/2) * (Parameter[,
1] + Parameter[, 2]/2))
hat_sigma_B = sqrt((exp(median(Parameter[,
3])) * (median((Parameter[, 5])^2) + 0.25 *
median((Parameter[, 4])^2))))
r = cor(exp(-Parameter[, 3]/2) * (-Parameter[,
1] + Parameter[, 2]/2), exp(Parameter[, 3]/2) *
(Parameter[, 1] + Parameter[, 2]/2))
hat_sigma_AB = -median((Parameter[, 5])^2) +
0.25 * median((Parameter[, 4])^2)
cos_fun_B = cos(range_t + acos(r))
probit_S_prediction = hat_mu_A + (sqrt(s_A^2 +
hat_sigma_A^2)) * bound_cte * cos_fun_A
probit_C_prediction = hat_mu_B + (sqrt(s_B^2 +
hat_sigma_B^2)) * bound_cte * cos(range_t +
acos(r))
S_prediction = pnorm(probit_S_prediction)
C_prediction = pnorm(probit_C_prediction)
}
pdf("Summary ROC curve.pdf")
default.x = range(1, 0)
default.y = range(0, 1)
plot(x = default.x, y = default.y, type = "n",
xlim = rev(range(default.x)), xlab = "", ylab = "")
title(xlab = "Specificity", ylab = "Sensitivity",
cex.lab = 1.5, main = "Summary ROC curve")
if (plot.ind.studies == TRUE) {
if (Gold_Std == TRUE) {
Scale_factor = 10
SENSi1 = data[[1]][, 1]/(data[[1]][, 1] +
data[[1]][, 3])
SPECi1 = data[[1]][, 4]/(data[[1]][, 4] +
data[[1]][, 2])
symbols(SPECi1, SENSi1, circles = rowSums(as.matrix(data[[1]])),
inches = 0.1 * Scale_factor/7, add = TRUE,
fg = 1)
}
else {
Scale_factor = 10
symbols(Speci1, Sensi1, circles = rowSums(as.matrix(data[[1]])),
inches = 0.1 * Scale_factor/7, add = TRUE,
fg = 1)
}
}
if (cred_region == TRUE) {
lines(C_credible, S_credible, lwd = 3, lty = lty.cred.region,
col = col.pooled.estimate)
}
if (predict_region == TRUE) {
lines(C_prediction, S_prediction, lwd = 3,
lty = lty.predict.region, col = col.predict.region)
}
lines(C_sroc, S_sroc, lwd = 3, col = "black",
lty = 1)
points(Ov_Sp, Ov_Se, pch = 19, cex = 2.5, col = col.pooled.estimate)
dev.off()
}
}
}
if (Gold_Std == FALSE) {
Results = list(parameter, parameters, refstd_parameters,
paste("See '", getwd(), "' for complete results",
sep = ""))
names(Results) = c("Between-study parameters", "Within-study parameters",
"Reference standard", "")
}
else {
Results = list(parameter, parameters, paste("See '",
getwd(), "' for complete results", sep = ""))
names(Results) = c("Between-study parameters", "Within-study parameters",
"")
}
return(Results)
cat(paste("See \"", getwd(), "\" for complete results", sep = ""))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.