R/utils.R
In muxuanliang/RSICF: Estimation and Inference for the Relative Contrast Function under a Monotonic Single-index Model Assumption
Defines functions
updateBeta
updateXi
fitLinkLinear
nsd
link
loss
screening
getPropensityModel
getOutcomeModel
model.gam
ks
rsfitSplit
rsfitSolver
predict.rsfitSplit
# utils
#
# This file codes each step in the algorithm proposed in the paper.
# In addition, it also provides defined link function can be used.
#
# predict.rsfit gets the predicted treatment for a splitted fit
predict.rsfitSplit <- function(fit, newx, derivative = FALSE,lossType = 'logistic', type = 'value'){
t <- fit
z <- newx %*% t$beta
spline <- splines2::bSpline(z, knots=t$knots, intercept = FALSE, Boundary.knots = t$boundaryPoint)
pre <- cbind(1,spline)%*%t$xi
if(sum(abs(t$xi))<1e-5){
pre <- z
}
if(derivative){
div <- nsd(z, knots = t$knots, intercept = FALSE, Boundary.knots = t$boundaryPoint)
pre <- cbind(0, div)%*%t$xi
if(sum(abs(t$xi))<1e-5){
pre <- 1
}
}
if (type=='eta'){
pre <- link(pre, lossType = lossType)
}
pre
}
# rsfitSolver obtains the optimizer of the proposed optimization given M_n
rsfitSolver <- function(covariate, response, treatment, estimatedNuisance, splitIndex = NULL, lossType = 'logistic', weights = NULL, tol = 10^(-3), m_0=0, constraint = TRUE, boundaryPoint=c(-1,1), numberKnots = 5){
# initialization
fit_earl <- ITRInference::ITRFitInfer(list(predictor = covariate, treatment = (treatment > 0), outcome = response), loss = 'logistic', sampleSplitIndex = splitIndex$nuisance, propensity = estimatedNuisance$pi, outcome = estimatedNuisance, test=FALSE)
#fit_init <- fitLinkLinear(covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex, lossType = lossType, weights = weights, tol = tol)
iter <- 0
# iteration start
fit_last <- NULL
tmp_beta<- (fit_earl$fit[[1]]$fit$beta[,fit_earl$fit[[1]]$fit$lambda==fit_earl$fit[[1]]$fit$lambda.min] + fit_earl$fit[[2]]$fit$beta[,fit_earl$fit[[2]]$fit$lambda==fit_earl$fit[[2]]$fit$lambda.min])/2
if (sum(tmp_beta^2)<10^-3){
fit_last$beta <- (tmp_beta+0.1)/sqrt(sum((tmp_beta+0.1)^2))
} else {
fit_last$beta <- tmp_beta/sqrt(sum(tmp_beta^2))
}
fit_last$xi <- 0
diff <- 1
step <- 0.1
while((diff > tol) & (iter < 100)){
if(any(is.na(fit_last$beta))){
print(fit_last)
print(fit_xi)
print(fit_beta)
print(iter)
}
# updateXi
fit_xi <- updateXi(fit_last, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = 1e-5, m_0=m_0, constraint=constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
# updateBeta
fit_beta <- updateBeta(fit_xi, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, withConstraint = TRUE, boundaryPoint = boundaryPoint, step=step)
# update no constraint
#fit_tmp <- updateBeta(fit_xi, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, withConstraint = FALSE, boundaryPoint = boundaryPoint)
#update
diff_coef <- max(abs(fit_beta$beta-fit_last$beta))
diff_xi <- max(c(abs(fit_beta$xi-fit_last$xi))) #abs(fit_beta$beta-fit_tmp$beta))
diff <- max(diff_coef, diff_xi)
iter <- iter+1
fit_last <- fit_beta
}
# final update
fit <- updateXi(fit_last, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, m_0=m_0, constraint=constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
fit_tmp <- updateBeta(fit, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, withConstraint = FALSE, boundaryPoint = boundaryPoint)
fit$solution <- fit_tmp$beta
fit$boundaryPoint <- boundaryPoint
fit
}
# rsfitSplit obtain the estimation of RSICF assuming a single index model on a slited data. The variance estimaor of the coefficiencts are provided.
rsfitSplit <- function(covariate, response, treatment, splitIndex = NULL, propensityModel = 'glmnet', estimatedPropensity = NULL, outcomeModel = 'kernel', estimatedOutcome = NULL, lossType = 'logistic', weights = NULL, tol = 1e-3, propensityFormula=NULL, outcomeFormula=NULL, constraint=TRUE, boundaryPoint=c(-1,1)){
# fit nuisance parameter
estimatedNuisance <- NULL
estimatedNuisance$pi <- estimatedPropensity
if (is.null(estimatedPropensity)){
data <- list(predictor = covariate, treatment = treatment)
predictedPropensityAll <- getPropensityModel(data, method = propensityModel, splitIndex = splitIndex, Formula = propensityFormula, predictAll = TRUE)
estimatedNuisance$pi <- predictedPropensityAll
}
estimatedNuisance$S <- estimatedOutcome$control+estimatedOutcome$treatment
estimatedNuisance$control <- estimatedOutcome$control
estimatedNuisance$treatment <- estimatedOutcome$treatment
if (is.null(estimatedOutcome)){
data <- list(predictor = covariate, treatment = treatment, outcome = response)
predictedOutcomeAll <- getOutcomeModel(data, method = outcomeModel, splitIndex = splitIndex, Formula = outcomeFormula, predictAll = TRUE)
estimatedNuisance$S <- predictedOutcomeAll$control+predictedOutcomeAll$treatment
estimatedNuisance$control <- predictedOutcomeAll$control
estimatedNuisance$treatment <- predictedOutcomeAll$treatment
}
# fit conditional variance
residual_square <- (response - (predictedOutcomeAll$control * (treatment==0) + predictedOutcomeAll$treatment * (treatment==1)))^2
data <- list(predictor = covariate, treatment = treatment, outcome = residual_square)
predictedVarAll <- getOutcomeModel(data, method = outcomeModel, splitIndex = splitIndex, Formula = outcomeFormula, predictAll = TRUE)
estimatedNuisance$var_control <- predictedVarAll$control
estimatedNuisance$var_treated <- predictedVarAll$treatment
# weights
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
# solve
m_0 <- 5
numberKnots <- 3
fit_hat <- rsfitSolver(covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex, lossType = lossType, weights =weights, tol = tol, m_0=m_0, constraint = constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
diff_min <- max(abs(fit_hat$solution-fit_hat$beta))
fit_hat_min <- fit_hat
dif_tmp <- diff_min
while((diff_min>10^(-3)) & (numberKnots <= 7) & (m_0 <= 2^5)){
early_step <- FALSE
if (sum(abs(fit_hat$xi[-1]))<m_0){
numberKnots <- numberKnots + 1
early_step <- TRUE
} else {
m_0 <- 2*m_0
}
fit_hat <- rsfitSolver(covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex, lossType = lossType, weights =weights, tol = tol, m_0=m_0, constraint = constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
dif_tmp <- max(abs(fit_hat$solution-fit_hat$beta))
if (dif_tmp < diff_min){
fit_hat_min <- fit_hat
diff_min <- dif_tmp
}
if ((dif_tmp > diff_min)&early_step&(sum(abs(fit_hat$xi[-1]))<m_0)){
break
}
if ((dif_tmp > diff_min)&(!early_step)&(sum(abs(fit_hat$xi[-1]))>=m_0)){
break
}
}
# return
list(fit=fit_hat_min, covariate=covariate, response=response, treatment=treatment, estimatedNuisance=estimatedNuisance, splitIndex=splitIndex, lossType = lossType, weights =weights)
}
# ks gets the kernel estimation
ks <- function(xx, yy, xx.test){
nobs <- nrow(as.matrix(xx))
nvars <- ncol(as.matrix(xx))
hopt <- (4/(nvars+2))^(1/(nvars+4)) * (nobs^(-1/(nvars+4)))
wm <- function(t){
if (ncol(as.matrix(xx))==1){
weight <- exp(-0.5 * (as.numeric(t)-xx)^2/(hopt^2)) * hopt
} else {
weight <- apply(xx,1,function(s){exp(-0.5 * sum((t-s)^2)/(hopt^2)) * hopt^(ncol(xx))})
}
weighted.mean(yy, weight)
}
if (nrow(as.matrix(xx.test))==1) {
yy.test <- wm(xx.test)
} else {
if (ncol((as.matrix(xx.test)))==1){
yy.test <- sapply(as.matrix(xx.test), function(t){
wm(t)
})
} else {
yy.test <- apply(as.matrix(xx.test),1,function(t){
wm(t)
})
}
}
yy.test
}
# model.gam gets the gam regression function given data
model.gam <- function(data){
p <- dim(data$predictor)[2]
expr <- "mgcv::gam(outcome~"
for (i in 1:(p-1)){
expr <- paste0(expr, "s(predictor[,",i,"])+")
}
expr <- paste0(expr, "s(predictor[,",p,"]), data = data,method ='REML')")
expr
}
# getOutcomeModel gets the outcome regression model
getOutcomeModel <- function(data, method=c('lm', 'glmnet', 'kernel', 'others'), splitIndex, Formula = NULL, predictAll = FALSE, screeningMethod="SIRS", outcomeScreeningFamily='Gaussian'){
sampleSplitIndex <- splitIndex$nuisance
p <- dim(data$predictor)[2]
size <- dim(data$predictor)[1]
fit <- NULL
supp <- NULL
dataPredict <- NULL
dataPredict$control=data$predictor[sampleSplitIndex,]
dataPredict$treatment=data$predictor[sampleSplitIndex,]
if (predictAll){
dataPredict$control=data$predictor
dataPredict$treatment=data$predictor
}
prediction <- NULL
dataControl <- list(predictor=data$predictor[(!sampleSplitIndex) & (data$treatment==FALSE),], outcome=data$outcome[(!sampleSplitIndex) & (data$treatment==FALSE)])
dataTreatment <- list(predictor=data$predictor[(!sampleSplitIndex) & (data$treatment==TRUE),], outcome=data$outcome[(!sampleSplitIndex) & (data$treatment==TRUE)])
if ((0.05*size >= p) | (p <= 5)){
supp$control <- supp$treatment <- rep(TRUE, times = p)
if ((method != 'lm')&&(method != 'glmnet')){
ans1 <- screening(data$predictor[data$treatment==FALSE,], data$outcome[data$treatment==FALSE], method = screeningMethod, family = outcomeScreeningFamily)
ans2 <- screening(data$predictor[data$treatment==TRUE,], data$outcome[data$treatment==TRUE], method = screeningMethod, family = outcomeScreeningFamily)
if (screeningMethod == 'glmnet'){
supp$control <- ans1
supp$treatment <- ans2
} else {
supp$control <- supp$treatment <- (ans1 <= floor(p/2))|(ans2 <= floor(p/2))
}
dataControl$predictor <- dataControl$predictor[,supp$control]
dataTreatment$predictor <- dataTreatment$predictor[,supp$treatment]
}
}
if ((0.05*size < p) || (method == 'glmnet')) {
fit$control <- glmnet::cv.glmnet(x = dataControl$predictor, y = dataControl$outcome)
fit$treatment <- glmnet::cv.glmnet(x = dataTreatment$predictor, y = dataTreatment$outcome)
supp$control <- abs(fit$control$glmnet.fit$beta[,fit$control$glmnet.fit$lambda==fit$control$lambda.min])>0
supp$treatment <- abs(fit$treatment$glmnet.fit$beta[,fit$treatment$glmnet.fit$lambda==fit$treatment$lambda.min])>0
if ((method != 'lm')&&(method != 'glmnet')){
ans1 <- screening(data$predictor[data$treatment==FALSE,], data$outcome[data$treatment==FALSE], method = screeningMethod, family = outcomeScreeningFamily)
ans2 <- screening(data$predictor[data$treatment==TRUE,], data$outcome[data$treatment==TRUE], method = screeningMethod, family = outcomeScreeningFamily)
if (screeningMethod == 'glmnet'){
supp$control <- ans1
supp$treatment <- ans2
} else {
supp$control <- supp$treatment <- (ans1 <= 5)|(ans2 <= 5)
}
}
dataControl$predictor <- dataControl$predictor[,supp$control]
dataTreatment$predictor <- dataTreatment$predictor[,supp$treatment]
prediction$control <- predict(fit$control, newx = dataPredict$control, s=fit$control$lambda.min)
prediction$treatment <- predict(fit$treatment, newx = dataPredict$treatment, s=fit$treatment$lambda.min)
}
if (is.null(Formula)){
Formula <- function(support){
expr <- (outcome ~ predictor)
if (sum(support)==0){
expr <- (outcome ~ 1)
}
expr
}
}
fit <- NULL
dataPredict <- NULL
dataPredict$control=data$predictor[sampleSplitIndex,supp$control]
dataPredict$treatment=data$predictor[sampleSplitIndex,supp$treatment]
if (predictAll){
dataPredict$control=data$predictor[,supp$control]
dataPredict$treatment=data$predictor[,supp$treatment]
}
if ((method == 'lm')||(method == 'glmnet')){
if (sum(supp$control) > 0){
fit$control <- lm(Formula(supp$control), data = dataControl)
prediction$control <- predict(fit$control, newdata = list(predictor=dataPredict$control))
}
if (sum(supp$treatment) > 0){
fit$treatment <- lm(Formula(supp$treatment), data = dataTreatment)
prediction$treatment <- predict(fit$treatment, newdata = list(predictor=dataPredict$treatment))
}
} else if (method == 'kernel') {
if (sum(supp$control) > 0){
prediction$control <- ks(dataControl$predictor, dataControl$outcome, dataPredict$control)
}
if (sum(supp$treatment) > 0){
prediction$treatment <- ks(dataTreatment$predictor, dataTreatment$outcome, dataPredict$treatment)
}
} else {
if (sum(supp$control) > 0){
fit$control <- eval(parse(text=model.gam(dataControl)))
prediction$control <- predict(fit$control, newdata = list(predictor=dataPredict$control))
}
if (sum(supp$treatment) > 0){
fit$treatment <- eval(parse(text=model.gam(dataTreatment)))
prediction$treatment <- predict(fit$treatment, newdata = list(predictor=dataPredict$treatment))
}
}
prediction
}
# getPropensityModel gets the estimation of the propensity model based on selected method
getPropensityModel <- function(data, method=c('lm', 'glmnet', 'kernel'), splitIndex, Formula = NULL, predictAll = FALSE, screeningMethod="SIRS"){
sampleSplitIndex <- splitIndex$nuisance
p <- dim(data$predictor)[2]
size <- dim(data$predictor)[1]
fit <- NULL
supp <- NULL
dataPredict <- NULL
dataPredict=data$predictor[sampleSplitIndex,]
if (predictAll){
dataPredict=data$predictor
}
prediction <- NULL
dataTrain <- list(predictor=data$predictor[(!sampleSplitIndex),], treatment=data$treatment[(!sampleSplitIndex)])
if (0.05*size >= p){
supp$control <- supp$treatment <- rep(TRUE, times = p)
if ((method != 'lm')&&(method != 'glmnet')){
ans <- screening(data$predictor, data$treatment, method = screeningMethod, family = 'binomial')
if (screeningMethod == 'glmnet'){
supp <- ans
} else {
supp <- (ans <= p/2)
}
}
dataTrain$predictor = dataTrain$predictor[,supp]
}
if ((0.05*size < p) || (method == 'glmnet')) {
fit <- glmnet::cv.glmnet(x = dataTrain$predictor, y = dataTrain$treatment, family='binomial')
supp <- abs(fit$glmnet.fit$beta[,fit$glmnet.fit$lambda==fit$lambda.min])>0
if ((method != 'lm')&&(method != 'glmnet')){
ans <- screening(data$predictor, data$treatment, method = screeningMethod, family = 'binomial')
if (screeningMethod == 'glmnet'){
supp <- ans
} else {
supp <- (ans <= 5)
}
}
dataTrain$predictor <- dataTrain$predictor[,supp]
prediction <- predict(fit, newx = dataPredict, type='response', s=fit$lambda.min)
}
if (is.null(Formula)){
Formula <- function(support){
expr <- (treatment ~ predictor)
if (sum(support)==0){
expr <- (treatment ~ 1)
}
expr
}
}
fit <- NULL
dataPredict <- NULL
dataPredict=data$predictor[sampleSplitIndex,supp]
if (predictAll){
dataPredict=data$predictor[,supp]
}
if ((method == 'lm')||(method == 'glmnet')){
if (sum(supp) > 0){
fit <- glm(Formula(supp), family=binomial, data = dataTrain)
prediction <- predict(fit, newdata = list(predictor=dataPredict), type="response")
}
} else if (method == 'kernel') {
if (sum(supp) > 0){
prediction <- ks(dataTrain$predictor, dataTrain$treatment, dataPredict)
prediction <- (prediction > 0.9) * 0.9 + (prediction < 0.1) * 0.1 + (prediction < 0.9) * (prediction > 0.1) * prediction
}
}
prediction
}
# screening gets the screened variable for high dimeniosnal data
screening <- function(x, y, method='glmnet', family='Gaussian'){
var <- apply(x, 2, sd)
supp <- order(var, decreasing = TRUE)
if (method=='glmnet'){
fit <- glmnet::cv.glmnet(x, y, family = family)
coef <- fit$glmnet.fit$beta[,fit$lambda==fit$lambda.min]
supp <- (abs(coef)>0)
} else {
fit <- VariableScreening::screenIID(x, y, method=method)
supp <- fit$rank
}
supp
}
# loss reutrn the value or the derivative of the chosen type of link
loss <- function(input, type = c('logistic', 'exponential'), order = 'original'){
return(switch(order, 'original'=switch(type, 'logistic'=exp(-input)/2, 'exponential'=exp(-2*input)/4),
'derivative'=switch(type, 'logistic'=-exp(-input)/2, 'exponential'=-exp(-2*input)/2),
'hessian'=switch(type, 'logistic'=exp(-input)/2, 'exponential'=exp(-2*input))))
}
# link returns the eta from eta_trans
link <- function(eta_trans, lossType='logistic'){
return(switch(lossType, 'logistic'=(exp(eta_trans/2)-exp(-eta_trans/2))/(exp(eta_trans/2)+exp(-eta_trans/2)),
'exponential'=(exp(eta_trans)-exp(-eta_trans))/(exp(eta_trans)+exp(-eta_trans))))
}
# nsd get the derivative of the ns
nsd <- function(x, knots = NULL, intercept = FALSE, Boundary.knots = range(x)){
div <- sapply(x,function(t){
if((t+1e-6)>=Boundary.knots[2]){
res <- (splines2::bSpline(t, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots)-splines2::bSpline(t-1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots))/(1e-6)
} else if ((t-1e-6)<=Boundary.knots[1]){
res <- (splines2::bSpline(t+1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots)-splines2::bSpline(t, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots))/(1e-6)
} else {
res <- (splines2::bSpline(t+1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots)-splines2::bSpline(t-1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots))/(2*1e-6)
}
res
})
t(div)
}
# fitLinkLinear impliments the Step 1 of the porposed algorithm.
# stimatedNuissance should contain two lists. One for propensity, the other for the estimated S.
# The predictions nad inputs is always for all the data.
# SplitIndex containes three list. $nuisance, $fit ,$efficient
# treatment is 0 or 1
# covariate does not contain intercept
fitLinkLinear <- function(covariate, response, treatment, estimatedNuisance, splitIndex = NULL, lossType = 'exponential', weights = NULL, tol = 1e-3, intercept = TRUE){
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
pi <- estimatedNuisance$pi[splitIndex$fit]
pi_T <- pi * treatment[splitIndex$fit] + (1-pi) * (1-treatment[splitIndex$fit])
S <- estimatedNuisance$S[splitIndex$fit]
trt_train <- 2*(treatment[splitIndex$fit]-0.5)
x_train <- covariate[splitIndex$fit,]
y_train <- response[splitIndex$fit]
w_train <- weights[splitIndex$fit]
if(intercept){
x_train <- cbind(1, covariate[splitIndex$fit,])
}
obj <- function(beta){
value <- weighted.mean(y_train/pi_T*loss(trt_train * x_train %*% beta, type = lossType)+trt_train/(2*pi_T)*(S-y_train)*x_train %*% beta, w_train)
grad <- (apply(apply(x_train,2,function(z){(1/pi_T*y_train*trt_train* loss(trt_train * x_train %*% beta, type = lossType, order = 'derivative')+trt_train/(2*pi_T)*(S-y_train))*z}),
2, function(t){weighted.mean(t, w_train)}))
list(objective = value, gradient = grad)
}
beta0 <- array(0,c(NCOL(x_train),1))
fit <- nloptr::nloptr(x0=beta0, eval_f = obj, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
if(intercept){
fit$solution <- fit$solution[-1]
}
fit$solution <- fit$solution/sqrt(sum(fit$solution^2))
fit <- list(beta=fit$solution, xi=NULL, knods = NULL, solution=NULL)
fit
}
# updateXi implements the Step 2 of the proposed algorithm
updateXi <- function(fit, covariate, response, treatment, estimatedNuisance, lossType='tanh', splitIndex = NULL, weights = NULL, tol = 1e-3, numberKnots = 5, m_0 =NULL, constraint = TRUE, boundaryPoint=c(-1,1)){
beta_last <- fit$beta
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
pi <- estimatedNuisance$pi[splitIndex$fit]
pi_T <- pi * treatment[splitIndex$fit] + (1-pi) * (1-treatment[splitIndex$fit])
S <- estimatedNuisance$S[splitIndex$fit]
trt_train <- 2*(treatment[splitIndex$fit]-0.5)
x_train <- covariate[splitIndex$fit,]
y_train <- response[splitIndex$fit]
w_train <- weights[splitIndex$fit]
z_train <- x_train %*% beta_last
knots <- quantile(z_train, probs=seq(0,1,length.out = numberKnots))
knots[1] <- knots[1]-0.3*abs(knots[2]-knots[1])
knots[numberKnots] <- knots[numberKnots]+0.3*abs(knots[numberKnots]-knots[numberKnots-1])
ns_train <- splines2::bSpline(z_train, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint)
ns_train <- cbind(1, ns_train)
# tuning start
obj <- function(xi){
value <- weighted.mean(y_train/pi_T*loss(trt_train * ns_train %*% xi, type = lossType)+trt_train/(2*pi_T)*(S-y_train)*ns_train %*% xi, w_train)
grad <- (apply(apply(ns_train,2,function(z){(1/pi_T*y_train*trt_train* loss(trt_train * ns_train %*% xi, type = lossType, order = 'derivative')+trt_train/(2*pi_T)*(S-y_train))*z}),
2, function(t){weighted.mean(t, w_train)}))
list(objective = value, gradient = grad)
}
#form contrast matrix
B <- matrix(0,NCOL(ns_train)-1, NCOL(ns_train))
for(i in 2:(NCOL(ns_train)-1)){
B[i,i] <- 1
B[i, i+1] <- -1
}
g_ineq <- function(xi){
value <- c(sum(abs(xi[-1]))-m_0, c(B%*% xi))
grad <- rbind(c(0,sign(xi[-1])), B)
list(constraints = value, jacobian = grad)
}
if(!constraint){
g_ineq <- function(xi){
value <- c(sum(abs(xi[-1]))-m_0)
grad <- c(0,sign(xi[-1]))
list(constraints = value, jacobian = grad)
}
}
fit_xi <- nloptr::nloptr(x0=rep(0, times=NCOL(ns_train)), eval_f = obj, eval_g_ineq = g_ineq, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
z_test <- seq(min(knots),max(knots), length.out=100)
ns_test <- splines2::bSpline(z_test, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint)
ns_test <- cbind(1, ns_test)
fit_test <- ns_test %*% fit_xi$solution
plot(z_test, fit_test)
fit <- list(beta=beta_last, xi=fit_xi$solution, knots = knots,solution=fit$solution)
fit
}
# updateBeta implements the Step 3 of the proposed algorithm
updateBeta <- function(fit, covariate, response, treatment, estimatedNuisance, lossType='tanh', splitIndex = NULL, weights = NULL, tol = 1e-3, withConstraint = TRUE, boundaryPoint = c(-1,1), step=0.1){
beta_last <- fit$beta
xi <- fit$xi
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
pi <- estimatedNuisance$pi[splitIndex$fit]
pi_T <- pi * treatment[splitIndex$fit] + (1-pi) * (1-treatment[splitIndex$fit])
S <- estimatedNuisance$S[splitIndex$fit]
trt_train <- 2*(treatment[splitIndex$fit]-0.5)
x_train <- covariate[splitIndex$fit,]
y_train <- response[splitIndex$fit]
w_train <- weights[splitIndex$fit]
knots <- fit$knots
# objective function
obj <- function(beta){
z_train <- x_train %*% beta
ns_train <- splines2::bSpline(z_train, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint)
ns_train <- cbind(1, ns_train)
nsd_train <- cbind(0,nsd(z_train, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint))
value <- weighted.mean(y_train/pi_T*loss(trt_train * ns_train %*% xi, type = lossType)+trt_train/(2*pi_T)*(S-y_train)*ns_train %*% xi, w_train)
grad <- (apply(apply(x_train,2,function(z){(1/pi_T*y_train*trt_train* loss(trt_train * ns_train %*% xi, type = lossType, order = 'derivative')+trt_train/(2*pi_T)*(S-y_train))*nsd_train%*% xi*z}),
2, function(t){weighted.mean(t, w_train, na.rm = TRUE)}))
list(objective = value, gradient = grad)
}
constraint <- function(beta){
return(list('constraints'=sum(beta^2)-1,
'jacobian'=2*beta))
}
if (withConstraint){
fit_beta <- nloptr::nloptr(x0=fit$beta, eval_f = obj, eval_g_eq = constraint, lb=beta_last-step, ub=beta_last+step, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
} else {
fit_beta <- nloptr::nloptr(x0=fit$beta, eval_f = obj, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
}
fit <- list(beta=fit_beta$solution, xi=fit$xi, knots = fit$knots)
fit
}
muxuanliang/RSICF documentation built on Feb. 1, 2022, 12:30 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions updateBeta updateXi fitLinkLinear nsd link loss screening getPropensityModel getOutcomeModel model.gam ks rsfitSplit rsfitSolver predict.rsfitSplit
# utils
#
# This file codes each step in the algorithm proposed in the paper.
# In addition, it also provides defined link function can be used.
#
# predict.rsfit gets the predicted treatment for a splitted fit
predict.rsfitSplit <- function(fit, newx, derivative = FALSE,lossType = 'logistic', type = 'value'){
t <- fit
z <- newx %*% t$beta
spline <- splines2::bSpline(z, knots=t$knots, intercept = FALSE, Boundary.knots = t$boundaryPoint)
pre <- cbind(1,spline)%*%t$xi
if(sum(abs(t$xi))<1e-5){
pre <- z
}
if(derivative){
div <- nsd(z, knots = t$knots, intercept = FALSE, Boundary.knots = t$boundaryPoint)
pre <- cbind(0, div)%*%t$xi
if(sum(abs(t$xi))<1e-5){
pre <- 1
}
}
if (type=='eta'){
pre <- link(pre, lossType = lossType)
}
pre
}
# rsfitSolver obtains the optimizer of the proposed optimization given M_n
rsfitSolver <- function(covariate, response, treatment, estimatedNuisance, splitIndex = NULL, lossType = 'logistic', weights = NULL, tol = 10^(-3), m_0=0, constraint = TRUE, boundaryPoint=c(-1,1), numberKnots = 5){
# initialization
fit_earl <- ITRInference::ITRFitInfer(list(predictor = covariate, treatment = (treatment > 0), outcome = response), loss = 'logistic', sampleSplitIndex = splitIndex$nuisance, propensity = estimatedNuisance$pi, outcome = estimatedNuisance, test=FALSE)
#fit_init <- fitLinkLinear(covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex, lossType = lossType, weights = weights, tol = tol)
iter <- 0
# iteration start
fit_last <- NULL
tmp_beta<- (fit_earl$fit[[1]]$fit$beta[,fit_earl$fit[[1]]$fit$lambda==fit_earl$fit[[1]]$fit$lambda.min] + fit_earl$fit[[2]]$fit$beta[,fit_earl$fit[[2]]$fit$lambda==fit_earl$fit[[2]]$fit$lambda.min])/2
if (sum(tmp_beta^2)<10^-3){
fit_last$beta <- (tmp_beta+0.1)/sqrt(sum((tmp_beta+0.1)^2))
} else {
fit_last$beta <- tmp_beta/sqrt(sum(tmp_beta^2))
}
fit_last$xi <- 0
diff <- 1
step <- 0.1
while((diff > tol) & (iter < 100)){
if(any(is.na(fit_last$beta))){
print(fit_last)
print(fit_xi)
print(fit_beta)
print(iter)
}
# updateXi
fit_xi <- updateXi(fit_last, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = 1e-5, m_0=m_0, constraint=constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
# updateBeta
fit_beta <- updateBeta(fit_xi, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, withConstraint = TRUE, boundaryPoint = boundaryPoint, step=step)
# update no constraint
#fit_tmp <- updateBeta(fit_xi, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, withConstraint = FALSE, boundaryPoint = boundaryPoint)
#update
diff_coef <- max(abs(fit_beta$beta-fit_last$beta))
diff_xi <- max(c(abs(fit_beta$xi-fit_last$xi))) #abs(fit_beta$beta-fit_tmp$beta))
diff <- max(diff_coef, diff_xi)
iter <- iter+1
fit_last <- fit_beta
}
# final update
fit <- updateXi(fit_last, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, m_0=m_0, constraint=constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
fit_tmp <- updateBeta(fit, covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex,lossType = lossType, weights = weights, tol = tol, withConstraint = FALSE, boundaryPoint = boundaryPoint)
fit$solution <- fit_tmp$beta
fit$boundaryPoint <- boundaryPoint
fit
}
# rsfitSplit obtain the estimation of RSICF assuming a single index model on a slited data. The variance estimaor of the coefficiencts are provided.
rsfitSplit <- function(covariate, response, treatment, splitIndex = NULL, propensityModel = 'glmnet', estimatedPropensity = NULL, outcomeModel = 'kernel', estimatedOutcome = NULL, lossType = 'logistic', weights = NULL, tol = 1e-3, propensityFormula=NULL, outcomeFormula=NULL, constraint=TRUE, boundaryPoint=c(-1,1)){
# fit nuisance parameter
estimatedNuisance <- NULL
estimatedNuisance$pi <- estimatedPropensity
if (is.null(estimatedPropensity)){
data <- list(predictor = covariate, treatment = treatment)
predictedPropensityAll <- getPropensityModel(data, method = propensityModel, splitIndex = splitIndex, Formula = propensityFormula, predictAll = TRUE)
estimatedNuisance$pi <- predictedPropensityAll
}
estimatedNuisance$S <- estimatedOutcome$control+estimatedOutcome$treatment
estimatedNuisance$control <- estimatedOutcome$control
estimatedNuisance$treatment <- estimatedOutcome$treatment
if (is.null(estimatedOutcome)){
data <- list(predictor = covariate, treatment = treatment, outcome = response)
predictedOutcomeAll <- getOutcomeModel(data, method = outcomeModel, splitIndex = splitIndex, Formula = outcomeFormula, predictAll = TRUE)
estimatedNuisance$S <- predictedOutcomeAll$control+predictedOutcomeAll$treatment
estimatedNuisance$control <- predictedOutcomeAll$control
estimatedNuisance$treatment <- predictedOutcomeAll$treatment
}
# fit conditional variance
residual_square <- (response - (predictedOutcomeAll$control * (treatment==0) + predictedOutcomeAll$treatment * (treatment==1)))^2
data <- list(predictor = covariate, treatment = treatment, outcome = residual_square)
predictedVarAll <- getOutcomeModel(data, method = outcomeModel, splitIndex = splitIndex, Formula = outcomeFormula, predictAll = TRUE)
estimatedNuisance$var_control <- predictedVarAll$control
estimatedNuisance$var_treated <- predictedVarAll$treatment
# weights
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
# solve
m_0 <- 5
numberKnots <- 3
fit_hat <- rsfitSolver(covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex, lossType = lossType, weights =weights, tol = tol, m_0=m_0, constraint = constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
diff_min <- max(abs(fit_hat$solution-fit_hat$beta))
fit_hat_min <- fit_hat
dif_tmp <- diff_min
while((diff_min>10^(-3)) & (numberKnots <= 7) & (m_0 <= 2^5)){
early_step <- FALSE
if (sum(abs(fit_hat$xi[-1]))<m_0){
numberKnots <- numberKnots + 1
early_step <- TRUE
} else {
m_0 <- 2*m_0
}
fit_hat <- rsfitSolver(covariate, response, treatment, estimatedNuisance, splitIndex = splitIndex, lossType = lossType, weights =weights, tol = tol, m_0=m_0, constraint = constraint, boundaryPoint = boundaryPoint, numberKnots = numberKnots)
dif_tmp <- max(abs(fit_hat$solution-fit_hat$beta))
if (dif_tmp < diff_min){
fit_hat_min <- fit_hat
diff_min <- dif_tmp
}
if ((dif_tmp > diff_min)&early_step&(sum(abs(fit_hat$xi[-1]))<m_0)){
break
}
if ((dif_tmp > diff_min)&(!early_step)&(sum(abs(fit_hat$xi[-1]))>=m_0)){
break
}
}
# return
list(fit=fit_hat_min, covariate=covariate, response=response, treatment=treatment, estimatedNuisance=estimatedNuisance, splitIndex=splitIndex, lossType = lossType, weights =weights)
}
# ks gets the kernel estimation
ks <- function(xx, yy, xx.test){
nobs <- nrow(as.matrix(xx))
nvars <- ncol(as.matrix(xx))
hopt <- (4/(nvars+2))^(1/(nvars+4)) * (nobs^(-1/(nvars+4)))
wm <- function(t){
if (ncol(as.matrix(xx))==1){
weight <- exp(-0.5 * (as.numeric(t)-xx)^2/(hopt^2)) * hopt
} else {
weight <- apply(xx,1,function(s){exp(-0.5 * sum((t-s)^2)/(hopt^2)) * hopt^(ncol(xx))})
}
weighted.mean(yy, weight)
}
if (nrow(as.matrix(xx.test))==1) {
yy.test <- wm(xx.test)
} else {
if (ncol((as.matrix(xx.test)))==1){
yy.test <- sapply(as.matrix(xx.test), function(t){
wm(t)
})
} else {
yy.test <- apply(as.matrix(xx.test),1,function(t){
wm(t)
})
}
}
yy.test
}
# model.gam gets the gam regression function given data
model.gam <- function(data){
p <- dim(data$predictor)[2]
expr <- "mgcv::gam(outcome~"
for (i in 1:(p-1)){
expr <- paste0(expr, "s(predictor[,",i,"])+")
}
expr <- paste0(expr, "s(predictor[,",p,"]), data = data,method ='REML')")
expr
}
# getOutcomeModel gets the outcome regression model
getOutcomeModel <- function(data, method=c('lm', 'glmnet', 'kernel', 'others'), splitIndex, Formula = NULL, predictAll = FALSE, screeningMethod="SIRS", outcomeScreeningFamily='Gaussian'){
sampleSplitIndex <- splitIndex$nuisance
p <- dim(data$predictor)[2]
size <- dim(data$predictor)[1]
fit <- NULL
supp <- NULL
dataPredict <- NULL
dataPredict$control=data$predictor[sampleSplitIndex,]
dataPredict$treatment=data$predictor[sampleSplitIndex,]
if (predictAll){
dataPredict$control=data$predictor
dataPredict$treatment=data$predictor
}
prediction <- NULL
dataControl <- list(predictor=data$predictor[(!sampleSplitIndex) & (data$treatment==FALSE),], outcome=data$outcome[(!sampleSplitIndex) & (data$treatment==FALSE)])
dataTreatment <- list(predictor=data$predictor[(!sampleSplitIndex) & (data$treatment==TRUE),], outcome=data$outcome[(!sampleSplitIndex) & (data$treatment==TRUE)])
if ((0.05*size >= p) | (p <= 5)){
supp$control <- supp$treatment <- rep(TRUE, times = p)
if ((method != 'lm')&&(method != 'glmnet')){
ans1 <- screening(data$predictor[data$treatment==FALSE,], data$outcome[data$treatment==FALSE], method = screeningMethod, family = outcomeScreeningFamily)
ans2 <- screening(data$predictor[data$treatment==TRUE,], data$outcome[data$treatment==TRUE], method = screeningMethod, family = outcomeScreeningFamily)
if (screeningMethod == 'glmnet'){
supp$control <- ans1
supp$treatment <- ans2
} else {
supp$control <- supp$treatment <- (ans1 <= floor(p/2))|(ans2 <= floor(p/2))
}
dataControl$predictor <- dataControl$predictor[,supp$control]
dataTreatment$predictor <- dataTreatment$predictor[,supp$treatment]
}
}
if ((0.05*size < p) || (method == 'glmnet')) {
fit$control <- glmnet::cv.glmnet(x = dataControl$predictor, y = dataControl$outcome)
fit$treatment <- glmnet::cv.glmnet(x = dataTreatment$predictor, y = dataTreatment$outcome)
supp$control <- abs(fit$control$glmnet.fit$beta[,fit$control$glmnet.fit$lambda==fit$control$lambda.min])>0
supp$treatment <- abs(fit$treatment$glmnet.fit$beta[,fit$treatment$glmnet.fit$lambda==fit$treatment$lambda.min])>0
if ((method != 'lm')&&(method != 'glmnet')){
ans1 <- screening(data$predictor[data$treatment==FALSE,], data$outcome[data$treatment==FALSE], method = screeningMethod, family = outcomeScreeningFamily)
ans2 <- screening(data$predictor[data$treatment==TRUE,], data$outcome[data$treatment==TRUE], method = screeningMethod, family = outcomeScreeningFamily)
if (screeningMethod == 'glmnet'){
supp$control <- ans1
supp$treatment <- ans2
} else {
supp$control <- supp$treatment <- (ans1 <= 5)|(ans2 <= 5)
}
}
dataControl$predictor <- dataControl$predictor[,supp$control]
dataTreatment$predictor <- dataTreatment$predictor[,supp$treatment]
prediction$control <- predict(fit$control, newx = dataPredict$control, s=fit$control$lambda.min)
prediction$treatment <- predict(fit$treatment, newx = dataPredict$treatment, s=fit$treatment$lambda.min)
}
if (is.null(Formula)){
Formula <- function(support){
expr <- (outcome ~ predictor)
if (sum(support)==0){
expr <- (outcome ~ 1)
}
expr
}
}
fit <- NULL
dataPredict <- NULL
dataPredict$control=data$predictor[sampleSplitIndex,supp$control]
dataPredict$treatment=data$predictor[sampleSplitIndex,supp$treatment]
if (predictAll){
dataPredict$control=data$predictor[,supp$control]
dataPredict$treatment=data$predictor[,supp$treatment]
}
if ((method == 'lm')||(method == 'glmnet')){
if (sum(supp$control) > 0){
fit$control <- lm(Formula(supp$control), data = dataControl)
prediction$control <- predict(fit$control, newdata = list(predictor=dataPredict$control))
}
if (sum(supp$treatment) > 0){
fit$treatment <- lm(Formula(supp$treatment), data = dataTreatment)
prediction$treatment <- predict(fit$treatment, newdata = list(predictor=dataPredict$treatment))
}
} else if (method == 'kernel') {
if (sum(supp$control) > 0){
prediction$control <- ks(dataControl$predictor, dataControl$outcome, dataPredict$control)
}
if (sum(supp$treatment) > 0){
prediction$treatment <- ks(dataTreatment$predictor, dataTreatment$outcome, dataPredict$treatment)
}
} else {
if (sum(supp$control) > 0){
fit$control <- eval(parse(text=model.gam(dataControl)))
prediction$control <- predict(fit$control, newdata = list(predictor=dataPredict$control))
}
if (sum(supp$treatment) > 0){
fit$treatment <- eval(parse(text=model.gam(dataTreatment)))
prediction$treatment <- predict(fit$treatment, newdata = list(predictor=dataPredict$treatment))
}
}
prediction
}
# getPropensityModel gets the estimation of the propensity model based on selected method
getPropensityModel <- function(data, method=c('lm', 'glmnet', 'kernel'), splitIndex, Formula = NULL, predictAll = FALSE, screeningMethod="SIRS"){
sampleSplitIndex <- splitIndex$nuisance
p <- dim(data$predictor)[2]
size <- dim(data$predictor)[1]
fit <- NULL
supp <- NULL
dataPredict <- NULL
dataPredict=data$predictor[sampleSplitIndex,]
if (predictAll){
dataPredict=data$predictor
}
prediction <- NULL
dataTrain <- list(predictor=data$predictor[(!sampleSplitIndex),], treatment=data$treatment[(!sampleSplitIndex)])
if (0.05*size >= p){
supp$control <- supp$treatment <- rep(TRUE, times = p)
if ((method != 'lm')&&(method != 'glmnet')){
ans <- screening(data$predictor, data$treatment, method = screeningMethod, family = 'binomial')
if (screeningMethod == 'glmnet'){
supp <- ans
} else {
supp <- (ans <= p/2)
}
}
dataTrain$predictor = dataTrain$predictor[,supp]
}
if ((0.05*size < p) || (method == 'glmnet')) {
fit <- glmnet::cv.glmnet(x = dataTrain$predictor, y = dataTrain$treatment, family='binomial')
supp <- abs(fit$glmnet.fit$beta[,fit$glmnet.fit$lambda==fit$lambda.min])>0
if ((method != 'lm')&&(method != 'glmnet')){
ans <- screening(data$predictor, data$treatment, method = screeningMethod, family = 'binomial')
if (screeningMethod == 'glmnet'){
supp <- ans
} else {
supp <- (ans <= 5)
}
}
dataTrain$predictor <- dataTrain$predictor[,supp]
prediction <- predict(fit, newx = dataPredict, type='response', s=fit$lambda.min)
}
if (is.null(Formula)){
Formula <- function(support){
expr <- (treatment ~ predictor)
if (sum(support)==0){
expr <- (treatment ~ 1)
}
expr
}
}
fit <- NULL
dataPredict <- NULL
dataPredict=data$predictor[sampleSplitIndex,supp]
if (predictAll){
dataPredict=data$predictor[,supp]
}
if ((method == 'lm')||(method == 'glmnet')){
if (sum(supp) > 0){
fit <- glm(Formula(supp), family=binomial, data = dataTrain)
prediction <- predict(fit, newdata = list(predictor=dataPredict), type="response")
}
} else if (method == 'kernel') {
if (sum(supp) > 0){
prediction <- ks(dataTrain$predictor, dataTrain$treatment, dataPredict)
prediction <- (prediction > 0.9) * 0.9 + (prediction < 0.1) * 0.1 + (prediction < 0.9) * (prediction > 0.1) * prediction
}
}
prediction
}
# screening gets the screened variable for high dimeniosnal data
screening <- function(x, y, method='glmnet', family='Gaussian'){
var <- apply(x, 2, sd)
supp <- order(var, decreasing = TRUE)
if (method=='glmnet'){
fit <- glmnet::cv.glmnet(x, y, family = family)
coef <- fit$glmnet.fit$beta[,fit$lambda==fit$lambda.min]
supp <- (abs(coef)>0)
} else {
fit <- VariableScreening::screenIID(x, y, method=method)
supp <- fit$rank
}
supp
}
# loss reutrn the value or the derivative of the chosen type of link
loss <- function(input, type = c('logistic', 'exponential'), order = 'original'){
return(switch(order, 'original'=switch(type, 'logistic'=exp(-input)/2, 'exponential'=exp(-2*input)/4),
'derivative'=switch(type, 'logistic'=-exp(-input)/2, 'exponential'=-exp(-2*input)/2),
'hessian'=switch(type, 'logistic'=exp(-input)/2, 'exponential'=exp(-2*input))))
}
# link returns the eta from eta_trans
link <- function(eta_trans, lossType='logistic'){
return(switch(lossType, 'logistic'=(exp(eta_trans/2)-exp(-eta_trans/2))/(exp(eta_trans/2)+exp(-eta_trans/2)),
'exponential'=(exp(eta_trans)-exp(-eta_trans))/(exp(eta_trans)+exp(-eta_trans))))
}
# nsd get the derivative of the ns
nsd <- function(x, knots = NULL, intercept = FALSE, Boundary.knots = range(x)){
div <- sapply(x,function(t){
if((t+1e-6)>=Boundary.knots[2]){
res <- (splines2::bSpline(t, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots)-splines2::bSpline(t-1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots))/(1e-6)
} else if ((t-1e-6)<=Boundary.knots[1]){
res <- (splines2::bSpline(t+1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots)-splines2::bSpline(t, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots))/(1e-6)
} else {
res <- (splines2::bSpline(t+1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots)-splines2::bSpline(t-1e-6, knots=knots, intercept=intercept, Boundary.knots =Boundary.knots))/(2*1e-6)
}
res
})
t(div)
}
# fitLinkLinear impliments the Step 1 of the porposed algorithm.
# stimatedNuissance should contain two lists. One for propensity, the other for the estimated S.
# The predictions nad inputs is always for all the data.
# SplitIndex containes three list. $nuisance, $fit ,$efficient
# treatment is 0 or 1
# covariate does not contain intercept
fitLinkLinear <- function(covariate, response, treatment, estimatedNuisance, splitIndex = NULL, lossType = 'exponential', weights = NULL, tol = 1e-3, intercept = TRUE){
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
pi <- estimatedNuisance$pi[splitIndex$fit]
pi_T <- pi * treatment[splitIndex$fit] + (1-pi) * (1-treatment[splitIndex$fit])
S <- estimatedNuisance$S[splitIndex$fit]
trt_train <- 2*(treatment[splitIndex$fit]-0.5)
x_train <- covariate[splitIndex$fit,]
y_train <- response[splitIndex$fit]
w_train <- weights[splitIndex$fit]
if(intercept){
x_train <- cbind(1, covariate[splitIndex$fit,])
}
obj <- function(beta){
value <- weighted.mean(y_train/pi_T*loss(trt_train * x_train %*% beta, type = lossType)+trt_train/(2*pi_T)*(S-y_train)*x_train %*% beta, w_train)
grad <- (apply(apply(x_train,2,function(z){(1/pi_T*y_train*trt_train* loss(trt_train * x_train %*% beta, type = lossType, order = 'derivative')+trt_train/(2*pi_T)*(S-y_train))*z}),
2, function(t){weighted.mean(t, w_train)}))
list(objective = value, gradient = grad)
}
beta0 <- array(0,c(NCOL(x_train),1))
fit <- nloptr::nloptr(x0=beta0, eval_f = obj, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
if(intercept){
fit$solution <- fit$solution[-1]
}
fit$solution <- fit$solution/sqrt(sum(fit$solution^2))
fit <- list(beta=fit$solution, xi=NULL, knods = NULL, solution=NULL)
fit
}
# updateXi implements the Step 2 of the proposed algorithm
updateXi <- function(fit, covariate, response, treatment, estimatedNuisance, lossType='tanh', splitIndex = NULL, weights = NULL, tol = 1e-3, numberKnots = 5, m_0 =NULL, constraint = TRUE, boundaryPoint=c(-1,1)){
beta_last <- fit$beta
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
pi <- estimatedNuisance$pi[splitIndex$fit]
pi_T <- pi * treatment[splitIndex$fit] + (1-pi) * (1-treatment[splitIndex$fit])
S <- estimatedNuisance$S[splitIndex$fit]
trt_train <- 2*(treatment[splitIndex$fit]-0.5)
x_train <- covariate[splitIndex$fit,]
y_train <- response[splitIndex$fit]
w_train <- weights[splitIndex$fit]
z_train <- x_train %*% beta_last
knots <- quantile(z_train, probs=seq(0,1,length.out = numberKnots))
knots[1] <- knots[1]-0.3*abs(knots[2]-knots[1])
knots[numberKnots] <- knots[numberKnots]+0.3*abs(knots[numberKnots]-knots[numberKnots-1])
ns_train <- splines2::bSpline(z_train, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint)
ns_train <- cbind(1, ns_train)
# tuning start
obj <- function(xi){
value <- weighted.mean(y_train/pi_T*loss(trt_train * ns_train %*% xi, type = lossType)+trt_train/(2*pi_T)*(S-y_train)*ns_train %*% xi, w_train)
grad <- (apply(apply(ns_train,2,function(z){(1/pi_T*y_train*trt_train* loss(trt_train * ns_train %*% xi, type = lossType, order = 'derivative')+trt_train/(2*pi_T)*(S-y_train))*z}),
2, function(t){weighted.mean(t, w_train)}))
list(objective = value, gradient = grad)
}
#form contrast matrix
B <- matrix(0,NCOL(ns_train)-1, NCOL(ns_train))
for(i in 2:(NCOL(ns_train)-1)){
B[i,i] <- 1
B[i, i+1] <- -1
}
g_ineq <- function(xi){
value <- c(sum(abs(xi[-1]))-m_0, c(B%*% xi))
grad <- rbind(c(0,sign(xi[-1])), B)
list(constraints = value, jacobian = grad)
}
if(!constraint){
g_ineq <- function(xi){
value <- c(sum(abs(xi[-1]))-m_0)
grad <- c(0,sign(xi[-1]))
list(constraints = value, jacobian = grad)
}
}
fit_xi <- nloptr::nloptr(x0=rep(0, times=NCOL(ns_train)), eval_f = obj, eval_g_ineq = g_ineq, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
z_test <- seq(min(knots),max(knots), length.out=100)
ns_test <- splines2::bSpline(z_test, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint)
ns_test <- cbind(1, ns_test)
fit_test <- ns_test %*% fit_xi$solution
plot(z_test, fit_test)
fit <- list(beta=beta_last, xi=fit_xi$solution, knots = knots,solution=fit$solution)
fit
}
# updateBeta implements the Step 3 of the proposed algorithm
updateBeta <- function(fit, covariate, response, treatment, estimatedNuisance, lossType='tanh', splitIndex = NULL, weights = NULL, tol = 1e-3, withConstraint = TRUE, boundaryPoint = c(-1,1), step=0.1){
beta_last <- fit$beta
xi <- fit$xi
if(is.null(weights)){
weights <- rep(1, times=NROW(covariate))
}
pi <- estimatedNuisance$pi[splitIndex$fit]
pi_T <- pi * treatment[splitIndex$fit] + (1-pi) * (1-treatment[splitIndex$fit])
S <- estimatedNuisance$S[splitIndex$fit]
trt_train <- 2*(treatment[splitIndex$fit]-0.5)
x_train <- covariate[splitIndex$fit,]
y_train <- response[splitIndex$fit]
w_train <- weights[splitIndex$fit]
knots <- fit$knots
# objective function
obj <- function(beta){
z_train <- x_train %*% beta
ns_train <- splines2::bSpline(z_train, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint)
ns_train <- cbind(1, ns_train)
nsd_train <- cbind(0,nsd(z_train, knots = knots, intercept = FALSE, Boundary.knots = boundaryPoint))
value <- weighted.mean(y_train/pi_T*loss(trt_train * ns_train %*% xi, type = lossType)+trt_train/(2*pi_T)*(S-y_train)*ns_train %*% xi, w_train)
grad <- (apply(apply(x_train,2,function(z){(1/pi_T*y_train*trt_train* loss(trt_train * ns_train %*% xi, type = lossType, order = 'derivative')+trt_train/(2*pi_T)*(S-y_train))*nsd_train%*% xi*z}),
2, function(t){weighted.mean(t, w_train, na.rm = TRUE)}))
list(objective = value, gradient = grad)
}
constraint <- function(beta){
return(list('constraints'=sum(beta^2)-1,
'jacobian'=2*beta))
}
if (withConstraint){
fit_beta <- nloptr::nloptr(x0=fit$beta, eval_f = obj, eval_g_eq = constraint, lb=beta_last-step, ub=beta_last+step, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
} else {
fit_beta <- nloptr::nloptr(x0=fit$beta, eval_f = obj, opts = list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=tol))
}
fit <- list(beta=fit_beta$solution, xi=fit$xi, knots = fit$knots)
fit
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.