#' Compute the L-optimal objective function assuming a logistic model (with or without weights)
#' @param D design matrix
#' @param cr matrix of contrasts
#' @param prior.scale prior scale parameter
#' @param wr matrix of weights, set to NULL if equal weights
#'
#' @return Value of the L-optimal objective function
#'
#' @export
calc.logit.wL <- function(D, cr, prior.scale=100, wr=NULL){
r <- nrow(cr)
var <- rep(0, r)
for (j in 1:r){
if(is.null(wr)){
var[j] <- t(cr[j,])%*% solve(D + diag(1/prior.scale ,ncol(Imat))) %*% cr[j,]
}else{
var[j] <- as.matrix(wr[j])%*%t(cr[j,])%*% solve(D + diag(1/prior.scale ,ncol(D))) %*% cr[j,]
}
}
opt <- sum(var)
return(opt)
}
#' Compute the L-optimal objective function assuming a logistic model (with or without weights) assuming continuous treatment,
#' given one new patient
#' @param t new treatment
#' @param D current deisgn matrix
#' @param z covariate values of new unit
#' @param beta current estimate of coefficients
#' @param cr matrix of contrasts
#' @param prior.scale prior scale parameter
#' @param wr matrix of weights, set to NULL if equal weights
#'
#' @return Value of the L-optimal objective function
#'
#' @export
wLopt.t <- function(t, D, z, beta, cr, prior.scale=100, wr=NULL){
if (is.null(dim(z))){
new.d <-c(1, z, t, z*t)
}else{
new.d <-cbind(1, z, t, t*z)
}
#rownames(new.d) <- rownames(D)
if (is.data.frame(new.d)){
colnames(new.d) <- colnames(D)
}
X <- as.matrix(rbind(D, new.d))
p <- apply(X, 1, probi, beta=beta) #calculate p(Y=1) for each row of the design matrix
W <- diag(p * (1 - p)) #weights
Imat <- t(X) %*% W %*% X
r <- nrow(cr)
var <- rep(0, r)
for (j in 1:r){
if(is.null(wr)){
var[j] <- t(cr[j,])%*% solve(Imat + diag(1/prior.scale ,ncol(Imat))) %*% cr[j,]
}else{
var[j] <- as.matrix(wr[j])%*%t(cr[j,])%*% solve(Imat + diag(1/prior.scale ,ncol(Imat))) %*% cr[j,]
}
}
opt <- sum(var)
return(opt)
}
#' Compute the L-optimal objective function assuming a logistic model (with or without weights) assuming continuous treatment,
#' given a vector of treatment values and a vector of covariates values for additional patients
#' @param t new treatment
#' @param z covariate values of new unit
#' @param beta current estimate of coefficients
#' @param cr matrix of contrasts
#' @param prior.scale prior scale parameter
#' @param wr matrix of weights, set to NULL if equal weights
#'
#' @return Value of the L-optimal objective function
#'
#' @export
wLopt.t.init <- function(t, z, beta, cr, prior.scale=100, wr=NULL){
X <-cbind(1, z, t, z*t)
X <- as.matrix(X)
#rownames(new.d) <- rownames(D)
#colnames(new.d) <- colnames(D)
p <- apply(X, 1, probi, beta=beta) #calculate p(Y=1) for each row of the design matrix
W <- diag(p * (1 - p)) #weights
Imat <- t(X) %*% W %*% X
r <- nrow(cr)
var <- rep(0, r)
for (j in 1:r){
if(is.null(wr)){
var[j] <- t(cr[j,])%*% solve(Imat + diag(1/prior.scale ,ncol(Imat))) %*% cr[j,]
}else{
var[j] <- as.matrix(wr[j])%*%t(cr[j,])%*% solve(Imat + diag(1/prior.scale ,ncol(Imat))) %*% cr[j,]
}
}
opt <- sum(var)
return(opt)
}
#' Compute the DA-optimal objective function assuming a logistic model (with or without weights)
#' @param D design matrix
#' @param cr matrix of contrasts
#' @param prior.scale prior scale parameter
#' @param wr matrix of weights, set to NULL if weights equal
#'
#' @return Value of the DA-optimal objective function
#'
#' @export
calc.logit.wDA <- function(D, cr, prior.scale=100, wr=NULL){
r <- nrow(cr)
if(is.null(wr)){
DA <- det( (cr)%*% solve(D+ diag(1/prior.scale ,ncol(D))) %*% t(cr))
}else{
var <-rep(0, r)
inv <- solve(D + diag(1/prior.scale, ncol(D)))
for (m in 1:r){
var[m] <- as.matrix( t(cr[m,]) %*% inv %*% as.matrix((cr[m,]))) ^ wr[m]
}
DA <- prod(var)
}
return(DA)
}
#' Assuming a logistic model for the response, allocate treatments using a coordinate exchange algorithm according
#' to an information matrix-based optimality criterion.
#' Responses are simulated assuming the true parameter values.
#'
#' @param covar a dataframe for the covariates
#' @param beta an estimate for the true values of beta
#' @param threshold the cut-off value for hypothesis tests
#' @param kappa the value of probability at which weights are set at zero
#' @param cr.lossfunc loss function appropriate for linear combinations of parameters
#' @param k the number of random starting designs in the coordinate exchange algorithm
#' @param wt set to T if the above lossfunction is weighted, NULL otherwise
#' @param ... further arguments to be passed to <lossfunc>
#'
#' @return design matrix D
#'
#'
#' @export
#'
cr.logit.coord <- function(covar, beta, threshold, kappa,cr.lossfunc=calc.logit.wL, k, wt=NULL, ... ){
covar <- as.data.frame(covar)
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
p <- length(beta)
cr1 <- matrix(0, 2^j, j+1)
cr2 <- expand.grid(rep(list(c(0,1)),j))
cr <- cbind(cr1, 1, cr2)
cr <- as.matrix(cr)
if(is.null(wt)){
wr <- rep(1, nrow(cr))
}else{
wr <- wt
}
Djs <- list() #list to store m potential designs
Ms <- list() #list to store the information matrices of the m potential designs
for (q in 1:k){
Dj <- as.matrix(cbind(rep(1, n), covar=covar, tmt=sample(c(0,1), n, replace=T))) #design matrix with random treatment assignment
Dj <- cbind(Dj, Dj[,2:(j+1)]*Dj[,(j+2)]) #append interaction columns
repeat{
D <- Dj
for (i in 1:n){
D[i, "tmt"] <- 1
D[i, (j+3):(2*j+2)] <- D[i, 2:(j+1)] #add interactino column
plusD <- cr.lossfunc(Imat.beta(D, beta), cr, prior.scale, wr) # ... needed in case lossfunc is DA
D[i, "tmt"] <- 0
D[i, (j+3):(2*j+2)] <- D[i, "tmt"]*D[i, 2:(j+1)]
minusD <- cr.lossfunc(Imat.beta(D, beta), cr, prior.scale, wr)
if (plusD < minusD){
D[i,"tmt"] <- 1
D[i, (j+3):(2*j+2)] <- D[i, 2:(j+1)]
} else {
D[i,"tmt"] <- 0
D[i, (j+3):(2*j+2)] <- D[i,"tmt"]* D[i, 2:(j+1)]
}
}
if (identical( cr.lossfunc(Imat.beta(D, beta), cr, prior.scale, wr) , cr.lossfunc(Imat.beta(Dj, beta), cr, prior.scale, wr) )) break #if no change to resulting design matrix, terminate
else Dj <- D
}
Djs[[q]] <- D
Ms[[q]] <- Imat.beta(D, beta)
}
mindet <- which.min(unlist(lapply(Ms, cr.lossfunc, cr, prior.scale, wr)))
D <- Djs[mindet][[1]]
return(D)
}
#' Assuming a logistic model for the response, allocate treatment sequentially based on an optimality criterion for
#' linear combinations of parameters. Responses are simulated assuming the true parameter values.
#'
#' @param covar a dataframe for the covariates
#' @param true.beta the true parameter values of the regression coefficients
#' @param threshold the cut-off value for hypothesis tests
#' @param kappa the value of probability at which weights are set at zero
#' @param init the number of units in the initial design
#' @param int set to T if you allow for treatment-covariate interactions in the model, NULL otherwise
#' @param cr.lossfunc loss function appropriate for olinear combinations of parameters
#' @param k the number of "outer loops" in the coordinate exchange algorithm
#' @param wt set to T if the above lossfunction is weighted, NULL otherwise
#' @param prior.scale the prior scale parameter
#' @param same.start the design matrix to be used for the initial design. If set to NULL, function generates initial design.
#' @param rand.start If set to T, function generates an initial design randomly. Else, coordinate exchange is used.
#' @param stoc set to T if treatments are allocated using a stochastic method where the probability is
#' determined by the optimality crtierion. Set to F if treatments are allocated deterministically.
#' @param bayes set to T if bayesglm is used instead of glm. Default prior assumed.
#' @param u vector of uniform random numbers for generating responses. If set to NULL, responses generated from the binomial distribution.
#' @param prior.default set to T if default priors for bayesglm is used. If set to False and bayes=T, normal priors used.
#' @param true.bvcov if set to T, use the true parameter values to compute obejctive function. If set to NULL, use estimated parameter values.
#' @param ... further arguments to be passed to <lossfunc>
#'
#'
#'
#'
#' @return design matrix D, responses y, all estimates of betas, final estimate of beta, all weights, all estimates of standard deviation,
#' beta, probabilities for treatment assignment, all values of the objective functions (weighted L, DA, weighted DA), proportion of treatment=1,
#' proportion of covariate in each group
#' Type 1 error, true value of power, empirical value of power.
#'
#' @export
cr.logit.des <- function(covar, true.beta, threshold, kappa, init, cr.lossfunc, k, wt, int=T, prior.scale=100,
same.start=NULL, rand.start=NULL, stoc=T, bayes=T, u=NULL, prior.default=T,
true.bvcov=NULL, ... ){
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
p <- length(true.beta)
cr1 <- matrix(0, 2^j, j+1)
cr2 <- expand.grid(rep(list(c(0,1)),j))
cr <- cbind(cr1, 1, cr2)
cr <- as.matrix(cr)
all.lopt <- all.dawopt <- c()
r <- nrow(cr)
loss.p <- loss.m <- c()
beta <- rep(0, p)
if (!is.null(same.start)) {
D <-same.start
}else if (!is.null(rand.start)) {
tmt <- sample(c(0,1), init, replace=T)
D <- cbind(rep(1, init), covar[1:init,], tmt, covar[1:init,]*tmt )
}else{
D <- as.matrix(logit.coord(as.data.frame(covar[1:init,]), beta, k, int=T, code=0,lossfunc=cr.lossfunc, cr, prior.scale))
}
pi <- apply(D, 1, probi, true.beta)
if (!is.null(u)){
y <- ifelse(u[1:init] < pi, 1, 0) #generate first observation based on true beta
}else{
y <- as.numeric(rbinom(init, 1, pi))
}
#find new estimate of beta by using logistic regression on the first init responses
if (bayes==TRUE){
if(prior.default==F){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit"), prior.df=Inf, prior.scale=prior.scale ))
}else{
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}
}else {
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- beta
crbeta <- cr %*% beta
pr <- pnorm(threshold, mean= crbeta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
true.pr <- pnorm(threshold, mean= cr %*% true.beta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
true.wr <- t(ifelse( true.pr >= kappa, true.pr, 0))
wr <- t(ifelse( pr >= kappa, pr, 0))
all.wr <- wr
yprop <- sum(y==1)/init
tmtprop <- sum(D[,"tmt"]==1)/init
if (!is.null(true.bvcov)){
all.lopt <- calc.logit.wL(Imat.beta(D, true.beta), cr, prior.scale, true.wr)
# all.dawopt <- calc.logit.wDA(Imat.beta(D,true.beta), cr, prior.scale, true.wr)
}else{
all.lopt <- calc.logit.wL(Imat.beta(D, beta), cr, prior.scale, wr)
#all.dawopt <- calc.logit.wDA(Imat.beta(D,beta), cr, prior.scale, wr)
}
for (i in (init+1):n){
#design matrix if tmt 1 is assigned - split into cases including/excluding interactions
if (wt!=T){
wr=NULL
}
if (!is.null(true.bvcov)){
beta <- true.beta
}
mat.plus <- rbind(D, c(1, as.numeric(covar[i, ]), 1,as.numeric(covar[i, ]) ))
d.plus <- cr.lossfunc(Imat.beta(mat.plus, beta), cr, prior.scale, wr)
mat.minus <- rbind(D, c(1, as.numeric(covar[i, ]), 0, rep(0, j) ))
d.minus <- cr.lossfunc(Imat.beta(mat.minus, beta), cr, prior.scale, wr)
loss.p <- c(loss.p, d.plus)
loss.m <- c(loss.m, d.minus)
probs <- (1/d.plus)/sum(1/d.plus+1/d.minus)
if (is.na(probs)){ #this is needed in case of separation occurring
probs <- 0
}
if (stoc==TRUE){
new.tmt <- sample(c(0,1), 1, prob=c(1-probs, probs)) #Assign treatments
}else{
if (probs > 0.5) {
new.tmt <- 1
} else {
new.tmt <- 0
}
}
#new row of design matrix
new.d <- as.numeric(c(1, as.numeric(covar[i, ]), new.tmt, as.numeric(covar[i, ])*new.tmt))
D <- as.matrix(rbind(D, as.numeric(new.d))) #Add the new row to the design matrix
row.names(D)<- NULL
pi <- probi(new.d, true.beta) #Compute new pi
if (!is.null(u)){
new.y <- ifelse(u[i] < pi, 1, 0)
}else{
new.y <- rbinom(1, 1, pi)
}
y <- c(y, new.y) #Simulate new observation
if(bayes==T){
if(prior.default==F){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit"), prior.df=Inf, prior.scale=prior.scale))
}else{
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit"))) #Update beta
}
all.beta <- rbind(all.beta, beta) #Store all betas
crbeta <- cr %*% beta
pr <- pnorm(threshold, mean= crbeta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
true.pr <- pnorm(threshold, mean= cr %*% true.beta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
true.wr <- t(ifelse( true.pr > kappa, true.pr, 0))
wr <- t(ifelse( pr >= kappa, pr, 0))
all.wr <- rbind(all.wr, wr)
if (!is.null(true.bvcov)){
lopt <- calc.logit.wL(Imat.beta(D, true.beta), cr, prior.scale, true.wr)
#dawopt <- calc.logit.wDA(Imat.beta(D,true.beta), cr, prior.scale, true.wr)
}else{
lopt <- calc.logit.wL(Imat.beta(D, beta), cr, prior.scale, wr)
#dawopt <- calc.logit.wDA(Imat.beta(D,beta), cr, prior.scale, wr)
}
all.lopt <- c(all.lopt,lopt)
# all.dawopt <- c(all.dawopt, dawopt)
yprop <- c( yprop, (sum(y==1)/i))
tmtprop <- c( tmtprop, (sum(D[,"tmt"]==1)/i))
}
D <- data.frame(D)
results <- list(D=D, y=y, all.beta=all.beta, all.wr= all.wr, beta = beta,
yprop=yprop, tmtprop=tmtprop, all.lopt=all.lopt,
#all.dawopt= all.dawopt,
loss.p=loss.p, loss.m=loss.m)
return(results)
}
#' Assuming a logistic model for the response, allocate treatment at random. Responses are simulated assuming the true parameter values.
#'
#' @param covar a dataframe for the covariates
#' @param true.beta the true parameter values of the regression coefficients
#' @param threshold the cut-off value for hypothesis tests
#' @param kappa the value of probability at which weights are set at zero
#' @param init the number of units in the initial design
#' @param int set to T if you allow for treatment-covariate interactions in the model, NULL otherwise
#' @param cr.lossfunc loss function appropriate for linear combinations of parameters
#' @param k the number of "outer loops" in the coordinate exchange algorithm
#' @param wt set to T if the above lossfunction is weighted, NULL otherwise
#' @param prior.scale the prior scale parameter
#' @param same.start the design matrix to be used for the initial design. If set to NULL, function generates initial design.
#' @param rand.start If set to T, function generates an initial design randomly. Else, coordinate exchange is used.
#' @param bayes set to T if bayesglm is used instead of glm. Default prior assumed.
#' @param u vector of uniform random numbers for generating responses. If set to NULL, responses generated from the binomial distribution.
#' @param prior.default set to T if default priors for bayesglm is used. If set to False and bayes=T, normal priors used.
#' @param true.bvcov if set to T, use the true parameter values to compute obejctive function. If set to NULL, use estimated parameter values.
#' @param ... further arguments to be passed to <lossfunc>
#'
#'
#'
#' @return design matrix D, responses y, all estimates of betas, final estimate of beta, all weights, all estimates of standard deviation,
#' beta, probabilities for treatment assignment, all values of the objective functions (weighted L, DA, weighted DA), proportion of treatment=1,
#' proportion of covariate in each group
#' Type 1 error, true value of power, empirical value of power.
#'
#' @export
cr.rand <- function(covar, true.beta, threshold, kappa, init, cr.lossfunc, k, wt, int=T, prior.scale=100,
same.start=NULL, rand.start=NULL, bayes=T, u=NULL,prior.default=T,
true.bvcov=NULL, ... ){
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
p <- length(true.beta)
cr1 <- matrix(0, 2^j, j+1)
cr2 <- expand.grid(rep(list(c(0,1)),j))
cr <- cbind(cr1, 1, cr2)
cr <- as.matrix(cr)
all.lopt <- all.dawopt <- c()
r <- nrow(cr)
beta <- rep(0, p)
if (!is.null(same.start)) {
D <-same.start
}else if (!is.null(rand.start)) {
tmt <- sample(c(0,1), init, replace=T)
D <- cbind(rep(1, init), covar[1:init,], tmt, covar[1:init,]*tmt )
}else{
D <- as.matrix(logit.coord(as.data.frame(covar[1:init,]), beta, k, int=T, code=0,lossfunc=cr.lossfunc, cr, prior.scale))
}
pi <- apply(D, 1, probi, true.beta)
if (!is.null(u)){
y <- ifelse(u[1:init] < pi, 1, 0) #generate first observation based on true beta
}else{
y <- as.numeric(rbinom(init, 1, pi))
}
#find new estimate of beta by using logistic regression on the first init responses
if (bayes==TRUE){
if(prior.default==F){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit"), prior.df=Inf, prior.scale=prior.scale ))
}else{
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}
}else {
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- beta
crbeta <- cr %*% beta
pr <- pnorm(threshold, mean= crbeta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
wr <- t(ifelse( pr >= kappa, pr, 0))
all.wr <- wr
yprop <- sum(y==1)/init
tmtprop <- sum(D[,"tmt"]==1)/init
if (!is.null(true.bvcov)){
lopt <- calc.logit.wL(Imat.beta(D, true.beta), cr, prior.scale, wr)
dawopt <- calc.logit.wDA(Imat.beta(D,true.beta), cr, prior.scale, wr)
}else{
lopt <- calc.logit.wL(Imat.beta(D, beta), cr, prior.scale, wr)
dawopt <- calc.logit.wDA(Imat.beta(D,beta), cr, prior.scale, wr)
}
all.lopt <- c(all.lopt,lopt)
all.dawopt <- c(all.dawopt, dawopt)
for (i in (init+1):n){
new.tmt <- sample(c(0,1), 1)
new.d <- as.numeric(c(1, as.numeric(covar[i, ]), new.tmt, as.numeric(covar[i, ])*new.tmt))
D <- as.matrix(rbind(D, as.numeric(new.d))) #Add the new row to the design matrix
row.names(D)<- NULL
pi <- probi(new.d, true.beta) #Compute new pi
if (!is.null(u)){
new.y <- ifelse(u[i] < pi, 1, 0)
}else{
new.y <- rbinom(1, 1, pi)
}
y <- c(y, new.y) #Simulate new observation
if(bayes==T){
if(prior.default==F){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit"), prior.df=Inf, prior.scale=prior.scale))
}else{
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit"))) #Update beta
}
all.beta <- rbind(all.beta, beta) #Store all betas
crbeta <- cr %*% beta
pr <- pnorm(threshold, mean= crbeta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
wr <- t(ifelse( pr > kappa, pr, 0))
all.wr <- rbind(all.wr, wr)
if (!is.null(true.bvcov)){
lopt <- calc.logit.wL(Imat.beta(D, true.beta), cr, prior.scale, wr)
dawopt <- calc.logit.wDA(Imat.beta(D,true.beta), cr, prior.scale, wr)
}else{
lopt <- calc.logit.wL(Imat.beta(D, beta), cr, prior.scale, wr)
dawopt <- calc.logit.wDA(Imat.beta(D,beta), cr, prior.scale, wr)
}
all.lopt <- c(all.lopt,lopt)
all.dawopt <- c(all.dawopt, dawopt)
yprop <- c( yprop, (sum(y==1)/i))
tmtprop <- c( tmtprop, (sum(D[,"tmt"]==1)/i))
}
D <- data.frame(D)
results <- list(D=D, y=y, all.beta=all.beta, all.wr= all.wr, beta = beta,
yprop=yprop, tmtprop=tmtprop, all.lopt=all.lopt,all.dawopt= all.dawopt)
return(results)
}
#' Assuming a logistic model for the response, allocate treatment sequentially based on an optimality criterion for
#' linear combinations of parameters. Responses are simulated assuming the true parameter values.
#'
#' @param covar a dataframe for the covariates
#' @param true.beta the true parameter values of the regression coefficients
#' @param threshold the cut-off value for hypothesis tests
#' @param kappa the value of probability at which weights are set at zero
#' @param init the number of units in the initial design
#' @param k the number of "outer loops" in the coordinate exchange algorithm
#' @param wt set to T if the above lossfunction is weighted, NULL otherwise
#' @param prior.scale the prior scale parameter
#' @param same.start the design matrix to be used for the initial design. If set to NULL, function generates initial design.
#' @param rand.start If set to T, function generates an initial design randomly. Else, coordinate exchange is used.
#' @param bayes set to T if bayesglm is used instead of glm. Default prior assumed.
#' @param u vector of uniform random numbers for generating responses. If set to NULL, responses generated from the binomial distribution.
#' @param prior.default set to T if default priors for bayesglm is used. If set to False and bayes=T, normal priors used.
#' @param true.bvcov if set to T, use the true parameter values to compute obejctive function. If set to NULL, use estimated parameter values.
#' @param ... further arguments to be passed to <lossfunc>
#'
#'
#'
#'
#' @return design matrix D, responses y, all estimates of betas, final estimate of beta, all weights, all estimates of standard deviation,
#' beta, probabilities for treatment assignment, all values of the objective functions (weighted L, DA, weighted DA), proportion of treatment=1,
#' proportion of covariate in each group
#' Type 1 error, true value of power, empirical value of power.
#'
#' @export
cr.logit.cont <- function(covar, true.beta, threshold, kappa, init, k, wt, prior.scale=100,
same.start=NULL, rand.start=NULL, bayes=T, u=NULL, prior.default=T,
true.bvcov=NULL, ... ){
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
p <- length(true.beta)
cr1 <- matrix(0, 2^j, j+1)
cr2 <- expand.grid(rep(list(c(0,1)),j))
cr <- cbind(cr1, 1, cr2)
cr <- as.matrix(cr)
all.lopt <- c()
allt.lopt <- betat.lopt <- not.lopt <- c()
r <- nrow(cr)
beta <- rep(0, p)
if (!is.null(same.start)) {
D <-same.start
}else if (!is.null(rand.start)) {
tmt <- sample(c(0,1), init, replace=T)
D <- cbind(rep(1, init), covar[1:init,], tmt, covar[1:init,]*tmt )
}else{
D <- as.matrix(logit.coord(as.data.frame(covar[1:init,]), beta, k, int=T, code=0,lossfunc=cr.lossfunc, cr, prior.scale))
}
pi <- apply(D, 1, probi, t(true.beta))
if (!is.null(u)){
y <- ifelse(u[1:init] < pi, 1, 0) #generate first observation based on true beta
}else{
y <- as.numeric(rbinom(init, 1, pi))
}
#find new estimate of beta by using logistic regression on the first init responses
if (bayes==TRUE){
if(prior.default==F){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit"), prior.df=Inf, prior.scale=prior.scale ))
}else{
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}
}else {
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- beta
crbeta <- cr %*% beta
pr <- pnorm(threshold, mean= crbeta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
wr <- t(ifelse( pr >= kappa, pr, 0))
all.wr <- wr
true.pr <- pnorm(threshold, mean= cr %*% (true.beta), sd=sqrt(diag( cr%*%solve( Imat.beta(D, t(true.beta)) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
true.wr <- t(ifelse( true.pr > kappa, true.pr, 0))
yprop <- sum(y==1)/init
tmtprop <- sum(D[,"tmt"]==1)/init
for (i in (init+1):n){
#design matrix if tmt 1 is assigned - split into cases including/excluding interactions
if (wt!=T){
wr=NULL
}
findopt <- optim(runif(1, min=0, max=1), fn=wLopt.t, D=D, z=as.numeric(covar[i, ]), beta =beta, cr=cr, prior.scale=prior.scale, wr=wr, method="L-BFGS-B",lower=0, upper=1)
new.tmt <- findopt$par
#new row of design matrix
new.d <- as.numeric(c(1, as.numeric(covar[i, ]), new.tmt, as.numeric(covar[i, ])*new.tmt))
D <- as.matrix(rbind(D, as.numeric(new.d))) #Add the new row to the design matrix
row.names(D)<- NULL
pi <- probi(new.d, t(true.beta)) #Compute new pi
if (!is.null(u)){
new.y <- ifelse(u[i] < pi, 1, 0)
}else{
new.y <- rbinom(1, 1, pi)
}
y <- c(y, new.y) #Simulate new observation
if(bayes==T){
if(prior.default==F){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit"), prior.df=Inf, prior.scale=prior.scale))
}else{
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit"))) #Update beta
}
all.beta <- rbind(all.beta, beta) #Store all betas
crbeta <- cr %*% beta
pr <- pnorm(threshold, mean= crbeta, sd=sqrt(diag( cr%*%solve( Imat.beta(D, beta) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
wr <- t(ifelse( pr > kappa, pr, 0))
all.wr <- rbind(all.wr, wr)
true.pr <- pnorm(threshold, mean= cr %*% (true.beta), sd=sqrt(diag( cr%*%solve( Imat.beta(D, t(true.beta)) +diag(1/prior.scale, ncol(D)))%*% t(cr))))
true.wr <- t(ifelse( true.pr > kappa, true.pr, 0))
#if (!is.null(true.bvcov)){
# new.opt <- wLopt.t(t=new.tmt, D=D, z=as.numeric(covar[i, ]), beta =t(true.beta), cr=cr, prior.scale=prior.scale, wr=true.wr)
#}else{
# new.opt <- findopt$value
#}
allt.lopt <- c( allt.lopt, wLopt.t(t=new.tmt, D=D, z=as.numeric(covar[i, ]), beta =t(true.beta), cr=cr, prior.scale=prior.scale, wr=true.wr))
not.lopt <- c(not.lopt, findopt$value)
yprop <- c( yprop, (sum(y==1)/i))
tmtprop <- c( tmtprop, (sum(D[,"tmt"]==1)/i))
}
D <- data.frame(D)
results <- list(D=D, y=y, all.beta=all.beta, all.wr= all.wr, beta = beta,
yprop=yprop, tmtprop=tmtprop, all.lopt=all.lopt,
allt.lopt= allt.lopt, not.lopt =not.lopt)
return(results)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.