R/selection.R In envlpaster: Enveloping the Aster Model

Documented in selection

```selection <- function(parm, index, model, data = NULL,
alpha = 0.05, type = c("canonical","mean-value"),
method = c("eigen","1d")){

# important initial quantities
modmat <- model\$modmat
offset <- as.vector(model\$origin)
nind <- dim(modmat)[1]
nnode <- dim(modmat)[2]
m <- length(index)
p <- length(parm)
modmat.mat <- matrix(modmat, nrow = nind * nnode)
x <- model\$x
root <- model\$root
pred <- model\$pred
fam <- model\$fam

# get the spectral structure of FI or inv FI
avar <- model\$fisher
if(type == "canonical") avar <- solve(avar, tol = 1e-20)
eig.avar.targ <- eigen(avar[index,index])

# get important quantities for 1D algorithm
avar.targ <- avar[index, index]
U <- parm %*% t(parm)
U.targ <- U[index, index]

# initialize internal quantities
beta.cand.test <- NULL
beta <- model\$coef
beta.foo <- rep(0, p)
tau.int <- parm
aic.env <- bic.env <- pval <- 0

# consider every possible envelope estimator constructed
# using eigenspaces
if(method == "eigen"){
beta.cand.test <- matrix(unlist(lapply(1:m, FUN = function(j) {
cand <- combs(1:m, j)
matrix(unlist(lapply(1:nrow(cand), FUN = function(k) {

# obtain the eigenspace of interest and the projection
# into that eigenspace
G <- eig.avar.targ\$vec[,cand[k,]]
P <-  projection(G)

# build the model matrix Menv
M <- t(modmat.mat)
#M2 <- M[(1:p %in% index),]; M2 <- P %*% M2
M[(1:p %in% index),] <- P %*% M[(1:p %in% index),]
modmat.mat.env <- t(M)
modmat.int <- array(modmat.mat.env, c(nind, nnode, p))

# obtain the corresponing beta vector
if(type == "mean-value"){
tau.env.int <- crossprod(P, parm[index])
tau.int[index] <- tau.env.int
ind.int <- max(which(!1:p %in% index)) + length(cand[k,])
beta.foo[1:(ind.int)] <- suppressWarnings(try(
aster(x, root, pred, fam, modmat = modmat.int,
origin = model\$origin, maxiter = 10000)\$coef,
silent = TRUE))
#beta.foo <- try(transformUnconditional(parm = tau.int,
#  modmat.mat.env, data, from = "tau", to = "beta",
#  offset = offset), silent = TRUE)
#print(beta.foo)
}

if(type == "canonical"){
beta.foo[index] <- P %*% beta.foo[index]
}

bic.env <- aic.env <- 1e9
pval <- 0
if(class(beta.foo) != "try-error"){
# obtain selection criteria values
env <-mlogl(beta.foo, pred, fam, x, root,
modmat = modmat.int, type = "unconditional")\$value
full <-mlogl(beta, pred, fam, x, root,
modmat = modmat, type = "unconditional")\$value
df <- j*(m-j) + j + j*(j+1)/2 + (m-j)*(m-j+1)/2
df.full <- m + m*(m+1)/2
LRT <- 2*(env - full)
if(j < m) pval <- pchisq(LRT, df = df.full - df, lower.tail = F)
if(j == m) pval <- 1
bic.env <-2*env + df*log(nind)
aic.env <- 2*df + 2*env
}

things <- c(j, aic.env, bic.env, pval)
return(things)
} )))
} )), ncol = 4, byrow = TRUE)
}

# consider every possible envelope estimator using
# the 1D algorithm
if(method == "1d"){
beta.cand.test <- matrix(unlist(lapply(1:m, FUN = function(j) {
if(j < m){

# obtain the eigenspace of interest and the projection
# into that eigenspace
G <- manifold1Dplus(M = avar.targ, U = U.targ, u = j)
P <- projection(G)

# build the model matrix Menv
M <- t(modmat.mat)
M2 <- M[(1:p %in% index),]; M2 <- P %*% M2
M[(1:p %in% index),] <- M2
modmat.mat.env <- t(M)
modmat.int <- array(modmat.mat.env, c(nind, nnode, p))

# obtain the corresponing beta vector
if(type == "mean-value"){
tau.env.int <- crossprod(P, parm[index])
tau.int[index] <- tau.env.int
beta.foo <- transformUnconditional(parm = tau.int,
modmat.mat.env, data, from = "tau", to = "beta",
offset = offset, tolerance = 1e-200)
}

if(type == "canonical"){
beta.foo[index] <- P %*% beta.foo[index]
}

# obtain selection criteria values
env <-mlogl(beta.foo, pred, fam, x, root,
modmat = modmat.int, type = "unconditional")\$value
full <-mlogl(beta, pred, fam, x, root,
modmat = modmat, type = "unconditional")\$value
k <- j*(m-j) + j + j*(j+1)/2 + (m-j)*(m-j+1)/2
k.full <- m + m*(m+1)/2
LRT <- 2*(env - full)
pval <- pchisq(LRT, df = k.full - k, lower.tail = F)
bic.env <-2*env + k*log(nind)
aic.env <- 2*k + 2*env
}

if(j == m){
full <- mlogl(beta, pred, fam, x, root,
modmat = modmat, type = "unconditional")\$value
k <- (m + m*(m+1)/2)
bic.env <-2*full + k*log(nind)
aic.env <- 2*k + 2*full
pval <- 1
}

things <- c(j, aic.env, bic.env, pval)
return(things)
})), ncol = 4, byrow = TRUE)
}

# select the dimension for each of the criteria
aic <- bic <- LRT <- m
if(method == "1d"){
aic <- which(beta.cand.test[,2] == min(beta.cand.test[,2]))
bic <- which(beta.cand.test[,3] == min(beta.cand.test[,3]))
LRT <- min(which(beta.cand.test[,4] > alpha))
}

# select the eigenspace for each of the criteria
if(method == "eigen"){
dims <- nrow(combs(1:m,1))
for(u in 2:m) dims[u] <- nrow(combs(1:m, u)) + dims[u-1]
aic <- which(beta.cand.test[,2] == min(beta.cand.test[,2]))
bic <- which(beta.cand.test[,3] == min(beta.cand.test[,3]))
LRT <- min(which(beta.cand.test[,4] > alpha))
u.aic <- min(which(dims >= aic))
u.bic <- min(which(dims >= bic))
u.LRT <- min(which(dims >= LRT))
if(u.aic == 1) aic <- combs(1:m, u.aic)[c(aic), ]
if(u.aic >= 2) aic <- combs(1:m, u.aic)[c(aic - dims[u.aic - 1]), ]
if(u.bic == 1) bic <- combs(1:m, u.bic)[c(bic), ]
if(u.bic >= 2) bic <- combs(1:m, u.bic)[c(bic - dims[u.bic - 1]), ]
if(u.LRT == 1) LRT <- combs(1:m, u.LRT)[c(LRT), ]
if(u.LRT >= 2) LRT <- combs(1:m, u.LRT)[c(LRT - dims[u.LRT - 1]), ]
}

beta.cand.test <- as.data.frame(beta.cand.test)
colnames(beta.cand.test) <- c("u","aic","bic","p-val")

out <- list(aic = aic, bic = bic, LRT = LRT, out = beta.cand.test)
return(out)
}
```

Try the envlpaster package in your browser

Any scripts or data that you put into this service are public.

envlpaster documentation built on May 2, 2019, 2:10 a.m.