Skip to content
Snippets Groups Projects
Commit ad905bfd authored by Christian Roever's avatar Christian Roever
Browse files

moved (prediction/shrinkage) moments into separate functions

parent accdeaed
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 1
− 1
View file @ ad905bfd
......@@ -2,7 +2,7 @@ Package: bayesmeta
Type: Package
Title: Bayesian Random-Effects Meta-Analysis and Meta-Regression
Version: 2.65
Date: 2021-04-19
Date: 2021-08-03
Authors@R: c(person(given="Christian", family="Roever", role=c("aut","cre"),
email="christian.roever@med.uni-goettingen.de",
comment=c(ORCID="0000-0002-6911-698X")),
......
R/bmr.R
+ 143
− 70
View file @ ad905bfd
......@@ -722,16 +722,17 @@ bmr.default <- function(y, sigma, labels = names(y),
stopifnot(!missing(theta), is.vector(theta), is.numeric(theta), all(is.finite(theta)),
!missing(x), is.vector(x), is.numeric(x),
length(x) == d, all(is.finite(x)))
meanvec <- as.vector(support$mean %*% x)
sdvec <- rep(NA_real_, length(support$weight))
for (i in 1:length(support$weight))
sdvec[i] <- as.vector(t(x) %*% support$covariance[i,,] %*% x)
if (!mean) sdvec <- sdvec + support$tau^2
sdvec <- sqrt(sdvec)
#meanvec <- as.vector(support$mean %*% x)
#sdvec <- rep(NA_real_, length(support$weight))
#for (i in 1:length(support$weight))
# sdvec[i] <- as.vector(t(x) %*% support$covariance[i,,] %*% x)
#if (!mean) sdvec <- sdvec + support$tau^2
#sdvec <- sqrt(sdvec)
predMom <- pred.moments(x=x, mean=mean)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
for (i in 1:length(theta))
result[i] <- sum(exp(logwts + dnorm(x=theta[i], mean=meanvec, sd=sdvec, log=TRUE)))
result[i] <- sum(exp(logwts + dnorm(x=theta[i], mean=predMom[,"mean"], sd=predMom[,"sd"], log=TRUE)))
if (log) result <- log(result)
return(result)
}
......@@ -741,16 +742,17 @@ bmr.default <- function(y, sigma, labels = names(y),
stopifnot(!missing(theta), is.vector(theta), is.numeric(theta), all(is.finite(theta)),
!missing(x), is.vector(x), is.numeric(x),
length(x) == d, all(is.finite(x)))
meanvec <- as.vector(support$mean %*% x)
sdvec <- rep(NA_real_, length(support$weight))
for (i in 1:length(support$weight))
sdvec[i] <- as.vector(t(x) %*% support$covariance[i,,] %*% x)
if (!mean) sdvec <- sdvec + support$tau^2
sdvec <- sqrt(sdvec)
#meanvec <- as.vector(support$mean %*% x)
#sdvec <- rep(NA_real_, length(support$weight))
#for (i in 1:length(support$weight))
# sdvec[i] <- as.vector(t(x) %*% support$covariance[i,,] %*% x)
#if (!mean) sdvec <- sdvec + support$tau^2
#sdvec <- sqrt(sdvec)
predMom <- pred.moments(x=x, mean=mean)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
for (i in 1:length(theta))
result[i] <- sum(exp(logwts + pnorm(q=theta[i], mean=meanvec, sd=sdvec, log.p=TRUE)))
result[i] <- sum(exp(logwts + pnorm(q=theta[i], mean=predMom[,"mean"], sd=predMom[,"sd"], log.p=TRUE)))
return(result)
}
......@@ -798,19 +800,20 @@ bmr.default <- function(y, sigma, labels = names(y),
is.finite(n), n>=1, n==round(n),
!missing(x), is.vector(x), is.numeric(x),
length(x) == d, all(is.finite(x)))
meanvec <- as.vector(support$mean %*% x)
sdvec <- rep(NA_real_, length(support$weight))
for (i in 1:length(support$weight))
sdvec[i] <- as.vector(t(x) %*% support$covariance[i,,] %*% x)
if (!mean) sdvec <- sdvec + support$tau^2
sdvec <- sqrt(sdvec)
#meanvec <- as.vector(support$mean %*% x)
#sdvec <- rep(NA_real_, length(support$weight))
#for (i in 1:length(support$weight))
# sdvec[i] <- as.vector(t(x) %*% support$covariance[i,,] %*% x)
#if (!mean) sdvec <- sdvec + support$tau^2
#sdvec <- sqrt(sdvec)
predMom <- pred.moments(x=x, mean=mean)
result <- NULL
# determine numbers of samples from each mixture component:
freq <- as.vector(stats::rmultinom(1, size=n, prob=support$weight))
# draw actual samples component-wise:
for (i in 1:length(support$weight)) {
if (freq[i] > 0) {
result <- c(result, rnorm(freq[i], mean=meanvec[i], sd=sdvec[i]))
result <- c(result, rnorm(freq[i], mean=predMom[i,"mean"], sd=predMom[i,"sd"]))
}
}
# re-shuffle rows:
......@@ -872,25 +875,26 @@ bmr.default <- function(y, sigma, labels = names(y),
stopifnot(is.element(which, labels))
idx <- which(which == labels)
}
# compute predictive moments:
meanvec <- as.vector(support$mean %*% X[idx,])
sdvec <- rep(NA_real_, length(support$weight))
for (i in 1:length(support$weight))
sdvec[i] <- as.vector(t(X[idx,]) %*% support$covariance[i,,] %*% X[idx,])
sdvec <- sqrt(sdvec)
# derive shrinkage moments:
# "inverse variance" weights:
ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + support$tau^-2)
ivw2 <- support$tau^-2 / (sigma[idx]^-2 + support$tau^-2)
ivw2[support$tau==0] <- 1.0
# shrinkage mean:
meanvec <- ivw1*y[idx] + ivw2*meanvec
# shrinkage stdev.:
sdvec <- sqrt((1 / (sigma[idx]^-2 + support$tau^-2)) + (ivw2*sdvec)^2)
## compute predictive moments:
#meanvec <- as.vector(support$mean %*% X[idx,])
#sdvec <- rep(NA_real_, length(support$weight))
#for (i in 1:length(support$weight))
# sdvec[i] <- as.vector(t(X[idx,]) %*% support$covariance[i,,] %*% X[idx,])
#sdvec <- sqrt(sdvec)
## derive shrinkage moments:
## "inverse variance" weights:
#ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + support$tau^-2)
#ivw2 <- support$tau^-2 / (sigma[idx]^-2 + support$tau^-2)
#ivw2[support$tau==0] <- 1.0
## shrinkage mean:
#meanvec <- ivw1*y[idx] + ivw2*meanvec
## shrinkage stdev.:
#sdvec <- sqrt((1 / (sigma[idx]^-2 + support$tau^-2)) + (ivw2*sdvec)^2)
shrinkMom <- shrink.moments(which=which)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
for (i in 1:length(theta))
result[i] <- sum(exp(logwts + dnorm(x=theta[i], mean=meanvec, sd=sdvec, log=TRUE)))
result[i] <- sum(exp(logwts + dnorm(x=theta[i], mean=shrinkMom[,"mean"], sd=shrinkMom[,"sd"], log=TRUE)))
if (log) result <- log(result)
return(result)
}
......@@ -907,25 +911,26 @@ bmr.default <- function(y, sigma, labels = names(y),
stopifnot(is.element(which, labels))
idx <- which(which == labels)
}
# compute predictive moments:
meanvec <- as.vector(support$mean %*% X[idx,])
sdvec <- rep(NA_real_, length(support$weight))
for (i in 1:length(support$weight))
sdvec[i] <- as.vector(t(X[idx,]) %*% support$covariance[i,,] %*% X[idx,])
sdvec <- sqrt(sdvec)
# derive shrinkage moments:
# "inverse variance" weights:
ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + support$tau^-2)
ivw2 <- support$tau^-2 / (sigma[idx]^-2 + support$tau^-2)
ivw2[support$tau==0] <- 1.0
# shrinkage mean:
meanvec <- ivw1*y[idx] + ivw2*meanvec
# shrinkage stdev.:
sdvec <- sqrt((1 / (sigma[idx]^-2 + support$tau^-2)) + (ivw2*sdvec)^2)
## compute predictive moments:
#meanvec <- as.vector(support$mean %*% X[idx,])
#sdvec <- rep(NA_real_, length(support$weight))
#for (i in 1:length(support$weight))
# sdvec[i] <- as.vector(t(X[idx,]) %*% support$covariance[i,,] %*% X[idx,])
#sdvec <- sqrt(sdvec)
## derive shrinkage moments:
## "inverse variance" weights:
#ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + support$tau^-2)
#ivw2 <- support$tau^-2 / (sigma[idx]^-2 + support$tau^-2)
#ivw2[support$tau==0] <- 1.0
## shrinkage mean:
#meanvec <- ivw1*y[idx] + ivw2*meanvec
## shrinkage stdev.:
#sdvec <- sqrt((1 / (sigma[idx]^-2 + support$tau^-2)) + (ivw2*sdvec)^2)
shrinkMom <- shrink.moments(which=which)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
for (i in 1:length(theta))
result[i] <- sum(exp(logwts + pnorm(q=theta[i], mean=meanvec, sd=sdvec, log.p=TRUE)))
result[i] <- sum(exp(logwts + pnorm(q=theta[i], mean=shrinkMom[,"mean"], sd=shrinkMom[,"sd"], log.p=TRUE)))
return(result)
}
......@@ -971,28 +976,29 @@ bmr.default <- function(y, sigma, labels = names(y),
stopifnot(is.element(which, labels))
idx <- which(which == labels)
}
# compute predictive moments:
meanvec <- as.vector(support$mean %*% X[idx,])
sdvec <- rep(NA_real_, length(support$weight))
for (i in 1:length(support$weight))
sdvec[i] <- as.vector(t(X[idx,]) %*% support$covariance[i,,] %*% X[idx,])
sdvec <- sqrt(sdvec)
# derive shrinkage moments:
# "inverse variance" weights:
ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + support$tau^-2)
ivw2 <- support$tau^-2 / (sigma[idx]^-2 + support$tau^-2)
ivw2[support$tau==0] <- 1.0
# shrinkage mean:
meanvec <- ivw1*y[idx] + ivw2*meanvec
# shrinkage stdev.:
sdvec <- sqrt((1 / (sigma[idx]^-2 + support$tau^-2)) + (ivw2*sdvec)^2)
## compute predictive moments:
#meanvec <- as.vector(support$mean %*% X[idx,])
#sdvec <- rep(NA_real_, length(support$weight))
#for (i in 1:length(support$weight))
# sdvec[i] <- as.vector(t(X[idx,]) %*% support$covariance[i,,] %*% X[idx,])
#sdvec <- sqrt(sdvec)
## derive shrinkage moments:
## "inverse variance" weights:
#ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + support$tau^-2)
#ivw2 <- support$tau^-2 / (sigma[idx]^-2 + support$tau^-2)
#ivw2[support$tau==0] <- 1.0
## shrinkage mean:
#meanvec <- ivw1*y[idx] + ivw2*meanvec
## shrinkage stdev.:
#sdvec <- sqrt((1 / (sigma[idx]^-2 + support$tau^-2)) + (ivw2*sdvec)^2)
shrinkMom <- shrink.moments(which=which)
result <- NULL
# determine numbers of samples from each mixture component:
freq <- as.vector(stats::rmultinom(1, size=n, prob=support$weight))
# draw actual samples component-wise:
for (i in 1:length(support$weight)) {
if (freq[i] > 0) {
result <- c(result, rnorm(freq[i], mean=meanvec[i], sd=sdvec[i]))
result <- c(result, rnorm(freq[i], mean=shrinkMom[i,"mean"], sd=shrinkMom[i,"sd"]))
}
}
# re-shuffle rows:
......@@ -1117,7 +1123,74 @@ bmr.default <- function(y, sigma, labels = names(y),
betaHat <- as.vector(Vbeta %*% t(rbind(X,Xp)) %*% sigmaTauInv %*% c(y,yp))
return(list("mean"=betaHat, "covariance"=Vbeta))
}
pred.moments <- function(tau, x, mean=TRUE)
# posterior predictive moments (mean and std. dev.)
# conditional on given heterogeneity ("tau") values
# and covariable vector "x".
{
stopifnot(!missing(x), is.vector(x), is.numeric(x),
length(x) == d, all(is.finite(x)))
if (missing(tau)) { # return moments based on "support" object
tau <- support$tau
meanMatrix <- support$mean
covArray <- support$covariance
} else { # return moments based on supplied "tau" argument
stopifnot(is.vector(tau), all(is.finite(tau)), all(tau>=0))
meanMatrix <- matrix(NA_real_, nrow=length(tau), ncol=d,
dimnames=list(NULL, betanames))
covArray <- array(NA_real_, dim=c(length(tau), d, d),
dimnames=list(NULL, betanames, betanames))
for (i in 1:length(tau)) {
cm <- conditionalmoments(tau=tau[i])
meanMatrix[i,] <- cm$mean
covArray[i,,] <- cm$covariance
}
}
# derive predictive moments:
meanvec <- as.vector(meanMatrix %*% x)
sdvec <- rep(NA_real_, nrow(meanMatrix))
for (i in 1:nrow(meanMatrix))
sdvec[i] <- as.vector(t(x) %*% covArray[i,,] %*% x)
if (!mean) sdvec <- sdvec + tau^2
sdvec <- sqrt(sdvec)
return(cbind("mean"=meanvec, "sd"=sdvec))
}
shrink.moments <- function(tau, which)
# shrinkage moments (mean and std. dev.)
# conditional on given heterogeneity ("tau") values.
{
stopifnot(!missing(which), is.vector(which), length(which) == 1,
(is.numeric(which) | is.character(which)))
# derive (numerical) study identifier "idx":
if (is.numeric(which)) {
stopifnot(is.element(which, 1:k))
idx <- which
} else if (is.character(which)) {
stopifnot(is.element(which, labels))
idx <- which(which == labels)
}
if (missing(tau)) { # return moments based on "support" object
tauVec <- support$tau
predMom <- pred.moments(x=X[idx,])
} else { # return moments based on supplied "tau" argument
stopifnot(is.vector(tau), all(is.finite(tau)), all(tau>=0))
tauVec <- tau
predMom <- pred.moments(tau=tau, x=X[idx,])
}
# derive shrinkage moments;
# "inverse variance" weights:
ivw1 <- sigma[idx]^-2 / (sigma[idx]^-2 + tauVec^-2)
ivw2 <- tauVec^-2 / (sigma[idx]^-2 + tauVec^-2)
ivw2[tauVec==0] <- 1.0
# shrinkage means:
meanvec <- ivw1*y[idx] + ivw2*predMom[,"mean"]
# shrinkage stdevs:
sdvec <- sqrt((1 / (sigma[idx]^-2 + tauVec^-2)) + (ivw2*predMom[,"sd"])^2)
return(cbind("mean"=meanvec, "sd"=sdvec))
}
discretize <- function(delta=0.01, epsilon=0.0001)
{
divergence <- function(tau1, tau2)
......
man/bayesmeta-package.Rd
+ 1
− 1
View file @ ad905bfd
......@@ -17,7 +17,7 @@
Package: \tab bayesmeta\cr
Type: \tab Package\cr
Version: \tab 2.65\cr
Date: \tab 2021-04-19\cr
Date: \tab 2021-08-03\cr
License: \tab GPL (>=2)
}
The main functionality is provided by the \code{\link{bayesmeta}()}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment