# R/metaLik.R In metaLik: Likelihood Inference in Meta-Analysis and Meta-Regression Models

```.likTerms <- function(y, X, sigma2){

nbeta <- ncol(X)
ans <- list()

ans\$lik <- function(beta, tau2){

if(!all(is.finite(beta)))
return(NA)
v <- tau2+sigma2
res <- y-X%*%beta
as.numeric( -0.5*sum(log(v))-0.5*t(res)%*%diag(1/v)%*%res )
}
ans\$prof.tau2 <- function(tau2){
v <- sigma2+tau2
beta <- coef( lm.wfit(x = X, y = y, w = 1/v) )
res <- y-X%*%beta
as.numeric( -0.5*sum(log(v))-0.5*t(res)%*%diag(1/v)%*%res )
}
ans\$Jmatrix <- function(beta, tau2){

npar <- nbeta+(tau2>0)
Info <- matrix(0, ncol=npar, nrow=npar)
v <- tau2+(sigma2)
res <- y-X%*%beta
Info[1:nbeta, 1:nbeta] <- t(X)%*%diag(1/v)%*%X
if(tau2 > 0){
Info[1:(npar-1), npar] <- t(X)%*%diag(1/v^2)%*%res
Info[npar, 1:(npar-1)] <- Info[1:(npar-1), npar]
Info[npar, npar] <- -0.5*sum(1/v^2)+t(res)%*%diag(1/v^3)%*%res
}

Info
}
ans\$Imatrix <- function(beta, tau2){

npar <- nbeta+(tau2>0)
Info <- matrix(0, ncol=npar, nrow=npar)
v <- tau2+sigma2
Info[1:nbeta, 1:nbeta] <- t(X)%*%diag(1/v)%*%X
if(tau2 > 0)
Info[npar, npar] <- 0.5*sum(1/v^2)

Info
}
ans\$Smatrix <- function(beta.mle, tau2.mle, beta.constr, tau2.constr){

npar <- nbeta+1
S <- matrix(0, ncol=npar, nrow=npar)
v <- tau2.constr+sigma2
S[1:(npar-1), 1:(npar-1)] <- t(X)%*%diag(1/v)%*%X
S[1:(npar-1), npar] <- t(X)%*%diag(1/v^2)%*%X%*%(beta.mle-beta.constr)
S[npar, npar] <- 0.5*sum(1/v^2)

S
}
ans\$qvect <- function(beta.mle, tau2.mle, beta.constr, tau2.constr){

npar <- nbeta+1
qvect <- rep(0.0, npar)
v.mle <- tau2.mle+sigma2
v.constr <- tau2.constr+sigma2
qvect[1:(npar-1)] <- t(X)%*%diag(1/v.constr)%*%(X)%*%(beta.mle-beta.constr)
qvect[npar] <- -0.5*sum(1/v.mle-1/v.constr)

qvect
}
class(ans) <- c("lik.metaLik")
ans
}

metaLik <- function(formula, data, subset, contrasts = NULL, offset, sigma2, weights=1/sigma2){

call <- match.call()
if (missing(data))
data <- environment(formula)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "offset", "sigma2", "weights"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf\$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "any")
N <- length(y)
sigma2 <- model.extract(mf, sigma2)
if(is.null(sigma2)){
w <- model.extract(mf, weights)
if(is.null(w))
stop("metaLik requires within-studies variances specified either by sigma2 or by weights\n")
sigma2 <- 1/w
}
if (is.empty.model(mt))
stop("metaLik requires model specification\n")
X <- model.matrix(mt, mf, contrasts)
offset <- as.vector(model.offset(mf))
if (!is.null(offset)) {
if (length(offset) != length(y))
stop(gettextf("number of offsets is %d should equal %d (number of observations)", length(offset), length(y)), domain = NA)
}
if (is.null(offset))
offset <- rep.int(0, length(y))
y <- y-offset

## compute DL estimates
D <- diag( 1/sigma2, ncol=NROW(X), nrow=NROW(X) )
P <- diag( NROW(X)) - X %*%  solve(t(X) %*% D %*% X) %*%t (X) %*% D
q.stat <- t(y) %*% t(P) %*% D %*% P %*%y
t.value <- ( q.stat-NROW(X)+NCOL(X) ) / ( sum( diag(D) ) - sum( diag( solve( t(X) %*% D %*%X ) %*% t(X) %*% D^2 %*%X ) ) )
tau2.dl <- max(0, t.value)
weights.dl <- diag(1/(tau2.dl+sigma2), ncol=NROW(X), nrow=NROW(X))
beta.dl <- solve( t(X) %*% weights.dl %*%X ) %*% t(X) %*% weights.dl %*%y
var.dl <- solve( t(X) %*% weights.dl %*%X )

## MLE random effects model
prof.tau2 <- .likTerms(y, X, sigma2)\$prof.tau2
tau2.mle <- optimize( prof.tau2, interval=c(0, 1e+4), maximum=TRUE)\$maximum
v <- sigma2+tau2.mle
beta.mle <- coef( lm.wfit(x = X, y = y, w = 1/v) )

## MLE fixed effects model
fit.fe <- lm(y ~ 0 + X, weights=1/sigma2)
beta.fe <- coef(fit.fe)
names(beta.fe) <- names(beta.mle)

terms <- .likTerms(y, X, sigma2)
lik <- terms\$lik
Jmatrix <- terms\$Jmatrix
Imatrix <- terms\$Imatrix

vcov.mle <- try( solve(Jmatrix(beta.mle, tau2.mle)) )
if(inherits(vcov.mle, "try-error") || any(diag(vcov.mle)<0)){
vcov.mle <- try( solve(Imatrix(beta.mle, tau2.mle)) )
if(inherits(vcov.mle, "try-error") || any(diag(vcov.mle)<0))
stop("convergence not reached, perhaps too few studies.\n")
}
fitted.values <- X%*%beta.mle

colnames(vcov.mle) <- rownames(vcov.mle) <- c(colnames(X), "tau^2")

## exit
m <- structure(list(y=y+offset, X=X, fitted.values=fitted.values, sigma2=sigma2, K=NROW(X), mle=c(beta.mle, tau2.mle), vcov=vcov.mle, max.lik=lik(beta.mle, tau2.mle), beta.mle=beta.mle, tau2.mle=tau2.mle, DL=beta.dl, tau2.DL=tau2.dl, vcov.DL=var.dl, call=call, formula=formula, terms=mt,  offset=offset, contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt, mf), model=mf), class="metaLik")
m
}

print.metaLik <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\nCall:\n",
paste(deparse(x\$call), sep="\n", collapse = "\n"), "\n\n", sep="")
if(length(coef(x))) {
cat("Coefficients:\n")
print.default(format(coef(x), digits=digits),
print.gap = 2, quote = FALSE)
} else cat("No coefficients\n")
cat("\n")
cat("Variance/covariance matrix:\n")
print.default(format(vcov(x), digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("Maximum likelihood estimate of tau^2:\n")
print.default(format(x\$tau2.mle, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("Maximized log-likelihood:\n")
print.default(format(as.numeric(x\$max.lik), digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("DL coefficients:\n")
print.default(format(x\$DL, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("DL Variance/covariance matrix:\n")
print.default(format(x\$vcov.DL, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("DL tau^2 estimate:\n")
print.default(format(x\$tau2.DL, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("Number of studies K:\n")
print.default(x\$K)
cat("\n")
cat("y:\n")
print.default(format(x\$y, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("X:\n")
print.default(format(x\$X, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("offset:\n")
print.default(format(x\$offset, digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("Within-study variances:\n")
print.default(format(x\$sigma2, digits=digits),
print.gap = 2, quote = FALSE)
invisible(x)
}

test.metaLik <- function(object, param=1, value=0, alternative=c("two.sided", "less", "greater"), print=TRUE){

digits <- max(3, getOption("digits") - 3)
if(class(object)!="metaLik")
stop("\nfunction designed for 'metaLik' objects")
beta.mle <- coef(object)
if(is.character(param)){
pnames <- names(beta.mle)
param <- match(param, pnames, nomatch=NA)
}
if(missing(param) || is.na(param)){
param <- 1L
warning("\nassumed test for intercept")
}
if(length(param)>1 || param>length(coef(object)) || param<=0 || !is.numeric(param))
stop("\n'param' must be one single fixed-effects component")
alternative <- match.arg(alternative)
if(!missing(value) && (length(value)!= 1 || is.na(value)))
stop("\n'value' must be a single number")
if(class(object)!="metaLik")
stop("\n'object' must be a metaLik object")

pl <- .profLik(object, param=param)
r <- predict(pl\$smooth.rs, x=value)\$y
r.skov <- predict(pl\$smooth.rskovs, x=value)\$y

if(alternative=="less"){
pval.r <- pnorm(r)
pval.rskov <- pnorm(r.skov)
}
else if(alternative=="greater"){
pval.r <- pnorm(r, lower.tail=FALSE)
pval.rskov <- pnorm(r.skov, lower.tail=FALSE)
}
else{
pval.r <- 2*pnorm(-abs(r))
pval.rskov <- 2*pnorm(-abs(r.skov))
}
if(print){
cat("\nSigned profile log-likelihood ratio test for parameter ", names(beta.mle[param]), sep="", "\n")
cat("\nFirst-order statistic")
cat("\nr:", formatC(r, digits), ", p-value:", formatC(pval.r, digits), sep="")
cat("\nSkovgaard's statistic")
cat("\nrSkov:", formatC(r.skov, digits), ", p-value:", formatC(pval.rskov, digits), sep="")
if(alternative=="two.sided")
cat("\nalternative hypothesis: parameter is different from ", round(value, digits), sep="", "\n")
else
cat("\nalternative hypothesis: parameter is ", alternative, " than ", round(value, digits), sep="", "\n")
}
## bye
ans <- list(r=r, pvalue.r=pval.r, rskov=r.skov, pvalue.rskov=pval.rskov)
invisible(ans)
}

coef.metaLik <- function(object, ...)
object\$beta.mle

vcov.metaLik <- function(object, ...){
nbeta <- length(object\$beta.mle)
object\$vcov[1:nbeta,1:nbeta]
}

logLik.metaLik <- function(object, ...) {
structure(object\$max.lik, df = sum(sapply(object\$mle, length))+(object\$tau2.mle>0), class = "logLik")
}

model.matrix.metaLik <- function(object, ...) {

rval <- if(!is.null(object\$X)) object\$X
else model.matrix(object\$terms, model.frame(object), contrasts = object\$contrasts)
return(rval)
}

residuals.metaLik <- function(object, type = c("pearson", "response", "working"), ...){

res <- object\$y-object\$fitted.values
type <- match.arg(type)
se <- sqrt(object\$sigma2+object\$tau2.mle)
switch(type, working=, response=res, pearson=res/se)
}

summary.metaLik <- function(object, ...){

digits <- max(3, getOption("digits") - 3)
if(class(object)!="metaLik")
stop("\n'summary.metaLik' designed for 'metaLik' objects\n")
beta <- coef(object)
nbeta <- length(beta)
beta.se <- sqrt(diag(as.matrix(vcov(object))))
r <- matrix(nrow=nbeta, ncol=4)
for(i in 1:nbeta)
r[i,] <- as.numeric( test.metaLik(object, param=i, print=FALSE) )
colnames(r) <- c("r", "pvalue.r", "rskov", "pvalue.rskov")##, "check.skov")
tau2 <- object\$tau2.mle

cat("\nLikelihood inference in random-effects meta-analysis models\n")
cat("\nCall:\n", paste(deparse(object\$call), sep = "\n", collapse = "\n"), "\n", sep = "")
cat("\nEstimated heterogeneity parameter tau2")
cat(" = ", formatC(unname(object\$tau2.mle), digits), sep="")

cat("\n\nFixed-effects:\n")

ans <- cbind( beta, beta.se, r[,"r"], r[,"pvalue.r"], r[,"rskov"], r[,"pvalue.rskov"]  )
dimnames(ans) <- list(names(beta), c("estimate", "std.err.", "signed logLRT", "p-value", "Skovgaard", "p-value"))
ans <- round(ans, digits=digits)

print.default(ans, print.gap = 2)
cat("\nLog-likelihood:", format(round(object\$max.lik, digits)), "\n")
invisible(object)
}

.profLik <- function(object, param=1){

if(class(object)!="metaLik")
stop("\ninternal function designed for 'metaLik' objects")
y <- object\$y
N <- length(y)
X <- object\$X
offset <- object\$offset
y <- y-offset
sigma2 <- object\$sigma2
nbeta <- length(object\$beta.mle)
terms <- .likTerms(y, X, sigma2)
lik <- terms\$lik
Jmatrix <- terms\$Jmatrix
Imatrix <- terms\$Imatrix
Smatrix <- terms\$Smatrix
qvect <- terms\$qvect

beta.mle <- object\$beta.mle
tau2.mle <- object\$tau2.mle
par.mle <- beta.mle[param]
se.mle <- sqrt(object\$vcov[param, param])
pts <- c(seq(-5, -1, by=.5), seq(1, 5, by=.5))
values <- par.mle+pts*se.mle
rs <- us <- rep(NA, length(values))

fixed <- rep(NA, nbeta)
beta.constr <- fixed
names(beta.constr) <- names(beta.mle)

.computeU <- function(beta.mle, tau2.mle, beta.constr, tau2.constr, id.par){
S <- Smatrix(beta.mle, tau2.mle, beta.constr, tau2.constr)
q <- qvect(beta.mle, tau2.mle, beta.constr, tau2.constr)
J.mle <- Jmatrix(beta.mle, tau2.mle)
J.constr <- Jmatrix(beta.constr, tau2.constr)
I.mle <- Imatrix(beta.mle, tau2.mle)

A <- abs( solve(S) %*%q )[id.par]
B <- abs( det( J.mle ) )^(.5)
C <- abs( det( solve(I.mle) ) )
D <- abs( det( S ) )
E <- abs( det( as.matrix(J.constr[-id.par,-id.par] ) ) )^(-.5)

return( sign(beta.mle[id.par]-beta.constr[id.par])*A*B*C*D*E )
}

for(i in 1:length(values)){
fixed[param] <- beta.constr[param] <- values[i]
ynew <- y-as.matrix(X[,param])%*%values[i]
Xnew <- as.matrix(X[,-param])
tau2.constr <- 0.0
if(tau2.mle>0){
prof.tau2 <- .likTerms(ynew, Xnew, sigma2)\$prof.tau2
tau2.constr <- optimize(prof.tau2, interval=c(0, 1e+4), maximum=TRUE)\$maximum
}
if(NCOL(Xnew)>0){
v <- sigma2+tau2.constr
beta.hat <- coef( lm.wfit(x = Xnew, y = ynew, w = 1/v) )
beta.constr[is.na(fixed)] <- beta.hat
}

## likelihood ratio test statistic (and Skovgaard correction)
lrt <- 2*(lik(beta.mle, tau2.mle)-lik(beta.constr, tau2.constr))

if(is.finite(lrt) && lrt>0){
rs[i] <- sign(beta.mle[param]-beta.constr[param])*sqrt(lrt)
us[i] <- try(.computeU(beta.mle, tau2.mle, beta.constr, tau2.constr, param), silent=TRUE)
if(inherits(us[i], "try-error"))
us[i] <- NA
}
}
smooth.rs <- smooth.rskovs <-  smooth.spline( x=values, y=rs )
smooth.rs.inv <- smooth.rskovs.inv <- smooth.spline( y=values, x=rs )
rskovs <- rs+1/rs*log(us/rs)
smooth.rskovs <- smooth.spline( x=values, y=rskovs )
smooth.rskovs.inv <- smooth.spline( y=values, x=rskovs )

return( list( smooth.rs=smooth.rs, smooth.rs.inv=smooth.rs.inv, smooth.rskovs=smooth.rskovs, smooth.rskovs.inv=smooth.rskovs.inv) )
}

simulate.metaLik <- function(object, nsim = 1, seed = NULL, ...){

if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
ftd <- fitted(object)
tau2 <- object\$tau2.mle
n <- length(ftd)
nm <- names(ftd)

val <- replicate(nsim, ftd + rnorm( n, sd = sqrt(object\$sigma2+tau2) ) )
val <- as.data.frame(val)
colnames(val) <- paste("sim", seq_len(nsim), sep = "_")
if (!is.null(nm))
rownames(val) <- nm
attr(val, "seed") <- RNGstate
val
}
```

## Try the metaLik package in your browser

Any scripts or data that you put into this service are public.

metaLik documentation built on May 1, 2019, 8:44 p.m.