cste_bin: Estimate the CSTE curve for binary outcome.

In CSTE: Covariate Specific Treatment Effect (CSTE) Curve

Description

Estimate covariate-specific treatment effect (CSTE) curve. Input data contains covariates X, treatment assignment Z and binary outcome Y. The working model is

logit(μ(X, Z)) = g_1(Xβ_1)Z + g_2(Xβ_2),

where μ(X, Z) = E(Y|X, Z). The model implies that CSTE(x) = g_1(xβ_1).

Usage

cste_bin(
  x,
  y,
  z,
  beta_ini = NULL,
  lam = 0,
  nknots = 1,
  max.iter = 200,
  eps = 0.001
)

Arguments

x	samples of covariates which is a n*p matrix.
y	samples of binary outcome which is a n*1 vector.
z	samples of treatment indicator which is a n*1 vector.
beta_ini	initial values for (β_1', β_2')', default value is NULL.
lam	value of the lasso penalty parameter λ for β_1 and β_2, default value is 0.
nknots	number of knots for the B-spline for estimating g_1 and g_2.
max.iter	maximum iteration for the algorithm.
eps	numeric scalar ≥q 0, the tolerance for the estimation of β_1 and β_2.

Value

A S3 class of cste, which includes:

• beta1: estimate of β_1.

• beta2: estimate of β_2.

• B1: the B-spline basis for estimating g_1.

• B2: the B-spline basis for estimating g_2.

• delta1: the coefficient of B-spline for estimating g_1.

• delta2: the coefficient for B-spline for estimating g_2.

• iter: number of iteration.

• g1: the estimate of g_1(Xβ_1).

• g2: the estimate of g_2(Xβ_2).

References

Guo W., Zhou X. and Ma S. (2021). Estimation of Optimal Individualized Treatment Rules Using a Covariate-Specific Treatment Effect Curve with High-dimensional Covariates, Journal of the American Statistical Association, 116(533), 309-321

See Also

cste_bin_SCB, predict_cste_bin, select_cste_bin

Examples

## Quick example for the cste

library(mvtnorm)
library(sigmoid)

# --------  Example 1: p = 20 ---------  #  
## generate data 
n <- 2000
p <- 20
set.seed(100)

# generate X
sigma <- outer(1:p, 1:p, function(i, j){ 2^(-abs(i-j)) } )
X <- rmvnorm(n, mean = rep(0,p), sigma = sigma)
X <- relu(X + 2) - 2
X <- 2 - relu(2 - X)

# generate Z
Z <- rbinom(n, 1, 0.5)

# generate Y
beta1 <- rep(0, p)
beta1[1:3] <- rep(1/sqrt(3), 3)
beta2 <- rep(0, p)
beta2[1:2] <- c(1, -2)/sqrt(5)
mu1 <- X %*% beta1
mu2 <- X %*% beta2
g1 <- mu1*(1 - mu1)
g2 <- exp(mu2)      
prob <- sigmoid(g1*Z + g2)
Y <- rbinom(n, 1, prob)

## estimate the CSTE curve
fit <- cste_bin(X, Y, Z)

## plot 
plot(mu1, g1, cex = 0.5, xlim = c(-2,2), ylim = c(-8, 3), 
     xlab = expression(X*beta), ylab = expression(g1(X*beta)))
     ord <- order(mu1)
     points(mu1[ord], fit$g1[ord], col = 'blue', cex = 0.5)
     
## compute 95% simultaneous confidence band (SCB)
res <- cste_bin_SCB(X, fit, alpha = 0.05)

## plot 
plot(res$or_x, res$fit_x, col = 'red', 
     type="l", lwd=2, lty = 3, ylim = c(-10,8),
     ylab=expression(g1(X*beta)), xlab = expression(X*beta), 
     main="Confidence Band")
lines(res$or_x, res$lower_bound, lwd=2.5, col = 'purple', lty=2)
lines(res$or_x, res$upper_bound, lwd=2.5, col = 'purple', lty=2)
abline(h=0, cex = 0.2, lty = 2)
legend("topleft", legend=c("Estimates", "SCB"), 
        lwd=c(2, 2.5), lty=c(3,2), col=c('red', 'purple'))
        
        
# --------  Example 2: p = 1 ---------  #  

## generate data 
set.seed(15)
p <- 1
n <- 2000
X <- runif(n)
Z <- rbinom(n, 1, 0.5)
g1 <- 2 * sin(5*X) 
g2 <- exp(X-3) * 2
prob <- sigmoid( Z*g1 + g2)
Y <- rbinom(n, 1, prob)

## estimate the CSTE curve
fit <- cste_bin(X, Y, Z)  

## simultaneous confidence band (SCB)
X <- as.matrix(X)
res <- cste_bin_SCB(X, fit)  

## plot 
plot(res$or_x, res$fit_x, col = 'red', type="l", lwd=2, 
     lty = 3, xlim = c(0, 1), ylim = c(-4, 4), 
     ylab=expression(g1(X)), xlab = expression(X), 
     main="Confidence Band")
lines(res$or_x, res$lower_bound, lwd=2.5, col = 'purple', lty=2)
lines(res$or_x, res$upper_bound, lwd=2.5, col = 'purple', lty=2)
abline(h=0, cex = 0.2)
lines(X[order(X)], g1[order(X)], col = 'blue', lwd = 1.5)
legend("topright", legend=c("Estimates", "SCB",'True CSTE Curve'), 
lwd=c(2, 2.5, 1.5), lty=c(3,2,1), col=c('red', 'purple','blue'))

CSTE documentation built on Dec. 16, 2021, 9:07 a.m.