Nothing
# R package: generalized meta-analysis (gmeta), file gmeta.r
# R interface to meta-analysis based on combine confidence distributions (CDs)
# *****************************************************************************
# main
# *****************************************************************************
## plot functions
plot.gmeta <- function(gmo, studies=NULL,
plot.option=c('confidence-density',
'confidence-curve', 'cv',
'confidence-distribution', 'cdf'),
type='l', xlab='x', ylab='density', xlim=NULL, ylim=NULL, ...) {
UseMethod("plot", gmo);
}
###plot functions - model based meta-anlaysis
plot.gmeta.m <- function(x, studies=NULL,
plot.option=c('confidence-density',
'confidence-curve', 'cv',
'confidence-distribution', 'cdf'),
type='l', xlab='x', ylab='density', xlim=NULL, ylim=NULL, ...) {
# match plot.option
mfplot <- match.call()
plot.option = match.arg(plot.option)
# plot.individual.studies
# if ( is.null(studies) ) {
# # take all studies
# gmi <- x$input
# nn1 <- ifelse(is.list(gmi), length(gmi), dim(gmi)[1])
# studies <- c(1:nn1)
# }
# plot via plot.option
if ( plot.option == 'confidence-density' ) {
gmeta.plot.cdd(x, studies, type, xlab, ylab, xlim, ylim, ...)
} else if ( plot.option == 'cv' || plot.option == 'confidence-curve' ) {
if ( ylab == 'density' ) {
ylab = 'confidence curve'
}
gmeta.plot.cvs(x, studies, type, xlab, ylab, xlim, ylim, ...)
} else if ( plot.option == 'cdf' || plot.option == 'confidence-distribution' ) {
if ( ylab == 'density' ) {
ylab = 'confidence distribution'
}
gmeta.plot.cdf(x, studies, type, xlab, ylab, xlim, ylim, ...)
} else {
stop('plot.option not recognize')
}
}
###plot functions - combine evidence from 2x2 tables
plot.gmeta.e <- function(x, studies=NULL,
plot.option=c('confidence-density',
'confidence-curve', 'cv',
'confidence-distribution', 'cdf'),
type='l', xlab='x', ylab='confidence density', xlim=NULL, ylim=NULL, ...) {
# match plot.option
mfplot <- match.call()
plot.option = match.arg(plot.option)
# plot.individual.studies
# if ( is.null(studies) ) {
# # take all studies
# gmi <- x$input
# nn1 <- ifelse(is.list(gmi), length(gmi), dim(gmi)[1])
# studies <- c(1:nn1)
# }
# take only non-zero/zero-event studies
idx <- studies[is.na(x$individual.mean)[studies]]
studies <- studies[!is.na(x$individual.mean)[studies]]
if ( length(idx) != 0 ) {
for (i in idx) {
cat('The confidence distribution of study', i, 'cannot be plot because it contains zero-zero events.','\n')
}
}
# plot via plot.option
if ( plot.option == 'confidence-density' ) {
gmeta.plot.cdd(x, studies, type, xlab, ylab, xlim, ylim, ...)
} else if ( plot.option == 'cv' || plot.option == 'confidence-curve' ) {
if ( ylab == 'confidence density' ) {
ylab = 'confidence curve'
}
gmeta.plot.cvs(x, studies, type, xlab, ylab, xlim, ylim, ...)
} else if ( plot.option == 'cdf' || plot.option == 'confidence-distribution' ) {
if ( ylab == 'confidence density' ) {
ylab = 'confidence distribution'
}
gmeta.plot.cdf(x, studies, type, xlab, ylab, xlim, ylim, ...)
} else {
stop('plot.option not recognize')
}
}
###plot functions - shared by model based and exact methods on 2x2 tables
#####
gmeta.plot.cdd <- function(x, studies, type, xlab, ylab, xlim, ylim, ...) {
#e <- studies
# number of layers
nn <- length(studies) + 1 # individual studies + combined CD
# extract study names
enames <- c(x$study.names[studies], 'combined.density')
# x.grids
x.grids <- x$x.grids
# individual CDs
ecds <- x$individual.cds[studies, ]
# individual densities
edns <- NULL
for (i in studies) {
edns <- rbind(edns, F2f(x.grids, x$individual.cds[i, ]))
}
# individual and combined density
edns <- rbind(edns, x$combined.density)
# unified y-range on all layers for visual comparison
if ( !is.null(studies) ) {
if ( max(edns, na.rm=T) > 1 ) {
edns <- edns / max(edns, na.rm=T)
}
}
# individual and combined medians
mdn <- c(x$individual.medians[studies], x$combined.median)
# individual and combined confidence intervals
emdncis <- rbind(x$individual.cis[studies,], x$combined.ci)
# set lower/upper bound
if ( !is.null(xlim) ) {
xl <- min(xlim)
xu <- max(xlim)
} else {
xl <- min(x.grids)
xu <- max(x.grids)
}
# set study.names position
xnames <- xl
# plot
if ( !is.null(studies) ) {
plot(0, type='n', xlab=xlab, ylab=ylab, xlim=c(xl,xu), ylim=c(0,nn), yaxt='n')
} else {
plot(0, type='n', xlab=xlab, ylab=ylab, xlim=c(xl,xu), ylim=c(0,ceiling(max(edns))), yaxt='n')
}
for ( i in 1:nn ) {
# supporting line
abline(h=i-1, lty=2)
# plot of density
lines(x.grids, edns[i, ]+i-1, lty=1)
# mark median, ci.lower, ci.upper
points(mdn[i], i-1, cex=1, col='dark red')
points(emdncis[i,1], i-1, pch='[', cex=1, col='dark red')
points(emdncis[i,2], i-1, pch=']', cex=1, col='dark red')
# legend study.names
legend(xnames, i-0.5, legend=paste('#', enames[i], sep=''), bty='n')
}
# overall median
abline(v=x$combined.median, lty=2)
}
#####
gmeta.plot.cvs <- function(gmo, studies, type, xlab, ylab, xlim, ylim, ...) {
# number of layers
nn <- length(studies) + 1 # individual studies + combined CD
# extract study names
enames <- c(gmo$study.names[studies], 'combined.cv')
# x.grids
x.grids <- gmo$x.grids
# construct individual and combined CVs
combined.cv <- 1 - 2 * abs( gmo$combined.cd - 0.5 )
individual.cvs <- 1 - 2 * abs( gmo$individual.cds - 0.5 )
# set cvs
ecvs <- rbind(individual.cvs[studies,], combined.cv)
# individual and combined medians
mdn <- c(gmo$individual.medians[studies], gmo$combined.median)
# individual and combined confidence intervals
emdncis <- rbind(gmo$individual.cis[studies,], gmo$combined.ci)
# set lower/upper bound
if ( !is.null(xlim) ) {
xl <- min(xlim)
xu <- max(xlim)
} else {
xl <- min(x.grids)
xu <- max(x.grids)
}
# set study.names position
xnames <- xl
# plot
if ( !is.null(studies) ) {
plot(0, type='n', xlab=xlab, ylab=ylab, xlim=c(xl,xu), ylim=c(0,nn), yaxt='n')
} else {
plot(0, type='n', xlab=xlab, ylab=ylab, xlim=c(xl,xu), ylim=c(0,ceiling(max(ecvs))), yaxt='n')
}
for ( i in 1:nn ) {
# supporting line
abline(h=i-1, lty=2)
# plot of density
lines(x.grids, ecvs[i, ]+i-1, lty=1)
# mark median, ci.lower, ci.upper
points(mdn[i], i-1, cex=1, col='dark red')
points(emdncis[i,1], i-1, pch='[', cex=1, col='dark red')
points(emdncis[i,2], i-1, pch=']', cex=1, col='dark red')
# legend study.names
legend(xnames, i-0.5, legend=paste('#', enames[i], sep=''), bty='n')
}
# overall median
abline(v=gmo$combined.median, lty=2)
}
#####
gmeta.plot.cdf <- function(gmo, studies, type, xlab, ylab, xlim, ylim, ...) {
#e <- studies
# number of layers
nn <- length(studies) + 1 # individual studies + combined CD
# extract study names
enames <- c(gmo$study.names[studies], 'combined.cdf')
# x.grids
x.grids <- gmo$x.grids
# individual CDs
ecds <- gmo$individual.cds[studies, ]
# individual and combined CDs
ecds <- rbind(ecds, gmo$combined.cd)
# unified y-range on all layers for visual comparison
if ( !is.null(studies) ) {
if ( max(ecds, na.rm=T) > 1 ) {
ecds <- ecds / max(ecds, na.rm=T)
}
}
# individual and combined medians
mdn <- c(gmo$individual.medians[studies], gmo$combined.median)
# individual and combined confidence intervals
emdncis <- rbind(gmo$individual.cis[studies,], gmo$combined.ci)
# set lower/upper bound
if ( !is.null(xlim) ) {
xl <- min(xlim)
xu <- max(xlim)
} else {
xl <- min(x.grids)
xu <- max(x.grids)
}
# set study.names position
xnames <- xl
# plot
if ( !is.null(studies) ) {
plot(0, type='n', xlab=xlab, ylab=ylab, xlim=c(xl,xu), ylim=c(0,nn), yaxt='n')
} else {
plot(0, type='n', xlab=xlab, ylab=ylab, xlim=c(xl,xu), ylim=c(0,ceiling(max(ecds))), yaxt='n')
}
for ( i in 1:nn ) {
# supporting line
abline(h=i-1, lty=2)
# plot of CDs
lines(x.grids, ecds[i, ]+i-1, lty=1)
# mark median, ci.lower, ci.upper
points(mdn[i], i-1, cex=1, col='dark red')
points(emdncis[i,1], i-1, pch='[', cex=1, col='dark red')
points(emdncis[i,2], i-1, pch=']', cex=1, col='dark red')
# legend study.names
legend(xnames, i-0.5, legend=paste('#', enames[i], sep=''), bty='n')
}
# overall median
abline(v=gmo$combined.median, lty=2)
}
# postprocessing [done]
# *****************************************************************************
# other functions used in this package
# *****************************************************************************
# other functions used in this package
### using CD to make inference
### get mean, median, sd, ci, given cd, used everywhere
##### get mean from cd
gmeta.cd.mean <- function(x, cd) {
x = x[!is.na(cd)]
cd = cd[!is.na(cd)]
if ( length(unique(cd)) == 1 ) {
mn <- NA
} else if (length(x) == length(cd)) {
dn <- F2f(x,cd) # density
mn <- sum(x*dn, na.rm=T) / sum(dn, na.rm=T) # mean
} else {
stop('length of x is not the same as length of cd.')
}
return(mn)
}
##### get median from cd
gmeta.cd.median <- function(x, cd) {
x = x[!is.na(cd)]
cd = cd[!is.na(cd)]
if ( length(unique(cd)) == 1 ) {
mdn <- NA
} else {
# same as approx in
# gmeta.cd.mdncis()
x1 <- rev(x[cd<0.5])[1]
x2 <- x[cd>0.5][1]
y1 <- rev(cd[cd<0.5])[1]
y2 <- cd[cd>0.5][1]
# linear interpolation
mdn <- x1+(x2-x1)/(y2-y1)*(0.5-y1)
}
return(mdn)
}
##### get stddev from cd
gmeta.cd.stddev <- function(x, cd) {
x = x[!is.na(cd)]
cd = cd[!is.na(cd)]
if ( length(unique(cd)) == 1 ) {
sdv <- NA
} else {
sdv <- diff(approx(x=cd, y=x, xout=c(0.25,0.75), ties='mean')$y) / (qnorm(0.75) - qnorm(0.25))
}
return(sdv)
}
##### get ci&median from cd
gmeta.cd.mdncis <- function(x, cd, alpha) {
x = x[!is.na(cd)]
cd = cd[!is.na(cd)]
if ( length(unique(cd)) == 1 ) {
mdncis <- c(NA, NA, NA)
} else {
mdncis <- approx(x=cd, y=x, xout=c(alpha/2, 0.5 ,1-alpha/2), ties='mean')$y
}
return(mdncis)
}
### computing density(pdf) from distribution(cdf), used everywhere
### mainly used for getting CD density from CD distribution, NA on two end to avoid sudden jump
F2f<-function (x, Fx) {
fx = diff(Fx, lag = 2)/diff(x, lag = 2)
fx = c(NA, fx, NA) # match the length of x-fx-cdf
# n.xgrids
n = length(fx)
# liner extrapolation
fx[1] = max(0, (fx[2]-fx[3])/(x[2]-x[3]) * (x[1]-x[2]) + fx[2])
fx[n] = max(0, (fx[n-1]-fx[n-2])/(x[n-1]-x[n-2]) * (x[n]-x[n-1]) + fx[n-1])
# return
return(fx)
#[update in v2.0: previously #return(c(NA, fx, NA)) # NA so that the edge of the plot of density is right.]
}
## [p-value combination tools]
## [update: p-value combination is implemented in C under this version]
## [p-value combination tools - done]
## [model based methods tools]
### computing tau2
### used in model based meta-analysis, random-effects models, estimate of tau2
##### Gtau2 - robust heterogeneity tau2 function
Gtau2 <- function(theta, sigma, method) {
# calculate tau2 by method
if ( method == 'random-mm' ) {
tau2 <- GtauDL2(theta, sigma)
}
else if ( method == 'random-reml' ) {
tau2 <- GtauREML2(theta, sigma)
}
else if ( method == 'random-robust1' ) {
tau2 <- GtauREML2(theta, sigma)
}
else if ( method == 'random-robust2' ) {
tau2 <- GtauREML2(theta, sigma)
}
else if ( method == 'random-robust2(sqrt12)' ) {
tau2 <- GtauREML2(theta, sigma)
}
else {
tau2 <- 0 # for fixed effects models
#cat('\nmethod not recongize\n') #commentOut
}
# return
return(tau2)
}
##### Gtau2m() methods for tau2
Gtau2m <- function(theta, sigma, tau2method) {
# calculate tau2 by method-name
if ( tau2method == 'DL' ) {
tau2 <- GtauDL2(theta, sigma)
}
else if ( tau2method == 'HS' ) {
tau2 <- GtauHS2(theta, sigma)
}
else if ( tau2method == 'SJ' ) {
tau2 <- GtauSJ2(theta, sigma)
}
else if ( tau2method == 'HE' ) {
tau2 <- GtauHE2(theta, sigma)
}
else if ( tau2method == 'ML' ) {
tau2 <- GtauML2(theta, sigma)
}
else if ( tau2method == 'REML' ) {
tau2 <- GtauREML2(theta, sigma)
}
else if ( tau2method == 'EB' ) {
tau2 <- GtauEB2(theta, sigma)
}
else {
stop('tau2method not recongize.')
}
# return
return(tau2)
}
##### Get tauDL2
# ?update: need robust version
GtauDL2 <- function(theta, sigma) {
nn1 = length(theta)
wis <- 1 / sigma^2
theta.mle <- sum(theta * wis) / sum(wis)
q.w <- sum(wis * (theta - theta.mle)^2)
tauDL2 <- max(0, (q.w - (nn1-1) ) / (sum(wis) - sum(wis^2)/sum(wis)) )
return(tauDL2)
}
##### Get tauHS2
GtauHS2 <- function(theta, sigma) {
nn1 = length(theta)
wis <- 1 / sigma^2
theta.mle <- sum(theta * wis) / sum(wis)
q.w <- sum(wis * (theta - theta.mle)^2)
tauHS2 <- max(0, (q.w - nn1) / sum(wis))
return(tauHS2)
}
##### Get tauSJ2
GtauSJ2 <- function(theta, sigma) {
nn1 = length(theta)
tauSJ2 <- var(theta) * (nn1 - 1)/nn1
wis <- 1 / (sigma^2 + tauSJ2)
theta.mle <- sum(theta * wis) / sum(wis)
q.w <- sum(wis * (theta - theta.mle)^2)
tauSJ2 <- max(0, tauSJ2 * q.w / (nn1 - 1))
return(tauSJ2)
}
##### Get tauHE2
GtauHE2 <- function(theta, sigma) {
nn1 = length(theta)
wis <- rep(1, nn1)
theta.mle <- sum(theta * wis) / sum(wis)
q.w <- sum(wis * (theta - theta.mle)^2)
adj.w <- sum(sigma^2) * (nn1 - 1) / nn1
tauHE2 <- max(0, (q.w - adj.w) / (nn1 - 1))
return(tauHE2)
}
##### Get tauML2
GtauML2 <- function(theta, sigma) {
nn1 = length(theta)
# initialization
tauML2 = max(0, GtauHE2(theta, sigma))
# update
nloop = 0
absch = 1 # absolute change in tauML2 value
while ( absch > 10^(-5) ) {
nloop = nloop + 1
if ( nloop > 10^5 ) {
stop("tauML2 via MLE method does not converge.")
}
else {
tauML2O <- tauML2 # tauML2Old
# update wML
wML <- 1 / (sigma^2 + tauML2)
WML <- diag(wML) - wML %o% wML / sum(wML)
# update tauML2
adj.tauML2 <- (t(theta) %*% WML %*% WML %*% theta - sum(wML)) / sum(wML^2)
while ( tauML2 + adj.tauML2 < 0 ) {
adj.tauML2 <- adj.tauML2 / 2
}
tauML2 <- tauML2 + adj.tauML2
absch <- abs(tauML2 - tauML2O)
}
}
# return
return(max(0, tauML2))
}
##### Get tauREML2
GtauREML2 <- function(theta, sigma) {
nn1 = length(theta)
# adopt tauHE2 as initial tauREML2
tauR2 <- GtauHE2(theta, sigma)
# update
nloop = 0
absch = 1 # absolute change in tauR2 value
while ( absch > 10^(-5) ) {
nloop = nloop + 1
if ( nloop > 10^5 ) {
stop("tauREML2 via REML method does not converge.")
}
else {
tauR2O <- tauR2 # tauR2Old
# update thetaR, wR
wR <- 1/(sigma^2 + tauR2O)
thetaR <- sum(wR*theta) / sum(wR)
# update tauR
tauR2 <- sum(wR^2*(nn1/(nn1-1)*(theta-thetaR)^2 - sigma^2)) / sum(wR^2)
absch <- abs(tauR2 - tauR2O)
}
}
# return
return(max(0, tauR2))
}
##### Get tauEB2
GtauEB2 <- function(theta, sigma) {
nn1 = length(theta)
# initialization
tauEB2 <- max(0, GtauHE2(theta, sigma))
# update
nloop = 0
absch = 1 # absolute change in tauEB2 value
while ( absch > 10^(-5) ) {
nloop = nloop + 1
if ( nloop > 10^5 ) {
stop("tauEB2 via EB method does not converge.")
}
else {
tauEB2O <- tauEB2 # tauEB2Old
# update thetaEB, wEB
wEB <- 1 / (sigma^2 + tauEB2)
thetaEB <- sum(theta * wEB) / sum(wEB)
adj.w <- sum(wEB * (theta - thetaEB)^2)
# update tauEB2
adj.tauEB2 <- 1/sum(wEB)*(adj.w*nn1/(nn1-1)-nn1)
while (tauEB2 + adj.tauEB2 < 0) {
adj.tauEB2 <- adj.tauEB2 / 2
}
tauEB2 <- tauEB2 + adj.tauEB2
absch <- abs(tauEB2 - tauEB2O)
}
}
# return
return(max(0, tauEB2))
}
### compute convolution distribution, used in double exp
convolution.sim <- function(rF0, weight, m=1e5) {
k = length(weight)
yy = numeric(m)
for (i in c(1:m)) {
yy[i] <- sum(rF0(k)*weight)
}
return(ecdf(yy))
}
### compute convolution distribution
### update in v2.0 - exact distribution functions
.convolution.de <- function(rF0, weight, verbose, m=1e5) {
k = length(weight)
if ( k <= 41 && all( weight == 1 ) ) {
# return
convolutionDEfunc <- .convolution.de.vfunc(k)
} else {
# verbose
if ( verbose ) {
if ( k > 41 ) {
cat('\nlarge number of study, convolution distribution obtained by simulation.\n') #commentOut
}
if ( ! all( weight == 1 ) ) {
cat('\nwarning: NOT all one weight, convolution distribution obtained by simulation.\n') #commentOut
}
}
# return
convolutionDEfunc <- convolution.sim(rF0, weight, m=m)
}
# return
return( convolutionDEfunc )
}
# exact DE convolution disribution function k=1-41
.convolution.de.vfunc <- function(k) {
# v1
v1 <- function(x) {
return( 1 )
}
# v2
v2 <- function(x) {
return( x/2 + 1 )
}
# v3
v3 <- function(x) {
return( x^2/8 + (5*x)/8 + 1 )
}
# v4
v4 <- function(x) {
return( x^3/48 + (3*x^2)/16 + (11*x)/16 + 1 )
}
# v5
v5 <- function(x) {
return( x^4/384 + (7*x^3)/192 + (29*x^2)/128 + (93*x)/128 + 1 )
}
# v6
v6 <- function(x) {
return( x^5/3840 + x^4/192 + (37*x^3)/768 + (65*x^2)/256 + (193*x)/256 + 1 )
}
# v7
v7 <- function(x) {
return( x^6/46080 + (3*x^5)/5120 + (23*x^4)/3072 + (11*x^3)/192 + (281*x^2)/1024 + (793*x)/1024 + 1 )
}
# v8
v8 <- function(x) {
return( x^7/645120 + x^6/18432 + (7*x^5)/7680 + (29*x^4)/3072 + (397*x^3)/6144 + (595*x^2)/2048 + (1619*x)/2048 + 1 )
}
# v9
v9 <- function(x) {
return( x^8/10321920 + (11*x^7)/2580480 + (67*x^6)/737280 + (299*x^5)/245760 + (1093*x^4)/98304 + (3473*x^3)/49152 + (9949*x^2)/32768 + (26333*x)/32768 + 1 )
}
# v10
v10 <- function(x) {
return( x^9/185794560 + x^8/3440640 + (79*x^7)/10321920 + (21*x^6)/163840 + (1471*x^5)/983040 + (103*x^4)/8192 + (14893*x^3)/196608 + (20613*x^2)/65536 + (53381*x)/65536 + 1 )
}
# v11
v11 <- function(x) {
return( x^10/3715891200 + (13*x^9)/743178240 + (23*x^8)/41287680 + (47*x^7)/4128768 + (647*x^6)/3932160 + (459*x^5)/262144 + (10889*x^4)/786432 + (15751*x^3)/196608 + (84883*x^2)/262144 + (215955*x)/262144 + 1 )
}
# v12
v12 <- function(x) {
return( x^11/81749606400 + x^10/1061683200 + (53*x^9)/1486356480 + x^8/1146880 + (839*x^7)/55050240 + (1567*x^6)/7864320 + (1559*x^5)/786432 + (11773*x^4)/786432 + (131975*x^3)/1572864 + (173965*x^2)/524288 + (436109*x)/524288 + 1 )
}
# v13
v13 <- function(x) {
return( (x^12 + 90*x^11 + 3993*x^10 + 115005*x^9 + 2386395*x^8 + 37469520*x^7 + 455259420*x^6 + 4302906300*x^5 + 31335467625*x^4 + 171172905750*x^3 + 664761133575*x^2 + 1645756410375*x + 1961990553600)/1961990553600 )
}
# v14
v14 <- function(x) {
return( (x^13 + 104*x^12 + 5343*x^11 + 178893*x^10 + 4341480*x^9 + 80424630*x^8 + 1167180300*x^7 + 13408094700*x^6 + 121696499925*x^5 + 860553193500*x^4 + 4601737965825*x^3 + 17600023616175*x^2 + 43105900812975*x + 51011754393600)/51011754393600 )
}
# v15
v15 <- function(x) {
return( (x^14 + 119*x^13 + 7007*x^12 + 269724*x^11 + 7561554*x^10 + 162912750*x^9 + 2775672900*x^8 + 37918881000*x^7 + 416674583325*x^6 + 3659572691775*x^5 + 25255014609825*x^4 + 132643472761800*x^3 + 500706514833525*x^2 + 1214871076343925*x + 1428329123020800)/1428329123020800 )
}
# v16
v16 <- function(x) {
return( (x^15 + 135*x^14 + 9030*x^13 + 395850*x^12 + 12686310*x^11 + 314143830*x^10 + 6196840650*x^9 + 98983684800*x^8 + 1288808846325*x^7 + 13659762691575*x^6 + 116744331904200*x^5 + 789273852617250*x^4 + 4082080279402125*x^3 + 15234653491682625*x^2 + 36659590336994625*x + 42849873690624000)/42849873690624000 )
}
# v17
v17 <- function(x) {
return( (x^16 + 152*x^15 + 11460*x^14 + 567420*x^13 + 20603310*x^12 + 580556340*x^11 + 13108004910*x^10 + 241511019750*x^9 + 3664417281525*x^8 + 45879983849700*x^7 + 471898161885150*x^6 + 3941370814030650*x^5 + 26181748152685125*x^4 + 133614981594344250*x^3 + 493699195087473375*x^2 + 1179297174137457375*x + 1371195958099968000)/1371195958099968000 )
}
# v18
v18 <- function(x) {
return( (x^17 + 170*x^16 + 14348*x^15 + 796620*x^14 + 32519130*x^13 + 1033829160*x^12 + 26460800730*x^11 + 556103137590*x^10 + 9702192775275*x^9 + 141154833169350*x^8 + 1710657725827050*x^7 + 17154519346814850*x^6 + 140481501759573975*x^5 + 919067426174898000*x^4 + 4635763624512145125*x^3 + 16977671416936605375*x^2 + 40288002704636061375*x + 46620662575398912000)/46620662575398912000 )
}
# v19
v19 <- function(x) {
return( (x^18 + 189*x^17 + 17748*x^16 + 1097928*x^15 + 50044770*x^14 + 1781769150*x^13 + 51272700570*x^12 + 1217623155840*x^11 + 24160874352615*x^10 + 403114038101775*x^9 + 5662993054568850*x^8 + 66763593395799300*x^7 + 655117082164019475*x^6 + 5273993980721691225*x^5 + 34045921262108881125*x^4 + 169957871025837394500*x^3 + 617528830880480644125*x^2 + 1456700757237661060125*x + 1678343852714360832000)/1678343852714360832000 )
}
# v20
v20 <- function(x) {
return( (x^19 + 209*x^18 + 21717*x^17 + 1488384*x^16 + 75297114*x^15 + 2982843630*x^14 + 95816929320*x^13 + 2550713370660*x^12 + 57036699560295*x^11 + 1079618519974995*x^10 + 17353300159520325*x^9 + 236653385032864800*x^8 + 2724788477433797775*x^7 + 26237740609970314425*x^6 + 208087722625924691550*x^5 + 1327519193937539352750*x^4 + 6566054316784789451625*x^3 + 23687738668934964248625*x^2 + 55576271870507820056625*x + 63777066403145711616000)/63777066403145711616000 )
}
# v21
v21 <- function(x) {
return( (x^20 + 230*x^19 + 26315*x^18 + 1987875*x^17 + 111018330*x^16 + 4865271480*x^15 + 173370863700*x^14 + 5137770462300*x^13 + 128456673938775*x^12 + 2733682807223550*x^11 + 49741855758770175*x^10 + 774605689977994875*x^9 + 10297696798485471375*x^8 + 116155760365285641000*x^7 + 1100170903364915382000*x^6 + 8610589485844903557000*x^5 + 54356745298536206150625*x^4 + 266631748389972173958750*x^3 + 955710341290036461504375*x^2 + 2231251669352950693824375*x + 2551082656125828464640000)/2551082656125828464640000 )
}
# v22
v22 <- function(x) {
return( (x^21 + 252*x^20 + 31605*x^19 + 2619435*x^18 + 160715205*x^17 + 7751748060*x^16 + 304733193660*x^15 + 9992154645900*x^14 + 277452017345475*x^13 + 6587383025386800*x^12 + 134486022782700225*x^11 + 2366345074258640475*x^10 + 35859684567759302250*x^9 + 466277451513791667750*x^8 + 5165622516149912817000*x^7 + 48216742006981857309000*x^6 + 372948556274797637759625*x^5 + 2332188069734348007682500*x^4 + 11354348528498951245895625*x^3 + 40459665320954409153999375*x^2 + 94032401099596806911439375*x + 107145471557284795514880000)/107145471557284795514880000 )
}
# v23
v23 <- function(x) {
return( (x^22 + 275*x^21 + 37653*x^20 + 3409560*x^19 + 228820515*x^18 + 12091058595*x^17 + 521782139340*x^16 + 18829417262040*x^15 + 577216656722475*x^14 + 15188395563096525*x^13 + 345282279595077825*x^12 + 6804383826087747900*x^11 + 116315417092553078400*x^10 + 1721366411385367246500*x^9 + 21951610770646412856000*x^8 + 239344775104528631538000*x^7 + 2205184752540108215501625*x^6 + 16877181764451455880307875*x^5 + 104641871317872871553195625*x^4 + 505987954989411410235720000*x^3 + 1793338344579681991379413125*x^2 + 4150538718839947492706773125*x + 4714400748520531002654720000)/4714400748520531002654720000 )
}
# v24
v24 <- function(x) {
return( (x^23 + 299*x^22 + 44528*x^21 + 4388538*x^20 + 320878635*x^19 + 18498033015*x^18 + 872422838595*x^17 + 34482881442240*x^16 + 1160928591845715*x^15 + 33659328578215725*x^14 + 846499333177263150*x^13 + 18543981332320393950*x^12 + 354468851005624254900*x^11 + 5908721426717278068900*x^10 + 85642167991905000976500*x^9 + 1073505984389092320066000*x^8 + 11539630981616724845483625*x^7 + 105084571866055784500372875*x^6 + 796606323660382562645818500*x^5 + 4900946550340072015469936250*x^4 + 23550820409124372631515373125*x^3 + 83057425880345955113400950625*x^2 + 191488643096318168174459510625*x + 216862434431944426122117120000)/216862434431944426122117120000 )
}
# v25
v25 <- function(x) {
return( (x^24 + 324*x^23 + 52302*x^22 + 5590794*x^21 + 443757699*x^20 + 27803513430*x^19 + 1427363829045*x^18 + 61527989438685*x^17 + 2264380797997395*x^16 + 71969972109124320*x^15 + 1990916504836597800*x^14 + 48171457993524604200*x^13 + 1022052178969158437100*x^12 + 19024068913925375500200*x^11 + 310173582207161567594700*x^10 + 4413550536073387358149500*x^9 + 54479870357180417648123625*x^8 + 578209442112341503165201500*x^7 + 5210158342034725511661479250*x^6 + 39155018467736522209240131750*x^5 + 239192468624087541312192568125*x^4 + 1142844344290942723531592741250*x^3 + 4012130233592232103390903239375*x^2 + 9216828659958898330321714119375*x + 10409396852733332453861621760000)/10409396852733332453861621760000 )
}
# v26
v26 <- function(x) {
return( (x^25 + 350*x^24 + 61050*x^23 + 7055250*x^22 + 605890725*x^21 + 41116244400*x^20 + 2289272745375*x^19 + 107203631968125*x^18 + 4294804449474000*x^17 + 148958919241035750*x^16 + 4509865528655949000*x^15 + 119844452167642125000*x^14 + 2804396124729568792500*x^13 + 57862051714753396110000*x^12 + 1052112269850251212102500*x^11 + 16820493824359850061937500*x^10 + 235435442336189299332253125*x^9 + 2866363997113919044386393750*x^8 + 30073164352865410147765143750*x^7 + 268401985517264444722345218750*x^6 + 2001168299672231040727998496875*x^5 + 12145697900998969623892450875000*x^4 + 57725814415266540109375762078125*x^3 + 201799079872386039293085069609375*x^2 + 462034001190719350639625613609375*x + 520469842636666622693081088000000)/520469842636666622693081088000000 )
}
# v27
v27 <- function(x) {
return( (x^26 + 377*x^25 + 70850*x^24 + 8825700*x^23 + 817548225*x^22 + 59898856875*x^21 + 3604992566175*x^20 + 182749632565500*x^19 + 7939727936390250*x^18 + 299277074972625750*x^17 + 9872386621333236000*x^16 + 286709476727912238000*x^15 + 7358485307099969542500*x^14 + 167233500579206579017500*x^13 + 3366594338440387056502500*x^12 + 59957096888220149758140000*x^11 + 941896182959303001933628125*x^10 + 12990088017570058915673278125*x^9 + 156193180225877848100766468750*x^8 + 1621694381396207901371776687500*x^7 + 14347659633466395497955878559375*x^6 + 106200607985593828538108380228125*x^5 + 640719313663217082056213404078125*x^4 + 3030363986220446504652497411437500*x^3 + 10551987994810021315293879094078125*x^2 + 24084203903363353505313987382078125*x + 27064431817106664380040216576000000)/27064431817106664380040216576000000 )
}
# v28
v28 <- function(x) {
return( (x^27 + 405*x^26 + 81783*x^25 + 10951200*x^24 + 1091144925*x^23 + 86060400075*x^22 + 5581654843050*x^21 + 305319379815450*x^20 + 14335965076182750*x^19 + 585107280682674750*x^18 + 20945638395320388750*x^17 + 661860168338575206000*x^16 + 18540154899488546824500*x^15 + 461572912863205360717500*x^14 + 10223167862187856796220000*x^13 + 201354059102716406131245000*x^12 + 3520051349152769441533648125*x^11 + 54433520067779391000752915625*x^10 + 740747141016530499306063984375*x^9 + 8806580671786588914007034250000*x^8 + 90567295559088166862429382871875*x^7 + 794888270391980812439990551078125*x^6 + 5844549104957314680423524035256250*x^5 + 35066329669381300607463167615343750*x^4 + 165100551292052793052571247077390625*x^3 + 572787579633484461900595700274140625*x^2 + 1303527238695364400161681547826140625*x + 1461479318123759876522171695104000000)/1461479318123759876522171695104000000 )
}
# v29
v29 <- function(x) {
return( (x^28 + 434*x^27 + 93933*x^26 + 13486473*x^25 + 1441583325*x^24 + 122068182600*x^23 + 8507708445600*x^22 + 500677299322200*x^21 + 25327462749538950*x^20 + 1115537988501436500*x^19 + 43179715061262029250*x^18 + 1478740065756070367250*x^17 + 45014561633031555064500*x^16 + 1221719263742235780522000*x^15 + 29609230202442481946355000*x^14 + 640950277176794248368705000*x^13 + 12379629949672291311308428125*x^12 + 212835830779654015869767081250*x^11 + 3244689064134382485340698103125*x^10 + 43621696299563522381392041515625*x^9 + 513283167804844434734767026871875*x^8 + 5232685752787300988699030311800000*x^7 + 45588962624556355302423051589162500*x^6 + 333138334022204349309062893413412500*x^5 + 1988549694099880424640655963075265625*x^4 + 9323116798112282493686871795375843750*x^3 + 32234056538903525342793849362629734375*x^2 + 73155477446368801885414656825541734375*x + 81842841814930553085241614925824000000)/81842841814930553085241614925824000000 )
}
# v30
v30 <- function(x) {
return( (x^29 + 464*x^28 + 107387*x^27 + 16492329*x^26 + 1886636934*x^25 + 171082015650*x^24 + 12780094836600*x^23 + 806954803363800*x^22 + 43852522824460350*x^21 + 2077981572983916600*x^20 + 86685696612818052750*x^19 + 3205928668206551537250*x^18 + 105642904329030440121750*x^17 + 3112330852329561093231000*x^16 + 82143158543358620508801000*x^15 + 1943756406084263454008325000*x^14 + 41222392422628032487900153125*x^13 + 782298808464579416189954775000*x^12 + 13247973110778121231219750921875*x^11 + 199366771378013881677745550465625*x^10 + 2650746286483457031422977061137500*x^9 + 30896844143029522725437381655393750*x^8 + 312455936016708705726073597490962500*x^7 + 2703764390499134825035061576049862500*x^6 + 19644881397276710938020989313986128125*x^5 + 116704800279505825424282293801440187500*x^4 + 545005480435079062495571798108301140625*x^3 + 1878262643624966221081870221132806859375*x^2 + 4251705056257952260553877053981702859375*x + 4746884825265972078944013665697792000000)/4746884825265972078944013665697792000000 )
}
# v31
v31 <- function(x) {
return( (x^30 + 495*x^29 + 122235*x^28 + 20036100*x^27 + 2447376120*x^26 + 237114308340*x^25 + 18939047400000*x^24 + 1279818312318000*x^23 + 74516805352284750*x^22 + 3788229963137870250*x^21 + 169804959532174716750*x^20 + 6760042229332091700000*x^19 + 240291908393705604686250*x^18 + 7654975738477870018466250*x^17 + 219085716045859308610965000*x^16 + 5640198540535401376904370000*x^15 + 130635187102504151372283103125*x^14 + 2719751252328096943121261971875*x^13 + 50798315917077933208337580121875*x^12 + 848517453806141822007513345637500*x^11 + 12619084855384151115310254584418750*x^10 + 166084904753685831328009211773406250*x^9 + 1919091831454243887448817443571437500*x^8 + 19263928999384696228516962243070875000*x^7 + 165648158484229991489914314420678703125*x^6 + 1197173277129724927015436706070677234375*x^5 + 7080474296087405286255380250988951640625*x^4 + 32943575028424472783329462713305971875000*x^3 + 113190938386505993083302349879684500703125*x^2 + 255597483144485155451622759850618260703125*x + 284813089515958324736640819941867520000000)/284813089515958324736640819941867520000000 )
}
# v32
v32 <- function(x) {
return( (x^31 + 527*x^30 + 138570*x^29 + 24192090*x^28 + 3148639620*x^27 + 325219848660*x^26 + 27712276808580*x^25 + 1999502113518000*x^24 + 124429719532686750*x^23 + 6768902177229260250*x^22 + 325122388020827397000*x^21 + 13891850529683429803500*x^20 + 530973724254985547786250*x^19 + 18227819707800916624661250*x^18 + 563559624277363459441946250*x^17 + 15718141478644929573008760000*x^16 + 395724518507668016086788493125*x^15 + 8990240233248296208990850921875*x^14 + 184066127281154683421279416743750*x^13 + 3388433249660038482424392351731250*x^12 + 55893474999497384037693435211931250*x^11 + 822277317233661689324142450163181250*x^10 + 10721591783399592947833305667561968750*x^9 + 122894887897913866150753104195928500000*x^8 + 1225164253450388284058347237789576828125*x^7 + 10473470152246604450450638313628684609375*x^6 + 75319351092481726126135272497017554000000*x^5 + 443611084201493979386141517270665167031250*x^4 + 2056861865063549887299740649964736841328125*x^3 + 7047053786334844740449763752631688302890625*x^2 + 15876259561329552807285629170829581422890625*x + 17658411549989416133671730836395786240000000)/17658411549989416133671730836395786240000000 )
}
# v33
v33 <- function(x) {
return( (x^32 + 560*x^31 + 156488*x^30 + 29042040*x^29 + 4019554860*x^28 + 441719514600*x^27 + 40070631057660*x^26 + 3080280909052620*x^25 + 204409804073406750*x^24 + 11870520678069417000*x^23 + 609416279464456327500*x^22 + 27872113214579007874500*x^21 + 1142215147561056459140250*x^20 + 42121637299275266275042500*x^19 + 1402039330836205624176363750*x^18 + 42205443819681012166780233750*x^17 + 1150195309482624635591208973125*x^16 + 28380741640124028997243487085000*x^15 + 633578138943569493870821962837500*x^14 + 12775805740998927336909642605662500*x^13 + 232092003981819385123761837501956250*x^12 + 3784631492207023043321894516395537500*x^11 + 55124566914017324171336997976373756250*x^10 + 712582435984891478281584915911836781250*x^9 + 8107277975733564788500521072761572828125*x^8 + 80307832598918736641776430867634563812500*x^7 + 682780619922784784252272294687481261343750*x^6 + 4887452798657915820828122594594700853031250*x^5 + 28673526917153188650468231686204646863203125*x^4 + 132515627555211387865733943400480635623906250*x^3 + 452793594314089926715170981833994256202109375*x^2 + 1017862763913751242992666368598659415882109375*x + 1130138339199322632554990773529330319360000000)/1130138339199322632554990773529330319360000000 )
}
# v34
v34 <- function(x) {
return( (x^33 + 594*x^32 + 176088*x^31 + 34675608*x^30 + 5094110340*x^29 + 594462599280*x^28 + 57297692127060*x^27 + 4683106151359020*x^26 + 330701321344564170*x^25 + 20455732449152500500*x^24 + 1119848668621441258500*x^23 + 54686429511015086284500*x^22 + 2396460242217111813492750*x^21 + 94663534087083863395494000*x^20 + 3381756283902143139103361250*x^19 + 109503331699818882127245693750*x^18 + 3218262056646994231763440426875*x^17 + 85890507114255260776803935741250*x^16 + 2080995962589894972730239804172500*x^15 + 45721868966064541018192384673212500*x^14 + 909209708254762533979972895602068750*x^13 + 16317599707225269840005033741501175000*x^12 + 263279298985403591554041196378128318750*x^11 + 3799558183169861631876456802588767131250*x^10 + 48724476826872379050550861279736269359375*x^9 + 550529955191465494374806653087805787843750*x^8 + 5420942743258990246117081715877920526281250*x^7 + 45854289994025002875964460275843576533656250*x^6 + 326808147635286053720983709956481398898109375*x^5 + 1910274296418709084194764307945168741142500000*x^4 + 8801278130292407362256409416064274300508203125*x^3 + 29996652800015506552763609205974291812817109375*x^2 + 67291217993593153427078304732442192351697109375*x + 74589130387155293748629391052935801077760000000)/74589130387155293748629391052935801077760000000 )
}
# v35
v35 <- function(x) {
return( (x^34 + 629*x^33 + 197472*x^32 + 41190864*x^31 + 6411783444*x^30 + 793132902540*x^29 + 81076196098260*x^28 + 7032311528568480*x^27 + 527391779701643010*x^26 + 34675889266968759810*x^25 + 2019900896384151280500*x^24 + 105079619598979942917000*x^23 + 4912035999723805782579750*x^22 + 207297165471288118629653250*x^21 + 7925605920082168582087073250*x^20 + 275209389611023895943310395000*x^19 + 8693428641637938338125725114375*x^18 + 250021872003251966596739397511875*x^17 + 6547302332531168533124044462717500*x^16 + 156014654983328974572895094294355000*x^15 + 3378190632422247748962361667955543750*x^14 + 66324133661237209208903542999614956250*x^13 + 1177064882590018702594323085902194118750*x^12 + 18806478225337866350456804996142081300000*x^11 + 269098851450724353699385355829884762971875*x^10 + 3425274087976935858357307468592245680046875*x^9 + 38452740054746919908605480901146267796906250*x^8 + 376531797332823407889106444944396728636812500*x^7 + 3169774127264836232030780247058783143984796875*x^6 + 22499187597441730468616738035203496871723390625*x^5 + 131058833101089788750721325124555073733628203125*x^4 + 602079731269021985099430221250152121345850312500*x^3 + 2047070302794616585909476512326745451997626328125*x^2 + 4583100735957896573362875808126562688641466328125*x + 5072060866326559974906798591599634473287680000000)/5072060866326559974906798591599634473287680000000 )
}
# v36
v36 <- function(x) {
return( (x^35 + 665*x^34 + 220745*x^33 + 48694800*x^32 + 8018227140*x^31 + 1049604240300*x^30 + 113594645102400*x^29 + 10437511764695400*x^28 + 829781175430087650*x^27 + 57881127573841052250*x^26 + 3580315913397745471950*x^25 + 197995060832650901820000*x^24 + 9850778120875863099678750*x^23 + 443074893458030796193481250*x^22 + 18083167028175286394940082500*x^21 + 671489685615132325047664057500*x^20 + 22729107511800157031234555259375*x^19 + 702080161368424760179277103459375*x^18 + 19798461548703522762751232530846875*x^17 + 509568980940012075361593495281100000*x^16 + 11958996656505341350471591854145068750*x^15 + 255502809076883083150795796405125406250*x^14 + 4957540501280539627501825036880246625000*x^13 + 87076241415558951100927543978469340187500*x^12 + 1378681066745658468376336850602267559671875*x^11 + 19571205175020397080320428818385331458359375*x^10 + 247394459421340007268401236485273254279765625*x^9 + 2760601680727132442222646260815465332045000000*x^8 + 26891277359232710929044486278555777048103984375*x^7 + 225364731742391249318586673846965618086750390625*x^6 + 1593506026934802269210809297412782370395648593750*x^5 + 9251962715940948042647037679470786987765311718750*x^4 + 42386412297819089587571301336086937582169597265625*x^3 + 143783881325991824415207278646345253424480056640625*x^2 + 321306011647421423536945229352332459989548856640625*x + 355044260642859198243475901411974413130137600000000)/355044260642859198243475901411974413130137600000000 )
}
# v37
v37 <- function(x) {
return( (x^36 + 702*x^35 + 246015*x^34 + 57303855*x^33 + 9966019140*x^32 + 1378351553040*x^31 + 157678023195000*x^30 + 15322081504098600*x^29 + 1289031693076685250*x^28 + 95221280468194996500*x^27 + 6242847781794433875450*x^26 + 366269908762344939001650*x^25 + 19354541040843106387038750*x^24 + 925763021380948088077740000*x^23 + 40236911701076826204614145000*x^22 + 1593731204052071931189608265000*x^21 + 57646571163787037933713086249375*x^20 + 1906722859493833082834708532206250*x^19 + 57710790262598459812432196117653125*x^18 + 1598484366118705827312911284477678125*x^17 + 40494859589146017570720827589886668750*x^16 + 937165618497687540127676532091394325000*x^15 + 19776703427739758450247981228377520187500*x^14 + 379571130991110789123535221154403891062500*x^13 + 6603255551679195534431989489310427970921875*x^12 + 103670750246505563651276976705123652598343750*x^11 + 1460809180272604626248267823225163804698046875*x^10 + 18346452624271552900131924159387680351670234375*x^9 + 203569784925769187231293846970745558298983984375*x^8 + 1973297760092517459706478281953091126213509375000*x^7 + 16467578321932624724237529771312513164121174375000*x^6 + 116016866520572700079771260606127205804273767500000*x^5 + 671519876981803556487569863540100452750210794140625*x^4 + 3068446329875509005957627070836343946594182267968750*x^3 + 10386177704466849132601454734596500199703152821484375*x^2 + 23167771087609780269366587185427579072388106421484375*x + 25563186766285862273530264901662157745369907200000000)/25563186766285862273530264901662157745369907200000000 )
}
# v38
v38 <- function(x) {
return( (x^37 + 740*x^36 + 273393*x^35 + 67144455*x^34 + 12315477195*x^33 + 1796924356920*x^32 + 216947139975720*x^31 + 22259914524678600*x^30 + 1978525360761122250*x^29 + 154516738349722518000*x^28 + 10718247963799598710950*x^27 + 665926602288477765023250*x^26 + 37301766570198008398119600*x^25 + 1893490073423103407429677500*x^24 + 87450825791505178696578885000*x^23 + 3686050612508066893829543805000*x^22 + 142114324978546850286762324294375*x^21 + 5019637702338333131255215189672500*x^20 + 162580547534759279368341165388996875*x^19 + 4830288620824219576809114267066253125*x^18 + 131608468457912249727556770805114321875*x^17 + 3285867474668156559229484532304821112500*x^16 + 75069474919103323317780896621769785362500*x^15 + 1566173336655496296138414094857055102312500*x^14 + 29757094542136953670967637979729065325734375*x^13 + 513072594450615686786060325042370186795500000*x^12 + 7992066445802455136278717355647953918393703125*x^11 + 111838240161718248980021064845798163852345234375*x^10 + 1396089205503420125739868917330944781157549218750*x^9 + 15408903565193283816971872182022297973946452343750*x^8 + 148678528248131294524571574662169618620426302500000*x^7 + 1235812192411991181327923285055378100964513625000000*x^6 + 8676683666901319861991786845362877653245100751640625*x^5 + 50074837718601757960087517518643375639108937857812500*x^4 + 228245867125627988555592976763976920407890093847265625*x^3 + 770972565809222917816671328076375593451015109568359375*x^2 + 1716810476161799821937291129437875430029701675968359375*x + 1891675820705153808241239602722999673157373132800000000)/1891675820705153808241239602722999673157373132800000000 )
}
# v39
v39 <- function(x) {
return( (x^38 + 779*x^37 + 302993*x^36 + 78353568*x^35 + 15135544305*x^34 + 2326489876305*x^33 + 296011811680200*x^32 + 32022535823586000*x^31 + 3002481428896337850*x^30 + 247507430305495263750*x^29 + 18135051404586279574950*x^28 + 1191120752514658101859800*x^27 + 70598096684621896649282100*x^26 + 3795880168503201835733777100*x^25 + 185912309609506555882922115000*x^24 + 8320944379457841364748224710000*x^23 + 341161058053982462553557689764375*x^22 + 12835925576158409897027143025863125*x^21 + 443680258865705934718633216010656875*x^20 + 14097379830305498500804560694239075000*x^19 + 411765208465716985182485398294957003125*x^18 + 11050767113684979293155334009806566103125*x^17 + 272221523399192716871549968560059052112500*x^16 + 6145411306720799018523048788785012009425000*x^15 + 126859841414777724031549311860866481669109375*x^14 + 2387769550605243768057021517765301302316765625*x^13 + 40828712490641859090586578158359645873305703125*x^12 + 631322559138877832838137692962915168547935937500*x^11 + 8777413056176824558421259197153623595412334687500*x^10 + 108946880333922432241543089643195598311502250000000*x^9 + 1196483970681274594078658883039685880848785061875000*x^8 + 11494630641533050503523361936745081785773928497500000*x^7 + 95183537135740702554946416799239344720761054501640625*x^6 + 666119378068595468161504383539407689019511091224296875*x^5 + 3833634182864954561681894238106299966423733619609765625*x^4 + 17432920865980066082374958631846690783226991960418750000*x^3 + 58768780346044295740370969407089669345404819784026953125*x^2 + 130652461532840140453538074310563656925384998830426953125*x + 143767362373591689426334209806947975159960358092800000000)/143767362373591689426334209806947975159960358092800000000 )
}
# v40
v40 <- function(x) {
return( (x^39 + 819*x^38 + 334932*x^37 + 91079274*x^36 + 18504747729*x^35 + 2992453825725*x^34 + 400703856113925*x^33 + 45639079160875200*x^32 + 4507540612604879850*x^31 + 391626538892519480550*x^30 + 30262915489555547498700*x^29 + 2097873322743972080607300*x^28 + 131345255062869459844131900*x^27 + 7466981196103540461496446300*x^26 + 387093958540176253176812301300*x^25 + 18360209098371195382426018920000*x^24 + 798812998924163737614710048814375*x^23 + 31940944827181427800026373652668125*x^22 + 1175328016706735298849180368484855000*x^21 + 39830834844516442714485287222857173750*x^20 + 1243510618453741396729954479255062428125*x^19 + 35756679621627998404104457907503986290625*x^18 + 946318317333976453754025343158259584403125*x^17 + 23023145757470747464559146839508673240400000*x^16 + 514020753738188062198501385554322238262884375*x^15 + 10506799401151018106076177476860756129139765625*x^14 + 196033733279982704014292976813104230523895468750*x^13 + 3326017583521240532816851851414651796186112343750*x^12 + 51076024518481639358576484625668939888124042500000*x^11 + 705810968153946502214188715049642002799541008750000*x^10 + 8713818713721922418745132068420182164342440311875000*x^9 + 95248508589222272089029483749523093445188882828750000*x^8 + 911302312684587288305105114308140151510709977824140625*x^7 + 7519334051841926052117646393084640508914307015342421875*x^6 + 52460348781872423737471714236483061264848043278983437500*x^5 + 301121850397986703646835132251712888298583279811541406250*x^4 + 1366237845294549251918492866795591478087429216815433203125*x^3 + 4597079767832206616721731749249378527176151302416475390625*x^2 + 10204006900402282504348765931720349558414605268035675390625*x + 11213854265140151775254068364941942062476907931238400000000)/11213854265140151775254068364941942062476907931238400000000 )
}
# v41
v41 <- function(x) {
return( (x^40 + 860*x^39 + 369330*x^38 + 105481350*x^37 + 22512235785*x^36 + 3825167473530*x^35 + 538356732097275*x^34 + 64472160398229675*x^33 + 6698216412326889450*x^32 + 612496028910158593200*x^31 + 49844242434181521526200*x^30 + 3641282012711305003041000*x^29 + 240434667845556008035711500*x^28 + 14428279714435621833235437000*x^27 + 790310943129767438097620401500*x^26 + 39650376818080889307150695491500*x^25 + 1826984708432950679030567108334375*x^24 + 77473285765858760844064846435087500*x^23 + 3027902816683258111250710040339606250*x^22 + 109175187830213825346586928963463618750*x^21 + 3633360709124727959599071712626492853125*x^20 + 111610827347306223604631681142062794406250*x^19 + 3163232453874912354808501733423506734421875*x^18 + 82641199749511264051062743458479027057796875*x^17 + 1987502082216315899930286783282877325648484375*x^16 + 43918148394133242148978767537891701616227250000*x^15 + 889482493755949899015320690285914135047120000000*x^14 + 16460277713280081701774481297892635241287108750000*x^13 + 277245220197925995590122410521865262028779681875000*x^12 + 4230056659929179617955966154220798855080099941250000*x^11 + 58120586484498177573738342121895122360310310924375000*x^10 + 713929637263478763819933860607356027113502144971875000*x^9 + 7769194931108590878715227944273802879564309541494140625*x^8 + 74044402877816798098390319737578871569196135396504687500*x^7 + 608891068618174951594177547324746458924506762414322656250*x^6 + 4235648009038418483957213699987942483162186525735299218750*x^5 + 24251498475541538729077962917925770988779758482492580078125*x^4 + 109797393855512499014445682492509922339908200997204832031250*x^3 + 368776228791314398608643842373171874678154406856520755859375*x^2 + 817330399396920469618806576970849557177230724106056755859375*x + 897108341211212142020325469195355364998152634499072000000000)/897108341211212142020325469195355364998152634499072000000000 )
}
# vk: k = 1 - 41
vv <- function(x) {
return( do.call(paste('v',k,sep=''),list(x=x)) )
}
# convolutionDEfunc
# x - single number
vvDE <- function(x) {
if ( x == 0 ) {
return( 0.5 )
} else if ( x > 0 ) {
return( 1 - 0.5 * vv(x) * exp(-x) )
} else if ( x < 0 ) {
return( 0.5 * vv(-x) * exp(x) )
} else {
stop('x not a number.')
}
}
# x - numeric vector
vvDEvec <- function(x) {
return( sapply(x,vvDE) )
}
# return
return( vvDEvec )
}
## [model based methods tools - done]
## [exact methods tools]
### compute quantile given cd function, used in exact method
.quantileCD <- function(CDF, prbblty) {
ff <- function(theta) {
CDF(theta) - prbblty
}
for ( power.l in 0:15 ) {
if ( ff(10^(-power.l)) < 0 ) {
break
}
}
for ( power.u in 0:15 ) {
if ( ff(10^( power.u)) > 0 ) {
break
}
}
if ( power.l == 15 ) {
qntl = 0
} else if ( power.u == 15 ) {
qntl = Inf
} else {
qntl = uniroot(f = function(theta) { CDF(theta) - prbblty }, interval=c(10^(-power.l), 10^(power.u)), tol=0.001)$root
}
return(qntl)
}
### compute cd function (based on p-value functions) using exact method
midp.oddsratio <- function(x, N, M, t, or){ # x and t can be vectors
mid.p = rep(NA, length(x))
for(i in min(t):max(t)) {
index = which( t == i )
if ( !identical(index, integer(0)) ) {
for (idx in index) {
mid.p[idx] = 1 - pFNCHypergeo.wrap(x[idx], N, M, i, or) + dFNCHypergeo(x[idx], N, M, i, or) / 2
}
} else {
#cat('\nindex null, mid.p all NA\n')
}
}
return(mid.p)
}
### conditional MLE, used in excat method
.conditionalMLE <- function(data_matrix){
# number of study
K = dim(data_matrix)[1]
# or.hat.CMLE
if ( sum(data_matrix[,1]) == 0 ) {
or.hat.CMLE = 0
} else if ( sum(data_matrix[,2]) == 0 ) {
or.hat.CMLE = Inf
} else {
.estimatingFunctionCMLE <- function(theta) {
summation = 0
for(j in 1:K) {
expectation = 0
for(x in max(0, data_matrix[j,1]+data_matrix[j,2]-data_matrix[j,4]):min(data_matrix[j,3], data_matrix[j,1]+data_matrix[j,2])) {
expectation = expectation + x * dFNCHypergeo(x, data_matrix[j,3], data_matrix[j,4], data_matrix[j,1]+data_matrix[j,2], theta)
}
summation = summation + expectation - data_matrix[j,1]
}
return(summation)
}
or.hat.CMLE = uniroot(.estimatingFunctionCMLE, c(0.0001,10000), tol=0.001)$root
}
return(or.hat.CMLE)
}
### calculate estimates of event rates and weights using empirical Bayes approach
### please refer to Efron 1996 JASA for details, functions are used in exact method
.Integrate <- function(f) { # input a function with range between 0 and 1
ss1 = seq(0.0001,0.1,0.0001)
ss2 = seq(0.101,0.999,0.001)
return( sum(f(ss1))*0.0001 + sum(f(ss2))*0.001 )
}
#####
.Estimates <- function(data_matrix) {
K = dim(data_matrix)[1]
x = data_matrix[,1]
y = data_matrix[,2]
n = data_matrix[,3]
m = data_matrix[,4]
if( sum(x)+sum(y) == 0 ) { # in the case of all zeros in two arms
p0.hat = NA
p1.hat = NA
psi.hat = NA
weight.hat = rep(1, K)
} else if ( sum(y) == 0 ) { # in the case of all zeros in the control arm
p0.hat = NA
p1.hat = rep(sum(x)/sum(n), K)
psi.hat = Inf
weight.hat = sqrt( m * p1.hat / (1 - p1.hat) )
} else if ( sum(x) == 0 ) { # in the case of all zeros in the treatment arm
p0.hat = .Estimates.1arm(y, m)
p1.hat = NA
psi.hat = 0
weight.hat = sqrt( n * p0.hat / (1 - p0.hat) )
} else {
# use moment estimates as initial values
mu = mean(y/m)
v.square = var(y/m)
beta1.initial = mu*(mu*(1-mu)/v.square-1)
beta2.initial = beta1.initial * (1-mu)/mu
psi.initial = .conditionalMLE(data_matrix)
# use optim() to obtain parameters estimates
parameter.initials = c(log(beta1.initial),log(beta2.initial),log(psi.initial))
.marginalLikelihood <- function(parameters) {
beta1 = exp(parameters[1])
beta2 = exp(parameters[2])
psi = exp(parameters[3])
ff = 0
for(i in 1:K) {
.individualLikelihood <- function(p0) {
p1 = psi * p0 / ( 1 - p0 + psi * p0 )
return(dbinom(x[i],n[i],p1)*dbinom(y[i],m[i],p0)*dbeta(p0,beta1,beta2))
}
ff = ff + log(.Integrate(.individualLikelihood))
}
return( -ff )
}
parameter.estimates = optim(parameter.initials, .marginalLikelihood, method='Nelder-Mead', lower=-Inf, upper=Inf)$par
# resolve parameters estimates
beta1.hat = exp(parameter.estimates[1])
beta2.hat = exp(parameter.estimates[2])
psi.hat = exp(parameter.estimates[3])
# update
p0.hat = rep(NA, K)
for(i in 1:K) {
.individualLikelihoodEmpirical <- function(p0) {
p1 = psi.hat * p0 / ( 1 - p0 + psi.hat * p0 )
return(dbinom(x[i],n[i],p1)*dbinom(y[i],m[i],p0)*dbeta(p0,beta1.hat,beta2.hat))
}
denominator <- .Integrate(.individualLikelihoodEmpirical)
nominator <- .Integrate( f = function(p0) { p0*.individualLikelihoodEmpirical(p0) } )
p0.hat[i] <- nominator / denominator
}
p1.hat = psi.hat * p0.hat / ( 1 - p0.hat + psi.hat * p0.hat )
weight.hat = 1 / sqrt( (n*p1.hat*(1-p1.hat))^(-1) + (m*p0.hat*(1-p0.hat))^(-1) )
}
# return
return( list(or.hat=psi.hat, pt.hat=p1.hat, pc.hat=p0.hat, weight.hat=weight.hat) )
}
#####
.Estimates.1arm <- function(y, m) {
K = length(m)
# initial
mu = mean(y/m)
v.square = var(y/m)
beta1.initial = mu*(mu*(1-mu)/v.square-1)
beta2.initial = beta1.initial*(1-mu)/mu
# optim()
parameter.initials = c(log(beta1.initial),log(beta2.initial))
.marginalLikelihood <- function(parameters) {
beta1 = exp(parameters[1])
beta2 = exp(parameters[2])
ff = 0
for(i in 1:K) {
ff = ff + lbeta(y[i]+beta1, m[i]-y[i]+beta2) - lbeta(beta1, beta2)
}
return(-ff)
}
parameter.estimates = optim(parameter.initials, .marginalLikelihood)$par
# update
beta1.hat = exp(parameter.estimates[1])
beta2.hat = exp(parameter.estimates[2])
p.hat = (y + beta1.hat) / (m + beta1.hat + beta2.hat)
# return
return(p.hat)
}
### mixed beta adjustment, used in exact method
adjust.beta<-function(x, n, pt, m, pc) {
ifelse(x<1/2, pbeta(x,1+1/(2*m*pc*(1-pc)),1+1/(2*m*pc*(1-pc))), pbeta(x,1+1/(2*n*pt*(1-pt)),1+1/(2*n*pt*(1-pt))))
}
### update on 2012.11.26
### pFNCHypergeo gives NaN when odds extremely large which results a probability very close to 0 or 1 and pFNCHypergeo can not handle it.
### wrap pFNCHypergeo in pFNCHypergeo.wrap to recursively adjust odds when it results NaN
pFNCHypergeo.wrap <- function(x, m1, m2, n, odds, precision=1E-7, lower.tail=TRUE) {
prbblty = pFNCHypergeo(x=x, m1=m1, m2=m2, n=n, odds=odds, precision=precision, lower.tail=lower.tail)
while ( is.na(prbblty) ) {
odds = odds / 10
prbblty = pFNCHypergeo(x=x, m1=m1, m2=m2, n=n, odds=odds, precision=precision, lower.tail=lower.tail)
}
return(prbblty)
}
## [exact methods tools - done]
# other functions used in this package [done]
# main function
gmeta <- function(gmi, gmi.type = c('pivot', 'cd', 'pvalue', '2x2'),
method = c('fixed-mle',
'fixed-robust1', 'fixed-robust2', 'fixed-robust2(sqrt12)',
'random-mm', 'random-reml', 'random-tau2',
'random-robust1', 'random-robust2', 'random-robust2(sqrt12)',
'fisher', 'normal', 'stouffer', 'min', 'tippett', 'max', 'sum',
'MH', 'Mantel-Haenszel', 'Peto', 'exact1', 'exact2'),
linkfunc = c('inverse-normal-cdf', 'inverse-laplace-cdf'),
weight = NULL, study.names = NULL, gmo.xgrid = NULL, ci.level = 0.95,
tau2 = NULL, mc.iteration = 10000, eta = 'Inf', verbose = FALSE,
report.error = FALSE)
{
UseMethod('gmeta')
}
# main function
gmeta.default <- function(gmi, # generalized meta-analysis input: vector of p-values, data.frame/matrix of mean/sd, list of cd-functions or matrix of 2x2 tables.
# parameters specifying input properties
gmi.type = c('pivot', # type of input (gmi), if model-based meta-analysis, input data.frame/matrix of mean/sd
'cd', # type of input (gmi), if model-based meta-analysis, input list of cd-funcitons (pointers to function)
'pvalue', # type of input (gmi), if combine p-values, input vector of p-values
'2x2'), # type of input (gmi), if combine 2x2 tables, input matrix of 2x2 tables in format c(rx-event_treatment_grounp, nx-number_observation_treatment_group, ry-event_control_group, ny-number_observation_treatment_group) for each row of the matrix
method = c('fixed-mle', # method of combination, if model-based meta-analysis, fixed-effects model with mle
'fixed-robust1', # method of combination, if model-based meta-analysis, fixed-effects model with robust method 1 - adaptive-weighting (xie2012section4)
'fixed-robust2', # method of combination, if model-based meta-analysis, fixed-effects model with robust method 2 - m-estimating-equation (xie2012section4)
'fixed-robust2(sqrt12)', # method of combination, if model-based meta-analysis, fixed-effects model with robust method 2 - m-estimator-simplified-variance (xie2012section4)
'random-mm', # method of combination, if model-based meta-analysis, random-effects model with moment estimator
'random-reml', # method of combination, if model-based meta-analysis, random-effects model with restricted-mle estimator
'random-tau2', # method of combination, if model-based meta-analysis, random-effects model with heterogeneity specified by user
'random-robust1', # method of combination, if model-based meta-analysis, random-effects model with robust method 1 - adaptive-weighting (xie2012section4)
'random-robust2', # method of combination, if model-based meta-analysis, random-effects model with robust method 2 - m-estimating-equation (xie2012section4)
'random-robust2(sqrt12)', # method of combination, if model-based meta-analysis, random-effects model with robust method 2 - m-estimator-simplified-variance (xie2012section4)
'fisher', # method of combination, if combine p-value vector, see fisher1932, bahadur optimal
'normal', # method of combination, if combine p-value vector, see stouffer(1949)
'stouffer', # method of combination, if combine p-value vector, see stouffer1949, alias 'normal'
'min', # method of combination, if combine p-value vector, see tippett(1931)
'tippett', # method of combination, if combine p-value vector, see tippett1931, alias 'min'
'max', # method of combination, if combine p-value vector, see marden(1991)
'sum', # method of combination, if combine p-value vector, see marden(1991)
'MH', # method of combination, if cobmine 2x2 tables, see Mantel-Haenszel method
'Mantel-Haenszel', # method of combination, if cobmine 2x2 tables, see Mantel-Haenszel method, alias 'MH'
'Peto', # method of combination, if cobmine 2x2 tables, see Peto's method
'exact1', # method of combination, if cobmine 2x2 tables, see liu2012exact
'exact2'), # method of combination, if cobmine 2x2 tables, see tian2009exact
linkfunc = c('inverse-normal-cdf', # 'link' function, F0^{-1}(\cdot), see xie2012section2
'inverse-laplace-cdf'), # 'link' function, F0^{-1}(\cdot), see xie2012section2
weight = NULL, # study-specific weight, in default inverse-of-adjusted-variance
study.names = NULL, # names of studies, in default, study1, study2, etc.
# parameters specifying output properties
#gmo.limit= NULL, # range of generalized meta-analysis output, in default, c(-1,1) with gmo.n.pts=2001, seq(from=-1,to=1,by=0.001)
#gmo.n.pts= 2001, # number of points in the range of generalized meta-analysis output, suppressed if gmo.limit is given as a vector
gmo.xgrid= NULL, # gridding positions where generalized meta-analysis output (fixed- or random-effects model-based) evaluating the combined and individual CDs. gmo.xgrids suppress gmo.limit and gmo.n.pts, in default, gmo.xgrid = seq(-1,1,by=0.001).
ci.level = 0.95, # confidence level
# parameters specifying model-based meta-analysis (random-effects model) heterogeneity
tau2 = NULL, # numeric value or method used to estimate heterogeneity, in default, REML. (only when method=='random-tau2')
# parameters specifying properties when use liu2012exact (exact1) method to combine 2x2 tables
mc.iteration = 10000, # monte-carlo iteration (see liu2012exact)
# parameters specifying properties when use tian2009exact (exact2) method to combine 2x2 tables
eta = 'Inf', # seq(0.05, 0.95, length=20), # number of confidence-level taken to make confidence interval for later combination (see tian2009exact)
# parameters specifying whether report detailed inforamtion
verbose = FALSE, # whether detailed inforamtion reported during meta-analysis
# parameters specifying whether report error of coverage probability when use liu2012exact/tian2009exact method to combine 2x2 tables
report.error = FALSE) # minor error message and exact error of coverage probability of the computed confidence interval for 2x2-"exact1"&"exact2" method
{
#check
mf <- match.call()
gmi.type <- match.arg(gmi.type)
gmi <- gmeta.check.gmi(gmi, gmi.type, verbose)
method <- match.arg(method)
linkfunc <- match.arg(linkfunc)
weight <- gmeta.check.weight(gmi, gmi.type, linkfunc, weight, verbose)
gmo.xgrid <- gmeta.check.gmo.xgrid(gmi, gmi.type, gmo.xgrid)
study.names <- gmeta.check.study.names(gmi, study.names)
#gmo.n.pts <- round(ifelse(is.numeric(gmo.n.pts)&&(gmo.n.pts>1), gmo.n.pts, 2001))
#if(! ((!is.null(gmo.limit))&&is.numeric(gmo.limit)&&(length(gmo.limit)==2)) ) { gmo.limit = c(-1,1) }
ci.level <- ifelse(is.numeric(ci.level)&&(ci.level>=0)&&(ci.level<=1), ci.level, 0.95)
mc.iteration <- round(ifelse(is.numeric(mc.iteration)&&(mc.iteration>1), mc.iteration, 10000))
#if(!(is.numeric(eta)&&(min(eta)>=0)&&(max(eta)<=1))) { eta = seq(0.05, 0.95, length=20) } # allow eta non-numeric as K\to\Infty
# meta-analysis
if (gmi.type == 'pvalue') {
gmo <- gmeta.p(gmi, method)
} else if (gmi.type == 'cd' || gmi.type == 'pivot') {
gmo <- gmeta.m(gmi, gmi.type, method, linkfunc, weight, gmo.xgrid, ci.level, tau2, verbose)
} else if (gmi.type == '2x2') { # combine 2x2 tables
gmo <- gmeta.e(gmi, method, weight, gmo.xgrid, ci.level, mc.iteration, eta, verbose, report.error)
} else {
stop("gmi.type must be 'cd', 'pivot', 'pvalue' or '2x2'.")
}
# post processing
gmo$call <- match.call()
gmo$input <- gmi
gmo$alpha <- 1 - ci.level
gmo$study.names <- study.names
# registr S3 class
if ( gmi.type=='pvalue' ) {
class(gmo) <- c('gmeta.p', 'gmeta')
} else if ( gmi.type == 'cd' || gmi.type == 'pivot' ) {
class(gmo) <- c('gmeta.m', 'gmeta')
} else if (gmi.type == '2x2') {
class(gmo) <- c('gmeta.e', 'gmeta')
} else {
stop("gmi.type must be 'cd', 'pivot', 'pvalue' or '2x2'.")
}
# return
return(gmo)
}
# preprocessing
## preprocessing - checks
# *****************************************************************************
# # checking - input, output, parameters formats, etc.
# *****************************************************************************
## main [gmeta.check.x()]
### check gmi - gmeta input formats
### update v2.0: 2x2 table order is changed - before x y M N, now x-M y-N.
gmeta.check.gmi <- function(gmi, gmi.type, verbose) {
# check by gmi.type
if ( gmi.type == 'cd' ) {
# verbose
if ( verbose ) {
cat('\ninput CDs must composite of a list.\n') #commentOut
}
# check
if ( !is.list(gmi) ) {
stop('input CDs must composite of a list.')
}
} else if ( gmi.type == 'pivot' ) {
# verbose
if ( verbose ) {
cat('\ninput pivots must be in form of of either a matrix or a data.frame of 2 columns.\n') #commentOut
}
# check
if ( ! (is.matrix(gmi) || is.data.frame(gmi)) ) {
stop('input pivots must be in form of either a matrix or a data.frame.')
}
if ( dim(gmi)[2] != 2 ) {
stop('input pivots must be in form of of either a matrix or a data.frame of 2 columns.')
}
} else if ( gmi.type == 'pvalue' ) {
# verbose
if ( verbose ) {
cat('\ninput p-values must be in form of a vector with all element in-between 0 and 1.\n') #commentOut
}
# check
if ( !is.vector(gmi) ) {
stop('input p-values must be in form of a vector.')
}
if ( any( gmi < 0 | gmi > 1 ) ) {
stop('input p-values must be between 0 and 1.')
}
} else if ( gmi.type == '2x2' ) {
# verbose
if ( verbose ) {
cat('\ninput 2x2 tables should be in format: r_group_1 n_group_1 r_group_2 n_group_2.\n') #commentOut
}
# check
if ( ! (is.matrix(gmi) || is.data.frame(gmi)) ) {
stop('input 2x2 tables must be in form of either a matrix or a data.frame.')
}
if ( dim(gmi)[2] != 4 ) {
stop('input 2x2 tables must composites of either a matrix or a data.frame of 4 columns.')
}
if ( any( gmi[,1] > gmi[,2] | gmi[,3] > gmi[,4] ) ) {
stop('input 2x2 tables should be in format: r_group_1 n_group_1 r_group_2 n_group_2 - second and fourth column should be larger than first and third column, respectively.')
}
if ( any( gmi < 0 ) ) {
stop('input data should be nonnegative.')
}
} else {
stop('gmi.type not recognize.')
}
# return
return(gmi)
}
### check weight - weight assigned to each study
### update v2.0: weights: need to distinguish whether use weight=1/s^2 (gaussian) or weight=1(laplacian, DE, warning if not), keep NULL set later.
gmeta.check.weight <- function(gmi, gmi.type, linkfunc, weight, verbose) {
if ( gmi.type == 'cd' || gmi.type == 'pivot' ) {
# number of studies
if ( gmi.type == 'cd' ) {
nn1 = length(gmi)
} else { # gmi.type == 'pivot'
nn1 = dim(gmi)[1]
}
if ( (!is.null(weight)) && (is.vector(weight)) ) {
nn2 = length(weight) # length of weight vector
if (nn2 != nn1) {
stop('weight must be assigned to each study.')
}
if (!match(typeof(weight), c('double', 'integer', 'logical'))) {
stop('weight must be a vector comprised by numeric numbers.')
}
} else if ( (!is.null(weight)) && (!is.vector(weight)) ) {
stop('weight must be a vector comprised by numeric numbers.')
} else {
# verbose
if ( verbose ) {
cat('\nweight is null - use default weight which depends on method.\n') #commentOut
} else {
#cat('\nweight is null - use default weight which depends on method.\n') #commentOut
}
}
} # note: weight=NULL is kept.
else if ( gmi.type == 'pvalue' ) {
if ( is.vector(gmi) ) {
if ( is.null(weight) ) {
# verbose
if ( verbose ) {
cat('\nweight is null - use default weight all one.\n') #commentOut
}
weight=rep(1, length(gmi))
} else if ( is.numeric(weight) ) {
if ( length(weight) != length(gmi) ) {
stop('weight must be assigned to each trial.')
}
} else {
stop('weight must be a vector comprised by numeric numbers.')
}
} else {
stop('input for p-value combination must be a vector.')
}
}
else if ( gmi.type == '2x2' ) {
#if (is.null(weight)) {
# weight = rep(1, dim(gmi)[1])
#}
if ( !is.null(weight) ) {
stopifnot( is.numeric(weight), length(weight) == dim(gmi)[1] )
} # note: weight=NULL is kept.
}
else {
stop('input type must be cd, pivot, pvalue or 2x2.')
}
# return
return(weight)
}
### checking gmo.xgrid: cd, pivot, exact
#####
gmeta.check.gmo.xgrid <- function(gmi, gmi.type, gmo.xgrid) {
if (gmi.type == 'cd') {
gmo.xgrid <- gmeta.check.gmo.xgrid.cd(gmi, gmo.xgrid)
}
else if (gmi.type == 'pivot') {
gmo.xgrid <- gmeta.check.gmo.xgrid.pivot(gmi, gmo.xgrid)
}
else if (gmi.type == 'pvalue') {
gmo.xgrid <- gmeta.check.gmo.xgrid.pvalue(gmi, gmo.xgrid)
}
else if (gmi.type == '2x2') {
gmo.xgrid <- gmeta.check.gmo.xgrid.2x2(gmi, gmo.xgrid)
}
else {
stop('input type must be cd, pivot, pvalue or 2x2!')
}
return(gmo.xgrid)
}
#####
gmeta.check.gmo.xgrid.cd <- function(gmi, gmo.xgrid) {
if( is.null(gmo.xgrid) ) {
nn1 <- length(gmi)
gmi.x <- list()
for (i in seq(1, nn1)) {
gmi.x[[i]] <- gmi[[i]][,1]
}
xl <- min(unlist(gmi.x))
xu <- max(unlist(gmi.x))
gmo.xgrid <- seq(xl, xu, by=0.001)
} else {
stopifnot( is.vector(gmo.xgrid), is.numeric(gmo.xgrid) )
}
return(gmo.xgrid)
}
#####
gmeta.check.gmo.xgrid.pivot <- function(gmi, gmo.xgrid) {
if ( is.null(gmo.xgrid) ) {
gmo.xgrid <- seq(from=-1, to=1, by=0.001)
} else {
stopifnot( is.vector(gmo.xgrid), is.numeric(gmo.xgrid) )
}
return(gmo.xgrid)
}
#####
gmeta.check.gmo.xgrid.pvalue <- function(gmi, gmo.xgrid) {
if ( !is.null(gmo.xgrid) ) {
warning('gmo.xgrid is not used in p-value combination.')
}
return(NULL)
}
#####
gmeta.check.gmo.xgrid.2x2 <- function(gmi, gmo.xgrid){
if( is.null(gmo.xgrid) ) {
gmo.xgrid = seq(from=-1, to=1, by=0.001)
} else {
stopifnot( is.vector(gmo.xgrid), is.numeric(gmo.xgrid), length(gmo.xgrid) > 2 )
}
return(gmo.xgrid)
}
### checking study.names
gmeta.check.study.names <- function(gmi, study.names) {
# nn1 - number of study
if ( is.matrix(gmi) || is.data.frame(gmi) ) {
nn1 = dim(gmi)[1]
} else if ( is.list(gmi) || is.vector(gmi) ) {
nn1 = length(gmi)
} else {
stop('impossible')
}
# nn2 - number of study.names
if ( is.null(study.names) ) {
# study.names
if ( is.list(gmi) && !is.data.frame(gmi) && !is.matrix(gmi) ) { # just a list
study.names <- names(gmi)
} else if ( is.data.frame(gmi) || is.matrix(gmi) ) { # R - is.list(data.frame)=TRUE
study.names <- rownames(gmi)
} else {
study.names = paste('study-', formatC(c(1:nn1),width=ceiling(log10(nn1)),format='d',flag='0'), sep='')
}
# study.names - again
if ( is.null(study.names) || length(unique(study.names) == nn1) ) {
study.names = paste('study-', formatC(c(1:nn1),width=ceiling(log10(nn1)),format='d',flag='0'), sep='')
}
} else if ( is.vector(study.names) ) {
nn2 = length(study.names)
if ( !(nn1 == nn2) ) {
stop('study.names must assign each study a name.')
}
} else {
stop('if not NULL, study.names must be a vector that assigns each study a name.')
}
# return
return(study.names)
}
# preprocessing[done]
# A unifying framework for meta-analysis
## meta-analysis - combine p-value vector
# *****************************************************************************
# gmeta.p() - combination of p-value vector
# *****************************************************************************
## main [gmeta.p()]
gmeta.p <- function(gmi, method) {
cmbd.pvalue <- Cpvaluecombine(gmi, method)
gmeta.cmbd <- list(individual.pvalues=gmi, method=method, combined.pvalue=cmbd.pvalue)
return(gmeta.cmbd)
}
### combine: p-value vector
### Cpvaluecombine: implemented in C
Cpvaluecombine <- function(RpVec, Rmethod) {
stopifnot(is.numeric(RpVec), is.character(Rmethod));
stopifnot(Rmethod %in% c('fisher', 'normal', 'stouffer', 'min', 'tippett', 'max', 'sum'));
.Call('pvaluecombine', # pvaluecombine() in C
RpVec, # generatlized meta-analysis input: vector of p-values.
Rmethod, # generatlized meta-analysis method: 'fisher', 'normal', 'stouffer', 'min', 'tippett', 'max', 'sum'.
package = 'pvaluecombine')
}
### combine: p-value vector
### R.pvalue.combine: implemented in R (deprecated)
R.pvalue.combine <- function(gmi, method) {
k = length(gmi)
if (method == 'fisher') {
cmbd.pvalue <- 1 - pchisq(-2*sum(log(gmi)), df=2*k)
}
else if (method == 'normal') {
cmbd.pvalue <- pnorm(sum(qnorm(gmi))/sqrt(k))
}
else if (method == 'stouffer') {
# alias for normal
cmbd.pvalue <- pnorm(sum(qnorm(gmi))/sqrt(k))
}
else if (method == 'tippett') {
cmbd.pvalue <- 1-(1-min(gmi))^k
}
else if (method == 'max') {
cmbd.pvalue <- max(gmi)^k
}
else if (method == 'min') {
# alias for 'tippett'
cmbd.pvalue <- 1-(1-min(gmi))^k
}
else if (method == 'sum') {
if (k <= 30) {
cmbd.pvalue <- psumunif(sum(gmi), k)
}
else {
cmbd.pvalue <- pnorm(sum(gmi), mean = k/2, sd = sqrt(k/12))
}
}
else {
print(method)
stop('method is not recognized')
}
# return
return(cmbd.pvalue)
# gmeta.cmbd <- list(indiv.pvalues=gmi, method=method, cmbd.pvalue=c0mbd.pvalue)
# return(gmeta.cmbd)
}
## [p-value combination tools]
### compute pdf of sum of n uniform(0,1), used in pvalue-based method
### input: q - quantile; n - number of uniform(0,1) r.v.s to sum.
### designed for n < 30, if n > 30 then use normal approximation.
psumunif <- function(q, n) {
psunif <- 0
for (i in 0:n) {
psunif = psunif + (-1)^i * choose(n,i) * max(0, q - i)^n
}
psunif = psunif / factorial(n)
psunif
}
## [p-value combination tools - done]
## meta-analysis - combine p-value vector [done]
## meta-analysis - model based meta-analysis
# *****************************************************************************
# gmeta.m() - model based meta-analysis methods
# *****************************************************************************
## main[gmeta.m]
## meta-analysis - model based meta-analysis
## main[gmeta.m]
gmeta.m <- function(gmi, gmi.type, method, linkfunc, weight, gmo.xgrid, ci.level, tau2, verbose) {
gmeta.data <- gmeta.cdpvt.dproc(gmi, gmi.type, method, gmo.xgrid, tau2)
gmeta.cmbd <- gmeta.cdpvt.combine(gmeta.data, method, linkfunc, weight, gmo.xgrid, ci.level, verbose)
#class(gmeta.cmbd) <- c('gmeta.m')
return(gmeta.cmbd)
}
### data processing
gmeta.cdpvt.dproc <- function(gmi, gmi.type, method, gmo.xgrid, tau2) {
if (gmi.type == 'cd') {
gmeta.data <- gmeta.cdpvt.dproc.cd(gmi, method, gmo.xgrid, tau2)
} else { # gmi.type == 'pivot'
gmeta.data <- gmeta.cdpvt.dproc.pivot(gmi, method, gmo.xgrid, tau2)
}
return(gmeta.data)
}
#####
gmeta.cdpvt.dproc.cd <- function(gmi, method, gmo.xgrid, tau2) {
# unified data structure
# x.grids and cdf values
gmi.x <- list()
gmi.cd <- list()
# maintain mean & stddev
nn1 <- length(gmi)
gmi.theta <- numeric(nn1)
gmi.sigma <- numeric(nn1)
# unifying data structure:
# x.grids, cdf-values-on-xgrids, mean, stddev for each study
for ( i in seq(1:nn1) ) {
gmi.x[[i]] <- gmi[[i]][,1]
gmi.cd[[i]] <- gmi[[i]][,2]
gmi.theta[i] <- gmeta.cd.mean(gmi.x[[i]], gmi.cd[[i]])
gmi.sigma[i] <- gmeta.cd.stddev(gmi.x[[i]], gmi.cd[[i]])
}
# calculate tau2
if ( is.element(method, c('fixed-mle',
'fixed-robust1',
'fixed-robust2',
'fixed-robust2(sqrt12)')) ) {
gmi.tau2 <- 0
} else if ( is.element(method, c('random-mm',
'random-reml',
#'random-tau2',
'random-robust1',
'random-robust2',
'random-robust2(sqrt12)')) ) {
#gmi.tau2 <- Gtau2(gmi.theta, gmi.sigma, method)
if ( is.numeric(tau2) ) {
gmi.tau2 <- tau2
} else if ( is.character(tau2) && is.element(tau2, c('DL', 'HS', 'SJ', 'HE', 'ML', 'REML', 'EB')) ) {
gmi.tau2 <- Gtau2m(gmi.theta, gmi.sigma, tau2method=tau2)
} else {
gmi.tau2 <- Gtau2(gmi.theta, gmi.sigma, method)
}
} else if ( method == 'random-tau2' ) {
if ( is.numeric(tau2) ) {
gmi.tau2 <- tau2
} else if ( is.character(tau2) && is.element(tau2, c('DL', 'HS', 'SJ', 'HE', 'ML', 'REML', 'EB')) ) {
gmi.tau2 <- Gtau2m(gmi.theta, gmi.sigma, tau2method=tau2)
} else {
stop("tau2 not recongize: must specify a numeric tau2 or a tau2method in c('DL', 'HS', 'SJ', 'HE', 'ML', 'REML', 'EB').")
}
} else {
stop('method not recongize.')
}
# unified x.grids
# xl <- min(unlist(gmi.x))
# xu <- max(unlist(gmi.x))
# gmi.unix <- seq(from=xl, to=xu, length.out=n)
gmi.unix <- gmo.xgrid # update in v2.0
# calculate cdf-value on unified x.grids
gmi.unix.cd <- NULL
for ( i in seq(1,nn1) ) {
gmi.unix.cd <- rbind(gmi.unix.cd,
approx(x=gmi.x[[i]], y=gmi.cd[[i]], xout=gmi.unix, rule=2)$y)
}
# return
gmeta.data <- list( # x.grids
gmi.x = gmi.unix, # [unified x.grids - gmo.xgrids]
# individual CDs
gmi.cd = gmi.unix.cd,
gmi.theta = gmi.theta,
gmi.sigma = gmi.sigma,
gmi.tau2 = gmi.tau2 )
# return
return(gmeta.data)
}
#####
gmeta.cdpvt.dproc.pivot <- function(gmi, method, gmo.xgrid, tau2) {
# unified data structure
# x.grids and cdf values
#gmi.x <- list()
#gmi.cd <- list()
# maintain mean & stddev
nn1 <- dim(gmi)[1]
gmi.theta <- gmi[,1]
gmi.sigma <- gmi[,2]
# calculate tau2
if ( is.element(method, c('fixed-mle',
'fixed-robust1',
'fixed-robust2',
'fixed-robust2(sqrt12)')) ) {
gmi.tau2 <- 0
} else if ( is.element(method, c('random-mm',
'random-reml',
#'random-tau2',
'random-robust1',
'random-robust2',
'random-robust2(sqrt12)')) ) {
#gmi.tau2 <- Gtau2(gmi.theta, gmi.sigma, method)
if ( is.numeric(tau2) ) {
gmi.tau2 <- tau2
} else if ( is.character(tau2) && is.element(tau2, c('DL', 'HS', 'SJ', 'HE', 'ML', 'REML', 'EB')) ) {
gmi.tau2 <- Gtau2m(gmi.theta, gmi.sigma, tau2method=tau2)
} else {
gmi.tau2 <- Gtau2(gmi.theta, gmi.sigma, method)
}
} else if ( method == 'random-tau2' ) {
if ( is.null(tau2) ) {
stop("tau2 must not NULL if method='random-tau2'")
}
if ( is.numeric(tau2) ) {
gmi.tau2 <- tau2
} else if ( is.character(tau2) && is.element(tau2, c('DL', 'HS', 'SJ', 'HE', 'ML', 'REML', 'EB')) ) {
gmi.tau2 <- Gtau2m(gmi.theta, gmi.sigma, tau2method=tau2)
} else {
stop("tau2 not recongize: must specify a numeric tau2 or a tau2method in c('DL', 'HS', 'SJ', 'HE', 'ML', 'REML', 'EB').")
}
} else {
stop('method not recongize.')
}
# unified x.grids
# xl <- min(unlist(gmi.x))
# xu <- max(unlist(gmi.x))
# gmi.unix <- seq(from=xl, to=xu, length.out=n)
gmi.unix <- gmo.xgrid # update in v2.0
# calculate cdf-value on unified x.grids
gmi.unix.cd <- NULL
if (is.element(method, c('fixed-mle',
'fixed-robust1',
'fixed-robust2',
'fixed-robust2(sqrt12)'))) {
for (i in seq(1,nn1)) {
gmi.unix.cd <- rbind(gmi.unix.cd,
pnorm(gmi.unix, gmi.theta[i], gmi.sigma[i]))
}
} else if (is.element(method, c('random-mm',
'random-reml',
'random-tau2',
'random-robust1',
'random-robust2',
'random-robust2(sqrt12)'))) {
for (i in seq(1,nn1)) {
gmi.unix.cd <- rbind(gmi.unix.cd,
pnorm(gmi.unix, gmi.theta[i], sqrt(gmi.sigma[i]^2+gmi.tau2)))
}
} else {
stop('method not recongize.')
}
# return
gmeta.data <- list( # x.grids
gmi.x = gmi.unix, # [unified x.grids - gmo.xgrids]
# individual CDs
gmi.cd = gmi.unix.cd,
gmi.theta = gmi.theta,
gmi.sigma = gmi.sigma,
gmi.tau2 = gmi.tau2 )
# return
return(gmeta.data)
}
### combine: gaussian & double exponential
gmeta.cdpvt.combine <- function(gmeta.data, method, linkfunc, weight, gmo.xgrid, ci.level, verbose) {
# use linkfunc combine
if ( linkfunc == 'inverse-normal-cdf' ) {
cmbdF <- gmeta.cdpvt.combine.gaussian(gmeta.data, method, weight, gmo.xgrid, verbose)
} else if ( linkfunc == 'inverse-laplace-cdf' ) {
cmbdF <- gmeta.cdpvt.combine.de(gmeta.data, method, weight, gmo.xgrid, verbose)
} else {
stop('linkfunc not recongize.')
}
# post-processing
cmbdf <- F2f(gmeta.data$gmi.x, cmbdF)
cmbdmn <- gmeta.cd.mean(gmeta.data$gmi.x, cmbdF)
cmbdmdn <- gmeta.cd.median(gmeta.data$gmi.x, cmbdF)
cmbdsdv <- gmeta.cd.stddev(gmeta.data$gmi.x, cmbdF)
# post-processing on individual CDs
# number of study
K = dim(gmeta.data$gmi.cd)[1]
# significance level
alpha = 1 - ci.level
# calculate individual and combined CIs
indiv.ci <- NULL
for (i in 1:K) {
indiv.ci <- rbind(indiv.ci,
gmeta.cd.mdncis(gmeta.data$gmi.x, gmeta.data$gmi.cd[i, ], alpha))
}
combined.ci <- gmeta.cd.mdncis(gmeta.data$gmi.x, cmbdF, alpha)[c(1,3)]
# return
gmeta.cmbd <- list( # x.grids
x.grids = gmeta.data$gmi.x,
# individual CDs
individual.cds = gmeta.data$gmi.cd,
individual.means = gmeta.data$gmi.theta,
individual.stddevs = gmeta.data$gmi.sigma,
individual.medians = indiv.ci[,2],
individual.cis = indiv.ci[,c(1,3)],
# combined CD
combined.cd = cmbdF,
combined.density = cmbdf,
combined.mean = cmbdmn,
combined.sd = cmbdsdv,
combined.median = cmbdmdn,
combined.ci = combined.ci,
# other information
method = method,
linkfunc = linkfunc,
weight = weight,
tau2 = gmeta.data$gmi.tau2,
ci.level = ci.level,
verbose = verbose )
# return
return(gmeta.cmbd)
}
#####
gmeta.cdpvt.combine.gaussian <- function(gmeta.data, method, weight, gmo.xgrid, verbose) {
# extract data structure
gmi.x <- gmeta.data$gmi.x
gmi.cd <- gmeta.data$gmi.cd
gmi.tau2 <- gmeta.data$gmi.tau2
gmi.theta <- gmeta.data$gmi.theta
gmi.sigma <- gmeta.data$gmi.sigma
# specify study-specific weight
if ( is.null(weight) ) {
if ( verbose ) {
cat('\nlinkfunc is inverse-normal-cdf - default weight is inverse-variance.\n') #commentOut
}
weight <- 1 / sqrt(gmi.sigma^2 + gmi.tau2) # tau2 is zero if fixed-effect model, method suppress tau2 if random-effects model
} else {
if ( verbose ) {
cat('\nuser-specified study-specific weights.\n') #commentOut
} else {
#cat('\nuser-specified study-specific weights.\n') #commentOut
}
}
# meta via CD-framework
if (method == 'fixed-mle') {
#adptws <- 1 / gmi.sigma
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(weight*qnorm(gmi.cd), 2, sum) / sqrt(weight.ss))
}
#else if (method == 'fixed-bayesian') {
# #adptws <- 1 / gmi.sigma
# #weight <- weight * adptws [no adptws in v2.0]
# weight.ss <- sum(weight^2)
# cmbdF = pnorm(apply(weight*qnorm(gmi.cd), 2, sum) / sqrt(weight.ss))
#}
else if (method == 'fixed-robust1') {
k <- length(weight)
#adptws <- 1 / gmi.sigma
#weight <- weight * adptws [no adptws in v2.0]
e.kernel <- function(u00, h) {
1/h * dnorm(u00/h)
}
data.iqr <- qnorm(0.75, gmi.theta, gmi.sigma) - qnorm(0.25, gmi.theta, gmi.sigma)
cmbdm.all <- NULL
cmbdF.all <- NULL
for (i in 1:k) {
# ith kernel weights
knwsi <- e.kernel(gmi.theta-gmi.theta[i], sqrt(data.iqr))
# adjusted weights [adjusted by kernel]
wts.adj <- weight * knwsi
weight.ss <- sum(wts.adj^2)
# ith kernel combined Fs
Fknl <- pnorm(apply(wts.adj*qnorm(gmi.cd), 2, sum)/sqrt(weight.ss))
# ith kernel combined F's median
mknl <- gmeta.cd.median(gmi.x, Fknl)
# save for later analysis
cmbdm.all <- c(cmbdm.all, mknl)
cmbdF.all <- rbind(cmbdF.all, Fknl)
}
# using the median one as robust one
srt.list <- sort.list(cmbdm.all)
if (k %% 2) {
# k is odd, median is (k+1)/2
cmbdF <- cmbdF.all[srt.list[(k+1)/2], ]
}
else {
# k is even, median is k/2 and k/2 + 1, combined with equal weight one
cmbdF <- pnorm( (qnorm(cmbdF.all[srt.list[k/2], ]) + qnorm(cmbdF.all[srt.list[k/2+1], ]))/sqrt(2) )
}
}
else if (method == 'fixed-robust2') {
cmbdFr0 <- gmeta.cdpvt.combine.gaussian(gmeta.data,
method='fixed-robust2(sqrt12)', weight, gmo.xgrid, verbose)
# latex format
# $\hat{theta^{(c)}) = H(c)^{-1}(1/2)$
x.median <- gmeta.cd.median(gmi.x, cmbdFr0)
# ith cd at x.median
nn1 <- length(weight)
mdn.cds <- numeric(nn1)
for (i in seq(1,nn1)) {
mdn.cds[i] <- approx(x=gmi.x, y=gmi.cd[i,], xout=x.median)$y
}
weight.ss <- sum(weight^2*(mdn.cds-0.5)^2)
cmbdF = pnorm(apply(weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else if (method == 'fixed-robust2(sqrt12)') {
#adptws <- 1 / gmi.sigma
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(sqrt(12)*weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-mm') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(weight*qnorm(gmi.cd), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-reml') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(weight*qnorm(gmi.cd), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-tau2') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(weight*qnorm(gmi.cd), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-robust1') {
k <- length(weight)
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
e.kernel <- function(u00, h) {
1/h * dnorm(u00/h)
}
data.iqr <- qnorm(0.75, gmi.theta, gmi.sigma) - qnorm(0.25, gmi.theta, gmi.sigma)
cmbdm.all <- NULL
cmbdF.all <- NULL
for (i in 1:k) {
# ith kernel weights
knwsi <- e.kernel(gmi.theta-gmi.theta[i], sqrt(data.iqr))
# adjusted weights
wts.adj <- weight * knwsi
weight.ss <- sum(wts.adj^2)
# ith kernel combined Fs
Fknl <- pnorm(apply(wts.adj*qnorm(gmi.cd), 2, sum)/sqrt(weight.ss))
# ith kernel combined F's median
mknl <- gmeta.cd.median(gmi.x, Fknl)
# save for later analysis
cmbdm.all <- c(cmbdm.all, mknl)
cmbdF.all <- rbind(cmbdF.all, Fknl)
}
# using the median one as robust one
srt.list <- sort.list(cmbdm.all)
if (k %% 2) {
# k is odd, median is (k+1)/2
cmbdF <- cmbdF.all[srt.list[(k+1)/2], ]
}
else {
# k is even, median is k/2 and k/2 + 1, combined with equal weight one
cmbdF <- pnorm( (qnorm(cmbdF.all[srt.list[k/2], ]) + qnorm(cmbdF.all[srt.list[k/2+1], ]))/sqrt(2) )
}
}
else if (method == 'random-robust2') {
cmbdFr0 <- gmeta.cdpvt.combine.gaussian(gmeta.data,
method='random-robust2(sqrt12)', weight, gmo.xgrid, verbose)
# latex format
# $\hat{theta^{(c)}) = H(c)^{-1}(1/2)$
x.median <- gmeta.cd.median(gmi.x, cmbdFr0)
# ith cd at x.median
nn1 <- length(weight)
mdn.cds <- numeric(nn1)
for (i in seq(1,nn1)) {
mdn.cds[i] <- approx(x=gmi.x, y=gmi.cd[i,], xout=x.median)$y
}
weight.ss <- sum(weight^2*(mdn.cds-0.5)^2)
cmbdF = pnorm(apply(weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-robust2(sqrt12)') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(sqrt(12)*weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else {
print(method)
stop('method is not recongized.')
}
# return
return(cmbdF)
}
#####
gmeta.cdpvt.combine.de <- function(gmeta.data, method, weight, gmo.xgrid, verbose) {
# extract data structure
gmi.x <- gmeta.data$gmi.x
gmi.cd <- gmeta.data$gmi.cd
gmi.tau2 <- gmeta.data$gmi.tau2
gmi.theta <- gmeta.data$gmi.theta
gmi.sigma <- gmeta.data$gmi.sigma
# specify study-specific weight
if ( is.null(weight) ) {
if ( is.element(method, c('fixed-mle',
'random-mm',
'random-reml',
'random-tau2')) ) {
if ( verbose ) {
cat('\nlinkfunc is inverse-laplace-cdf - default weight is all one (Bahadur Efficiency).\n') #commentOut
}
weight <- rep(1, length(gmi.theta)) # tau2 is zero if fixed-effect model, method suppress tau2 if random-effects model
} else if ( is.element(method, c('fixed-robust1',
'fixed-robust2',
'fixed-robust2(sqrt12)',
'random-robust1',
'random-robust2',
'random-robust2(sqrt12)')) ) {
if ( verbose ) {
cat('\nif use robust method, linkfunc is forced to use inverse-normal-cdf and default weight is inverse-variance.\n') #commentOut
}
weight <- 1 / sqrt(gmi.sigma^2 + gmi.tau2) # tau2 is zero if fixed-effect model, method suppress tau2 if random-effects model
} else {
stop('method not recongize.')
}
} else {
if ( ( verbose ) && !( weight == 1 ) ) {
cat("\nwarning: user-specified study-specific weights are not recommend if take linkfunc='inverse-laplace-cdf', the default weights with all ones give Bahadur Efficiency.\n") #commentOut
} else {
#cat("\nwarning: use user-specified weight which is not recommend.\n") #commentOut
}
}
# meta via CD-framework
# F0: double.exponential
# linkfunc: inverse-laplace-cdf
glog <- function(x) {
ifelse(x<=0, -5000, log(x))
}
qde <- function(x) {
ifelse(x>0.5, -glog(2*(1-x)), glog(2*x))
}
rde <- function(k) {
rexp(k)*(2*rbinom(k,1,0.5)-1)
}
# meta via CD-framework
if (method == 'fixed-mle') {
#adptws <- 1 / gmi.sigma
#weight <- weight * adptws [no adptws in v2.0]
#p.cmbd.de <- convolution.sim(rde, weight)
p.cmbd.de <- .convolution.de(rde, weight, verbose) # [update in v2.0: .convolution.de()]
cmbd.q = apply(weight*qde(gmi.cd), 2, sum)
cmbdF <- p.cmbd.de(cmbd.q)
}
#else if (method == 'fixed-bayesian') {
# #adptws <- 1 / gmi.sigma
# #weight <- weight * adptws [no adptws in v2.0]
# #p.cmbd.de <- convolution.sim(rde, weight)
# p.cmbd.de <- .convolution.de(rde, weight, verbose) # [update in v2.0: .convolution.de()]
# cmbd.q = apply(weight*qde(gmi.cd), 2, sum)
# cmbdF <- p.cmbd.de(cmbd.q)
#}
else if (method == 'fixed-robust1') {
k <- length(weight)
#adptws <- 1 / gmi.sigma
#weight <- weight * adptws [no adptws in v2.0]
e.kernel <- function(u00, h) {
1/h * dnorm(u00/h)
}
data.iqr <- qnorm(0.75, gmi.theta, gmi.sigma) - qnorm(0.25, gmi.theta, gmi.sigma)
cmbdm.all <- NULL
cmbdF.all <- NULL
for (i in 1:k) {
# ith kernel weights
knwsi <- e.kernel(gmi.theta-gmi.theta[i], sqrt(data.iqr))
# adjusted weights
wts.adj <- weight * knwsi
# ith kernel combined Fs
p.cmbd.de <- convolution.sim(rde, wts.adj)
#p.cmbd.de <- .convolution.de(rde, wts.adj, verbose=FALSE) # [update in v2.0: .convolution.de() - no need exact convolution since weight not all one]
cmbd.q = apply(wts.adj*qde(gmi.cd), 2, sum)
Fknl <- p.cmbd.de(cmbd.q)
# ith kernel combined F's median
mknl <- gmeta.cd.median(gmi.x, Fknl)
# save for later analysis
cmbdm.all <- c(cmbdm.all, mknl)
cmbdF.all <- rbind(cmbdF.all, Fknl)
}
# using the median one as robust one
srt.list <- sort.list(cmbdm.all)
if (k %% 2) {
# k is odd, median is (k+1)/2
cmbdF <- cmbdF.all[srt.list[(k+1)/2], ]
}
else {
# k is even, median is k/2 and k/2 + 1, combined with equal weight one
cmbdF <- pnorm( (qnorm(cmbdF.all[srt.list[k/2], ]) + qnorm(cmbdF.all[srt.list[k/2+1], ]))/sqrt(2) )
}
}
else if (method == 'fixed-robust2') {
cmbdFr0 <- gmeta.cdpvt.combine.de(gmeta.data,
method='fixed-robust2(sqrt12)', weight, gmo.xgrid, verbose)
# latex format
# $\hat{theta^{(c)}) = H(c)^{-1}(1/2)$
x.median <- gmeta.cd.median(gmi.x, cmbdFr0)
# ith cd at x.median
nn1 <- length(weight)
mdn.cds <- numeric(nn1)
for (i in seq(1,nn1)) {
mdn.cds[i] <- approx(x=gmi.x, y=gmi.cd[i,], xout=x.median)$y
}
weight.ss <- sum(weight^2*(mdn.cds-0.5)^2)
cmbdF = pnorm(apply(weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else if (method == 'fixed-robust2(sqrt12)') {
#adptws <- 1 / gmi.sigma
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(sqrt(12)*weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-mm') {
#adptws <- 1 / sqrt((gmi.sigma^2 + gmi.tau2))
#weight <- weight * adptws [no adptws in v2.0]
#p.cmbd.de <- convolution.sim(rde, weight)
p.cmbd.de <- .convolution.de(rde, weight, verbose) # [update in v2.0: .convolution.de()]
cmbd.q = apply(weight*qde(gmi.cd), 2, sum)
cmbdF <- p.cmbd.de(cmbd.q)
}
else if (method == 'random-reml') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
#p.cmbd.de <- convolution.sim(rde, weight)
p.cmbd.de <- .convolution.de(rde, weight, verbose) # [update in v2.0: .convolution.de()]
cmbd.q = apply(weight*qde(gmi.cd), 2, sum)
cmbdF <- p.cmbd.de(cmbd.q)
}
else if (method == 'random-tau2') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
#p.cmbd.de <- convolution.sim(rde, weight)
p.cmbd.de <- .convolution.de(rde, weight, verbose) # [update in v2.0: .convolution.de()]
cmbd.q = apply(weight*qde(gmi.cd), 2, sum)
cmbdF <- p.cmbd.de(cmbd.q)
}
else if (method == 'random-robust1') {
k <- length(weight)
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
e.kernel <- function(u00, h) {
1/h * dnorm(u00/h)
}
data.iqr <- qnorm(0.75, gmi.theta, gmi.sigma) - qnorm(0.25, gmi.theta, gmi.sigma)
cmbdm.all <- cmbdF.all <- NULL
for (i in 1:k) {
# ith kernel weights
knwsi <- e.kernel(gmi.theta-gmi.theta[i], sqrt(data.iqr))
# adjusted weights
wts.adj <- weight * knwsi
# ith kernel combined Fs
p.cmbd.de <- convolution.sim(rde, wts.adj)
#p.cmbd.de <- .convolution.de(rde, wts.adj, verbose=FALSE) # [update in v2.0: .convolution.de() - no need exact convolution since weight not all one]
cmbd.q = apply(wts.adj*qde(gmi.cd), 2, sum)
Fknl <- p.cmbd.de(cmbd.q)
# ith kernel combined F's median
mknl <- gmeta.cd.median(gmi.x, Fknl)
# save for later analysis
cmbdm.all <- c(cmbdm.all, mknl)
cmbdF.all <- rbind(cmbdF.all, Fknl)
}
# using the median one as robust one
srt.list <- sort.list(cmbdm.all)
if (k %% 2) {
# k is odd, median is (k+1)/2
cmbdF <- cmbdF.all[srt.list[(k+1)/2], ]
}
else {
# k is even, median is k/2 and k/2 + 1, combined with equal weight one
cmbdF <- pnorm( (qnorm(cmbdF.all[srt.list[k/2], ]) + qnorm(cmbdF.all[srt.list[k/2+1], ]))/sqrt(2) )
}
}
else if (method == 'random-robust2') {
cmbdFr0 <- gmeta.cdpvt.combine.de(gmeta.data,
method='random-robust2(sqrt12)', weight, gmo.xgrid, verbose)
# latex format
# $\hat{theta^{(c)}) = H(c)^{-1}(1/2)$
x.median <- gmeta.cd.median(gmi.x, cmbdFr0)
# ith cd at x.median
nn1 <- length(weight)
mdn.cds <- numeric(nn1)
for (i in seq(1,nn1)) {
mdn.cds[i] <- approx(x=gmi.x, y=gmi.cd[i,], xout=x.median)$y
}
weight.ss <- sum(weight^2*(mdn.cds-0.5)^2)
cmbdF = pnorm(apply(weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else if (method == 'random-robust2(sqrt12)') {
#adptws <- 1 / sqrt(gmi.sigma^2 + gmi.tau2)
#weight <- weight * adptws [no adptws in v2.0]
weight.ss <- sum(weight^2)
cmbdF = pnorm(apply(sqrt(12)*weight*(gmi.cd-0.5), 2, sum)/sqrt(weight.ss))
}
else {
print(method)
stop('method is not recongized.')
}
# return
return(cmbdF)
}
## meta-analysis - model based meta-analysis[done]
## meta-analysis - model based meta-analysis[done]
## meta-analysis - combine evidence from 2x2 tables
# *****************************************************************************
# gmeta.e() - combine evidence from 2x2 tables (exact methods)
# *****************************************************************************
## main[gmeta.e]
gmeta.e <- function(gmi, method, weight, gmo.xgrid, ci.level, mc.iteration, eta, verbose, report.error) {
data_matrix = as.matrix(cbind(gmi[,1], gmi[,3], gmi[,2], gmi[,4], gmi[,1]+gmi[,3]))
colnames(data_matrix) = c("x","y","N","M","x+y")
if (method == 'exact1') {
gmeta.data <- gmeta.exact.indiv(data_matrix, gmo.xgrid, ci.level)
gmeta.cmbd <- gmeta.exact.combine(gmeta.data, weight, gmo.xgrid, ci.level, mc.iteration, report.error)
} else if (method == 'exact2') {
gmeta.cmbd <- gmeta.exact.LT(data_matrix, weight, gmo.xgrid, ci.level, mc.iteration, eta, verbose, report.error)
} else if (method == 'MH' || method == 'Mantel-Haenszel') {
gmeta.cmbd <- gmeta.MH(data_matrix, weight, gmo.xgrid, ci.level)
} else if (method == 'Peto') {
gmeta.cmbd <- gmeta.peto(data_matrix, weight, gmo.xgrid, ci.level)
} else {
stop('gmi.type 2x2 only match methods MH, Mantel-Haenszel, Peto, exact1, exact2.')
}
#gmeta.cmbd$method <- method
return(gmeta.cmbd)
}
### 2x2 with exact1(LLX)
##### data processing
gmeta.exact.indiv <- function(data_matrix, gmo.xgrid, ci.level) {
# number of study
K = nrow(data_matrix)
# input - individual study mean and standard deviation
gmeta.theta = log( (data_matrix[,1]*(data_matrix[,4]-data_matrix[,2])) / (data_matrix[,2]*(data_matrix[,3]-data_matrix[,1])) )
gmeta.sigma = sqrt(1/data_matrix[,1] + 1/data_matrix[,2] + 1/(data_matrix[,3]-data_matrix[,1]) + 1/(data_matrix[,4]-data_matrix[,2]))
index = is.finite(gmeta.theta) & is.finite(gmeta.sigma)
# output - individual study CDs
indiv.cds = matrix(NA, K, length(gmo.xgrid))
indiv.cis = matrix(NA, K, 2)
indiv.medians = rep(NA, K)
indiv.means = rep(NA, K)
indiv.stddevs = rep(NA, K)
# construct individual CDs
alpha <- 1 - ci.level
pfunc <- function(theta) {
return(1 - pFNCHypergeo.wrap(data_matrix[i,1],data_matrix[i,3],data_matrix[i,4],data_matrix[i,5], theta) +
dFNCHypergeo(data_matrix[i,1],data_matrix[i,3],data_matrix[i,4],data_matrix[i,5], theta) / 2)
}
for (i in 1:K) {
indiv.cds[i, ] <- sapply(exp(gmo.xgrid), pfunc)
}
indiv.cds <- ifelse(indiv.cds < 1-0.1^12, indiv.cds, 1-0.1^12) # otherwise qnorm() will resolve Inf.
# individual study - CI median
for (i in 1:K) {
indiv.cis[i, ] <- log(c(.quantileCD(pfunc, alpha/2), .quantileCD(pfunc, 1-alpha/2)))
indiv.medians[i] <- log(.quantileCD(pfunc, 0.5))
}
# individual study - mean and standard deviation
na.moment <- is.na(gmeta.theta*0)
for (i in c(1:K)[!na.moment]) {
indiv.means[i] <- gmeta.cd.mean(gmo.xgrid, indiv.cds[i,])
indiv.stddevs[i] <- gmeta.cd.stddev(gmo.xgrid, indiv.cds[i,])
}
indiv.means[na.moment] <- gmeta.theta[na.moment]
indiv.stddevs[na.moment] <- gmeta.sigma[na.moment]
# return
gmeta.data <- list( # input
data_matrix = data_matrix,
# individual CDs
indiv.cds = indiv.cds,
indiv.cis = indiv.cis,
indiv.medians = indiv.medians,
indiv.means = indiv.means,
indiv.stddevs = indiv.stddevs,
# output gridding points
x.grids = gmo.xgrid )
# return
return(gmeta.data)
}
##### combine: exact1 method
gmeta.exact.combine <- function(gmeta.data, weight, gmo.xgrid, ci.level, mc.iteration, report.error) {
#gmeta.cmbd <- gmeta.data
#x.grids <- gmeta.data$x.grids
# data_matrix
data_matrix = gmeta.data$data_matrix
# number of study
K = nrow(data_matrix)
# combine individual CDs
alpha = 1 - ci.level
xstmts <- .Estimates(data_matrix[,1:4])
# obtain weight
if ( is.null(weight) ) {
weight <- xstmts$weight.hat
} else {
warning('combine 2x2 with method="exact1", default weight=NULL is strongly recommended.')
}
indiv.cds <- gmeta.data$indiv.cds
# combined cdf
combined.cd <- weight %*% qnorm( indiv.cds )
combined.cd <- combined.cd / sqrt( sum( weight^2 ) )
combined.cd <- as.vector( pnorm(combined.cd) ) # combined CD with values on limited points
# combined CD function
combinedCDF <- function(theta) { # combined CD function, used when searching for quantiles
ff = sapply(1:K, FUN = function(j) { midp.oddsratio(data_matrix[j,1], data_matrix[j,3], data_matrix[j,4], data_matrix[j,5], theta) })
ff = ifelse(ff<1-0.1^12, ff, 1-0.1^12)
ff = pnorm( weight %*% qnorm(ff) / sqrt(sum(weight^2)) )
return(ff)
}
# inference derived from combined CD function
mn.xstmt = log(.quantileCD(combinedCDF, 1/2)) # mean
lower.ci = log(.quantileCD(combinedCDF, alpha/2)) # ci.lower
upper.ci = log(.quantileCD(combinedCDF, 1 - alpha/2)) # ci.upper
o.pvalue = 2 * min( combinedCDF(1), 1 - combinedCDF(1) ) # null hypothesis: odd-ratio==1
# adjust individual CDs and combined CD function to reduce test.size.error
if( is.na(xstmts$or.hat)==TRUE || xstmts$or.hat==0 || xstmts$or.hat==Inf ) {
mn.xstmt.adj = xstmts$or.hat
lower.ci.adj = NA
upper.ci.adj = NA
o.pvalue.adj = NA
coverage.prbblty.error = NA
coverage.prbblty.error.adj = NA
} else {
indiv.cds.adj <- t( sapply(1:K, FUN = function(j) { adjust.beta(gmeta.data$indiv.cds[j,], data_matrix[j,3], xstmts$pt.hat[j], data_matrix[j,4], xstmts$pc.hat[j]) }) )
combined.cd.adj <- weight %*% qnorm( indiv.cds.adj )
combined.cd.adj <- combined.cd.adj / sqrt( sum( weight^2 ) )
combined.cd.adj <- as.vector( pnorm(combined.cd.adj) ) # combined CD with values on limited points
# combined CD function - adjusted
combinedCDF.adj <- function(theta) { # combined CD function, used when searching for quantiles
ff = sapply(1:K, FUN = function(j) { midp.oddsratio(data_matrix[j,1], data_matrix[j,3], data_matrix[j,4], data_matrix[j,5], theta) })
ff = sapply(1:K, FUN = function(j) { adjust.beta(ff[j], data_matrix[j,3], xstmts$pt.hat[j], data_matrix[j,4], xstmts$pc.hat[j]) })
ff = ifelse(ff<1-0.1^12, ff, 1-0.1^12)
ff = pnorm( weight %*% qnorm(ff) / sqrt(sum(weight^2)) )
return(ff)
}
# inference derived from combined CD function - adjusted
mn.xstmt.adj = log(.quantileCD(combinedCDF.adj,1/2)) # mean
lower.ci.adj = log(.quantileCD(combinedCDF.adj,alpha/2)) # ci.lower
upper.ci.adj = log(.quantileCD(combinedCDF.adj,1-alpha/2)) # ci.upper
o.pvalue.adj = 2 * min(combinedCDF.adj(1), 1 - combinedCDF.adj(1)) # null hypothesis: odd-ratio==1
# calculate coverage probability error
coverage.prbblty.error <- 'Not Requested'
coverage.prbblty.error.adj <- 'Not Requested'
if ( report.error ) {
ee = test.size.error(alpha/2, data_matrix[,3], xstmts$pt.hat, data_matrix[,4], xstmts$pc.hat, weight, result.1minusS=TRUE, mc.iteration)
coverage.prbblty.error = ee$Test.Size.Error.1minusS - ee$Test.Size.Error
coverage.prbblty.error.adj = ee$Test.Size.Error.Adjusted.1minusS - ee$Test.Size.Error.Adjusted
}
}
# return
gmeta.cmbd <- list( # input
data_matrix = gmeta.data$data_matrix[,1:4],
# individual CDs
individual.cds = gmeta.data$indiv.cds,
individual.cis = gmeta.data$indiv.cis,
individual.medians = gmeta.data$indiv.medians,
individual.means = gmeta.data$indiv.means,
individual.stddevs = gmeta.data$indiv.stddevs,
# combined CD function
combined.cd = combined.cd,
combined.density = F2f(gmo.xgrid, combined.cd),
combined.mean = mn.xstmt,
combined.median = gmeta.cd.median(gmo.xgrid, combined.cd),
combined.sd = gmeta.cd.stddev(gmo.xgrid, combined.cd),
combined.ci = c(lower.ci, upper.ci),
# p-value for null hypothesis: odd-ratio==1
pvalue = o.pvalue,
# combined CD function - adjusted
combined.cd.adjusted = combined.cd.adj,
combined.density.adjusted = F2f(gmo.xgrid, combined.cd.adj),
combined.mean.adjusted = mn.xstmt.adj,
combined.median.adjusted = gmeta.cd.median(gmo.xgrid, combined.cd.adj),
combined.sd.adjusted = gmeta.cd.stddev(gmo.xgrid, combined.cd.adj),
combined.ci.adjusted = c(lower.ci.adj, upper.ci.adj),
# p-value for null hypothesis: odd-ratio==1 - adjusted
pvalue.adj = o.pvalue.adj,
# coverage probability error - check liu2012exact
coverage.prbblty.error = coverage.prbblty.error,
coverage.prbblty.error.adj = coverage.prbblty.error.adj,
# other information
method = 'exact1',
linkfunc = 'inverse-fisher-exact-test-function', #[? adjusted]
weight = weight,
tau2 = NULL,
ci.level = ci.level,
verbose = report.error,
mc.iteration = mc.iteration,
report.error = report.error,
# output gridding points
x.grids = gmo.xgrid )
# return
return(gmeta.cmbd)
}
##### overall test size error for exact1 method
test.size.error <- function(s, n.vec, p1.vec, m.vec, p0.vec, weight.vec, result.1minusS=FALSE, mc.iteration=1000000) {
# number of study
K = length(n.vec)
psi = (p1.vec[1]/(1-p1.vec[1])) / (p0.vec[1]/(1-p0.vec[1]))
# generate p.i(\psi), i=1,2,...K.
mid.p.sample = matrix(NA, nrow=mc.iteration, ncol=K) # each column is for an i
mid.p.sample.adj = matrix(NA, nrow=mc.iteration, ncol=K) # each column is for an i
for(i in 1:K){
x = rbinom(mc.iteration, n.vec[i], p1.vec[i])
y = rbinom(mc.iteration,m.vec[i],p0.vec[i])
mid.p.sample[,i] = midp.oddsratio(x, n.vec[i], m.vec[i], x + y, or=psi)
mid.p.sample.adj[,i] = adjust.beta(mid.p.sample[,i], n=n.vec[i], pt=p1.vec[i], m=m.vec[i], pc=p0.vec[i])
}
# uniform sample
uniform.sample = matrix(runif(mc.iteration*K,0,1), nrow=mc.iteration, ncol=K)
# quantile
qnorm.mid.p.sample = qnorm(mid.p.sample)
qnorm.mid.p.sample.adj = qnorm(mid.p.sample.adj)
qnorm.uniform.sample = qnorm(uniform.sample)
#
bb1 = sapply(1:K, function(i) { ww<-weight.vec/weight.vec[i]; ww[1:K<=i]<-0; ww })
bb2 = sapply(1:K, function(i) { ww<-weight.vec/weight.vec[i]; ww[1:K>=i]<-0; ww })
#
xoffset = sqrt(sum(weight.vec^2))/weight.vec*qnorm(s)
#
main = qnorm.mid.p.sample %*% bb1 + qnorm.uniform.sample %*% bb2
main.adj = qnorm.mid.p.sample.adj %*% bb1 + qnorm.uniform.sample %*% bb2
#
aa = t(xoffset-t(main))
aa.adj = t(xoffset-t(main.adj))
#
pnorm.aa = pnorm(aa)
pnorm.aa.adj = pnorm(aa.adj)
#
dd = sapply(1:K, FUN = function(i) { findInterval(pnorm.aa[,i], sort(mid.p.sample[,i])) / mc.iteration - pnorm.aa[,i] })
dd.adj = sapply(1:K, FUN = function(i) { findInterval(pnorm.aa.adj[,i], sort(mid.p.sample.adj[,i])) / mc.iteration - pnorm.aa.adj[,i] })
# calculate test.size.error
test.size.error = sum(colMeans(dd))
test.size.error.adj = sum(colMeans(dd.adj))
# calculate test.size.error.1minusS
test.size.error.1minusS = 'Not requested'
test.size.error.adj.1minusS = 'Not requested'
if( result.1minusS ) {
xoffset = sqrt(sum(weight.vec^2))/weight.vec*qnorm(1-s) # use 1 - s instead of s
#
aa = t(xoffset-t(main))
aa.adj = t(xoffset-t(main.adj))
#
pnorm.aa = pnorm(aa)
pnorm.aa.adj = pnorm(aa.adj)
#
dd = sapply(1:K, FUN = function(i) { findInterval(pnorm.aa[,i], sort(mid.p.sample[,i])) / mc.iteration - pnorm.aa[,i] })
dd.adj = sapply(1:K, FUN = function(i) { findInterval(pnorm.aa.adj[,i], sort(mid.p.sample.adj[,i])) / mc.iteration - pnorm.aa.adj[,i] })
# calculate test.size.error.1minusS
test.size.error.1minusS = sum(colMeans(dd))
test.size.error.adj.1minusS = sum(colMeans(dd.adj))
}
# return
return( list( Test.Size.Error = test.size.error,
Test.Size.Error.Adjusted = test.size.error.adj,
Test.Size.Error.1minusS = test.size.error.1minusS,
Test.Size.Error.Adjusted.1minusS = test.size.error.adj.1minusS ) )
}
### 2x2 with exact2(LT's method) - risk difference[delta = p1 - p2]
gmeta.exact.LT <- function(data_matrix, weight, gmo.xgrid, ci.level, mc.iteration, eta, verbose, report.error) {
# n=length(gmo.xgrid) - number of gridding delta [output]
# mc.iteration - number of simulation for null distribution
#if (!is.null(gmo.xgrid)) { warning('input gmo.xgrid will not be utilized, output range and griddings will be automatically set.') }#
# expand data_matrix and set gmo.xgrid
#data.mi = data_matrix[,1:4]
data.mi = data_matrix[,c(2,1,4,3)] # [update in v2.0, target p1-p2, following program is on p2-p1, so change case-ctrl order]
# number of trial
nstudy = dim(data.mi)[1]
# delta - individual studies
n1 = data.mi[,3]
n2 = data.mi[,4]
p1 = data.mi[,1]/data.mi[,3]
p2 = data.mi[,2]/data.mi[,4]
deltap = p2 - p1
# var weight - non-zero-event studies
id = (1:nstudy)[p1*p2==0]
n1[id] = data.mi[id,3]+1
n2[id] = data.mi[id,4]+1
p1[id] = (data.mi[id,1]+0.5)/(data.mi[id,3]+1)
p2[id] = (data.mi[id,2]+0.5)/(data.mi[id,4]+1)
varp = p1*(1-p1)/n1+p2*(1-p2)/n2
wght = (n1*n2/(n1+n2))/sum(n1*n2/(n1+n2))
# Mantel-Haenszel's method
mu.MH = sum(deltap*wght)
sd.MH = sqrt(sum(wght^2*varp))
ci.MH = c(mu.MH-1.96*sd.MH, mu.MH+1.96*sd.MH)
p.MH = 1 - pchisq(mu.MH^2/sd.MH^2,1)
# gridding detla
d0 = max(abs(ci.MH))
delta.grd = sort(c(0, seq(from=max(-1,-d0*15), to=min(1,d0*15), length=length(gmo.xgrid)-1)))
# exact p-values [for observed data] given true delta[risk difference]
diff.exact = function(x1, x2, n1, n2, delta.grd, n.grd=15, midp=TRUE) {
# fit
fit = binom.confint(x1, n1, 0.9995, methods='exact')
# l,u
l = fit$lower
u = fit$upper
# grd
p1.grd = seq(l, u, length=n.grd)
pnull1.tot = matrix(0, n1 + 1, n.grd)
for( b in 1:n.grd ) {
p1 = p1.grd[b]
pnull1.tot[,b] = dbinom(c(0:n1), n1, p1)
}
# df & sd
dfnull = matrix(0, n1 + 1, n2 + 1)
sdnull = matrix(0, n1 + 1, n2 + 1)
for(i in 0:n1) {
p1 = (i + 0.5) / (n1 + 1)
p2 = (c(0:n2) + 0.5) / (n2 + 1)
dfnull[i+1, ] = c(0:n2) / n2 - i / n1
sdnull[i+1, ] = sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
}
# pv1 & pv2
pv1 = numeric(0)
pv2 = numeric(0)
for( theta in delta.grd ) {
p1 = (x1 + 0.5) / (n1 + 1)
p2 = (x2 + 0.5) / (n2 + 1)
tt = (x2/n2-x1/n1-theta) / sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
# null hypothesis detla==theta
tnull = (dfnull-theta) / sdnull
pvalue1 = rep(0, n.grd)
pvalue2 = rep(0, n.grd)
error = 1e-6
for ( b in 1:n.grd ) {
p1 = p1.grd[b]
p2 = p1 + theta
if ( p2 >= 0 && p2 <= 1) {
pnull1 = pnull1.tot[,b]
pnull2 = dbinom(c(0:n2), n2, p2)
n1.adj = n1+1-max(c(1,(1:(n1+1))[cumsum(sort(pnull1))<error]))
n2.adj = n2+1-max(c(1,(1:(n2+1))[cumsum(sort(pnull2))<error]))
id1 = order(pnull1)
id2 = order(pnull2)
id1 = (id1[(n1+1):1])[1:n1.adj]
id2 = (id2[(n2+1):1])[1:n2.adj]
pnull = pnull1[id1] %*% t(pnull2[id2])
if ( midp==TRUE ) {
pvalue1[b] = sum(pnull[tnull[id1,id2]>tt])+sum(pnull[tnull[id1,id2]==tt])*0.5
pvalue2[b] = sum(pnull[tnull[id1,id2]<tt])+sum(pnull[tnull[id1,id2]==tt])*0.5
} else {
pvalue1[b] = sum(pnull[tnull[id1,id2]>tt])+sum(pnull[tnull[id1,id2]==tt])
pvalue2[b] = sum(pnull[tnull[id1,id2]<tt])+sum(pnull[tnull[id1,id2]==tt])
}
}
}
pv1=c(pv1, max(pvalue1)+(1-0.9995))
pv2=c(pv2, max(pvalue2)+(1-0.9995))
}
return(list(pv1=pv1, pv2=pv2))
}
# pv1.pool&pv2.pool are the CDs of individual studies.
pv1.pool = numeric(0)
pv2.pool = numeric(0)
for ( kk in 1:nstudy ) {
# verbose
if ( verbose ) {
cat('2x2 exact2 processing trial:', kk, '\n') # use report.error as a surrogate as verbose [use parameter verbose - update v2.0]
} # update status - due to slow speed of diff.exact().
# resolve
x1 = data.mi[kk,1]
x2 = data.mi[kk,2]
n1 = data.mi[kk,3]
n2 = data.mi[kk,4]
# fit
fit = diff.exact(x1,x2,n1,n2,delta.grd,n.grd=15, midp=TRUE)
pv1.pool = rbind(pv1.pool, fit$pv1)
pv2.pool = rbind(pv2.pool, fit$pv2)
}
n = length(gmo.xgrid)
for ( i in 1:nstudy ) {
for ( j in 1:length(gmo.xgrid)) {
pv1.pool[i,(n-j+1)] = max(pv1.pool[i,1:(n-j+1)])
pv2.pool[i,j] = max(pv2.pool[i,j:n])
}
}
# pv1.pool&pv2.pool are the CDs of individual studies.
# weight for individual study
if ( is.null(weight) ) {
weight = data.mi[,3] + data.mi[,4]
} else {
if ( verbose ) { # use report.error as a surrogate as verbose [use parameter verbose - update v2.0]
cat('\nwarning: use user specified weights, instead of the default weight determined by study size.\n')
}
}
# combine individual CDs
if ( is.numeric(eta) ) {
# eta in [0,1]
if ( ! (min(eta)>=0 && max(eta)<=1) ) {
eta = seq(0.05, 0.95, length=20)
}
# F0^{-1}()
F0inv <- function(ui) {
sum( ((ui > (1-eta)) - eta) / ( eta*(1-eta) ) )
}
} else { # eta = 'Inf'
# for the case K\to\inf:
# pv1.pool&pv2.pool - adjust to make sense ln(ui/(1-ui))
pv1.pool[pv1.pool<=0] = 1e-6
pv1.pool[pv1.pool>=1] = 1 - 1e-6
pv2.pool[pv2.pool<=0] = 1e-6
pv2.pool[pv2.pool>=1] = 1 - 1e-6
# pv2.pool = 1+(1-0.9995)*2 - pv2.pool
# F0^{-1}()
F0inv <- function(ui) { log(ui/(1-ui)) }
}
# combine recipe g_c()
gc <- function(u) { sum( sapply(u, F0inv) * weight ) }
# simulation for G_c()
T = numeric(mc.iteration)
for (b in 1:mc.iteration) {
u = runif(nstudy); T[b] = gc(u)
}
Gc = ecdf(T)
# combined CDs
Hc1 = Gc( apply(pv1.pool, 2, gc) ) # HC1 non-decreasing
Hc2 = Gc( apply(pv2.pool, 2, gc) ) # HC2 non-increasing
combined.ci = c(max(delta.grd[Hc1<0.025]), min(delta.grd[Hc2<0.025]))
# post processing
integrated2CDs <- function(cd1, cd2, delta.grd) {
mdnpt = delta.grd[ round(median(which( abs(cd1-cd2) == min(abs(cd1-cd2)) ))) ]
intgrtdcd = ifelse(delta.grd<mdnpt, cd1, 1+(1-0.9995)*2-cd2) # 1+(1-0.9995)*2-cd2
for(j in 1:n) {
intgrtdcd[j] = max(intgrtdcd[1:j])
}
intgrtdcd[intgrtdcd<=0] = 1e-6
intgrtdcd[intgrtdcd>=1] = 1 - 1e-6;
return(intgrtdcd)
}
# individual CDs
indiv.cds <- NULL
for (i in 1:nstudy) {
indiv.cds = rbind(indiv.cds, integrated2CDs(pv1.pool[i,],pv2.pool[i,],delta.grd))
}
indiv.means = numeric(0)
indiv.stddevs = numeric(0)
indiv.ci = numeric(0)
for (i in 1:nstudy) {
indiv.means[i] = gmeta.cd.mean(delta.grd, indiv.cds[i,])
indiv.stddevs[i] = gmeta.cd.stddev(delta.grd, indiv.cds[i,])
indiv.ci = rbind(indiv.ci, gmeta.cd.mdncis(delta.grd, indiv.cds[i,], 1-ci.level))
}
# combined CD
cmbdF = integrated2CDs(Hc1, Hc2, delta.grd)
cmbdf = F2f(delta.grd, cmbdF)
cmbdmn = gmeta.cd.mean(delta.grd, cmbdF)
cmbdmdn = gmeta.cd.median(delta.grd, cmbdF)
cmbdsdv = gmeta.cd.stddev(delta.grd, cmbdF)
# report.error
if ( min(gmo.xgrid) < -1 || max(gmo.xgrid) > 1 ) {
#gmo.xgrid = seq(from=-1, to=1, length=length(gmo.xgrid))
gmo.xgrid = delta.grd
if ( report.error ) {
warning('\nuse exact2method combine 2x2 tables, parameter is risk-difference, only evaluate gmo.xgrid within [-1,1]\n')
}
}
# return
gmeta.cmbd <- list( # input
data_matrix = data.mi,
# individual CDs
individual.cds = indiv.cds,
individual.cis = indiv.ci[,c(1,3)],
individual.medians = indiv.ci[,2],
individual.means = indiv.means,
individual.stddevs = indiv.stddevs,
# combined CD
combined.cd = cmbdF,
combined.density = cmbdf,
combined.mean = cmbdmn,
combined.median = cmbdmdn,
combined.sd = cmbdsdv,
combined.ci = combined.ci,
# other information
method = 'exact2',
linkfunc = 'log(u/(1-u))',
weight = weight,
tau2 = NULL,
ci.level = ci.level,
mc.iteration = mc.iteration,
eta = eta,
verbose = verbose,
report.error = report.error,
# output gridding points
x.grids = delta.grd )
# return
return(gmeta.cmbd)
}
### Mantel-Haenszel's method
gmeta.MH <- function(data_matrix, weight, gmo.xgrid, ci.level) {
# data matrix
dm = data_matrix
# number of study
K = nrow(data_matrix)
alpha = 1 - ci.level
Ni = dm[,3]+dm[,4]
R = dm[,1]*(dm[,4]-dm[,2])/Ni
S = dm[,2]*(dm[,3]-dm[,1])/Ni
P = (dm[,1]+dm[,4]-dm[,2])/Ni
Q = (dm[,2]+dm[,3]-dm[,1])/Ni
gmeta.theta = R/S
gmeta.sigma = sqrt(1/dm[,1]+1/(dm[,3]-dm[,1])+1/dm[,2]+1/(dm[,4]-dm[,2])) * gmeta.theta
# this is the same as var_U and var_US in Robins et.al(1986)'s paper when K=1.
index = is.finite(gmeta.theta) & is.finite(gmeta.sigma)
# individual CDs
indiv.cds = matrix(NA,K,length(gmo.xgrid))
indiv.cis = matrix(NA,K,2)
for (ii in 1:K) {
indiv.cds[ii,] = pnorm( ( gmo.xgrid - gmeta.theta[ii] ) / gmeta.sigma[ii] )
indiv.cis[ii,] = c(gmeta.theta[ii]+qnorm(alpha/2)*gmeta.sigma[ii], gmeta.theta[ii]-qnorm(alpha/2)*gmeta.sigma[ii])
}
# combined CD
cmbd.mean = sum(R) / sum(S)
cmbd.stddev = sqrt(sum(P*R)/2/sum(R)^2+sum(P*S+Q*R)/2/sum(R)/sum(S)+sum(Q*S)/2/sum(S)^2) * cmbd.mean
cmbd.cd = pnorm( ( gmo.xgrid - cmbd.mean ) / cmbd.stddev )
cmbd.density = dnorm( ( gmo.xgrid - cmbd.mean ) / cmbd.stddev ) / cmbd.stddev
cmbd.ci = c(cmbd.mean+qnorm(alpha/2)*cmbd.stddev, cmbd.mean-qnorm(alpha/2)*cmbd.stddev)
# weight - MH default
if (!is.null(weight)) { warning('user specific weight will be surrogated by MH default weight.') }
weight = gmeta.sigma*S / (cmbd.stddev*sum(S))
# return
gmeta.cmbd <- list( # input
data_matrix = data_matrix,
# individual CDs
individual.cds = indiv.cds,
individual.cis = indiv.cis,
individual.medians = gmeta.theta,
individual.means = gmeta.theta,
individual.stddevs = gmeta.sigma,
# combined CD
combined.cd = cmbd.cd,
combined.density = cmbd.density,
combined.mean = cmbd.mean,
combined.median = cmbd.mean,
combined.sd = cmbd.stddev,
combined.ci = cmbd.ci,
# other information
method = 'Mantel-Haenszel',
linkfunc = 'inverse-normal-cdf',
tau2 = NULL,
weight = weight,
ci.level = ci.level,
verbose = FALSE,
# output gridding points
x.grids = gmo.xgrid )
# return
return(gmeta.cmbd)
}
### Peto's method
gmeta.peto <- function(data_matrix, weight, gmo.xgrid, ci.level) {
# data matrix
dm = data_matrix
# number of study
K = nrow(data_matrix)
alpha = 1 - ci.level
Ni = dm[,3]+dm[,4]
gmeta.theta = (dm[,1]-dm[,5]*dm[,3]/Ni)*Ni^2*(Ni-1) / (dm[,5]*dm[,3]*dm[,4]*(Ni-dm[,5]))
gmeta.sigma = 1 / sqrt(dm[,5]*dm[,3]*dm[,4]*(Ni-dm[,5])/Ni^2/(Ni-1))
index = is.finite(gmeta.theta) & is.finite(gmeta.sigma)
# individual CDs
indiv.cds = matrix(NA,K,length(gmo.xgrid))
indiv.cis = matrix(NA,K,2)
for (ii in 1:K) {
indiv.cds[ii,] = pnorm((gmo.xgrid-gmeta.theta[ii])/gmeta.sigma[ii])
indiv.cis[ii,] = c(gmeta.theta[ii]+qnorm(alpha/2)*gmeta.sigma[ii], gmeta.theta[ii]-qnorm(alpha/2)*gmeta.sigma[ii])
}
# combined CD
cmbd.mean = (sum(dm[,1])-sum(dm[,5]*dm[,3]/Ni))/sum(dm[,5]*dm[,3]*dm[,4]*(Ni-dm[,5])/Ni^2/(Ni-1))
cmbd.stddev = 1 / sqrt(sum(dm[,5]*dm[,3]*dm[,4]*(Ni-dm[,5])/Ni^2/(Ni-1)))
cmbd.cd = pnorm( ( gmo.xgrid - cmbd.mean ) / cmbd.stddev )
cmbd.density = dnorm( ( gmo.xgrid - cmbd.mean ) / cmbd.stddev ) / cmbd.stddev
cmbd.ci = c(cmbd.mean+qnorm(alpha/2)*cmbd.stddev, cmbd.mean-qnorm(alpha/2)*cmbd.stddev)
# weight - Peto's method default
if (!is.null(weight)) { warning('user specific weight will be surrogated by Peto-s method default weight.') }
weight = cmbd.stddev / gmeta.sigma
# return
# return
gmeta.cmbd <- list( # input
data_matrix = data_matrix,
# individual CDs
individual.cds = indiv.cds,
individual.cis = indiv.cis,
individual.medians = gmeta.theta,
individual.means = gmeta.theta,
individual.stddevs = gmeta.sigma,
# combined CD
combined.cd = cmbd.cd,
combined.density = cmbd.density,
combined.mean = cmbd.mean,
combined.median = cmbd.mean,
combined.sd = cmbd.stddev,
combined.ci = cmbd.ci,
# other information
method = 'Peto',
linkfunc = 'inverse-normal-cdf',
weight = weight,
tau2 = NULL,
ci.level = ci.level,
verbose = FALSE,
# output gridding points
x.grids = gmo.xgrid )
# return
return(gmeta.cmbd)
}
## meta-analysis - combine evidence from 2x2 tables[done]
# unified meta-analysis[done]
# postprocessing
# *****************************************************************************
# # postprocessing - print, summary, plot, etc..
# *****************************************************************************
# print, summary, plot, etc.
## print&summary [gmeta.p]
###print
print.gmeta.p <- function(x, ...) {
# Title
cat('\t\tP-value combination through CD-Framework\n')
# Call
cat('\nCall:\n')
print(x$call) #[*$call]: type of language
# Results
cat('\nCombine Method: ', x$method, '\n')
cat('\nCombined p-value: ', x$combined.pvalue, '\n')
# Details
cat('\nIndividual p-values:\n')
print(x$individual.pvalues) #[nicer visual display]
}
###summary
summary.gmeta.p <- function(object, ...) {
# set
object.sry <- object
# set class
class(object.sry) <- "summary.gmeta.p"
# return
return(object.sry)
}
###print of summary
print.summary.gmeta.p <- function(x, ...) {
# Title
cat('\t\tP-value combination through CD-Framework\n')
# Call
cat('\nCall:\n')
print(x$call) #[*$call]: type of language
# Results
cat('\nCombine Method: ', x$method, '\n')
cat('\nCombined p-value: ', x$combined.pvalue, '\n')
# Details
cat('\nIndividual p-values:\n')
print(x$individual.pvalues) #[nicer visual display]
}
## print&summary [gmeta.m]
###print
print.gmeta.m <- function(x, ...) {
# Title
cat('\t\tModel-Based Meta-Analysis through CD-Framework\n')
# Call
cat('\nCall:\n')
print(x$call) #[*$call]: type of language
# Results
cat('\nSummary of Combined CD:\n')
cmbd.cd.summary = data.frame( mean = format(round(x$combined.mean,4),nsmall=4),
median = format(round(x$combined.median,4),nsmall=4),
standard.deviation = format(round(x$combined.sd,4),nsmall=4) )
row.names(cmbd.cd.summary) <- 'Combined CD'
print(cmbd.cd.summary)
# Details
cat('\nCombined Confidence Distribution:\n')
# construct combined CD
cmbd.cd <- data.frame( x = x$x.grids,
density = x$combined.density, probability = x$combined.cd )
# set lower/upper bound
xl <- min(x$x.grids)
xu <- max(x$x.grids)
# count
n.grids <- sum(cmbd.cd$x >= xl & cmbd.cd$x <= xu)
# print head/tail
if (n.grids <= 10) {
print(cmbd.cd) # simple print all x-cd-points
} else {
cmbd.cdhead <- cmbd.cd[1:5, ]
cmbd.cdtail <- cmbd.cd[(n.grids-5):n.grids, ]
names(cmbd.cdtail) <- c(' ', ' ', ' ') # avoid duplicate titles
# print head/tail 5 x-cd-points
# output
print(cmbd.cdhead)
cat('\t...\n\t...\n')
print(cmbd.cdtail)
#cat('\n \n')
}
}
###summary
summary.gmeta.m <- function(object, ...) {
# set
object.sry <- list()
# set Call
object.sry$call <- object$call
# set combined CD
object.sry$cmbd <- data.frame( mean = object$combined.mean,
median = object$combined.median,
stddev = object$combined.sd,
ci.lower = object$combined.ci[1],
ci.upper = object$combined.ci[2])
row.names(object.sry$cmbd) <- 'Combined CD'
# set individual CDs
object.sry$idiv <- data.frame( mean = object$individual.means,
median = object$individual.medians,
stddev = object$individual.stddevs,
ci.lower = object$individual.cis[,1],
ci.upper = object$individual.cis[,2])
# set individual study names
row.names(object.sry$idiv) <- object$study.names
# set ci.level
object.sry$ci.level <- object$ci.level
# set number of study
object.sry$n.study <- dim(object.sry$idiv)[1]
# set class
class(object.sry) <- 'summary.gmeta.m'
# return
return(object.sry)
}
###print of summary
print.summary.gmeta.m <- function(x, ...) {
# Title
cat('\t\tModel-Based Meta-Analysis through CD-Framework\n')
# Call
cat('\nCall:\n')
print(x$call) #[*$call]: type of language
# Results
cat('\nSummary of Combined CD:\n')
print(x$cmbd)
cat('\nConfidence level =', x$ci.level, '\n')
# Details
cat('\nSummary of Individual CDs:\n')
print(x$idiv)
cat('\nConfidence level =', x$ci.level, '\n')
}
## print&summary [gmeta.e]
###print
print.gmeta.e <- function(x, ...) {
# Title
cat('\t\tExact Meta-Analysis Approach through CD-Framework\n')
# Call
cat('\nCall:\n')
print(x$call) #[*$call]: type of language
# Results
cat('\nSummary of Combined CD:\n')
cmbd.cd.summary = data.frame( mean = format(round(x$combined.mean,4),nsmall=4),
median = format(round(x$combined.median,4),nsmall=4),
standard.deviation = format(round(x$combined.sd,4),nsmall=4) )
row.names(cmbd.cd.summary) <- 'Combined CD'
print(cmbd.cd.summary)
# Details
cat('\nCombined Confidence Distribution:\n')
# construct combined CD
cmbd.cd <- data.frame( x = x$x.grids,
density = x$combined.density, probability = x$combined.cd )
# set lower/upper bound
xl <- min(x$x.grids)
xu <- max(x$x.grids)
# count
n.grids <- sum(cmbd.cd$x >= xl & cmbd.cd$x <= xu)
# print head/tail
if (n.grids <= 10) {
print(cmbd.cd) # simple print all x-cd-points
} else {
cmbd.cdhead <- cmbd.cd[1:5, ]
cmbd.cdtail <- cmbd.cd[(n.grids-5):n.grids, ]
names(cmbd.cdtail) <- c(' ', ' ', ' ') # avoid duplicate titles
# print head/tail 5 x-cd-points
# output
print(cmbd.cdhead)
cat('\t...\n\t...\n')
print(cmbd.cdtail)
#cat('\n \n')
}
}
###summary
summary.gmeta.e <- function(object, ...) {
# set
object.sry <- list()
# set Call
object.sry$call <- object$call
# set combined CD
object.sry$cmbd <- data.frame( mean = object$combined.mean,
median = object$combined.median,
stddev = object$combined.sd,
ci.lower = object$combined.ci[1],
ci.upper = object$combined.ci[2])
row.names(object.sry$cmbd) <- 'Combined CD'
# set individual CDs
object.sry$idiv <- data.frame( mean = object$individual.means,
median = object$individual.medians,
stddev = object$individual.stddevs,
ci.lower = object$individual.cis[,1],
ci.upper = object$individual.cis[,2])
# set individual study names
row.names(object.sry$idiv) <- object$study.names
# set ci.level
object.sry$ci.level <- object$ci.level
# set number of study
object.sry$n.study <- dim(object.sry$idiv)[1]
# set class
class(object.sry) <- 'summary.gmeta.e'
# return
return(object.sry)
}
###print of summary
print.summary.gmeta.e <- function(x, ...) {
# Title
cat('\t\tExact Meta-Analysis Approach through CD-Framework\n')
# Call
cat('\nCall:\n')
print(x$call) #[*$call]: type of language
# Results
cat('\nSummary of Combined CD:\n')
print(x$cmbd)
cat('\nConfidence level =', x$ci.level, '\n')
# Details
cat('\nSummary of Individual CDs:\n')
print(x$idiv)
cat('\nConfidence level =', x$ci.level, '\n')
}
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.