From ff67426aae4e4268c9d4a73a387417c4c1b67834 Mon Sep 17 00:00:00 2001
From: Christian Roever <christian.roever@med.uni-goettingen.de>
Date: Tue, 3 Aug 2021 21:00:58 +0200
Subject: [PATCH] finalized moments implementation
---
R/bmr.R | 80 ++----------------------------------------------------
man/bmr.Rd | 8 +++++-
2 files changed, 10 insertions(+), 78 deletions(-)
diff --git a/R/bmr.R b/R/bmr.R
index 67d29fb..3d84496 100644
--- a/R/bmr.R
+++ b/R/bmr.R
@@ -722,12 +722,6 @@ 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)
predMom <- pred.moments(x=x, mean=mean)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
@@ -742,12 +736,6 @@ 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)
predMom <- pred.moments(x=x, mean=mean)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
@@ -758,10 +746,6 @@ bmr.default <- function(y, sigma, labels = names(y),
qpredict <- function(p, x, mean=TRUE)
{
- #stopifnot(!missing(p), is.vector(p), is.numeric(p),
- # all(is.finite(p)), all(p >= 0), all(p <= 1),
- # !missing(x), is.vector(x), is.numeric(x),
- # length(x) == d, all(is.finite(x)))
stopifnot(!missing(p), is.vector(p), is.numeric(p),
all(is.finite(p)), all(p >= 0), all(p <= 1),
!missing(x), is.numeric(x), all(is.finite(x)),
@@ -800,12 +784,6 @@ 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)
predMom <- pred.moments(x=x, mean=mean)
result <- NULL
# determine numbers of samples from each mixture component:
@@ -824,10 +802,6 @@ bmr.default <- function(y, sigma, labels = names(y),
pred.interval <- function(level=0.95, x, mean=TRUE,
method=interval.type)
{
- #stopifnot(!missing(level), is.vector(level), is.numeric(level), length(level)==1,
- # all(is.finite(level)), all(level > 0), all(level < 1),
- # !missing(x), is.vector(x), is.numeric(x),
- # length(x) == d, all(is.finite(x)))
stopifnot(is.vector(level), is.numeric(level), length(level)==1,
all(is.finite(level)), all(level > 0), all(level < 1),
!missing(x), is.numeric(x), all(is.finite(x)),
@@ -875,21 +849,6 @@ 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)
shrinkMom <- shrink.moments(which=which)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
@@ -911,21 +870,6 @@ 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)
shrinkMom <- shrink.moments(which=which)
result <- rep(NA_real_, length(theta))
logwts <- log(support$weight)
@@ -976,21 +920,6 @@ 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)
shrinkMom <- shrink.moments(which=which)
result <- NULL
# determine numbers of samples from each mixture component:
@@ -1105,11 +1034,6 @@ bmr.default <- function(y, sigma, labels = names(y),
# posterior moments (mean and covariance) of regression parameters (beta)
# conditional on a given heterogeneity (tau) value.
{
- ##sigmaTau <- diag(sigma^2 + tau^2, nrow=k, ncol=k)
- #sigmaTauInv <- diag(1 / (sigma^2 + tau^2), nrow=k, ncol=k)
- #VbetaInv <- t(X) %*% sigmaTauInv %*% X
- #Vbeta <- solve(VbetaInv)
- #betaHat <- as.vector(Vbeta %*% t(X) %*% sigmaTauInv %*% y)
stopifnot(length(tau)==1, is.finite(tau), tau>=0)
if (!is.null(Sigmap)) { # set up block diagonal matrix:
sigmaTauInv <- cbind(rbind(diag(1 / (sigma^2 + tau^2), nrow=k, ncol=k),
@@ -1551,7 +1475,9 @@ bmr.default <- function(y, sigma, labels = names(y),
"qshrink" = qshrink,
"rshrink" = rshrink,
"shrink.interval" = shrink.interval,
- "cond.moment" = conditionalmoments,
+ "post.moments" = conditionalmoments,
+ "pred.moments" = pred.moments,
+ "shrink.moments" = shrink.moments,
"summary" = sumstats,
"interval.type" = interval.type,
"ML" = ml.estimate,
diff --git a/man/bmr.Rd b/man/bmr.Rd
index 2aaf4af..342bd38 100644
--- a/man/bmr.Rd
+++ b/man/bmr.Rd
@@ -277,9 +277,15 @@
(\eqn{i=1,\ldots,k}{i=1,...,k}). May be identified using the
\sQuote{\code{which}} argument via its index (\eqn{i}) or a character
string giving the corresponding study label.}
- \item{cond.moment}{a \code{function(tau)}
+ \item{post.moments}{a \code{function(tau)}
returning conditional posterior moments (mean and covariance) of
\eqn{\beta} as a function of \eqn{\tau}.}
+ \item{pred.moments}{a \code{function(tau, x, mean=TRUE)}
+ returning conditional posterior predictive moments (means and
+ standard deviations) as a function of \eqn{\tau}.}
+ \item{shrink.moments}{a \code{function(tau, which)}
+ returning conditional moments (means and standard deviations of
+ shrinkage distributions) as a function of \eqn{\tau}.}
\item{summary}{a \code{matrix} listing some summary statistics, namely
marginal posterior mode, median, mean, standard deviation
and a (shortest) 95\% credible intervals,
--
GitLab