X_tilde constructed such that element-wise X_tilde * theta is not equal to X * Beta + z

[svteichman]

I may be missing another transformation elsewhere, but I've been working with the radEmu code base and I'm confused about the construction of X_tilde. As a toy example, take the following:

n <- 4
J <- 3
p <- 2
X <- matrix(c(rep(1, n), rnorm(n)), nrow = n)
B <- matrix(rnorm(J*p), ncol = J)
z <- rnorm(n)

Now, apply the code from the Firth penalty in emuFit lines 138 through 153 (ignoring weights and making slight updates so the internal functions can be run). Here my understanding is that D_tilde corresponds with X_tilde in the paper (incorporating columns of 0's to account for the z's in the theta vector).

X_tilde <- radEmu:::X_to_X_tilde(X,J)
S <- Matrix::sparseMatrix(i = 1:(n*J),
                          j = rep(1:n,each = J),
                          x = rep(1, n*J))
D_tilde <- cbind(X_tilde,S)
rownames(D_tilde) <- 1:(n*J)
B_tilde <- radEmu:::B_to_B_tilde(B)
theta <- rbind(B_tilde,Matrix::Matrix(z,ncol = 1))

Now, compare X * Beta + z to X_tilde * theta

X %*% B+ z
D_tilde %*% theta

On my end these don't match. I think that the issue is with the ordering of the elements of X_tilde. Construct an alternative version of X_tilde:

X_tilde_alt <- cbind(rbind(cbind(diag(X[1, 1], 3), diag(X[1, 2], 3)),
                cbind(diag(X[2, 1], 3), diag(X[2, 2], 3)),
                cbind(diag(X[3, 1], 3), diag(X[3, 2], 3)),
                cbind(diag(X[4, 1], 3), diag(X[4, 2], 3))),
                c(rep(1, 3), rep(0, 9)), 
                c(rep(0, 3), rep(1, 3), rep(0, 6)),
                c(rep(0, 6), rep(1, 3), rep(0, 3)),
                c(rep(0, 9), rep(1, 3)))

Now again check. Here on my end they are now equal element-wise.

X %*% B + z
X_tilde_alt %*% theta

[svteichman]

Closing this issue because it is no longer a problem for the new version of radEmu