Revision 2ce4f3930a435f1dd3cbbfd95164c4ed3805340f authored by Unknown author on 08 August 1977, 00:00:00 UTC, committed by Gabor Csardi on 08 August 1977, 00:00:00 UTC
1 parent f9b9cc1
meta.R
meta.MH <- function ( ntrt, nctrl, ptrt, pctrl, conf.level = .95,
names = NULL, data = NULL,
subset = NULL, na.action = na.fail ) {
if (conf.level>1 & conf.level<100)
conf.level<-conf.level/100
if ( is.null( data ) )
data <- sys.frame( sys.parent() )
mf <- match.call()
mf$data <- NULL
mf$subset <- NULL
mf[[1]] <- as.name( "data.frame" )
mf <- eval( mf,data )
if ( !is.null( subset ) )
mf <- mf[subset,]
mf <- na.action( mf )
A <- mf$ptrt
B <- mf$pctrl
C <- mf$ntrt - mf$ptrt
D <- mf$nctrl - mf$pctrl
logORs <- log( ( A * D ) / ( B * C ) )
vars <- 1/A + 1/B + 1/C + 1/D
Ti <- A + B + C + D
P <- ( A + D ) / Ti
Q <- ( B + C ) / Ti
R <- A * D / Ti
S <- B * C / Ti
logMH <- log( sum( R ) / sum( S ) )
varMH <- sum( P * R ) / ( 2 * sum( R )^2 ) +
sum( P * S + Q * R ) / ( 2 * sum( R ) * sum( S ) ) +
sum( Q * S ) / ( 2 * sum( S ) ^ 2 )
MHchisq <- ( sum( A - ( A + B ) * ( A + C ) / Ti )^2 ) /
sum( ( A + B ) * ( C + D ) * ( A + C ) * ( B + D ) /
( Ti * Ti * ( Ti - 1 )
)
)
ok <- is.finite( logORs )
heterog <- sum( ( ( logORs[ok] - logMH ) ^ 2 ) / vars[ok] )
df <- sum( ok ) - 1
rval <- list( logOR = logORs, selogOR = sqrt( vars ),
logMH = logMH, selogMH = sqrt( varMH ),
MHtest = c( MHchisq, 1 - pchisq( MHchisq, 1 ) ),
het = c( heterog, df, 1 - pchisq( heterog, df ) ),
call = match.call(), names = mf$names,
conf.level = conf.level )
class( rval ) <- "meta.MH"
rval
}
# Print out the output of meta.MH function.
print.meta.MH <- function( obj ) {
cat( "Fixed effects ( Mantel-Haenszel ) Meta-Analysis\n" )
cat( "Call: " )
print( obj$call )
conf.level <- obj$conf.level
ci.value <- -qnorm( ( 1 - conf.level ) / 2 )
#lower <- mn - ci.value * se
#upper <- mn + ci.value * se
ci <- exp( obj$logMH + c( -ci.value, 0, ci.value ) * obj$selogMH )
cat( paste( "Mantel-Haensel OR = ", round( ci[2],2 ), " ",
conf.level * 100, "% CI ( ", round( ci[1],2 ),
", ", round( ci[3],2 ), " )\n", sep = ""
)
)
cat( paste( "Test for heterogeneity: X^2( ",obj$het[2]," ) = ",
round( obj$het[1],2 ), " ( p-value ",
round( obj$het[3],4 ), " )\n", sep="" ) )
}
# Summarise meta.MH function, same with print.meta.MH, but with odds
# ratios for each variable.
summary.meta.MH <- function( obj ,conf.level=NULL) {
if (is.null(conf.level))
conf.level <- obj$conf.level
if (conf.level>1 & conf.level<100)
conf.level<-conf.level/100
ci.value <- -qnorm( ( 1 - conf.level ) / 2 )
m <- exp( outer( obj$selogOR, c( 0, -ci.value, ci.value ), "*" ) +
cbind( obj$logOR, obj$logOR, obj$logOR )
)
dimnames( m ) <- list( as.character( obj$names ),
c( "OR", "(lower ", paste(100*conf.level,"% upper)",sep="")) )
rval <- list( ors = m, call = obj$call,
MHci = exp( obj$logMH + c( -ci.value, 0, ci.value ) *
obj$selogMH
),
het = obj$het,
conf.level=conf.level
)
class( rval ) <- "summary.meta.MH"
rval
}
# Print out the summarised meta.MH function.
print.summary.meta.MH <- function( obj ) {
cat( "Fixed effects ( Mantel-Haenszel ) meta-analysis\n" )
cat( "Call: " )
print( obj$call )
cat( "------------------------------------\n" )
print( round( obj$ors,2 ) )
cat( "------------------------------------\n" )
conf.level <- obj$conf.level
cat( paste( "Mantel-Haensel OR = ", round( obj$MHci[2], 2 ), " ",
conf.level*100, "% CI ( ", round( obj$MHci[1],2 ), ",",
round( obj$MHci[3],2 ), " )\n", sep=""
)
)
cat( paste( "Test for heterogeneity: X^2( ",obj$het[2]," ) = ",
round( obj$het[1],2 ), " ( p-value ",
round( obj$het[3],4 ), " )\n", sep = ""
)
)
}
# Plot the Odds Ratio of meta.MH
plot.meta.MH <- function( obj, summary = T, summlabel = "Summary", conf.level=NULL,colors=list(box="black",lines="gray",summary="black",zero="lightgray"),... ){
if (is.null(conf.level))
conf.level <- obj$conf.level
if (conf.level>1 & conf.level<100)
conf.level<-conf.level/100
if (is.null(colors))
colors<-par("fg")
if (!is.list(colors))
colors<-list(box=colors,lines=colors,summary=colors,zero=colors)
if ( summary )
metaplot( obj$logOR, obj$selogOR, labels=obj$names,
summn = obj$logMH, sumse = obj$selogMH,
conf.level = conf.level, sumnn = obj$selogMH ^ ( -2 ),
summlabel = summlabel, logeffect = T,colors=colors,... )
else
metaplot( obj$logOR, obj$selogOR, labels=obj$names, logeffect = T,colors=colors,... )
}
# Meta Analysis using DerSimonian-Laird Method.
meta.DSL <- function( ntrt, nctrl, ptrt, pctrl, conf.level = .95,
names = NULL, data = NULL,
subset = NULL, na.action = na.fail ) {
if (conf.level>1 & conf.level<100)
conf.level<-conf.level/100
if ( is.null( data ) )
data <- sys.frame( sys.parent() )
mf <- match.call()
mf$data <- NULL
mf$subset <- NULL
mf[[1]] <- as.name( "data.frame" )
mf <- eval( mf,data )
if ( !is.null( subset ) )
mf <- mf[subset,]
mf <- na.action( mf )
A <- mf$ptrt
B <- mf$pctrl
C <- mf$ntrt - mf$ptrt
D <- mf$nctrl - mf$pctrl
logORs <- log( ( A * D ) / ( B * C ) )
ok <- is.finite( logORs )
if ( any( !ok ) )
warning( "Studies with 0/Inf odds ratio omitted" )
vars <- 1/A + 1/B + 1/C + 1/D
tau2 <- max( 0, var( logORs[ok] ) - sum( vars[ok] ) / ( length( ok )-1 ) )
wts <- 1 / ( vars + tau2 )
logDSL <- sum( wts[ok] * logORs[ok] ) / sum( wts[ok] )
varDSL <- 1 / sum( wts[ok] )
summary.test <- logDSL / sqrt( varDSL )
heterog <- sum( ( ( logORs[ok] - logDSL )^2 ) / vars[ok] )
df <- sum( ok ) - 1
rval <- list( logOR = logORs, selogOR = sqrt( vars ),
logDSL = logDSL, selogDSL = sqrt( varDSL ),
test = c( summary.test,
1 - pchisq( summary.test ^ 2, 1 ) ),
het = c( heterog, df, 1 - pchisq( heterog, df ) ),
call = match.call(),
names = mf$names, conf.level = conf.level,
omitted = !ok, tau2=tau2 )
class( rval ) <- "meta.DSL"
rval
}
# Print out the output of meta.DSL function.
print.meta.DSL <- function( obj ){
conf.level <- obj$conf.level
ci.value <- -qnorm( ( 1 - conf.level ) / 2 )
cat( "Random effects ( DerSimonian-Laird ) meta-analysis\n" )
cat( "Call: " )
print( obj$call )
ci <- exp( obj$logDSL + c( -ci.value, 0, ci.value ) * obj$selogDSL )
cat( paste( "Summary OR = ", round( ci[2],2 ), " ",
conf.level * 100, "% CI ( ",round( ci[1],2 ), ", ",
round( ci[3],2 )," )\n", sep="" ) )
cat( paste( "Estimated random effects variance:",round( obj$tau2,2 ),"\n" ) )
}
# Summarise meta.DSL function, same with print.meta.DSL, but with odds
# ratios for each variable.
summary.meta.DSL <- function( obj,conf.level=NULL ){
if (is.null(conf.level))
conf.level <- obj$conf.level
if (conf.level>1 & conf.level<100)
conf.level<-conf.level/100
ci.value <- -qnorm( ( 1 - conf.level ) / 2 )
m <- exp( outer( obj$selogOR, c( 0,-ci.value,ci.value ), "*" ) +
cbind( obj$logOR, obj$logOR, obj$logOR ) )
dimnames( m ) <- list( as.character( obj$names ),
c( "OR", "(lower ", paste(100*conf.level,"% upper)",sep="")) )
rval <- list( ors = m, call = obj$call,
ci = exp( obj$logDSL + c( -ci.value, 0, ci.value ) *
obj$selogDSL ),
het = obj$het,
omitted = as.character( obj$names )[obj$omitted],
conf.level=conf.level,
tau2 = obj$tau2 )
class( rval ) <- "summary.meta.DSL"
rval
}
# Print out the summarised meta.MH function.
print.summary.meta.DSL <- function( obj ) {
conf.level <- obj$conf.level
cat( "Random effects ( DerSimonian-Laird ) meta-analysis\n" )
cat( "Call: " )
print( obj$call )
cat( "------------------------------------\n" )
print( round( obj$ors,2 ) )
cat( "------------------------------------\n" )
cat( paste( "Summary OR = ",round( obj$ci[2],2 ), " ",
conf.level * 100, "% CI ( ",round( obj$ci[1],2 ), ",",
round( obj$ci[3],2 )," )\n",sep="" ) )
cat( paste( "Test for heterogeneity: X^2( ",obj$het[2]," ) = ",
round( obj$het[1],2 ),
" ( p-value ", round( obj$het[3],4 ), " )\n", sep="" ) )
cat( paste( "Estimated random effects variance:",
round( obj$tau2,2 ), "\n" ) )
if ( length( obj$omitted )>0 ){
cat( paste( "( ", length( obj$omitted ),
"studies with zero or infinite odds ratio omitted )\n" ) )
}
}
# Plot the Odds Ratio of meta.DSL
plot.meta.DSL <- function( obj,summary=T,summlabel="Summary",conf.level=NULL,colors=list(box="black",lines="gray",summary="black",zero="lightgray"),... ){
if (is.null(conf.level))
conf.level <- obj$conf.level
if (conf.level>1 & conf.level<100)
conf.level<-conf.level/100
if (is.null(colors))
colors<-par("fg")
if (!is.list(colors))
colors<-list(box=colors,lines=colors,summary=colors,zero=colors)
if ( summary )
metaplot( obj$logOR, obj$selogOR, conf.level = conf.level,
labels = obj$names, summn = obj$logDSL,
sumse = obj$selogDSL, sumnn = obj$selogDSL^( -2 ),
summlabel = summlabel, logeffect = T,colors=colors,... )
else
metaplot( obj$logOR, obj$selogOR, labels = obj$names,
logeffect = T,colors=colors,... )
}
metaplot <- function( mn, se, nn=NULL, labels=NULL, conf.level = .95,
xlab = "Odds ratio", ylab = "Study Reference",
boxwidth = NULL, xlim = NULL, summn = NULL,
sumse = NULL, sumnn = NULL,
summlabel = "Summary", logeffect = F,
lwd = 2, boxsize = 1,
zero = if ( logeffect )
1
else
0,
colors=NULL,
... ) {
ci.value <- -qnorm( ( 1 - conf.level ) / 2 )
ok <- is.finite( mn + se )
if ( is.null( xlim ) )
xlim <- c( min( mn[ok] - ci.value * se[ok], na.rm = T ),
max( mn[ok] + ci.value * se[ok], na.rm = T ) )
##par( pty="s" )
n <- length( mn )
if ( logeffect ) {
xlog <- "x"
nxlim <- exp( xlim )
}
else {
xlog <- ""
nxlim <- xlim
}
if ( !is.null( labels ) ) {
if ( logeffect )
nxlim[1] <- nxlim[1] / sqrt( nxlim[2] / nxlim[1] )
else
nxlim[1] <- nxlim[1] - 0.5 * ( nxlim[2] - nxlim[1] )
}
par( xaxt = "n",yaxt = "n" )
plot( nxlim,c( 1,-n-2-3 * !is.null( summn ) ),
type = "n", bty = "n", xaxt = "n", yaxt = "n",
log = xlog, xlab=xlab,ylab=ylab,... )
par( xaxt = "s" )
axis( 1,at = round( 10 ^ pretty( log( exp( xlim ),10 ), 6 ), 2 ) )
if ( !is.null( zero ) )
abline( v = zero, lty = 2, lwd = 2 ,col=colors$zero)
ci.value <- -qnorm( ( 1 - conf.level ) / 2 )
lower <- mn - ci.value * se
upper <- mn + ci.value * se
if ( logeffect ){
lower <- exp( lower )
upper <- exp( upper )
}
for ( i in 1:n ){
if ( is.na( lower[i]+upper[i] ) )
next
lines( c( lower[i], upper[i] ), c( -i, -i ), lwd = lwd, col=colors$lines,... )
}
if ( !is.null( labels ) )
text( rep( nxlim[1], n ), -( 1:n ), labels,..., adj=0 )
if ( is.null( nn ) )
nn <- se ^ -2
yscale <- 0.3 * boxsize / max( sqrt( nn ), na.rm = T )
if ( logeffect ) {
scale <- ( nxlim[2] / nxlim[1] ) ^ ( yscale / ( 4 + n ) )
xl <- exp( mn ) * ( scale ^ -sqrt( nn ) )
xr <- exp( mn ) * ( scale ^ sqrt( nn ) )
}
else {
scale <- yscale * ( nxlim[2] - nxlim[1] ) / ( 4 + n )
xl <- mn - scale * sqrt( nn )
xr <- mn + scale * sqrt( nn )
}
yb <- ( 1:n ) - yscale * sqrt( nn )
yt <- ( 1:n ) + yscale * sqrt( nn )
for ( i in 1:n ) {
if ( !is.finite( mn[i] ) )
next
rect( xl[i], -yb[i], xr[i], -yt[i], col = colors$box,border= colors$box)
}
if ( !is.null( summn ) ) {
if ( logeffect ) {
x0 <- exp( summn )
xl <- exp( summn - ci.value * sumse )
xr <- exp( summn + ci.value * sumse )
}
else{
x0 <- summn
xl <- summn - ci.value * sumse
xr <- summn + ci.value * sumse
}
y0 <- n + 3
yb <- n + 3 - sqrt( sumnn ) * yscale
yt <- n + 3 + sqrt( sumnn ) * yscale
polygon( c( xl, x0, xr, x0 ), -c( y0, yt, y0, yb ),
col = colors$summary, border = colors$summary )
text( nxlim[1], -y0, labels = summlabel, adj = 0 )
}
}
Computing file changes ...