Nothing
#' A wrapper for the different qris fit
#'
#' This version is more flexible for future extension
#'
#' @noRd
qris.fit <- function(info, method) {
switch(method,
smooth = qris.smooth(info),
iterative = qris.iter(info),
nonsmooth = qris.nonsmooth(info))
}
qris.nonsmooth <- function(info) {
## 1. : L1-minimization : Estimating equation for estimating beta (using rq)
out <- with(info, {
M <- 1e6
pseudo1 <- -colSums(X*I*W)
pseudo2 <- 2 * colSums(X*I*Q)
Y.reg <- c(data[,2], M, M)
X.reg <- rbind(X, rbind(pseudo1, pseudo2))
W.reg <- c(I * W, 1, 1)
Li.fit <- rq.wfit(X.reg, Y.reg, weights = W.reg)
coefficient <- as.vector(Li.fit$coefficients)
uTime <- sort(unique(data$Z))
if (all(Li.fit$coefficients <= 10)){
if (nB <= 1) {
se <- rep(NA_real_, nc)
vcov <- matrix(NA_real_, nc, nc)
return(list(coefficient = coefficient, stderr = se, vcov = vcov))
}
if (se == "fmb"){
fmb.result <- matrix(NA_real_, nc, nB)
for (j in 1:nB){
## generating perturbation variable
eta <- rexp(n)
if (all(data$delta == 1)) {
W_star <- rep(1, n)
} else {
survp <- drop(ghatC(data$Z, 1 - data$delta, eta))
ghatstart0 <- ifelse(t0 > 0, survp[findInterval(t0, uTime)], 1)
W_star <- data$delta / survp[findInterval(data$Z, uTime)] * ghatstart0
W_star[is.na(W_star)] <- max(W_star, na.rm = TRUE)
}
rev_pseudo1 <- -colSums(X * W_star * eta)
rev_pseudo2 <- 2 * colSums(X * I * eta * Q)
rev_Y.reg <- c(data[,2], M, M)
rev_X.reg <- rbind(X, rbind(rev_pseudo1, rev_pseudo2))
rev_W.reg <- c(W_star * I * eta, 1, 1)
fmb.fit <- rq.wfit(rev_X.reg, rev_Y.reg, weights = rev_W.reg)
if (all(fmb.fit$coef <= 10)){
fmb.result[, j] <- as.vector(fmb.fit$coef)
}
}
fmb.sigma <- try(cov(t(fmb.result), use = "complete.obs"), silent = T)
if(class(fmb.sigma)[1] == "try-error" | sum(!is.na(colSums(fmb.result))) <= 1){
stop("More resampling iterations are necessary")
} else {
fmb.se <- sqrt(diag(fmb.sigma))
out <- list(coefficient = coefficient, stderr = fmb.se, vcov = fmb.sigma)
}
} else {
stop("Only full multiplier bootstrapping (fmb) is available for nonsmooth estimating equation approach")
}
} else {
coefficient <- se <- rep(NA, nc)
vcov <- matrix(NA, nc, nc)
out <- list(coefficient = coefficient, stderr = se, vcov = vcov)
}})
out
}
qris.iter <- function(info) {
out <- with(info, {
new_h <- old_h <- H
new_sigma <- old_sigma <- H*n
iter_beta_result <- iter_SE_result <- matrix(NA_real_, control$maxiter, nc)
iter_beta_result[1,] <- old_beta <- new_beta <- betastart
iter_SE_result[1,] <- sqrt(diag(new_sigma))
iter_norm_result <- rep(NA_real_, control$maxiter)
uTime <- sort(unique(data$Z))
if (se == "fmb") {
for (k in 1:control$maxiter){
old_beta <- new_beta
old_sigma <- new_sigma
old_h <- new_h
slope_a <- Amat(old_beta, X, W, old_h, I, logZ, Q) / n
## Step 1 : Update beta()
## Singular matrix 'a' error message and break
if (qr(slope_a)$rank < min(dim(slope_a))) {
warning("'A' matrix is singular during iteration. Please try the non-iterative method.")
break
} else {
new_beta <- old_beta + qr.solve(slope_a) %*% (isObj(old_beta, X, W, old_h, I, logZ, Q)/n)
iter_beta_result[k,] <- t(new_beta)
## Step 2 : Update Sigma()
result.fmb <- matrix(NA_real_, nc, nB)
for (j in 1:nB){
## generating perturbation variable
eta <- rexp(n)
if (all(data$delta == 1)) {
W_star <- rep(1, n)
} else {
survp <- drop(ghatC(data$Z, 1 - data$delta, eta))
ghatstart0 <- ifelse(t0 > 0, survp[findInterval(t0, uTime)], 1)
W_star <- data$delta / survp[findInterval(data$Z, uTime)] * ghatstart0
W_star[is.na(W_star)] <- max(W_star, na.rm = TRUE)
}
fmb.fit <- nleqslv(old_beta, function(b)
rev_isObj(b, X, W_star, old_h, eta, I, logZ, Q)/n)
if (fmb.fit$termcd == 1 | fmb.fit$termcd == 2) {
result.fmb[, j] <- fmb.fit$x
}
}
new_sigma <- try(cov(t(result.fmb), use = "complete.obs"), silent = T)
## Trace the result
if (verbose) {
cat("\n Step:", k)
cat("\n beta:", as.numeric(new_beta))
cat("\n se:", as.numeric(sqrt(diag(new_sigma))), "\n")
}
if (qr(new_sigma)$rank < min(dim(new_sigma))) {
warning("Futher fmb method is inapplicable to this dataset. Please try other estimation methods")
new_sigma <- old_sigma
break
} else {
new_h <- new_sigma
iter_SE_result[k,] <- sqrt(diag(new_sigma))
iter_norm_result[k] <- norm(new_beta-old_beta, "F")
if(iter_norm_result[k]>=100*iter_norm_result[1]) {
warning("Point estimation failed to converge.")
break
}
if(norm(new_beta-old_beta, "i") < control$tol) break
}
}
} ## end for loop
} else {
for (k in 1:control$maxiter){
old_beta <- new_beta
old_sigma <- new_sigma
old_h <- new_h
slope_a <- Amat(old_beta, X, W, old_h, I, logZ, Q)/n
## Step 1 : Update beta()
## Singular matrix 'a' error message and break
if (qr(slope_a)$rank < min(dim(slope_a))) {
warning("'A' matrix is singular during iteration. Please try the non-iterative method.")
break
} else {
new_beta <- old_beta + qr.solve(slope_a) %*% (isObj(old_beta, X, W, old_h, I, logZ, Q)/n)
iter_beta_result[k,] <- t(new_beta)
## Step 2 : Update Sigma()
result.pmb <- matrix(NA, nc, nB)
m2 <- isObjL(old_beta, X, H, logZ)
for (j in 1:nB){
## generating perturbation variable
eta <- rexp(n)
if (all(data$delta == 1)) {
W_star <- rep(1, n)
} else {
survp <- drop(ghatC(data$Z, 1 - data$delta, eta))
ghatstart0 <- ifelse(t0 > 0, survp[findInterval(t0, uTime)], 1)
W_star <- data$delta / survp[findInterval(data$Z, uTime)] * ghatstart0
W_star[is.na(W_star)] <- max(W_star, na.rm = TRUE)
}
result.pmb[,j] <- t(X * I * eta) %*% (m2 * W_star - Q) / n
## rev_isObj(old_beta, X, W_star, old_h, eta, I, logZ, Q)/n
}
new_V <- cov(t(result.pmb), use = "complete.obs")
new_sigma <- t(qr.solve(slope_a)) %*% new_V %*% qr.solve(slope_a)
new_h <- new_sigma
## Trace the result
if (verbose) {
cat("\n beta:", as.numeric(new_beta), "\n")
cat("\n se:", as.numeric(sqrt(diag(new_sigma))), "\n")
}
iter_SE_result[k,] <- sqrt(diag(new_sigma))
iter_norm_result[k] <- norm(new_beta-old_beta, "F")
if(iter_norm_result[k]>=100*iter_norm_result[1]) {
warning("Point estimation result is diverging")
break
}
if(norm(new_beta-old_beta, "i") < control$tol) break
}
} ## end for loop
}
iter_beta_result <- iter_beta_result[!is.na(iter_beta_result[,1]),,drop = FALSE]
iter_SE_result <- iter_SE_result[!is.na(iter_SE_result[,1]),,drop = FALSE]
iter_norm_result <- iter_norm_result[!is.na(iter_norm_result)]
out <- list(coefficient = as.numeric(tail(iter_beta_result, n = 1)),
trace.coefficient = iter_beta_result,
stderr = as.numeric(tail(iter_SE_result, n = 1)),
trace.stderr = iter_SE_result,
vcov = new_sigma, iterno = k,
norm = iter_norm_result)
})
out
}
qris.smooth <- function(info) {
out <- with(info, {
uTime <- sort(unique(data$Z))
smooth.fit <- nleqslv(betastart, function(b) isObj(b, X, W, H, I, logZ, Q))
if (smooth.fit$termcd == 1 | smooth.fit$termcd == 2) {
coefficient <- smooth.fit$x
if (nB <= 1) {
se <- rep(NA, nc)
vcov <- matrix(NA, nc, nc)
return(list(coefficient = coefficient, stderr = se, vcov = vcov))
}
if (se == "fmb"){
## Full multiplier bootstrap
smooth.fmb.result <- matrix(NA, nc, nB)
for (j in 1:nB){
## generating perturbation variable
eta <- rexp(n)
if (all(data$delta == 1)) {
W_star <- rep(1, n)
} else {
survp <- drop(ghatC(data$Z, 1 - data$delta, eta))
ghatstart0 <- ifelse(t0 > 0, survp[findInterval(t0, uTime)], 1)
W_star <- data$delta / survp[findInterval(data$Z, uTime)] * ghatstart0
W_star[is.na(W_star)] <- max(W_star, na.rm = TRUE)
}
fmb.fit <- nleqslv(coefficient, function(b) rev_isObj(b, X, W_star, H, eta, I, logZ, Q))
if (fmb.fit$termcd == 1 | fmb.fit$termcd == 2) {
smooth.fmb.result[,j] <- fmb.fit$x
}
}
fmb.sigma <- try(cov(t(smooth.fmb.result), use = "complete.obs"), silent = T)
if(class(fmb.sigma)[1] == "try-error" | sum(!is.na(colSums(smooth.fmb.result))) <= 1){
stop("More resampling iterations are needed")
} else {
fmb.se <- sqrt(diag(fmb.sigma))
out <- list(coefficient = coefficient, stderr = fmb.se, vcov = fmb.sigma)
}
} else if (se == "pmb") {
## Partial Multiplier Bootstrap
smooth.pmb.result <- isObjE(coefficient, X, H, I, logZ, data$delta, t0, Q, nB)
pmb.v <- try(cov(t(smooth.pmb.result), use = "complete.obs"), silent = T)
pmb.a <- Amat(coefficient, X, W, H, I, logZ, Q) / n
## Singular matrix 'a' error message and break
if (qr(pmb.a)$rank < min(dim(pmb.a))) {
pmb.se <- rep(NA,nc+1)
stop("'A' matrix is singular during iteration. Please try the nonsmooth method.")
} else {
pmb.inva <- qr.solve(pmb.a)
pmb.sigma <- t(pmb.inva) %*% pmb.v %*% pmb.inva
pmb.se <- sqrt(diag(pmb.sigma))
}
out <- list(coefficient=coefficient, stderr = pmb.se, vcov = pmb.sigma)
} else {
stop("Select either 'fmb' or 'pmb'")
}
} else {
coefficient <- se <- rep(NA, nc)
vcov <- matrix(NA, nc, nc)
out <- list(coefficient = coefficient, stderr = se, vcov = vcov)
}})
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.