Skip to content

Commit

Permalink
Merge pull request #6 from zackfisher/ImportantCorrectionOn-A_MBB-Cal…
Browse files Browse the repository at this point in the history 
…culation

Important correction on a mbb calculation.
  • Loading branch information
@zackfisher
zackfisher committed Feb 20, 2017
2 parents 603b951 + 152d86a commit 7e76317
Showing 1 changed file with 29 additions and 86 deletions.
115 changes: 29 additions & 86 deletions R/robu.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,102 +314,45 @@ robu <- function(formula, data, studynum,var.eff.size, userweights,
function(x){matrix(x, ncol =M)})
#ImHj <- by(data.full$ImH, data.full$study,
# function(x) as.matrix(x))

diag_one <- by(rep(1, M), X.full$study,
function(x) diag(x, nrow = length(x)))
ImHii <- Map(function(X, Q, W, D)
D - X %*% Q %*% t(X) %*% W,
X, Q.list, W.r, diag_one)

if (!user_weighting){

switch(modelweights,

HIER = { # Begin HIER

# inside = Wj^(-1/2) * (I-Hjj) * Wj^(-3/2)
inside <- Map(function(W, I)
solve(sqrt(W)) %*% I %*% solve(sqrt(W)^3),
W.r, ImHii)
I <- inside
eigenvec <- lapply(inside, function(x) eigen(x)$vectors)
eigenval <- lapply(inside, function(x) eigen(x)$values)

}, # End HIER
if (!user_weighting){

CORR = { # Begin CORR
Working_Matrx_E_j <- W.r
Working_Matrx_E <- W.r.big

} else { # Begin userweights

eigenvec <- lapply(ImHii, function(x) eigen(x)$vectors)
eigenval <- lapply(ImHii, function(x) eigen(x)$values)
I <- ImHii

} # End CORR
)
V.big <- diag(c(1), dim(Xreg)[1], dim(Xreg)[1]) %*%
diag(data.full$avg.var.eff.size)
v.j <- by(data.full$avg.var.eff.size, data.full$study,
function(x) diag(x, nrow = length(x)))
Working_Matrx_E_j <- v.j
Working_Matrx_E <- V.big
} # End userweights

} else { # Begin userweights

V.big <- diag(c(1), dim(Xreg)[1], dim(Xreg)[1]) %*%
diag(data.full$avg.var.eff.size)
V.big.list <- rep(list(V.big), N)
v.j <- by(data.full$avg.var.eff.size, data.full$study,
function(x) diag(x, nrow = length(x)))
v.j.sqrt <- lapply(v.j, function (x) sqrt(x))
inside <- Map(function(V, I)
I %*% V %*% t(I),
V.big.list, ImHj)
eigenvec <- lapply(inside, function(x) eigen(x)$vectors)
eigenval <- lapply(inside, function(x) eigen(x)$values)
I <- inside

} # End userweights

A.MBB <- Map(function (eigenvec, eigenval, k_list)
eigenvec %*%
diag(1/sqrt(eigenval), k_list, k_list) %*% t(eigenvec),
eigenvec, eigenval, k_list)
A.MBB1 <- Map(function(K, A, I)
if (K > 1) A else matrix(sqrt(solve(I))),
k_list, A.MBB, I)
A.MBB_inv_square <- Map(
function (W_E, ImH_j) {
tcrossprod(sqrt(W_E) %*% ImH_j %*%sqrt(Working_Matrx_E))
},
Working_Matrx_E_j, ImHj)

This comment has been minimized.

Sign in to view

Copy link
@jepusto

jepusto Feb 21, 2017

Contributor

This formula is incorrect. From Equation (4) in Tipton & Pustejovsky (2015), let A_j = D_j B_j^{-1/2} D_j, where D_j = W_j^{-1/2} and B_j = D_j (W_j^{-1} - X_j M X_j') D_j. Lines 334-338 should be calculating B_j. The sqrt(W_E) needs to be replaced by its inverse. Also, this calculation is quite inefficient because it involves multiplication of m matrices, each of which is n_j X N. Using the reduced form from the paper should improve efficiency.

This comment has been minimized.

Sign in to view

Copy link
@windshield999

windshield999 via email Feb 21, 2017

Collaborator

if (!user_weighting){

switch(modelweights,

HIER = { # Begin HIER

A.MBB2 <- Map(function(W, A)
solve(sqrt(W)) %*% A %*% solve(sqrt(W)),
W.r, A.MBB1)
sumXWA.MBBeeA.MBBWX.r <- Map(function(X,W,A,S)
t(X) %*% W %*% A %*% S %*% A %*% W %*%X,
X, W.r, A.MBB2, sigma.hat.r)
}, # End HIER

CORR = { # Begin CORR

A.MBB2 <- A.MBB1
sumXWA.MBBeeA.MBBWX.r <- Map(function(X,W,A,S)
t(X) %*% W %*% A %*% S %*% A %*% W %*%X,
X, W.r, A.MBB2, sigma.hat.r)
} # End CORR

)
eigenvec <- lapply(A.MBB_inv_square, function(x) eigen(x)$vectors)
eigenval <- lapply(A.MBB_inv_square, function(x) eigen(x)$values)

This comment has been minimized.

Sign in to view

Copy link
@jepusto

jepusto Feb 21, 2017

Contributor

These two lines have the effect of taking the eigen-decomposition TWICE. Better to just do

eigens <- lapply(A.MBB_inv_square, eigen)

and then refer to

eigens[[j]]$vectors
eigens[[j]]$values

as needed.

A.MBB <- Map(function (eigenvec, eigenval, k_list, W_E)
sqrt(W_E) %*% eigenvec %*% diag(1/sqrt(eigenval), k_list, k_list) %*% t(eigenvec) %*%sqrt(W_E),
eigenvec, eigenval, k_list,Working_Matrx_E_j)

This comment has been minimized.

Sign in to view

Copy link
@jepusto

jepusto Feb 21, 2017

Contributor

For consistency with Pustejovsky & Tipton (2016), the square-root of the matrix inverse here should use the Moore-Penrose generalized inverse. If

B = vec %*% diag(val) %*% t(vec)

Then set

B^{-1/2} = vec %*% ifelse(val > tol, val^{-1/2}, 0) %*% t(vec)

for some small tolerance value tol, say tol = 10^(-12). Without using the Moore-Penrose inverse, you may get errors if the matrices B_j are not of full rank (as can occur if there is a meta-regression coefficient identified in exactly one cluster).

sumXWA.MBBeeA.MBBWX.r <- Map(function(X,W,A,S)
t(X) %*% W %*% A %*% S %*% A %*% W %*%X,
X, W.r, A.MBB, sigma.hat.r)

} else { # Begin userweights

A.MBB2 <- Map(function(V, A)
V %*% A,
v.j.sqrt, A.MBB1)
sumXWA.MBBeeA.MBBWX.r <- Map(function(X,W,A,S)
t(X) %*% W %*% A %*% S %*% A %*% W %*%X,
X, W.r, A.MBB2, sigma.hat.r)
} # End userweights

sumXWA.MBBeeA.MBBWX.r <- Reduce("+", sumXWA.MBBeeA.MBBWX.r)
giTemp <- Map(function(I, A, W, X, Q)
t(I) %*% A %*% W %*% X %*% Q,
ImHj, A.MBB2, W.r, X, Q.list)
ImHj, A.MBB, W.r, X, Q.list)



Expand Down Expand Up @@ -448,7 +391,7 @@ robu <- function(formula, data, studynum,var.eff.size, userweights,
} # End small = TRUE

reg_table <- data.frame(cbind(b.r, SE, t, dfs, prob, CI.L, CI.U))
names(X.full)[2] <- "intercept"
#names(X.full)[2] <- "intercept"
labels <- c(colnames(X.full[2:length(X.full)]))
sig <- ifelse(prob < .01, "***",
ifelse(prob > .01 & prob < .05, "**",
Expand Down

0 comments on commit 7e76317

Please sign in to comment.