# R/BSDesparseLasso.R In WaverlyWei/debiased.subgroup: Sharp Inference on Selected Subgroup in Observational Studies

#### Documented in BSDesparseLasso

#' Bootstrap-calibrated Desparsified Lasso
#'
#'This method first constructs the debiased estimator of \eqn{\beta} via the
#'desparsified Lasso procedure. Then it calculates the calibration term
#' \eqn{\hat{b}_{max} =(1-n^{r-0.5})(\hat{\beta}_{max}-\hat{\beta}_{j,lasso})}. Through B bootstrap iterations,
#' it recalibrates the bootstrap statistic \eqn{T_b}. The bias-reduced estimate
#' is computed as: \eqn{\hat{b}_{max}-\frac{1}{B}\sum_{b=1}^BT_b}.
#'
#'
#' @param y response
#' @param x design matrix
#' @param r tuning parameter
#' @param G subgroup indicator
#' @param B bootstrap iterations
#' @param alpha level of CI
#' @return
#' \item{LowerBound}{lower confidence bound}
#' \item{UpperBound}{upper confidence bound}
#' \item{betaMax}{bias-reduced maximum beta estimate}
#' \item{betaEst}{debiased beta estimate for each subgroup}
#' \item{op}{optimal tuning}
#' @export
BSDesparseLasso <- function(y, x,
r = NULL,
G = NULL,
B = NULL,
alpha = 0.95,
fold = 3){

if(is.null(r)){

stop("Tuning parameter is missing.")

}else if(is.null(G)){

stop("Number of subgroups is missing.")

}else if(is.null(B)){

stop("Specify bootstrap iterations.")
}

p <- length(x[1,])

n <- length(y)

k <- length(G) # number of subgroups

cc <- length(r) # number of candidate tuning parameters

y <- y - mean(y)

for(i in 1:p)
x[,i] <- x[,i]-mean(x[,i])

# lasso estimates
fit.lasso <- cv.glmnet(x = x, y = y)

lambda <- fit.lasso$lambda.1se lambda <- lambda *1.1 gamma.lasso <- coef(fit.lasso, s = lambda) beta.lasso <- gamma.lasso[G+1] # calculate residual pred <- gamma.lasso[1] + x %*% gamma.lasso[-1] residual <- y - pred epsilion <- residual-mean(residual) #Desparsified lasso estimate beta.Dlasso <- 0 Z <- Zmatrix(x,G) for(i in 1:k){ index <- G[i] beta.Dlasso[i] <- beta.lasso[i] + sum(Z[,i]*residual)/sum(Z[,i]*x[,index]) } #calculate the correction term TB <- matrix(0, B, cc+1) correction <- matrix(0, cc+1, k) for(i in 1:cc) { r0 <- r[i] correction[i,] = (1-n^(r0-0.5))*(max(beta.lasso)-beta.lasso) } #the simulataneous one correction[cc+1,] <- (max(beta.lasso)-beta.lasso) TB_op <- matrix(0, B, cc) c_op <- matrix(0, cc, k) for(i in 1:cc) { r0 <- r[i] rp <- r0/sqrt(k/2) c_op[i,] <- (1-n^(rp-0.5))*(max(beta.lasso)-beta.lasso) } for(i in 1:B){ #generate bootstrap sample Bepsilion <- rnorm(n)*epsilion By <- pred + Bepsilion #calculate bootstrap desparsified lasso Bfit.lasso <- cv.glmnet(x = x, y = By) Blambda <- Bfit.lasso$lambda.1se

Blambda <- Blambda * 1.1

Bgamma.lasso <- coef(Bfit.lasso, s = Blambda)

Bbeta.lasso <- Bgamma.lasso[G+1]

Bpred <- Bgamma.lasso[1] + x %*% Bgamma.lasso[-1]

Bresidual <- By - Bpred

Bbeta.Dlasso <- 0

for(j in 1:k){

Bindex <- G[j]
Bbeta.Dlasso[j] <- Bbeta.lasso[j] + sum(Z[,j]*Bresidual)/sum(Z[,j]*x[,Bindex])
}

#correct the maximum quantity
for(j in 1:(cc+1)){
TB[i,j] <- max(Bbeta.Dlasso+correction[j,])-max(beta.lasso)
}

for(j in 1:(cc)){

TB_op[i,j] <- max(Bbeta.Dlasso+c_op[j,])-max(beta.lasso)

}
}

#choose the optimal tuning parameter

op <- cvDesparse(y, x, r, G, B, fold) # index of the correction term

# collect results
result <- list()

for(j in 1:(cc+1)) {

result[j] = list(c(BSciCoverfun(beta.Dlasso, TB[,j], G, alpha),
betaEst = list(beta.Dlasso),
op = r[op]))
}

if(is.integer(op)){
result[j+1] = list(c(BSciCoverfun(beta.Dlasso, TB_op[,op], G, alpha),
betaEst = list(beta.Dlasso),
op = r[op]))
}else{
result[j+1] = list(c(BSciCoverfun(beta.Dlasso, TB[,cc], G, alpha),
betaEst = list(beta.Dlasso),
op = r[op]))
}
return(result[[12]])
}

WaverlyWei/debiased.subgroup documentation built on Jan. 27, 2021, 12:15 a.m.