tests/testthat/helper.R
In miheerdew/bmdupdate: This package implements the pvalue calculation for the ChiSq test
m <- 400
nY <- 10
rho <- 0.5
ndata <- 300
Beta <- 0.2
s2 <- 1
mvrnormR <- function(n, mu, sigma) {
ncols <- ncol(sigma)
mu <- rep(mu, each = n) ## not obliged to use a matrix (recycling)
mu + matrix(rnorm(n * ncols), ncol = ncols) %*% chol(sigma)
}
generate_data <- function(m, nY, rho, ndata, beta, s2) {
Sig <- diag(1 - rho, m) + matrix(rep(rho, m^2), ncol = m)
X <- mvrnormR(ndata, rep(0, m), Sig)
Y <- Beta * rowSums(X) + matrix(rnorm(ndata * nY, sd = sqrt(s2)), ncol = nY)
list(X = X, Y = Y, dx=m, dy=nY, m=ndata)
}
precomputations <- function(X, Y) {
X <- scale(X)
Y <- scale(Y)
ndata <- nrow(X)
X2 <- X^2
X3 <- X^3
X4 <- X^4
Y2 <- Y^2
Y3 <- Y^3
Y4 <- Y^4
tX <- t(X)
tY <- t(Y)
tX2 <- t(X2)
tX3 <- t(X3)
XXt <- tcrossprod(X)
XXt2 <- XXt^2
Y4ColSums <- colSums(Y4)
X4ColSums <- colSums(X4)
X1RowSums <- rowSums(X)
X2RowSums <- rowSums(X2)
M <- ndata - 1
allr <- crossprod(Y, X) / M
allrSums <- rowSums(allr)
allr22 <- crossprod(Y2, X2) / M
allr31 <- crossprod(Y3, X) / M
return(environment())
}
traceXR <- function (i, env) {
y <- env$Y[ , i]
e <- matrix(1, nrow=env$ndata, ncol=1)
y2 <- env$Y2[ , i] # ndata x 1 ## GET THIS FROM Y2
y3 <- env$Y3[ , i] ## GET THIS FROM Y3
M <- env$ndata - 1
r <- env$allr[i, , drop = FALSE]# 1 x dx ## GET THIS FROM allr
r2 <- env$allr22[i, , drop = FALSE]# 1 x dx
r3 <- env$allr31[i, , drop = FALSE]# 1 x dx
rX <- tcrossprod(env$X, r) # m x 1
r3X <- tcrossprod(env$X, r3) # m x 1
rX3 <- tcrossprod(env$X3, r) # m x 1
rSq <- r^2
rSqX2 <- tcrossprod(env$X2,rSq) # m x 1
r2rSqX2 <- tcrossprod(env$X2, r2 * rSq) # m x 1
r3rX2 <- tcrossprod(env$X2, r3*r) # m x 1
tXr <- env$tX * as.vector(r) # dx x m
y4Sum <- env$Y4ColSums[i] # scalar ## GET THIS FROM Y4ColSums
rSqSum <- sum(rSq) #scalar
Z1 <- env$X2 %*% tXr # m x m
Z2 <- env$X2 %*% tXr^2 # m x m
#All of these are scalars
#TODO: Normalize them by (m-1)^2 after they are working correctly.
AA <- emulator::quad.form(env$XXt2, y2)
AB1 <- y4Sum * crossprod(y2, rX^2)
AB2 <- emulator::quad.3form(env$XXt, y, y2 * tcrossprod(env$X, r2*r))
AB3 <- AB2
AB4 <- emulator::quad.3form(Z1^2, e, y2)
AC1 <- crossprod(rX * r3X, y2) * M
AC2 <- emulator::quad.3form(Z1 * env$XXt, y, y2)
AC3 <- AC1
AC4 <- AC2
B1B1 <- (y4Sum*rSqSum)^2
B1B2 <- rSqSum * y4Sum * tcrossprod(r2,rSq) * M
B1B3 <- B1B2
B1B4 <- y4Sum * sum(rSqX2^2)
B1C1 <- y4Sum * rSqSum * tcrossprod(r,r3) * M
B1C2 <- y4Sum * crossprod(y, rX * rSqX2)
B1C3 <- B1C1
B1C4 <- B1C2
B2B2 <- tcrossprod(rSq, r2^2) * rSqSum * M^2
B2B3 <- tcrossprod(rSq, r2)^2 * M^2
B2B4 <- sum(rSqX2 * r2rSqX2) * M
B2C1 <- tcrossprod(rSq, r2) * tcrossprod(r,r3) * M^2
B2C2 <- crossprod(y, r2rSqX2 * rX) * M
B2C3 <- sum(r * r2 * r3) * rSqSum * M^2
B2C4 <- crossprod(y, tcrossprod(env$X, r * r2) * rSqX2) * M
B3B3 <- B2B2
B3B4 <- B2B4
B3C1 <- B2C3
B3C2 <- B2C4
B3C3 <- B2C1
B3C4 <- B2C2
B4B4 <- sum(Z2^2)
B4C1 <- sum(rSqX2 * r3rX2) * M
B4C2 <- emulator::quad.3form(Z2 * Z1, e, y)
B4C3 <- B4C1
B4C4 <- B4C2
C1C1 <- rSqSum * sum(r3^2) * M^2
C1C2 <- crossprod(y, rSqX2 * r3X) * M
C1C3 <- sum(r * r3)^2 * M^2
C1C4 <- crossprod(y, rX * r3rX2) * M
C2C2 <- emulator::quad.form(Z2 * env$XXt, y)
C2C3 <- C1C4
C2C4 <- emulator::quad.form(Z1*t(Z1), y)
C3C3 <- C1C1
C3C4 <- C1C2
C4C4 <- C2C2
resvec <- c(AA, c(AB1, AB2, AB3, AB4)/4, -1/2*c( AC1, AC2, AC3, AC4),
c(B1B1, B1B2, B1B3, B1B4)/16, -1/8*c(B1C1, B1C2, B1C3, B1C4),
c(B2B2, B2B3, B2B4)/16, -1/8*c(B2C1, B2C2, B2C3, B2C4),
c(B3B3, B3B4)/16, -1/8*c(B3C1, B3C2, B3C3, B3C4),
B4B4/16, -1/8*c(B4C1, B4C2, B4C3, B4C4),
c(C1C1, C1C2, C1C3, C1C4,
C2C2, C2C3, C2C4,
C3C3, C3C4,
C4C4)/4)
resmat <- matrix(0, 9, 9)
resmat[row(resmat) >= col(resmat)] <- resvec
trprod <- sum(diag(resmat)) + 2 * sum(resmat[row(resmat) > col(resmat)])
#Trace of straight up varmat
star1 <- crossprod(y2, env$X2RowSums) #O(m)
star2 <- y4Sum * rSqSum
star3 <- crossprod(y2, rSqX2) #O(m)
star4 <- tcrossprod(rSq, env$X4ColSums) #O(n)
dagger1 <- crossprod(y3, rX) #O(m)
dagger2 <- crossprod(y, rX3) #O(m)
trmat <- star1 + (star2 + 2 * star3 + star4) / 4 - dagger1 - dagger2
return(c(trmat/M, trprod/M^2))
}
traceMR <- function(i, env){
r <- env$allr[i, ]
r2 <- env$allr22[i, ]
r3 <- env$allr31[i, ]
yX <- env$Y[ , i] * env$X
tX2r <- env$tX2 * r
# star 1
S1 <- crossprod(yX)
# star 2
S2 <- env$Y4ColSums[i] * tcrossprod(r)
# star 3
S3 <- tcrossprod(r, r * r2) * env$M
# star 4
S4 <- tcrossprod(tX2r)
# dagger 1
D1 <- tcrossprod(r, r3) * env$M
# dagger 2
D2 <- tX2r %*% yX
bigmat <- (S1 + (S2 + S3 + t(S3) + S4) / 4 - (D1 + t(D1) + D2 + t(D2)) / 2) / env$M
trSig <- sum(diag(bigmat))
trSig2 <- sum(bigmat^2)
return(c(trSig, trSig2))
}
traceXR_fast <- function (i, env) {
y <- env$Y[ , i]
e <- matrix(1, nrow=env$ndata, ncol=1)
y2 <- env$Y2[ , i] # ndata x 1 ## GET THIS FROM Y2
y3 <- env$Y3[ , i] ## GET THIS FROM Y3
M <- env$ndata - 1
r <- env$allr[i, , drop = FALSE]# 1 x dx ## GET THIS FROM allr
r2 <- env$allr22[i, , drop = FALSE]# 1 x dx
r3 <- env$allr31[i, , drop = FALSE]# 1 x dx
rX <- tcrossprod(env$X, r) # m x 1
r3X <- tcrossprod(env$X, r3) # m x 1
rX3 <- tcrossprod(env$X3, r) # m x 1
rSq <- r^2
rSqX2 <- tcrossprod(env$X2,rSq) # m x 1
r2rSqX2 <- tcrossprod(env$X2, r2 * rSq) # m x 1
r3rX2 <- tcrossprod(env$X2, r3*r) # m x 1
tXr <- env$tX * as.vector(r) # dx x m
y4Sum <- env$Y4ColSums[i] # scalar ## GET THIS FROM Y4ColSums
rSqSum <- sum(rSq) #scalar
Z1 <- env$X2 %*% tXr # m x m
Z2 <- env$X2 %*% tXr^2 # m x m
Z1 <- diag(diag(Z1))
Z2[row(Z1) != col(Z1)] <- rSqSum
#All of these are scalars
#TODO: Normalize them by (m-1)^2 after they are working correctly.
AA <- emulator::quad.form(env$XXt2, y2)
AB1 <- y4Sum * crossprod(y2, rX^2)
AB2 <- emulator::quad.3form(env$XXt, y, y2 * tcrossprod(env$X, r2*r))
AB3 <- AB2
AB4 <- emulator::quad.3form(Z1^2, e, y2)
AC1 <- crossprod(rX * r3X, y2) * M
AC2 <- emulator::quad.3form(Z1 * env$XXt, y, y2)
AC3 <- AC1
AC4 <- AC2
B1B1 <- (y4Sum*rSqSum)^2
B1B2 <- rSqSum * y4Sum * tcrossprod(r2,rSq) * M
B1B3 <- B1B2
B1B4 <- y4Sum * sum(rSqX2^2)
B1C1 <- y4Sum * rSqSum * tcrossprod(r,r3) * M
B1C2 <- y4Sum * crossprod(y, rX * rSqX2)
B1C3 <- B1C1
B1C4 <- B1C2
B2B2 <- tcrossprod(rSq, r2^2) * rSqSum * M^2
B2B3 <- tcrossprod(rSq, r2)^2 * M^2
B2B4 <- sum(rSqX2 * r2rSqX2) * M
B2C1 <- tcrossprod(rSq, r2) * tcrossprod(r,r3) * M^2
B2C2 <- crossprod(y, r2rSqX2 * rX) * M
B2C3 <- sum(r * r2 * r3) * rSqSum * M^2
B2C4 <- crossprod(y, tcrossprod(env$X, r * r2) * rSqX2) * M
B3B3 <- B2B2
B3B4 <- B2B4
B3C1 <- B2C3
B3C2 <- B2C4
B3C3 <- B2C1
B3C4 <- B2C2
B4B4 <- sum(Z2^2)
B4C1 <- sum(rSqX2 * r3rX2) * M
B4C2 <- emulator::quad.3form(Z2 * Z1, e, y)
B4C3 <- B4C1
B4C4 <- B4C2
C1C1 <- rSqSum * sum(r3^2) * M^2
C1C2 <- crossprod(y, rSqX2 * r3X) * M
C1C3 <- sum(r * r3)^2 * M^2
C1C4 <- crossprod(y, rX * r3rX2) * M
C2C2 <- emulator::quad.form(Z2 * env$XXt, y)
C2C3 <- C1C4
C2C4 <- emulator::quad.form(Z1*t(Z1), y)
C3C3 <- C1C1
C3C4 <- C1C2
C4C4 <- C2C2
resvec <- c(AA, c(AB1, AB2, AB3, AB4)/4, -1/2*c( AC1, AC2, AC3, AC4),
c(B1B1, B1B2, B1B3, B1B4)/16, -1/8*c(B1C1, B1C2, B1C3, B1C4),
c(B2B2, B2B3, B2B4)/16, -1/8*c(B2C1, B2C2, B2C3, B2C4),
c(B3B3, B3B4)/16, -1/8*c(B3C1, B3C2, B3C3, B3C4),
B4B4/16, -1/8*c(B4C1, B4C2, B4C3, B4C4),
c(C1C1, C1C2, C1C3, C1C4,
C2C2, C2C3, C2C4,
C3C3, C3C4,
C4C4)/4)
resmat <- matrix(0, 9, 9)
resmat[row(resmat) >= col(resmat)] <- resvec
trprod <- sum(diag(resmat)) + 2 * sum(resmat[row(resmat) > col(resmat)])
#Trace of straight up varmat
star1 <- crossprod(y2, env$X2RowSums) #O(m)
star2 <- y4Sum * rSqSum
star3 <- crossprod(y2, rSqX2) #O(m)
star4 <- tcrossprod(rSq, env$X4ColSums) #O(n)
dagger1 <- crossprod(y3, rX) #O(m)
dagger2 <- crossprod(y, rX3) #O(m)
trmat <- star1 + (star2 + 2 * star3 + star4) / 4 - dagger1 - dagger2
return(c(trmat/M, trprod/M^2))
}
miheerdew/bmdupdate documentation built on May 17, 2019, 1:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
m <- 400
nY <- 10
rho <- 0.5
ndata <- 300
Beta <- 0.2
s2 <- 1
mvrnormR <- function(n, mu, sigma) {
ncols <- ncol(sigma)
mu <- rep(mu, each = n) ## not obliged to use a matrix (recycling)
mu + matrix(rnorm(n * ncols), ncol = ncols) %*% chol(sigma)
}
generate_data <- function(m, nY, rho, ndata, beta, s2) {
Sig <- diag(1 - rho, m) + matrix(rep(rho, m^2), ncol = m)
X <- mvrnormR(ndata, rep(0, m), Sig)
Y <- Beta * rowSums(X) + matrix(rnorm(ndata * nY, sd = sqrt(s2)), ncol = nY)
list(X = X, Y = Y, dx=m, dy=nY, m=ndata)
}
precomputations <- function(X, Y) {
X <- scale(X)
Y <- scale(Y)
ndata <- nrow(X)
X2 <- X^2
X3 <- X^3
X4 <- X^4
Y2 <- Y^2
Y3 <- Y^3
Y4 <- Y^4
tX <- t(X)
tY <- t(Y)
tX2 <- t(X2)
tX3 <- t(X3)
XXt <- tcrossprod(X)
XXt2 <- XXt^2
Y4ColSums <- colSums(Y4)
X4ColSums <- colSums(X4)
X1RowSums <- rowSums(X)
X2RowSums <- rowSums(X2)
M <- ndata - 1
allr <- crossprod(Y, X) / M
allrSums <- rowSums(allr)
allr22 <- crossprod(Y2, X2) / M
allr31 <- crossprod(Y3, X) / M
return(environment())
}
traceXR <- function (i, env) {
y <- env$Y[ , i]
e <- matrix(1, nrow=env$ndata, ncol=1)
y2 <- env$Y2[ , i] # ndata x 1 ## GET THIS FROM Y2
y3 <- env$Y3[ , i] ## GET THIS FROM Y3
M <- env$ndata - 1
r <- env$allr[i, , drop = FALSE]# 1 x dx ## GET THIS FROM allr
r2 <- env$allr22[i, , drop = FALSE]# 1 x dx
r3 <- env$allr31[i, , drop = FALSE]# 1 x dx
rX <- tcrossprod(env$X, r) # m x 1
r3X <- tcrossprod(env$X, r3) # m x 1
rX3 <- tcrossprod(env$X3, r) # m x 1
rSq <- r^2
rSqX2 <- tcrossprod(env$X2,rSq) # m x 1
r2rSqX2 <- tcrossprod(env$X2, r2 * rSq) # m x 1
r3rX2 <- tcrossprod(env$X2, r3*r) # m x 1
tXr <- env$tX * as.vector(r) # dx x m
y4Sum <- env$Y4ColSums[i] # scalar ## GET THIS FROM Y4ColSums
rSqSum <- sum(rSq) #scalar
Z1 <- env$X2 %*% tXr # m x m
Z2 <- env$X2 %*% tXr^2 # m x m
#All of these are scalars
#TODO: Normalize them by (m-1)^2 after they are working correctly.
AA <- emulator::quad.form(env$XXt2, y2)
AB1 <- y4Sum * crossprod(y2, rX^2)
AB2 <- emulator::quad.3form(env$XXt, y, y2 * tcrossprod(env$X, r2*r))
AB3 <- AB2
AB4 <- emulator::quad.3form(Z1^2, e, y2)
AC1 <- crossprod(rX * r3X, y2) * M
AC2 <- emulator::quad.3form(Z1 * env$XXt, y, y2)
AC3 <- AC1
AC4 <- AC2
B1B1 <- (y4Sum*rSqSum)^2
B1B2 <- rSqSum * y4Sum * tcrossprod(r2,rSq) * M
B1B3 <- B1B2
B1B4 <- y4Sum * sum(rSqX2^2)
B1C1 <- y4Sum * rSqSum * tcrossprod(r,r3) * M
B1C2 <- y4Sum * crossprod(y, rX * rSqX2)
B1C3 <- B1C1
B1C4 <- B1C2
B2B2 <- tcrossprod(rSq, r2^2) * rSqSum * M^2
B2B3 <- tcrossprod(rSq, r2)^2 * M^2
B2B4 <- sum(rSqX2 * r2rSqX2) * M
B2C1 <- tcrossprod(rSq, r2) * tcrossprod(r,r3) * M^2
B2C2 <- crossprod(y, r2rSqX2 * rX) * M
B2C3 <- sum(r * r2 * r3) * rSqSum * M^2
B2C4 <- crossprod(y, tcrossprod(env$X, r * r2) * rSqX2) * M
B3B3 <- B2B2
B3B4 <- B2B4
B3C1 <- B2C3
B3C2 <- B2C4
B3C3 <- B2C1
B3C4 <- B2C2
B4B4 <- sum(Z2^2)
B4C1 <- sum(rSqX2 * r3rX2) * M
B4C2 <- emulator::quad.3form(Z2 * Z1, e, y)
B4C3 <- B4C1
B4C4 <- B4C2
C1C1 <- rSqSum * sum(r3^2) * M^2
C1C2 <- crossprod(y, rSqX2 * r3X) * M
C1C3 <- sum(r * r3)^2 * M^2
C1C4 <- crossprod(y, rX * r3rX2) * M
C2C2 <- emulator::quad.form(Z2 * env$XXt, y)
C2C3 <- C1C4
C2C4 <- emulator::quad.form(Z1*t(Z1), y)
C3C3 <- C1C1
C3C4 <- C1C2
C4C4 <- C2C2
resvec <- c(AA, c(AB1, AB2, AB3, AB4)/4, -1/2*c( AC1, AC2, AC3, AC4),
c(B1B1, B1B2, B1B3, B1B4)/16, -1/8*c(B1C1, B1C2, B1C3, B1C4),
c(B2B2, B2B3, B2B4)/16, -1/8*c(B2C1, B2C2, B2C3, B2C4),
c(B3B3, B3B4)/16, -1/8*c(B3C1, B3C2, B3C3, B3C4),
B4B4/16, -1/8*c(B4C1, B4C2, B4C3, B4C4),
c(C1C1, C1C2, C1C3, C1C4,
C2C2, C2C3, C2C4,
C3C3, C3C4,
C4C4)/4)
resmat <- matrix(0, 9, 9)
resmat[row(resmat) >= col(resmat)] <- resvec
trprod <- sum(diag(resmat)) + 2 * sum(resmat[row(resmat) > col(resmat)])
#Trace of straight up varmat
star1 <- crossprod(y2, env$X2RowSums) #O(m)
star2 <- y4Sum * rSqSum
star3 <- crossprod(y2, rSqX2) #O(m)
star4 <- tcrossprod(rSq, env$X4ColSums) #O(n)
dagger1 <- crossprod(y3, rX) #O(m)
dagger2 <- crossprod(y, rX3) #O(m)
trmat <- star1 + (star2 + 2 * star3 + star4) / 4 - dagger1 - dagger2
return(c(trmat/M, trprod/M^2))
}
traceMR <- function(i, env){
r <- env$allr[i, ]
r2 <- env$allr22[i, ]
r3 <- env$allr31[i, ]
yX <- env$Y[ , i] * env$X
tX2r <- env$tX2 * r
# star 1
S1 <- crossprod(yX)
# star 2
S2 <- env$Y4ColSums[i] * tcrossprod(r)
# star 3
S3 <- tcrossprod(r, r * r2) * env$M
# star 4
S4 <- tcrossprod(tX2r)
# dagger 1
D1 <- tcrossprod(r, r3) * env$M
# dagger 2
D2 <- tX2r %*% yX
bigmat <- (S1 + (S2 + S3 + t(S3) + S4) / 4 - (D1 + t(D1) + D2 + t(D2)) / 2) / env$M
trSig <- sum(diag(bigmat))
trSig2 <- sum(bigmat^2)
return(c(trSig, trSig2))
}
traceXR_fast <- function (i, env) {
y <- env$Y[ , i]
e <- matrix(1, nrow=env$ndata, ncol=1)
y2 <- env$Y2[ , i] # ndata x 1 ## GET THIS FROM Y2
y3 <- env$Y3[ , i] ## GET THIS FROM Y3
M <- env$ndata - 1
r <- env$allr[i, , drop = FALSE]# 1 x dx ## GET THIS FROM allr
r2 <- env$allr22[i, , drop = FALSE]# 1 x dx
r3 <- env$allr31[i, , drop = FALSE]# 1 x dx
rX <- tcrossprod(env$X, r) # m x 1
r3X <- tcrossprod(env$X, r3) # m x 1
rX3 <- tcrossprod(env$X3, r) # m x 1
rSq <- r^2
rSqX2 <- tcrossprod(env$X2,rSq) # m x 1
r2rSqX2 <- tcrossprod(env$X2, r2 * rSq) # m x 1
r3rX2 <- tcrossprod(env$X2, r3*r) # m x 1
tXr <- env$tX * as.vector(r) # dx x m
y4Sum <- env$Y4ColSums[i] # scalar ## GET THIS FROM Y4ColSums
rSqSum <- sum(rSq) #scalar
Z1 <- env$X2 %*% tXr # m x m
Z2 <- env$X2 %*% tXr^2 # m x m
Z1 <- diag(diag(Z1))
Z2[row(Z1) != col(Z1)] <- rSqSum
#All of these are scalars
#TODO: Normalize them by (m-1)^2 after they are working correctly.
AA <- emulator::quad.form(env$XXt2, y2)
AB1 <- y4Sum * crossprod(y2, rX^2)
AB2 <- emulator::quad.3form(env$XXt, y, y2 * tcrossprod(env$X, r2*r))
AB3 <- AB2
AB4 <- emulator::quad.3form(Z1^2, e, y2)
AC1 <- crossprod(rX * r3X, y2) * M
AC2 <- emulator::quad.3form(Z1 * env$XXt, y, y2)
AC3 <- AC1
AC4 <- AC2
B1B1 <- (y4Sum*rSqSum)^2
B1B2 <- rSqSum * y4Sum * tcrossprod(r2,rSq) * M
B1B3 <- B1B2
B1B4 <- y4Sum * sum(rSqX2^2)
B1C1 <- y4Sum * rSqSum * tcrossprod(r,r3) * M
B1C2 <- y4Sum * crossprod(y, rX * rSqX2)
B1C3 <- B1C1
B1C4 <- B1C2
B2B2 <- tcrossprod(rSq, r2^2) * rSqSum * M^2
B2B3 <- tcrossprod(rSq, r2)^2 * M^2
B2B4 <- sum(rSqX2 * r2rSqX2) * M
B2C1 <- tcrossprod(rSq, r2) * tcrossprod(r,r3) * M^2
B2C2 <- crossprod(y, r2rSqX2 * rX) * M
B2C3 <- sum(r * r2 * r3) * rSqSum * M^2
B2C4 <- crossprod(y, tcrossprod(env$X, r * r2) * rSqX2) * M
B3B3 <- B2B2
B3B4 <- B2B4
B3C1 <- B2C3
B3C2 <- B2C4
B3C3 <- B2C1
B3C4 <- B2C2
B4B4 <- sum(Z2^2)
B4C1 <- sum(rSqX2 * r3rX2) * M
B4C2 <- emulator::quad.3form(Z2 * Z1, e, y)
B4C3 <- B4C1
B4C4 <- B4C2
C1C1 <- rSqSum * sum(r3^2) * M^2
C1C2 <- crossprod(y, rSqX2 * r3X) * M
C1C3 <- sum(r * r3)^2 * M^2
C1C4 <- crossprod(y, rX * r3rX2) * M
C2C2 <- emulator::quad.form(Z2 * env$XXt, y)
C2C3 <- C1C4
C2C4 <- emulator::quad.form(Z1*t(Z1), y)
C3C3 <- C1C1
C3C4 <- C1C2
C4C4 <- C2C2
resvec <- c(AA, c(AB1, AB2, AB3, AB4)/4, -1/2*c( AC1, AC2, AC3, AC4),
c(B1B1, B1B2, B1B3, B1B4)/16, -1/8*c(B1C1, B1C2, B1C3, B1C4),
c(B2B2, B2B3, B2B4)/16, -1/8*c(B2C1, B2C2, B2C3, B2C4),
c(B3B3, B3B4)/16, -1/8*c(B3C1, B3C2, B3C3, B3C4),
B4B4/16, -1/8*c(B4C1, B4C2, B4C3, B4C4),
c(C1C1, C1C2, C1C3, C1C4,
C2C2, C2C3, C2C4,
C3C3, C3C4,
C4C4)/4)
resmat <- matrix(0, 9, 9)
resmat[row(resmat) >= col(resmat)] <- resvec
trprod <- sum(diag(resmat)) + 2 * sum(resmat[row(resmat) > col(resmat)])
#Trace of straight up varmat
star1 <- crossprod(y2, env$X2RowSums) #O(m)
star2 <- y4Sum * rSqSum
star3 <- crossprod(y2, rSqX2) #O(m)
star4 <- tcrossprod(rSq, env$X4ColSums) #O(n)
dagger1 <- crossprod(y3, rX) #O(m)
dagger2 <- crossprod(y, rX3) #O(m)
trmat <- star1 + (star2 + 2 * star3 + star4) / 4 - dagger1 - dagger2
return(c(trmat/M, trprod/M^2))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.