R/doseFinding_quadprog.R
In lixiaomao/personalized-dosing-interval: Determining Personalized Dosing Interval
#library(CVXR)
#library(LowRankQP)
#library(quadprog)
#' @export
doseFind<-function(train, test=NULL,eps=NULL,aL=NULL,aU=NULL,step=1,thres=NULL,maxiter=20,lambda=0){
if (is.null(thres)) thres=quantile(train$Y,0.8)
if (is.null(aL)) aL=min(train$A)-0.5*sd(train$A)
if (is.null(aU)) aU=max(train$A)+0.5*sd(train$A)
p=dim(train$X)[2]
n=dim(train$X)[1]
index=train$Y>=thres
beta <- Variable(p)
alpha<-Variable(1)
obj <- Minimize(sum_entries(abs(train$A[index] -alpha- train$X[index,] %*% beta)))
prob<- Problem(obj)
solv <- solve(prob)
beta0= solv$getValue(beta)
alpha0=solv$getValue(alpha)
pred0=alpha0+train$X%*% beta0
for (i in 1:maxiter){
solv=fit.dose(train,pred0=pred0,eps=eps)
if (sum((solv$beta-beta0)^2)<1e-5){
cat('Converged at iteration: ',i,'-th\n')
break;
} else {
beta0=solv$beta
alpha0=solv$alpha
}
pred0=alpha0+train$X%*% beta0
dosedif=pred0-train$A
index=(abs(dosedif)<0.1)
cat('iteration:',i,' sample size:',sum(index),'\n')
print(quantile(pred0))
}
pred_dose=predict.doseCVX(test$X,solv)
pred_dose=pmin(pmax(pred_dose,3),12)
return(list(pred=pred_dose,solv=solv,beta=solv$beta,alpha=solv$alpha))
}
#' @export
fit.dose<-function(train,pred0=NULL,beta0=NULL,alpha0=NULL,eps,aL=NULL,aU=NULL){
if (is.null(pred0)) pred0=alpha0+train$X%*% beta0
p=dim(train$X)[2]
Q=(abs(train$A-pred0)<eps)
Q1=train$A<pred0-eps
Q2=train$A>pred0+eps
beta <- Variable(p)
alpha<-Variable(1)
pred=alpha+train$X %*% beta
res=train$A - pred
part1=sum_entries(Q*abs(res))
part2=sum_entries((Q1)*2*pos(train$A-pred))
part3=sum_entries((Q2)*2*pos(pred-train$A))
partAll=sum_entries(1/train$propensity*train$Y*(part1+part2+part3))
obj <- Minimize(partAll)
prob<- Problem(obj)
solv <- solve(prob)
class(solv)='doseCVX'
solv$beta=solv$getValue(beta)
solv$alpha=solv$getValue(alpha)
return(solv)
}
#' @export
predict.doseCVX<-function(X,solv){
return(solv$alpha+X%*%solv$beta)
}
##############################################################################
## one-sided dosing interval
##############################################################################
#' @export
lineOne<-function(b,s1,s2,s3,s4,a,fit){
res=a-fit-b
return(sum(s1*pmax(-res,0)+s2*pmax(res,0)))
}
#' @export
fit.Int1<-function(train,eps,XX=NULL,CC=NULL,pred0=NULL,rbfkernel=NULL,
Intcep0=NULL,beta0=NULL,lambda=0.0001,type='PDI',lower=TRUE,method='completeLinear',...){
# if (FALSE){
# alpha=c(0.82,0.18)
# eps=0.3
# lambda=0.0001
# n=length(train$A)
# pred0=6
#
# }
n=length(train$A)
p=dim(train$X)[2]
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (is.null(CC)){
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = sigma)
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
CC=cbind(rbind(XX,-XX),rbind(-XX,XX))
}
d=matrix(c(-train$A,train$A),ncol = 1)
Aeq=matrix(c(rep(1,n),rep(-1,n)),ncol = 1)
if (type=='PDI'){
H1=(train$Y>train$S)*(train$alpha[,2])*train$propensity/lambda
H2=(train$Y<train$S)*(train$alpha[,1])*train$propensity/lambda
} else{
H1=pmax(train$Y-train$S,0)*(train$alpha[,2])*train$propensity/lambda
H2=pmax(train$S-train$Y,0)*(train$alpha[,1])*train$propensity/lambda
}
Q1=train$A<(pred0-eps)
Q2=train$A>(pred0+eps)
W1=train$A>(pred0-eps)
W2=train$A<(pred0+eps)
#Q1=Q1*0
#Q2=Q2*0
if (lower){
ub=c(H1*W1+H2*Q2,H1*Q1+H2*W2)
} else{
ub=c(H2*W1+H1*Q2,H2*Q1+H1*W2)
}
zero=(ub==0)
ubMat=diag(2*n)
Amat <- cbind(Aeq[!zero,],-ubMat[!zero,!zero],ubMat[!zero,!zero])
bvec=matrix(c(0,-ub[!zero],rep(0,sum(!zero))),ncol=1)
solv=solve.QP(Dmat=CC[!zero,!zero]+diag(0.1,sum(!zero)), dvec=d[!zero,], Amat=Amat, bvec=bvec, meq=1, factorized=FALSE)
alfa=rep(0,2*n)
alfa[!zero]=solv$solution
if (method=='completeLinear'){
beta=t((alfa[-c(1:n)]-alfa[1:n])%*%X)
pred=X%*%beta
} else {
beta=(alfa[-c(1:n)]-alfa[1:n])
pred=XX%*%beta
}
opt=optimize(interval = c(min(train$A),max(train$A)),f=lineOne,s1=H1*W1+H2*Q2,s2=H1*Q1+H2*W2,a=train$A,fit=pred)
intcept=opt$minimum
output=list(beta=beta,X=X,intcept=intcept,pred=pred+intcept,solv=solv,method=method,rbfkernel=rbfkernel)
class(output)='Quad1'
return(output)
}
#' @export
predict.Quad1<-function(fit,X){
X=apply(X,2,function(x) (x-mean(x))/sd(x))
if (fit$method=='completeLinear'){
pred=X%*%fit$beta+fit$intcept
} else {
XX=as.matrix(as.data.frame(kernelMatrix(fit$rbfkernel,X,fit$X)))
pred=XX%*%fit$beta+fit$intcept
}
return(pred)
}
#' @export
predict.RDCDI1<-function(obj,test){
fit.model=obj$fit.model
lower=obj$lower
type=obj$type
n1=dim(test$X)[1]
pred0=obj$pred0
misclass<-NULL->region->pred
if (isTRUE(all.equal(fit.model,NA))){
fit.model=pred0$fit
pred=predict.init(fit.model,test)
} else{
pred=predict.Quad1(fit.model,test$X)
}
if (!is.null(test$A)){
if (lower){
region=pred<test$A
} else {
region=pred>test$A
}
if ((!is.null(test$S))&(!is.null(test$Y))){
if (type=="EDI"){
misclass=(sum(test$alpha[,1]*(region)*pmax((test$S-test$Y),0))+sum(test$alpha[,2]*(!region)*pmax((test$Y-test$S),0)))/n1
} else if (type=="PDI"){
misclass=(sum(test$alpha[,1]*(region)*(test$Y<test$S))+sum(test$alpha[,2]*(!region)*(test$Y>test$S)))/n1
}
}
}
return(list(misclass=misclass,region=region,pred=pred))
}
#' @export
RDCDI1<-function(train,type='PDI',method=NULL,pred0=NULL,
maxiter=25,lower=TRUE,eps=NULL,aL=NULL,aU=NULL,
lambda=0.0001,...){
n=length(train$A)
p=dim(train$X)[2]
sigma=1/dim(train$X)[2]
if (is.null(aL)) aL<-min(train$A)-0.5*sd(train$A)
if (is.null(aU)) aU<-max(train$A)+0.5*sd(train$A)
if (is.null(eps)) eps=0.1*(aU-aL)
if (is.null(pred0)){
if (length(table(train$A))>4){
init0=initiate.continuous(train,train,side=1)
pred0=init0$pred
Intcep0=drop(coef(init0$fit))[1]
beta0=drop(coef(init0$fit))[-1]
} else {
init0=initiate.binary(train,train,side=1)
pred0=init0$pred
Intcep0=drop(coef(init0$fit))[1]
beta0=drop(coef(init0$fit))[-1]
}
} else {
Intcep0=median(pred0)
if (length(pred0)==1) pred0=rep(pred0,n)
if (method=='completeLinear'){
beta0=rep(0,p)
} else {
beta0=rep(0,n)
}
}
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = sigma)
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
CC=cbind(rbind(XX,-XX),rbind(-XX,XX))
for (i in 1:maxiter){
before<-Sys.time()
fit=fit.Int1(train,XX=XX,CC=CC,eps=eps,pred0=pred0,lambda=lambda,type=type,lower=lower,method=method,rbfkernel=rbfkernel)
after<-Sys.time()
after-before
fit0=fit
Intcep0=fit$intcept
pred0=drop(predict(fit,train$X))
if (sum(abs(fit$beta-beta0))<1e-15){
beta0=fit$beta
cat('Converged at iteration: ',i-1,'-th\n')
break;
}
beta0=fit$beta
print(quantile(pred0))
}
#1-sum((pred0-train$opt)^2)/sum(train$opt^2)
if (lower){
region=pred0<train$A
} else{
region=pred0>train$A
}
if (type=="EDI"){
if (lower){
misclass=(sum(train$alpha[,1]*(region)*pmax((train$S-train$Y),0))+sum(train$alpha[,2]*(!region)*pmax((train$Y-train$S),0)))/n
} else{
misclass=(sum(train$alpha[,1]*(!regopm)*pmax((train$S-train$Y),0))+sum(train$alpha[,2]*(region)*pmax((train$Y-train$S),0)))/n
}
} else if (type=="PDI"){
if (lower){
misclass=(sum(train$alpha[,1]*(region)*(train$Y<train$S))+sum(train$alpha[,2]*(!region)*(train$Y>train$S)))/n
} else{
misclass=(sum(train$alpha[,1]*(!region)*(train$Y<train$S))+sum(train$alpha[,2]*(region)*(train$Y>train$S)))/n
}
}
output=list(beta=fit$beta,intcept=fit$intcept,pred=pred0,solv=fit$solv,region=region,type=type,lower=lower,numiter=i,lambda=lambda,fit.model=fit,misclass=misclass)
class(output)='RDCDI1'
return(output)
}
#
# fit_L=RDCDI1(train=train,alpha=c(0.82,0.18),type='PDI',pred0=6,
# maxiter=5,lower=TRUE,eps=0.2,aL=NULL,aU=NULL,
# lambda=0.000001)
#
# pred_L=predict(fit_L,test$X)
#
# fit_R=RDCDI1(train=train,alpha=c(0.82,0.18),type='PDI',pred0=8,
# maxiter=5,lower=FALSE,eps=0.2,aL=NULL,aU=NULL,
# lambda=0.000001)
#
#
# pred_R=predict(fit_R,test$X)
#
# plotdata=data.frame(observed_A1c=test$A,recommended_A1c=pred_dose,
# lowerBound=pred_L,upperBound=pred_R)
# colnames(plotdata)[1]='observed_A1c'
# plotdata<- melt(plotdata)
# p<-ggplot(plotdata,aes(x=value, fill=variable,colour=variable)) +
# scale_x_continuous(breaks = sort(c(seq(2,12,0.5))))+
# geom_density(alpha=0.25)+
# xlab('A1c')
# p
#####################################################################
## two-sided version
#####################################################################
# lineTwo<-function(b,s1,s2,s3,s4,s5,s6,s7,s8,a,fitL,fitR,eps){
# resL1=a-fitL-b[1]
# resL2=a-fitL-b[1]+eps
# resR1=a-fitR-b[2]
# resR2=a-fitR-b[2]-eps
# return(sum(s1*pmax(-resL1,0)+s2*pmax(resL1,0)+s3*pmax(-resL2,0)+s4*pmax(resL2,0)+s5*pmax(resR1,0)+s6*pmax(-resR1,0)+s7*pmax(resR2,0)+s8*pmax(-resR2,0)))
# }
# fit.Int2<-function(train,eps,lambda=0.0001,type='PDI',CC=NULL,XX=NULL,lower=TRUE,
# pred0L=NULL,Intcep0L=NULL,beta0L=NULL,pred0R=NULL,Intcep0R=NULL,beta0R=NULL,...){
# if (is.null(pred0L)) pred0L=Intcep0L+train$X %*% beta0L
# if (is.null(pred0R)) pred0R=Intcep0R+train$X %*% beta0R
# # library(quadprog)
# # train=subsetData(train,sample(1:length(train$A),1000))
# # alpha=c(0.82,0.18)
# # eps=0.3
# # lambda=0.000001
# # pred0L=6.5
# # pred0R=8
# n=length(train$A)
# p=dim(train$X)[2]
#
# if (is.null(X)) X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
# if (is.null(CC)){
# XX=X%*%t(X)
# col1=rbind(XX,-XX,XX,-XX,XX,matrix(0,nrow=4*n,ncol=n))
# col2=rbind(matrix(0,nrow=4*n,ncol=n),XX,-XX,XX,-XX,XX)
# CC=cbind(col1,-col1,col1,-col1,col1+col2,-col2,col2,-col2,col2)
# }
#
# Q1L=train$A<(pred0L-eps)
# Q2L=train$A<(pred0L)
# W1L=train$A>(pred0L-eps)
# W2L=train$A>(pred0L)
#
# Q1R=train$A>(pred0R+eps)
# Q2R=train$A>(pred0R)
# W1R=train$A<(pred0R+eps)
# W2R=train$A<(pred0R)
#
# #Q1L=Q1L*0;Q2L=Q2L*0;Q1R=Q1R*0;Q2R=Q2R*0
# if (type=='PDI'){
# H1=(train$Y>train$S)*(train$alpha[,2])*train$propensity/lambda
# H2=(train$Y<train$S)*(train$alpha[,1])*train$propensity/lambda
# } else {
# H1=pmax(train$Y-train$S,0)*(train$alpha[,2])*train$propensity/lambda
# H2=pmax(train$S-train$Y,0)*(train$alpha[,1])*train$propensity/lambda
# }
#
#
#
# ub=c(H1*W1L,H1*Q1L,H2*W2L,H2*Q2L,rep(100000,n),H1*W1R,H1*Q1R,H2*W2R,H2*Q2R)
# zero0<-(ub==0)->zero1
# zero1[(4*n+1):(5*n)]=TRUE
#
#
#
# d=matrix(c(-train$A,train$A,-train$A-rep(eps,n),train$A+rep(eps,n),rep(0,n),-train$A,train$A,-train$A-rep(eps,n),train$A+rep(eps,n)),ncol = 1)
# Aeq=cbind(c(rep(1,n),rep(-1,n),rep(1,n),rep(-1,n),rep(1,n),rep(0,4*n)),c(rep(0,4*n),rep(-1,n),rep(1,n),rep(-1,n),rep(1,n),rep(-1,n)))
#
# ubMat=diag(9*n)
#
# Amat <- cbind(Aeq[!zero0,],-ubMat[!zero0,!zero1],ubMat[!zero0,!zero0])
# bvec=matrix(c(rep(0,2),-ub[!zero1],rep(0,sum(!zero0))),ncol=1)
# solv=solve.QP(Dmat=CC[!zero0,!zero0]+diag(0.00000001,sum(!zero0)), dvec=d[!zero0,], Amat=Amat, bvec=bvec, meq=2, factorized=FALSE)
# alfa=rep(0,9*n)
# alfa[!zero0]=solv$solution
# betaL=t((alfa[(n+1):(2*n)]-alfa[(1:n)]+alfa[(3*n+1):(4*n)]-alfa[(2*n+1):(3*n)]-alfa[(4*n+1):(5*n)])%*%X)
# betaR=t((alfa[(6*n+1):(7*n)]-alfa[(5*n+1):(6*n)]+alfa[(8*n+1):(9*n)]-alfa[(7*n+1):(8*n)]+alfa[(4*n+1):(5*n)])%*%X)
#
# #lineTwo(c(6,8),s1=H1*W1L+H2*Q2L,s2=H1*Q1L+H2*W2L,s3=H1*W1R+H2*Q2R,s4=H1*Q1R+H2*W2R,a=train$A,fitL=X%*%betaL,fitR=X%*%betaR)
# fit=optim(par=c(quantile(train$A,0.2),quantile(train$A,0.8)),fn=lineTwo,eps=eps,s1=H1*W1L,s2=H1*Q1L,s3=H2*W2L,s4=H2*Q2L,s5=H1*W1R,s6=H1*Q1R,s7=H2*W2R,s8=H2*Q2R,a=train$A,fitL=X%*%betaL,fitR=X%*%betaR)
#
#
# # mu=ub[-c((2*n+1):(3*n))]-alfa[-c((2*n+1):(3*n))]
# # aa=c(train$A-X%*%betaL,X%*%betaL-train$A,train$A-X%*%betaR,X%*%betaR-train$A)
# # aa[which((mu!=0)&(abs(alfa)[-c((2*n+1):(3*n))]>0.0000001))]
#
# output=list()
# output$betaL=betaL
# output$betaR=betaR
# output$intcepL=fit$par[1]
# output$intcepR=fit$par[2]
# output$pred_L=X%*%betaL+fit$par[1]
# output$pred_R=X%*%betaR+fit$par[2]
# output$fit=fit
# #output$solv=solv
# class(output)='Quad2'
# return(output)
# }
#' @export
lineTwo<-function(b,s1,s2,s3,s4,a,fitL,fitR){
resL=a-fitL-b[1]
resR=a-fitR-b[2]
return(sum(s1*pmax(-resL,0)+s2*pmax(resL,0)+s3*pmax(-resR,0)+s4*pmax(resR,0)+10*pmax(fitL+b[1]-fitR-b[2],0)))
}
#' @export
fit.Int2<-function(train,eps,lambda=0.0001,type='PDI',CC=NULL,XX=NULL,method='completeLinear',rbfkernel=NULL,
pred0L=NULL,Intcep0L=NULL,beta0L=NULL,pred0R=NULL,Intcep0R=NULL,beta0R=NULL,...){
if (is.null(pred0L)) pred0L=Intcep0L+train$X %*% beta0L
if (is.null(pred0R)) pred0R=Intcep0R+train$X %*% beta0R
#library(quadprog)
#train=subsetData(train,sample(1:length(train$A),1000))
#alpha=c(0.82,0.18)
#eps=0.3
#lambda=0.000001
n=length(train$A)
sigma=1/dim(train$X)[2]
#pred0L=6.5
#pred0R=8
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (is.null(CC)){
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = sigma)
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
col1=rbind(XX,-XX,XX,matrix(0,nrow=2*n,ncol=n))
col2=rbind(matrix(0,nrow=2*n,ncol=n),XX,-XX,XX)
CC=cbind(col1,-col1,col1+col2,-col2,col2)
}
Q1L=train$A<(pred0L-eps)
Q2L=train$A<(pred0L+eps)
W1L=train$A>(pred0L-eps)
W2L=train$A>(pred0L+eps)
Q1R=train$A>(pred0R+eps)
Q2R=train$A>(pred0R-eps)
W1R=train$A<(pred0R+eps)
W2R=train$A<(pred0R-eps)
#Q1L=Q1L*0;Q2L=Q2L*0;Q1R=Q1R*0;Q2R=Q2R*0
if (type=='PDI'){
H1=(train$Y>train$S)*(train$alpha[,2])*train$propensity/lambda
H2=(train$Y<train$S)*(train$alpha[,1])*train$propensity/lambda
} else{
H1=pmax(train$Y-train$S,0)*(train$alpha[,2])*train$propensity/lambda
H2=pmax(train$S-train$Y,0)*(train$alpha[,1])*train$propensity/lambda
}
ub=c(H1*W1L+H2*W2L,H1*Q1L+H2*Q2L,rep(1e05,n),H1*Q1R+H2*Q2R,H1*W1R+H2*W2R)
zero0<-(ub==0)->zero1
zero1[(2*n+1):(3*n)]=TRUE
d=matrix(c(-train$A,train$A,rep(2*eps,n),-train$A,train$A),ncol = 1)
Aeq=cbind(c(rep(1,n),rep(-1,n),rep(1,n),rep(0,2*n)),c(rep(0,2*n),rep(-1,n),rep(1,n),rep(-1,n)))
ubMat=diag(5*n)
Amat <- cbind(Aeq[!zero0,],-ubMat[!zero0,!zero1],ubMat[!zero0,!zero0])
bvec=matrix(c(rep(0,2),-ub[!zero1],rep(0,sum(!zero0))),ncol=1)
solv=solve.QP(Dmat=CC[!zero0,!zero0]+diag(0.1,sum(!zero0)), dvec=d[!zero0,], Amat=Amat, bvec=bvec, meq=2, factorized=FALSE)
alfa=rep(0,5*n)
alfa[!zero0]=solv$solution
if (method=='completeLinear'){
betaL=t((alfa[(n+1):(2*n)]-alfa[(1:n)]-alfa[(2*n+1):(3*n)])%*%X)
betaR=t((alfa[(4*n+1):(5*n)]-alfa[(3*n+1):(4*n)]+alfa[(2*n+1):(3*n)])%*%X)
pred_L=X%*%betaL
pred_R=X%*%betaR
} else {
betaL=(alfa[(n+1):(2*n)]-alfa[(1:n)]-alfa[(2*n+1):(3*n)])
betaR=(alfa[(4*n+1):(5*n)]-alfa[(3*n+1):(4*n)]+alfa[(2*n+1):(3*n)])
pred_L=XX%*%betaL
pred_R=XX%*%betaR
}
#lineTwo(c(6,8),s1=H1*W1L+H2*Q2L,s2=H1*Q1L+H2*W2L,s3=H1*W1R+H2*Q2R,s4=H1*Q1R+H2*W2R,a=train$A,fitL=X%*%betaL,fitR=X%*%betaR)
fit=optim(par=c(quantile(train$A,0.2),quantile(train$A,0.8)),fn=lineTwo,s1=H1*W1L+H2*W2L,s2=H1*Q1L+H2*Q2L,s3=H1*Q1R+H2*Q2R,s4=H1*W1R+H2*W2R,a=train$A,fitL=pred_L,fitR=pred_R)
# mu1L=ub[c((1):(n))]-alfa[c((1):(n))]
# mu2L=ub[c((n+1):(2*n))]-alfa[c((n+1):(2*n))]
# index1L=!zero0[c((1):(n))]
# bLub=(-X%*%betaL-train$A)[index1L&(abs(mu1L)<0.00000001)]
# index2L=!zero0[c((1+n):(2*n))]
# bLlb=(train$A-X%*%betaL)[index2L&(abs(mu2L)<0.00000001)]
#
# ##
# ## # this is supposed to be the right one
# bLlb=(train$A-X%*%betaL)[((alfa[c((1):(n))]>1e5)&index1L)|((mu2L>1e5)&index2L)]
# quantile(bLlb)
# bLub=(X%*%betaL-train$A)[((alfa[c((n+1):(2*n))]>1e5)&index2L)|((mu1L>1e5)&index1L)]
# quantile(bLub)
#
# mu1L=ub[c((1):(n))]-alfa[c((1):(n))]
# mu2L=ub[c((n+1):(2*n))]-alfa[c((n+1):(2*n))]
# index1L=!zero0[c((1):(n))]
# bLub=(-X%*%betaL-train$A)[index1L&abs(mu2L)==0]
# bLub=(-X%*%betaL-train$A)[index1L&(abs(alfa[c((1):(n))])<1e-10)]
# bLub=(-X%*%betaL-train$A)[index1L&(alfa[c((1):(n))]>1e-2)&(mu1L>1e-2)]
#
# index2L=!zero0[c((1+n):(2*n))]
# bLlb=(train$A-X%*%betaL)[index2L&abs(mu2L)==0]
# index2L=!zero0[c((1+n):(2*n))]
# bLlb=(train$A-X%*%betaL)[index2L&(alfa[c((n+1):(2*n))]>1e-2)&(mu2L>1e-2)]
# bLlb=(train$A-X%*%betaL)[index2L&(alfa[c((n+1):(2*n))]>1e-2)&(mu2L>1e-2)]
#
# index=(abs(mu)<0.00001)&(ub[-c((2*n+1):(3*n))]!=0)
# aa=c(train$A-X%*%betaL,X%*%betaL-train$A,train$A-X%*%betaR,X%*%betaR-train$A)
# aa[index]
# aa[which((mu!=0)&(abs(alfa)[-c((2*n+1):(3*n))]>0.0000001))]
output=list()
output$betaL=betaL
output$betaR=betaR
output$solv=solv
output$intcepL=fit$par[1]
output$intcepR=fit$par[2]
output$pred_L=pred_L+fit$par[1]
output$pred_R=pred_R+fit$par[2]
output$rbfkernel=rbfkernel
output$method=method
output$X=X
#output$solv=solv
class(output)='Quad2'
return(output)
}
#' @export
RDCDI2<-function(train,method='completeLinear',type='PDI',pred0L=NULL,pred0R=NULL,
maxiter=25,eps=NULL,aL=NULL,aU=NULL,lower=TRUE,
lambda=0.0001,...){
n=length(train$A)
p=dim(train$X)[2]
rbfkernel=NULL
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = 1/dim(train$X)[2])
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
col1=rbind(XX,-XX,XX,matrix(0,nrow=2*n,ncol=n))
col2=rbind(matrix(0,nrow=2*n,ncol=n),XX,-XX,XX)
CC=cbind(col1,-col1,col1+col2,-col2,col2)
if (is.null(aL)) aL<-min(train$A)-0.5*sd(train$A)
if (is.null(aU)) aU<-max(train$A)+0.5*sd(train$A)
if (is.null(eps)) eps=0.2*sd(train$A)
if (is.null(pred0L)|is.null(pred0R)){
if (length(table(train$A))>4){
init0=initiate.continuous(train,train,side=2)
} else {
init0=initiate.binary(train,train,side=2)
}
pred0L=fit$pred_L
pred0R=fit$pred_R
Intcep0L=drop(coef(init0$fit$fit_L))[1]
beta0L=drop(coef(init0$fit$fit_L))[-1]
Intcep0R=drop(coef(init0$fit$fit_R))[1]
beta0R=drop(coef(init0$fit$fit_R))[-1]
} else {
if (method=='completeLinear'){
beta0L=rep(0,p)
beta0R=rep(0,p)
} else {
beta0L=rep(0,n)
beta0R=rep(0,n)
}
Intcep0L=median(pred0L)
Intcep0R=median(pred0R)
if (length(pred0L)==1) pred0L=rep(pred0L,n)
if (length(pred0R)==1) pred0R=rep(pred0R,n)
}
#pred0L=rep(-0.5,n);pred0R=rep(0.5,n)
for (i in 1:maxiter){
before=Sys.time()
fit=fit.Int2(train,eps=eps,lambda=lambda,type=type,pred0L=pred0L,pred0R=pred0R,rbfkernel=rbfkernel,CC=CC,XX=XX,method=method)
after=Sys.time()
print(after-before)
pred0L=fit$pred_L
pred0R=fit$pred_R
if (sum(abs(fit$betaL-beta0L)+abs(fit$betaR-beta0R))<1e-15){
cat('Converged at iteration: ',i-1,'-th\n')
beta0L=fit$betaL
beta0R=fit$betaR
break;
}
beta0L=fit$betaL
beta0R=fit$betaR
cat('two-sided iteration: ',i,'\n')
print(quantile(pred0L))
print(quantile(pred0R))
}
region=(pred0L<train$A)&(pred0R>train$A)
if (type=="EDI"){
misclass=(sum(train$alpha[,1]*(region)*pmax((train$S-train$Y),0))+sum(train$alpha[,2]*(!region)*pmax((train$Y-train$S),0)))/n
} else if (type=="PDI"){
misclass=(sum(train$alpha[,1]*(region)*(train$Y<train$S))+sum(train$alpha[,2]*(!region)*(train$Y>train$S)))/n
}
output=list(method=method,betaL=fit$betaL,betaR=fit$betaR,intcepL=fit$intcepL,intcepR=fit$intcepR,pred_L=fit$pred_L,pred_R=fit$pred_R,region=region,type=type,lower=lower,numiter=i,lambda=lambda,fit.model=fit,misclass=misclass)
class(output)='RDCDI2'
return(output)
}
#' @export
predict.Quad2<-function(fit,X){
X=apply(X,2,function(x) (x-mean(x))/sd(x))
if (fit$method=='completeLinear'){
pred_L=fit$intcepL+X %*% fit$betaL
pred_R=fit$intcepR+X %*% fit$betaR
} else {
XX=as.matrix(as.data.frame(kernelMatrix(fit$rbfkernel,X,fit$X)))
pred_L=fit$intcepL+XX %*% fit$betaL
pred_R=fit$intcepR+XX %*% fit$betaR
}
return(list(pred_L=pred_L,pred_R=pred_R))
}
#' @export
predict.RDCDI2<-function(obj,test){
fit.model=obj$fit.model
lower=obj$lower
type=obj$type
n1=dim(test$X)[1]
pred0=obj$pred0
misclass<-NULL->region->pred
if (isTRUE(all.equal(fit.model,NA))){
fit.model=pred0$fit
pred=predict.init(fit.model,test)
} else{
pred=predict.Quad2(fit.model,test$X)
pred_L=pred$pred_L
pred_R=pred$pred_R
}
if (!is.null(test$A)){
region=(test$A<pred_R)&(test$A>pred_L)
if ((!is.null(test$S))&(!is.null(test$Y))){
if (type=="EDI"){
misclass=(sum(test$alpha[,1]*(region)*pmax((test$S-test$Y),0))+sum(test$alpha[,2]*(!region)*pmax((test$Y-test$S),0)))/n1
} else if (type=="PDI"){
misclass=(sum(test$alpha[,1]*(region)*(test$Y<test$S))+sum(test$alpha[,2]*(!region)*(test$Y>test$S)))/n1
}
}
}
return(list(misclass=misclass,region=region,pred_L=pred_L,pred_R=pred_R))
}
# fit_two=RDCDI2(train=train,alpha=c(0.85,0.15),type='PDI',pred0=c(5.5,7.6),
# maxiter=5,eps=0.2,aL=NULL,aU=NULL,
# lambda=0.000000001)
#
# pred_L=predict(fit_two,test$X)$pred_L
# pred_R=predict(fit_two,test$X)$pred_R
#
# plotdata=data.frame(observed_A1c=test$A,recommended_A1c=pred_dose,
# lowerBound=pred_L,upperBound=pred_R)
# colnames(plotdata)[1]='observed_A1c'
# plotdata<- melt(plotdata)
# p<-ggplot(plotdata,aes(x=value, fill=variable,colour=variable)) +
# scale_x_continuous(breaks = sort(c(seq(2,12,0.5))))+
# geom_density(alpha=0.25)+
# xlab('A1c')
# p
###############################################################
lixiaomao/personalized-dosing-interval documentation built on May 16, 2019, 9:14 p.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.
#library(CVXR)
#library(LowRankQP)
#library(quadprog)
#' @export
doseFind<-function(train, test=NULL,eps=NULL,aL=NULL,aU=NULL,step=1,thres=NULL,maxiter=20,lambda=0){
if (is.null(thres)) thres=quantile(train$Y,0.8)
if (is.null(aL)) aL=min(train$A)-0.5*sd(train$A)
if (is.null(aU)) aU=max(train$A)+0.5*sd(train$A)
p=dim(train$X)[2]
n=dim(train$X)[1]
index=train$Y>=thres
beta <- Variable(p)
alpha<-Variable(1)
obj <- Minimize(sum_entries(abs(train$A[index] -alpha- train$X[index,] %*% beta)))
prob<- Problem(obj)
solv <- solve(prob)
beta0= solv$getValue(beta)
alpha0=solv$getValue(alpha)
pred0=alpha0+train$X%*% beta0
for (i in 1:maxiter){
solv=fit.dose(train,pred0=pred0,eps=eps)
if (sum((solv$beta-beta0)^2)<1e-5){
cat('Converged at iteration: ',i,'-th\n')
break;
} else {
beta0=solv$beta
alpha0=solv$alpha
}
pred0=alpha0+train$X%*% beta0
dosedif=pred0-train$A
index=(abs(dosedif)<0.1)
cat('iteration:',i,' sample size:',sum(index),'\n')
print(quantile(pred0))
}
pred_dose=predict.doseCVX(test$X,solv)
pred_dose=pmin(pmax(pred_dose,3),12)
return(list(pred=pred_dose,solv=solv,beta=solv$beta,alpha=solv$alpha))
}
#' @export
fit.dose<-function(train,pred0=NULL,beta0=NULL,alpha0=NULL,eps,aL=NULL,aU=NULL){
if (is.null(pred0)) pred0=alpha0+train$X%*% beta0
p=dim(train$X)[2]
Q=(abs(train$A-pred0)<eps)
Q1=train$A<pred0-eps
Q2=train$A>pred0+eps
beta <- Variable(p)
alpha<-Variable(1)
pred=alpha+train$X %*% beta
res=train$A - pred
part1=sum_entries(Q*abs(res))
part2=sum_entries((Q1)*2*pos(train$A-pred))
part3=sum_entries((Q2)*2*pos(pred-train$A))
partAll=sum_entries(1/train$propensity*train$Y*(part1+part2+part3))
obj <- Minimize(partAll)
prob<- Problem(obj)
solv <- solve(prob)
class(solv)='doseCVX'
solv$beta=solv$getValue(beta)
solv$alpha=solv$getValue(alpha)
return(solv)
}
#' @export
predict.doseCVX<-function(X,solv){
return(solv$alpha+X%*%solv$beta)
}
##############################################################################
## one-sided dosing interval
##############################################################################
#' @export
lineOne<-function(b,s1,s2,s3,s4,a,fit){
res=a-fit-b
return(sum(s1*pmax(-res,0)+s2*pmax(res,0)))
}
#' @export
fit.Int1<-function(train,eps,XX=NULL,CC=NULL,pred0=NULL,rbfkernel=NULL,
Intcep0=NULL,beta0=NULL,lambda=0.0001,type='PDI',lower=TRUE,method='completeLinear',...){
# if (FALSE){
# alpha=c(0.82,0.18)
# eps=0.3
# lambda=0.0001
# n=length(train$A)
# pred0=6
#
# }
n=length(train$A)
p=dim(train$X)[2]
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (is.null(CC)){
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = sigma)
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
CC=cbind(rbind(XX,-XX),rbind(-XX,XX))
}
d=matrix(c(-train$A,train$A),ncol = 1)
Aeq=matrix(c(rep(1,n),rep(-1,n)),ncol = 1)
if (type=='PDI'){
H1=(train$Y>train$S)*(train$alpha[,2])*train$propensity/lambda
H2=(train$Y<train$S)*(train$alpha[,1])*train$propensity/lambda
} else{
H1=pmax(train$Y-train$S,0)*(train$alpha[,2])*train$propensity/lambda
H2=pmax(train$S-train$Y,0)*(train$alpha[,1])*train$propensity/lambda
}
Q1=train$A<(pred0-eps)
Q2=train$A>(pred0+eps)
W1=train$A>(pred0-eps)
W2=train$A<(pred0+eps)
#Q1=Q1*0
#Q2=Q2*0
if (lower){
ub=c(H1*W1+H2*Q2,H1*Q1+H2*W2)
} else{
ub=c(H2*W1+H1*Q2,H2*Q1+H1*W2)
}
zero=(ub==0)
ubMat=diag(2*n)
Amat <- cbind(Aeq[!zero,],-ubMat[!zero,!zero],ubMat[!zero,!zero])
bvec=matrix(c(0,-ub[!zero],rep(0,sum(!zero))),ncol=1)
solv=solve.QP(Dmat=CC[!zero,!zero]+diag(0.1,sum(!zero)), dvec=d[!zero,], Amat=Amat, bvec=bvec, meq=1, factorized=FALSE)
alfa=rep(0,2*n)
alfa[!zero]=solv$solution
if (method=='completeLinear'){
beta=t((alfa[-c(1:n)]-alfa[1:n])%*%X)
pred=X%*%beta
} else {
beta=(alfa[-c(1:n)]-alfa[1:n])
pred=XX%*%beta
}
opt=optimize(interval = c(min(train$A),max(train$A)),f=lineOne,s1=H1*W1+H2*Q2,s2=H1*Q1+H2*W2,a=train$A,fit=pred)
intcept=opt$minimum
output=list(beta=beta,X=X,intcept=intcept,pred=pred+intcept,solv=solv,method=method,rbfkernel=rbfkernel)
class(output)='Quad1'
return(output)
}
#' @export
predict.Quad1<-function(fit,X){
X=apply(X,2,function(x) (x-mean(x))/sd(x))
if (fit$method=='completeLinear'){
pred=X%*%fit$beta+fit$intcept
} else {
XX=as.matrix(as.data.frame(kernelMatrix(fit$rbfkernel,X,fit$X)))
pred=XX%*%fit$beta+fit$intcept
}
return(pred)
}
#' @export
predict.RDCDI1<-function(obj,test){
fit.model=obj$fit.model
lower=obj$lower
type=obj$type
n1=dim(test$X)[1]
pred0=obj$pred0
misclass<-NULL->region->pred
if (isTRUE(all.equal(fit.model,NA))){
fit.model=pred0$fit
pred=predict.init(fit.model,test)
} else{
pred=predict.Quad1(fit.model,test$X)
}
if (!is.null(test$A)){
if (lower){
region=pred<test$A
} else {
region=pred>test$A
}
if ((!is.null(test$S))&(!is.null(test$Y))){
if (type=="EDI"){
misclass=(sum(test$alpha[,1]*(region)*pmax((test$S-test$Y),0))+sum(test$alpha[,2]*(!region)*pmax((test$Y-test$S),0)))/n1
} else if (type=="PDI"){
misclass=(sum(test$alpha[,1]*(region)*(test$Y<test$S))+sum(test$alpha[,2]*(!region)*(test$Y>test$S)))/n1
}
}
}
return(list(misclass=misclass,region=region,pred=pred))
}
#' @export
RDCDI1<-function(train,type='PDI',method=NULL,pred0=NULL,
maxiter=25,lower=TRUE,eps=NULL,aL=NULL,aU=NULL,
lambda=0.0001,...){
n=length(train$A)
p=dim(train$X)[2]
sigma=1/dim(train$X)[2]
if (is.null(aL)) aL<-min(train$A)-0.5*sd(train$A)
if (is.null(aU)) aU<-max(train$A)+0.5*sd(train$A)
if (is.null(eps)) eps=0.1*(aU-aL)
if (is.null(pred0)){
if (length(table(train$A))>4){
init0=initiate.continuous(train,train,side=1)
pred0=init0$pred
Intcep0=drop(coef(init0$fit))[1]
beta0=drop(coef(init0$fit))[-1]
} else {
init0=initiate.binary(train,train,side=1)
pred0=init0$pred
Intcep0=drop(coef(init0$fit))[1]
beta0=drop(coef(init0$fit))[-1]
}
} else {
Intcep0=median(pred0)
if (length(pred0)==1) pred0=rep(pred0,n)
if (method=='completeLinear'){
beta0=rep(0,p)
} else {
beta0=rep(0,n)
}
}
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = sigma)
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
CC=cbind(rbind(XX,-XX),rbind(-XX,XX))
for (i in 1:maxiter){
before<-Sys.time()
fit=fit.Int1(train,XX=XX,CC=CC,eps=eps,pred0=pred0,lambda=lambda,type=type,lower=lower,method=method,rbfkernel=rbfkernel)
after<-Sys.time()
after-before
fit0=fit
Intcep0=fit$intcept
pred0=drop(predict(fit,train$X))
if (sum(abs(fit$beta-beta0))<1e-15){
beta0=fit$beta
cat('Converged at iteration: ',i-1,'-th\n')
break;
}
beta0=fit$beta
print(quantile(pred0))
}
#1-sum((pred0-train$opt)^2)/sum(train$opt^2)
if (lower){
region=pred0<train$A
} else{
region=pred0>train$A
}
if (type=="EDI"){
if (lower){
misclass=(sum(train$alpha[,1]*(region)*pmax((train$S-train$Y),0))+sum(train$alpha[,2]*(!region)*pmax((train$Y-train$S),0)))/n
} else{
misclass=(sum(train$alpha[,1]*(!regopm)*pmax((train$S-train$Y),0))+sum(train$alpha[,2]*(region)*pmax((train$Y-train$S),0)))/n
}
} else if (type=="PDI"){
if (lower){
misclass=(sum(train$alpha[,1]*(region)*(train$Y<train$S))+sum(train$alpha[,2]*(!region)*(train$Y>train$S)))/n
} else{
misclass=(sum(train$alpha[,1]*(!region)*(train$Y<train$S))+sum(train$alpha[,2]*(region)*(train$Y>train$S)))/n
}
}
output=list(beta=fit$beta,intcept=fit$intcept,pred=pred0,solv=fit$solv,region=region,type=type,lower=lower,numiter=i,lambda=lambda,fit.model=fit,misclass=misclass)
class(output)='RDCDI1'
return(output)
}
#
# fit_L=RDCDI1(train=train,alpha=c(0.82,0.18),type='PDI',pred0=6,
# maxiter=5,lower=TRUE,eps=0.2,aL=NULL,aU=NULL,
# lambda=0.000001)
#
# pred_L=predict(fit_L,test$X)
#
# fit_R=RDCDI1(train=train,alpha=c(0.82,0.18),type='PDI',pred0=8,
# maxiter=5,lower=FALSE,eps=0.2,aL=NULL,aU=NULL,
# lambda=0.000001)
#
#
# pred_R=predict(fit_R,test$X)
#
# plotdata=data.frame(observed_A1c=test$A,recommended_A1c=pred_dose,
# lowerBound=pred_L,upperBound=pred_R)
# colnames(plotdata)[1]='observed_A1c'
# plotdata<- melt(plotdata)
# p<-ggplot(plotdata,aes(x=value, fill=variable,colour=variable)) +
# scale_x_continuous(breaks = sort(c(seq(2,12,0.5))))+
# geom_density(alpha=0.25)+
# xlab('A1c')
# p
#####################################################################
## two-sided version
#####################################################################
# lineTwo<-function(b,s1,s2,s3,s4,s5,s6,s7,s8,a,fitL,fitR,eps){
# resL1=a-fitL-b[1]
# resL2=a-fitL-b[1]+eps
# resR1=a-fitR-b[2]
# resR2=a-fitR-b[2]-eps
# return(sum(s1*pmax(-resL1,0)+s2*pmax(resL1,0)+s3*pmax(-resL2,0)+s4*pmax(resL2,0)+s5*pmax(resR1,0)+s6*pmax(-resR1,0)+s7*pmax(resR2,0)+s8*pmax(-resR2,0)))
# }
# fit.Int2<-function(train,eps,lambda=0.0001,type='PDI',CC=NULL,XX=NULL,lower=TRUE,
# pred0L=NULL,Intcep0L=NULL,beta0L=NULL,pred0R=NULL,Intcep0R=NULL,beta0R=NULL,...){
# if (is.null(pred0L)) pred0L=Intcep0L+train$X %*% beta0L
# if (is.null(pred0R)) pred0R=Intcep0R+train$X %*% beta0R
# # library(quadprog)
# # train=subsetData(train,sample(1:length(train$A),1000))
# # alpha=c(0.82,0.18)
# # eps=0.3
# # lambda=0.000001
# # pred0L=6.5
# # pred0R=8
# n=length(train$A)
# p=dim(train$X)[2]
#
# if (is.null(X)) X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
# if (is.null(CC)){
# XX=X%*%t(X)
# col1=rbind(XX,-XX,XX,-XX,XX,matrix(0,nrow=4*n,ncol=n))
# col2=rbind(matrix(0,nrow=4*n,ncol=n),XX,-XX,XX,-XX,XX)
# CC=cbind(col1,-col1,col1,-col1,col1+col2,-col2,col2,-col2,col2)
# }
#
# Q1L=train$A<(pred0L-eps)
# Q2L=train$A<(pred0L)
# W1L=train$A>(pred0L-eps)
# W2L=train$A>(pred0L)
#
# Q1R=train$A>(pred0R+eps)
# Q2R=train$A>(pred0R)
# W1R=train$A<(pred0R+eps)
# W2R=train$A<(pred0R)
#
# #Q1L=Q1L*0;Q2L=Q2L*0;Q1R=Q1R*0;Q2R=Q2R*0
# if (type=='PDI'){
# H1=(train$Y>train$S)*(train$alpha[,2])*train$propensity/lambda
# H2=(train$Y<train$S)*(train$alpha[,1])*train$propensity/lambda
# } else {
# H1=pmax(train$Y-train$S,0)*(train$alpha[,2])*train$propensity/lambda
# H2=pmax(train$S-train$Y,0)*(train$alpha[,1])*train$propensity/lambda
# }
#
#
#
# ub=c(H1*W1L,H1*Q1L,H2*W2L,H2*Q2L,rep(100000,n),H1*W1R,H1*Q1R,H2*W2R,H2*Q2R)
# zero0<-(ub==0)->zero1
# zero1[(4*n+1):(5*n)]=TRUE
#
#
#
# d=matrix(c(-train$A,train$A,-train$A-rep(eps,n),train$A+rep(eps,n),rep(0,n),-train$A,train$A,-train$A-rep(eps,n),train$A+rep(eps,n)),ncol = 1)
# Aeq=cbind(c(rep(1,n),rep(-1,n),rep(1,n),rep(-1,n),rep(1,n),rep(0,4*n)),c(rep(0,4*n),rep(-1,n),rep(1,n),rep(-1,n),rep(1,n),rep(-1,n)))
#
# ubMat=diag(9*n)
#
# Amat <- cbind(Aeq[!zero0,],-ubMat[!zero0,!zero1],ubMat[!zero0,!zero0])
# bvec=matrix(c(rep(0,2),-ub[!zero1],rep(0,sum(!zero0))),ncol=1)
# solv=solve.QP(Dmat=CC[!zero0,!zero0]+diag(0.00000001,sum(!zero0)), dvec=d[!zero0,], Amat=Amat, bvec=bvec, meq=2, factorized=FALSE)
# alfa=rep(0,9*n)
# alfa[!zero0]=solv$solution
# betaL=t((alfa[(n+1):(2*n)]-alfa[(1:n)]+alfa[(3*n+1):(4*n)]-alfa[(2*n+1):(3*n)]-alfa[(4*n+1):(5*n)])%*%X)
# betaR=t((alfa[(6*n+1):(7*n)]-alfa[(5*n+1):(6*n)]+alfa[(8*n+1):(9*n)]-alfa[(7*n+1):(8*n)]+alfa[(4*n+1):(5*n)])%*%X)
#
# #lineTwo(c(6,8),s1=H1*W1L+H2*Q2L,s2=H1*Q1L+H2*W2L,s3=H1*W1R+H2*Q2R,s4=H1*Q1R+H2*W2R,a=train$A,fitL=X%*%betaL,fitR=X%*%betaR)
# fit=optim(par=c(quantile(train$A,0.2),quantile(train$A,0.8)),fn=lineTwo,eps=eps,s1=H1*W1L,s2=H1*Q1L,s3=H2*W2L,s4=H2*Q2L,s5=H1*W1R,s6=H1*Q1R,s7=H2*W2R,s8=H2*Q2R,a=train$A,fitL=X%*%betaL,fitR=X%*%betaR)
#
#
# # mu=ub[-c((2*n+1):(3*n))]-alfa[-c((2*n+1):(3*n))]
# # aa=c(train$A-X%*%betaL,X%*%betaL-train$A,train$A-X%*%betaR,X%*%betaR-train$A)
# # aa[which((mu!=0)&(abs(alfa)[-c((2*n+1):(3*n))]>0.0000001))]
#
# output=list()
# output$betaL=betaL
# output$betaR=betaR
# output$intcepL=fit$par[1]
# output$intcepR=fit$par[2]
# output$pred_L=X%*%betaL+fit$par[1]
# output$pred_R=X%*%betaR+fit$par[2]
# output$fit=fit
# #output$solv=solv
# class(output)='Quad2'
# return(output)
# }
#' @export
lineTwo<-function(b,s1,s2,s3,s4,a,fitL,fitR){
resL=a-fitL-b[1]
resR=a-fitR-b[2]
return(sum(s1*pmax(-resL,0)+s2*pmax(resL,0)+s3*pmax(-resR,0)+s4*pmax(resR,0)+10*pmax(fitL+b[1]-fitR-b[2],0)))
}
#' @export
fit.Int2<-function(train,eps,lambda=0.0001,type='PDI',CC=NULL,XX=NULL,method='completeLinear',rbfkernel=NULL,
pred0L=NULL,Intcep0L=NULL,beta0L=NULL,pred0R=NULL,Intcep0R=NULL,beta0R=NULL,...){
if (is.null(pred0L)) pred0L=Intcep0L+train$X %*% beta0L
if (is.null(pred0R)) pred0R=Intcep0R+train$X %*% beta0R
#library(quadprog)
#train=subsetData(train,sample(1:length(train$A),1000))
#alpha=c(0.82,0.18)
#eps=0.3
#lambda=0.000001
n=length(train$A)
sigma=1/dim(train$X)[2]
#pred0L=6.5
#pred0R=8
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (is.null(CC)){
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = sigma)
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
col1=rbind(XX,-XX,XX,matrix(0,nrow=2*n,ncol=n))
col2=rbind(matrix(0,nrow=2*n,ncol=n),XX,-XX,XX)
CC=cbind(col1,-col1,col1+col2,-col2,col2)
}
Q1L=train$A<(pred0L-eps)
Q2L=train$A<(pred0L+eps)
W1L=train$A>(pred0L-eps)
W2L=train$A>(pred0L+eps)
Q1R=train$A>(pred0R+eps)
Q2R=train$A>(pred0R-eps)
W1R=train$A<(pred0R+eps)
W2R=train$A<(pred0R-eps)
#Q1L=Q1L*0;Q2L=Q2L*0;Q1R=Q1R*0;Q2R=Q2R*0
if (type=='PDI'){
H1=(train$Y>train$S)*(train$alpha[,2])*train$propensity/lambda
H2=(train$Y<train$S)*(train$alpha[,1])*train$propensity/lambda
} else{
H1=pmax(train$Y-train$S,0)*(train$alpha[,2])*train$propensity/lambda
H2=pmax(train$S-train$Y,0)*(train$alpha[,1])*train$propensity/lambda
}
ub=c(H1*W1L+H2*W2L,H1*Q1L+H2*Q2L,rep(1e05,n),H1*Q1R+H2*Q2R,H1*W1R+H2*W2R)
zero0<-(ub==0)->zero1
zero1[(2*n+1):(3*n)]=TRUE
d=matrix(c(-train$A,train$A,rep(2*eps,n),-train$A,train$A),ncol = 1)
Aeq=cbind(c(rep(1,n),rep(-1,n),rep(1,n),rep(0,2*n)),c(rep(0,2*n),rep(-1,n),rep(1,n),rep(-1,n)))
ubMat=diag(5*n)
Amat <- cbind(Aeq[!zero0,],-ubMat[!zero0,!zero1],ubMat[!zero0,!zero0])
bvec=matrix(c(rep(0,2),-ub[!zero1],rep(0,sum(!zero0))),ncol=1)
solv=solve.QP(Dmat=CC[!zero0,!zero0]+diag(0.1,sum(!zero0)), dvec=d[!zero0,], Amat=Amat, bvec=bvec, meq=2, factorized=FALSE)
alfa=rep(0,5*n)
alfa[!zero0]=solv$solution
if (method=='completeLinear'){
betaL=t((alfa[(n+1):(2*n)]-alfa[(1:n)]-alfa[(2*n+1):(3*n)])%*%X)
betaR=t((alfa[(4*n+1):(5*n)]-alfa[(3*n+1):(4*n)]+alfa[(2*n+1):(3*n)])%*%X)
pred_L=X%*%betaL
pred_R=X%*%betaR
} else {
betaL=(alfa[(n+1):(2*n)]-alfa[(1:n)]-alfa[(2*n+1):(3*n)])
betaR=(alfa[(4*n+1):(5*n)]-alfa[(3*n+1):(4*n)]+alfa[(2*n+1):(3*n)])
pred_L=XX%*%betaL
pred_R=XX%*%betaR
}
#lineTwo(c(6,8),s1=H1*W1L+H2*Q2L,s2=H1*Q1L+H2*W2L,s3=H1*W1R+H2*Q2R,s4=H1*Q1R+H2*W2R,a=train$A,fitL=X%*%betaL,fitR=X%*%betaR)
fit=optim(par=c(quantile(train$A,0.2),quantile(train$A,0.8)),fn=lineTwo,s1=H1*W1L+H2*W2L,s2=H1*Q1L+H2*Q2L,s3=H1*Q1R+H2*Q2R,s4=H1*W1R+H2*W2R,a=train$A,fitL=pred_L,fitR=pred_R)
# mu1L=ub[c((1):(n))]-alfa[c((1):(n))]
# mu2L=ub[c((n+1):(2*n))]-alfa[c((n+1):(2*n))]
# index1L=!zero0[c((1):(n))]
# bLub=(-X%*%betaL-train$A)[index1L&(abs(mu1L)<0.00000001)]
# index2L=!zero0[c((1+n):(2*n))]
# bLlb=(train$A-X%*%betaL)[index2L&(abs(mu2L)<0.00000001)]
#
# ##
# ## # this is supposed to be the right one
# bLlb=(train$A-X%*%betaL)[((alfa[c((1):(n))]>1e5)&index1L)|((mu2L>1e5)&index2L)]
# quantile(bLlb)
# bLub=(X%*%betaL-train$A)[((alfa[c((n+1):(2*n))]>1e5)&index2L)|((mu1L>1e5)&index1L)]
# quantile(bLub)
#
# mu1L=ub[c((1):(n))]-alfa[c((1):(n))]
# mu2L=ub[c((n+1):(2*n))]-alfa[c((n+1):(2*n))]
# index1L=!zero0[c((1):(n))]
# bLub=(-X%*%betaL-train$A)[index1L&abs(mu2L)==0]
# bLub=(-X%*%betaL-train$A)[index1L&(abs(alfa[c((1):(n))])<1e-10)]
# bLub=(-X%*%betaL-train$A)[index1L&(alfa[c((1):(n))]>1e-2)&(mu1L>1e-2)]
#
# index2L=!zero0[c((1+n):(2*n))]
# bLlb=(train$A-X%*%betaL)[index2L&abs(mu2L)==0]
# index2L=!zero0[c((1+n):(2*n))]
# bLlb=(train$A-X%*%betaL)[index2L&(alfa[c((n+1):(2*n))]>1e-2)&(mu2L>1e-2)]
# bLlb=(train$A-X%*%betaL)[index2L&(alfa[c((n+1):(2*n))]>1e-2)&(mu2L>1e-2)]
#
# index=(abs(mu)<0.00001)&(ub[-c((2*n+1):(3*n))]!=0)
# aa=c(train$A-X%*%betaL,X%*%betaL-train$A,train$A-X%*%betaR,X%*%betaR-train$A)
# aa[index]
# aa[which((mu!=0)&(abs(alfa)[-c((2*n+1):(3*n))]>0.0000001))]
output=list()
output$betaL=betaL
output$betaR=betaR
output$solv=solv
output$intcepL=fit$par[1]
output$intcepR=fit$par[2]
output$pred_L=pred_L+fit$par[1]
output$pred_R=pred_R+fit$par[2]
output$rbfkernel=rbfkernel
output$method=method
output$X=X
#output$solv=solv
class(output)='Quad2'
return(output)
}
#' @export
RDCDI2<-function(train,method='completeLinear',type='PDI',pred0L=NULL,pred0R=NULL,
maxiter=25,eps=NULL,aL=NULL,aU=NULL,lower=TRUE,
lambda=0.0001,...){
n=length(train$A)
p=dim(train$X)[2]
rbfkernel=NULL
X=apply(train$X,2,function(x) (x-mean(x))/sd(x))
if (method=='completeLinear'){
XX=X%*%t(X)
} else {
rbfkernel <- rbfdot(sigma = 1/dim(train$X)[2])
XX=as.matrix(as.data.frame(kernelMatrix(rbfkernel,X)))
}
col1=rbind(XX,-XX,XX,matrix(0,nrow=2*n,ncol=n))
col2=rbind(matrix(0,nrow=2*n,ncol=n),XX,-XX,XX)
CC=cbind(col1,-col1,col1+col2,-col2,col2)
if (is.null(aL)) aL<-min(train$A)-0.5*sd(train$A)
if (is.null(aU)) aU<-max(train$A)+0.5*sd(train$A)
if (is.null(eps)) eps=0.2*sd(train$A)
if (is.null(pred0L)|is.null(pred0R)){
if (length(table(train$A))>4){
init0=initiate.continuous(train,train,side=2)
} else {
init0=initiate.binary(train,train,side=2)
}
pred0L=fit$pred_L
pred0R=fit$pred_R
Intcep0L=drop(coef(init0$fit$fit_L))[1]
beta0L=drop(coef(init0$fit$fit_L))[-1]
Intcep0R=drop(coef(init0$fit$fit_R))[1]
beta0R=drop(coef(init0$fit$fit_R))[-1]
} else {
if (method=='completeLinear'){
beta0L=rep(0,p)
beta0R=rep(0,p)
} else {
beta0L=rep(0,n)
beta0R=rep(0,n)
}
Intcep0L=median(pred0L)
Intcep0R=median(pred0R)
if (length(pred0L)==1) pred0L=rep(pred0L,n)
if (length(pred0R)==1) pred0R=rep(pred0R,n)
}
#pred0L=rep(-0.5,n);pred0R=rep(0.5,n)
for (i in 1:maxiter){
before=Sys.time()
fit=fit.Int2(train,eps=eps,lambda=lambda,type=type,pred0L=pred0L,pred0R=pred0R,rbfkernel=rbfkernel,CC=CC,XX=XX,method=method)
after=Sys.time()
print(after-before)
pred0L=fit$pred_L
pred0R=fit$pred_R
if (sum(abs(fit$betaL-beta0L)+abs(fit$betaR-beta0R))<1e-15){
cat('Converged at iteration: ',i-1,'-th\n')
beta0L=fit$betaL
beta0R=fit$betaR
break;
}
beta0L=fit$betaL
beta0R=fit$betaR
cat('two-sided iteration: ',i,'\n')
print(quantile(pred0L))
print(quantile(pred0R))
}
region=(pred0L<train$A)&(pred0R>train$A)
if (type=="EDI"){
misclass=(sum(train$alpha[,1]*(region)*pmax((train$S-train$Y),0))+sum(train$alpha[,2]*(!region)*pmax((train$Y-train$S),0)))/n
} else if (type=="PDI"){
misclass=(sum(train$alpha[,1]*(region)*(train$Y<train$S))+sum(train$alpha[,2]*(!region)*(train$Y>train$S)))/n
}
output=list(method=method,betaL=fit$betaL,betaR=fit$betaR,intcepL=fit$intcepL,intcepR=fit$intcepR,pred_L=fit$pred_L,pred_R=fit$pred_R,region=region,type=type,lower=lower,numiter=i,lambda=lambda,fit.model=fit,misclass=misclass)
class(output)='RDCDI2'
return(output)
}
#' @export
predict.Quad2<-function(fit,X){
X=apply(X,2,function(x) (x-mean(x))/sd(x))
if (fit$method=='completeLinear'){
pred_L=fit$intcepL+X %*% fit$betaL
pred_R=fit$intcepR+X %*% fit$betaR
} else {
XX=as.matrix(as.data.frame(kernelMatrix(fit$rbfkernel,X,fit$X)))
pred_L=fit$intcepL+XX %*% fit$betaL
pred_R=fit$intcepR+XX %*% fit$betaR
}
return(list(pred_L=pred_L,pred_R=pred_R))
}
#' @export
predict.RDCDI2<-function(obj,test){
fit.model=obj$fit.model
lower=obj$lower
type=obj$type
n1=dim(test$X)[1]
pred0=obj$pred0
misclass<-NULL->region->pred
if (isTRUE(all.equal(fit.model,NA))){
fit.model=pred0$fit
pred=predict.init(fit.model,test)
} else{
pred=predict.Quad2(fit.model,test$X)
pred_L=pred$pred_L
pred_R=pred$pred_R
}
if (!is.null(test$A)){
region=(test$A<pred_R)&(test$A>pred_L)
if ((!is.null(test$S))&(!is.null(test$Y))){
if (type=="EDI"){
misclass=(sum(test$alpha[,1]*(region)*pmax((test$S-test$Y),0))+sum(test$alpha[,2]*(!region)*pmax((test$Y-test$S),0)))/n1
} else if (type=="PDI"){
misclass=(sum(test$alpha[,1]*(region)*(test$Y<test$S))+sum(test$alpha[,2]*(!region)*(test$Y>test$S)))/n1
}
}
}
return(list(misclass=misclass,region=region,pred_L=pred_L,pred_R=pred_R))
}
# fit_two=RDCDI2(train=train,alpha=c(0.85,0.15),type='PDI',pred0=c(5.5,7.6),
# maxiter=5,eps=0.2,aL=NULL,aU=NULL,
# lambda=0.000000001)
#
# pred_L=predict(fit_two,test$X)$pred_L
# pred_R=predict(fit_two,test$X)$pred_R
#
# plotdata=data.frame(observed_A1c=test$A,recommended_A1c=pred_dose,
# lowerBound=pred_L,upperBound=pred_R)
# colnames(plotdata)[1]='observed_A1c'
# plotdata<- melt(plotdata)
# p<-ggplot(plotdata,aes(x=value, fill=variable,colour=variable)) +
# scale_x_continuous(breaks = sort(c(seq(2,12,0.5))))+
# geom_density(alpha=0.25)+
# xlab('A1c')
# p
###############################################################
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.