# R/LF_logistic.R In SIHR: Statistical Inference in High Dimensional Regression

#### Documented in ITE_LogisticLF_logistic

#' Inference for the case probability or a linear combination of regression coefficients in high dimensional logistic regression.
#'
#' @description
#' Computes the bias corrected estimator of the case probability or the linear combination of regression coefficients in the high dimensional logistic regression model and the corresponding standard error.
#' It also constructs the confidence interval for the case probability or the linear combination and tests whether the case probability is above \eqn{0.5} or not.
#' Here case probability refers to the conditional probability of the binary response variable taking value 1 given the predictors take value \code{loading}.
#'
#' @param X Design matrix, of dimension \eqn{n} x \eqn{p}
#' @param y Outcome vector, of length \eqn{n}
#' @param weight The weight vector used for bias correction, of length \eqn{n}; if set to \code{NULL}, the weight is the inverse of the first derivative of the logit function (default = \code{NULL})
#' @param trans Should results for the case probability (\code{TRUE}) or the linear combination (\code{FALSE}) be reported (default = \code{TRUE})
#' @param intercept Should intercept be fitted for the initial estimator (default = \code{TRUE})
#' @param init.coef Initial estimator of the regression vector (default = \code{NULL})
#' @param lambda The tuning parameter used in the construction of \code{init.coef} (default = \code{NULL})
#' @param mu The dual tuning parameter used in the construction of the projection direction (default = \code{NULL})
#' @param step The step size used to compute \code{mu}; if set to \code{NULL} it is
#' computed to be the number of steps (< \code{maxiter}) to obtain the smallest \code{mu}
#' such that the dual optimization problem for constructing the projection direction converges (default = \code{NULL})
#' @param resol The factor by which \code{mu} is increased/decreased to obtain the smallest \code{mu}
#' such that the dual optimization problem for constructing the projection direction converges (default = 1.5)
#' @param maxiter Maximum number of steps along which \code{mu} is increased/decreased to obtain the smallest \code{mu}
#' such that the dual optimization problem for constructing the projection direction converges (default = 6)
#' @param alpha Level of significance to test the null hypothesis that the case probability is less than or equal to 0.5 (default = 0.05)
#' @param verbose Should inetrmediate message(s) be printed (default = \code{TRUE})
#'
#' @return
#' \item{prop.est}{The bias corrected estimator of the case probability or the linear combination of regression coefficients}
#' \item{se}{The standard error of the bias-corrected estimator}
#' \item{CI}{The confidence interval for the case probability or the linear combination}
#' \item{decision}{\code{decision}\eqn{=1} implies the case probability is above 0.5 \eqn{\newline}
#' \code{decision}\eqn{=0} implies the case probability is not above 0.5}
#' \item{proj}{The projection direction, of length \eqn{p}}
#' \item{plug.in}{The plug-in estimator of the case probability or the linear combination}
#'
#' @export
#'
#' @importFrom Rdpack reprompt
#' @importFrom stats coef na.omit
#' @import CVXR Matrix glmnet
#'
#' @references
#'
#' \insertRef{linlog}{SIHR}
#' @examples
#' \donttest{
#' A1gen <- function(rho,p){
#'  A1 <- matrix(0,p,p)
#'  for(i in 1:p){
#'    for(j in 1:p){
#'      A1[i,j] <- rho^(abs(i-j))
#'    }
#'  }
#'  A1
#' }
#' n <- 100
#' p <- 400
#' mu <- rep(0,p)
#' rho <- 0.5
#' Cov <- (A1gen(rho,p))/2
#' beta <- rep(0,p)
#' beta[1:10] <- 0.5*c(1:10)/10
#' a0 <- 0
#' X <- MASS::mvrnorm(n,mu,Cov)
#' exp_val <- X%*%beta+a0
#' prob <- exp(exp_val)/(1+exp(exp_val))
#' y <- rbinom(n,1,prob)
#' }
LF_logistic <- function(X, y, loading, weight = NULL, trans = TRUE, intercept.loading = TRUE, intercept = TRUE, init.coef = NULL, lambda = NULL, mu = NULL, step = NULL, resol = 1.5, maxiter = 6, alpha = 0.05, verbose = TRUE){
p <- ncol(X)
n <- nrow(X)
n_y <- length(y)

if(n_y!=n)
{
stop("Error: Check dimensions of X and y")
} else {
data <- na.omit(data.frame(y,X))
X <- as.matrix(data[,-1])
y <- as.vector(data[,1])
p <- ncol(X)
n <- nrow(X)
mean = colMeans(X)
M = matrix(rep(mean,nrow(X)),byrow = T, nrow = nrow(X), ncol = ncol(X))
X = X - M
xnew = xnew - mean
}
col.norm <- 1/sqrt((1/n)*diag(t(X)%*%X)+0.0001)
Xnor <- X %*% diag(col.norm)
if(is.null(init.coef))
{
init.coef <-  Initialization.step_log(X,y,lambda,intercept)
htheta <- init.coef$lasso.est } else { htheta <- init.coef } if (intercept==TRUE){ Xb <- cbind(rep(1,n),Xnor) Xc <- cbind(rep(1,n),X) col.norm <- c(1,col.norm) pp <- (p+1) } else { Xb <- Xnor Xc <- X pp <- p } sparsity <- sum(abs(htheta) > 0.001) sd.est <- sqrt(sum((y - Xb %*% htheta)^2) / max(0.9*n, n - sparsity)) if(intercept == TRUE){ loading <- rep(0,pp) if(intercept.loading == TRUE){ loading[1] <- 1 } if(intercept.loading == FALSE){ loading[1] <- 0 } loading[-1] <- xnew } else { if(intercept.loading == TRUE){ print(paste("Setting intercept = FALSE and intercept.loading = FALSE")) } loading <- xnew } loading.norm <- sqrt(sum(loading^2)) lasso.plugin <- sum(loading*htheta) deriv.vec <- exp(Xc%*%htheta)/(1+exp(Xc%*%htheta))^2 if(is.null(weight)){ weight <- 1/deriv.vec } X.weight = diag(c(sqrt(deriv.vec*weight))) %*% Xc ####modified count=0 for(i in 1:ncol(X)){ if(length(unique(X[,i])) == 1){ count = count+1 } } if(count!=0 && intercept==TRUE) { stop("Data is singular") }else{ if ((n >= 6*p)){ gamma.hat <- (1/n)*(t(X.weight)%*%X.weight) tmp <- eigen(gamma.hat) tmp <- min(tmp$values)/max(tmp$values) }else{ tmp <- 0 } if ((n >= 6*p) && (tmp >= 1e-4)){ direction <- solve(gamma.hat)%*%loading/loading.norm }else{ if(is.null(step)){ step.vec<-rep(NA,3) for(t in 1:3){ index.sel <- sample(1:n,size=ceiling(0.5*min(n,p)), replace=FALSE) Direction.Est.temp <- Direction_searchtuning_glm(Xc[index.sel,], loading, mu = NULL, weight = weight[index.sel], deriv.vec = deriv.vec[index.sel], resol, maxiter) step.vec[t] <- Direction.Est.temp$step
}
step<- getmode_log(step.vec)
}
Direction.Est <-  Direction_fixedtuning_glm(Xc, loading, mu = sqrt(2.01*log(pp)/n)*resol^{-(step-1)}, weight = weight, deriv.vec = deriv.vec)

while(is.na(Direction.Est) || length(Direction.Est$proj)==0){ step <- step-1 Direction.Est <- Direction_fixedtuning_glm(Xc, loading, mu = sqrt(2.01*log(pp)/n)*resol^{-(step-1)}, weight = weight, deriv.vec = deriv.vec) } if(verbose == TRUE){ print(paste("step is", step)) } direction <- Direction.Est$proj
}
exp_pred <- Xc%*%(htheta)

weighed.residual <- (y - exp(exp_pred)/(1+ exp(exp_pred)))*weight

correction <- sum((Xc%*%direction)*weighed.residual)/n
rho_hat <- exp(debias.est)/(1+exp(debias.est))^2
CI_linear <- c(debias.est - qnorm(1-alpha/2)*se_linear, debias.est + qnorm(1-alpha/2)*se_linear)
if(debias.est - qnorm(1-alpha)*se_linear > 0){
dec <- 1
}else{
dec <- 0
}
if(trans == TRUE){
prop.est = expo(debias.est)
se = rho_hat*se_linear
CI = c(expo(CI_linear[1]),expo(CI_linear[2]))
plug.in = expo(lasso.plugin)
}else{
prop.est = debias.est
se = se_linear
CI = CI_linear
plug.in = lasso.plugin
}
returnList <- list("prop.est" = prop.est,
"se" = se,
"CI" = CI,
"decision" = dec,
"proj"=direction,
"plug.in"= plug.in
)
return(returnList)
}
}
}

#' Inference for difference of case probabilities in high dimensional logistic regressions
#'
#' @description
#' Computes the bias corrected estimator of the difference between case probabilities or a linear combination of the difference between two regression vectors with respect to two high dimensional logistic regression models
#' and the corresponding standard error. It also constructs the confidence interval for the difference of case probabilities or a linear combination of the difference between the regression vectors and test
#' whether it is above zero or not. Here the case probability refers to the conditional probability of the binary response variable taking value 1 given the predictors are assigned to \code{loading}.
#'
#' @param X1 Design matrix for the first sample, of dimension \eqn{n_1} x \eqn{p}
#' @param y1 Outcome vector for the first sample, of length \eqn{n_1}
#' @param X2 Design matrix for the second sample, of dimension \eqn{n_2} x \eqn{p}
#' @param y2 Outcome vector for the second sample, of length \eqn{n_2}
#' @param weight The weight vector used for bias correction, of length \eqn{n}; if set to \code{NULL}, the weight is
#' the inverse of the first derivative of the logit function (default = \code{NULL})
#' @param trans Should results for the case probability (\code{TRUE}) or the linear combination (\code{FALSE}) be reported (default = \code{TRUE})
#' @param intercept Should intercept(s) be fitted for the initial estimators (default = \code{TRUE})
#' @param init.coef1 Initial estimator of the first regression vector (default = \code{NULL})
#' @param init.coef2 Initial estimator of the second regression vector (default = \code{NULL})
#' @param lambda1 The tuning parameter in the construction of \code{init.coef1} (default = \code{NULL})
#' @param lambda2 The tuning parameter in the construction of \code{init.coef2} (default = \code{NULL})
#' @param mu1 The dual tuning parameter used in the construction of the first projection direction (default = \code{NULL})
#' @param mu2 The dual tuning parameter used in the construction of the second projection direction (default = \code{NULL})
#' @param step1 The step size used to compute \code{mu1}; if set to \code{NULL} it is
#' computed to be the number of steps (< \code{maxiter}) to obtain the smallest \code{mu1}
#' such that the dual optimization problem for constructing the projection direction converges (default = \code{NULL})
#' @param step2 The step size used to compute \code{mu2}; if set to \code{NULL} it is
#' computed to be the number of steps (< \code{maxiter}) to obtain the smallest \code{mu2}
#' such that the dual optimization problem for constructing the second projection direction converges (default = \code{NULL})
#' @param resol The factor by which \code{mu1} (and \code{mu2}) is increased/decreased to obtain the smallest \code{mu1} (and \code{mu2})
#' such that the dual optimization problem for constructing the first (and the second) projection direction converges (default = 1.5)
#' @param maxiter Maximum number of steps along which \code{mu1} (and \code{mu2}) is increased/decreased to obtain the smallest \code{mu} (and \code{mu2})
#' such that the dual optimization problem for constructing the first (and the second) projection direction converges (default = 6)
#' @param alpha Level ofsignificance to test the null hypothesis which claims that the first case probability is not greater than the second case probability (default = 0.05)
#' @param verbose Should inetrmediate message(s) be printed (default = \code{TRUE})
#'
#' @return
#' \item{prop.est}{The bias-corrected estimator for the difference between case probabilities or the linear combination of the difference between two regression vectors}
#' \item{se}{The standard error for the bias-corrected estimator}
#' \item{CI}{The confidence interval for the difference between case probabilities or the linear combination of the difference between two regression vectors}
#' \item{decision}{\code{decision}\eqn{=1} implies the first case probability or linear combination is greater than the second one\eqn{\newline}
#' \code{decision}\eqn{=0} implies the first case probability or linear combination is less than the second one}
#' @export
#'
#' @importFrom Rdpack reprompt
#' @importFrom stats coef na.omit
#' @import CVXR Matrix glmnet
#' @examples
#' \donttest{
#' A1gen <- function(rho,p){
#' A1 <- matrix(0,p,p)
#' for(i in 1:p){
#'   for(j in 1:p){
#'     A1[i,j] <- rho^(abs(i-j))
#'   }
#' }
#' A1
#' }
#' n1 <- 100
#' n2 <- 100
#' p <- 400
#' mu <- rep(0,p)
#' rho <- 0.5
#' Cov <- (A1gen(rho,p))/2
#' beta1 <- rep(0,p)
#' beta1[1:10] <- c(1:10)/5
#' beta2 <- rep(0,p)
#' beta2[1:5] <- c(1:5)/10
#' X1 <- MASS::mvrnorm(n1,mu,Cov)
#' X2 <- MASS::mvrnorm(n2,mu,Cov)
#' exp_val1 <- X1%*%beta1
#' exp_val2 <- X2%*%beta2
#' prob1 <- exp(exp_val1)/(1+exp(exp_val1))
#' prob2 <- exp(exp_val2)/(1+exp(exp_val2))
#' y1 <- rbinom(n1,1,prob1)
#' y2 <- rbinom(n2,1,prob2)
#' }
ITE_Logistic <- function(X1, y1, X2, y2, loading, weight = NULL, trans = TRUE, intercept = TRUE, intercept.loading = TRUE, init.coef1 = NULL, init.coef2 = NULL, lambda1 = NULL, lambda2 = NULL, mu1 = NULL, mu2 = NULL, step1 = NULL, step2 = NULL, resol = 1.5, maxiter = 6, alpha = 0.05, verbose = TRUE){
if(trans == TRUE){
prop.est = expo(Est1$prop.est) - expo(Est2$prop.est)
se <- sqrt(((exp(Est1$prop.est)/(1+exp(Est1$prop.est))^2)*Est1$se)^2 + ((exp(Est2$prop.est)/(1+exp(Est2$prop.est))^2)*Est2$se)^2)
}
else{
prop.est <- Est1$prop.est - Est2$prop.est
se <- sqrt((Est1$se)^2 + (Est2$se)^2)
}
CI <- c((Est1$prop.est - Est2$prop.est) - qnorm(1-alpha/2)*sqrt((Est1$se)^2 + (Est2$se)^2), (Est1$prop.est - Est2$prop.est) + qnorm(1-alpha/2)*sqrt((Est1$se)^2 + (Est2$se)^2))
if(trans == TRUE){
if(CI[1]<-1){
CI[1] = -1
}
if(CI[2] > 1){
CI[2] = 1
}
}else{
CI = CI
}
if((Est1$prop.est - Est2$prop.est) - qnorm(1-alpha)*sqrt((Est1$se)^2 + (Est2$se)^2) > 0){
dec <- 1
}else{
dec <- 0
}
returnList <- list("prop.est" = prop.est,
"se"=se,
"CI"=CI,
"decision" = dec
)
return(returnList)
}


## Try the SIHR package in your browser

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

SIHR documentation built on Oct. 7, 2021, 9:08 a.m.