# smrt.stepdown
#'
#' Perform stepdown testing on parameter estimates for each predictor
#'
#' @param beta Vector of estimates
#' @param betastar Matrix of perturbed estimates
#' @param inv.info Full Fisher information matrix
#' @param inv.info.beta Sub-matrix of Fisher information corresponding to \code{beta}.
#' @param z.th Threshold for rejection, either scalar or p-length vector.
#' @param get.u Function to obtain the influence function of beta based on \code{psd}.
#' @details
#'
#' Estimates should be arranged such that all estimates for first outcome
#' come first, then all estimates for second outcomes, and so on. Only one of
#' inv.info and inv.info.beta need to be supplied.
#'
#'@return Indices of beta corresponding to rejected tests.
#'
#'@export
smrt.stepdown <- function(beta, betastar, inv.info = NULL, inv.info.beta = NULL, z.th, p, m) {
if (is.null(inv.info) & is.null(inv.info.beta)) {
stop('please provide inverse information.')
} else if (!is.null(inv.info.beta)) {
inv.info <- inv.info.beta
} else {
inv.info <- lapply(inv.info, function(x) {
p.prime <- dim(x)[1]
x[(p.prime - p + 1):p.prime, (p.prime - p + 1):p.prime]
})
}
if (length(z.th) == 1) {
z.th <- rep(z.th, p)
} else if (length(z.th) != p) {
stop('number of thresholds specified must be equal to p')
}
rej.set <- c()
sigma <- sqrt(unlist(lapply(inv.info, diag)))
z <- abs(beta/sigma)
z.star <- abs(t((t(betastar) - beta)/sigma))
for (i in 1:p) {
this.snp <- matrix(z, ncol = m)[i,]
this.snpstar <- t(apply(z.star, 1, function(x) matrix(x, ncol = m)[i,]))
this.snp.rej <- stepdown(this.snp, this.snpstar, z.th[i])
if (!is.null(this.snp.rej)) {
rej.set <- c(rej.set, get.x.null(i, p, m)[this.snp.rej])
}
}
rej.set
}
#' Basic stepdown function
#'
#' @param z Index vector of estimates
#' @param zstar Matrix of perturbed estimates
#' @param thr Threshold for rejection
#' @return Indices corresponding to rejected tests.
#'
#'@export
stepdown <- function(z, zstar, thr) {
rej.set <- c()
names(z) <- 1:length(z)
for (i in 1:length(z)) {
max.ind <- which.max(z)
if (length(z) > 1) {
max.z <- max(z)
max.zstar <- apply(zstar, 1, max)
} else {
max.z <- z
max.zstar <- zstar
}
if (max.z <= quantile(max.zstar, thr, na.rm = TRUE)) {
return(rej.set)
} else {
rej.set <- c(rej.set, as.numeric(names(z)[z == max.z]))
if (length(z) > 1) {
z <- z[-max.ind]
zstar <- zstar[,-max.ind]
}
}
}
rej.set
}
get.x.null <- function(inds, p, m) {
if (all(inds == 0)) {
NULL
} else if (any(inds == 0)) {
stop('if 0 is an index, it must be the only index.')
} else {
sort(c(sapply(inds, function(x) x + 0:(m-1)*p)))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.