#' Compute P(y=1) for logistic regression
#'
#' @param x a row of the design matrix
#' @param beta coefficients of the model
#' @return Probability P(y=1)'
#'
#' @export
probi <- function(x, beta){
x <- as.matrix(x)
eta <- beta %*% x
p <- (exp(eta)/(1+exp(eta))) #logistic function
return(p)
}
#' Compute the information matrix, given beta, for logistic regression
#'
#' @param X the design matrix
#' @param beta coefficients of the model
#' @return Information matrix
#'
#' @export
Imat.beta <- function(X, beta){
X <- as.matrix(X)
p <- apply(X, 1, probi, beta=beta) #calculate p(Y=1) for each row of the design matrix
W <- diag(p * (1 - p)) #weights
I.mat <- t(X) %*% W %*% X
return(I.mat)
}
#' Given an information matrix and assuming logistic regression, compute the D-optimal objective function
#'
#' @param I the information matrix
#' @param epsilon a small real number used for regularization. If set to zero,
#' no regularization takes place
#' @return value of the D-optimal objective function
#'
#' @export
#'
calc.y.D <- function(I, epsilon=0.00001){
k <- ncol(I)
I <- as.matrix(I)
M <- I + diag(k)*epsilon
D <- 1/det(M)
return(D)
}
#' Given an information matrix and assuming logistic regression, compute the A-optimal objective function
#'
#' @param I the information matrix
#' @param epsilon a small real number used for regularization. If set to zero,
#' no regularization takes place
#' @return value of the A-optimal objective function
#'
#' @export
#'
calc.y.A <- function(I, epsilon=0.00001){
k <- ncol(I)
I <- as.matrix(I)
M <- solve(I + diag(k)*epsilon)
D <- sum(diag(M))
return(D)
}
#' Given a design matrix with an additional row, and assuming assuming logistic regression, compute the D-optimal objective function
#' (to be used in conjunction with optim)
#'
#' @param t new treatment
#' @param D current deisgn matrix
#' @param z covariate values of new unit
#' @param int set to TRUE if treatment-covariate interactions are included in the model
#' @param beta current estimate of regression parameters
#' @param epsilon a small real number used for regularization. If set to zero,
#' no regularization takes place
#' @return value of the D-optimal objective function
#'
#' @export
#'
Dopt.y.t <- function(t, D, z, int, beta, epsilon=0.00001){
k <- ncol(D)
if (!is.null(int)&length(z)==1){
new.d <-c(1, z, t, z*t)
}else if (!is.null(int)&length(z)>1){
new.d <-cbind(1, z, t, t%*%as.matrix(z))
}else {
new.d <-as.matrix(cbind(1, z, t))
}
row.names(new.d) <- rownames(D) <- NULL
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
I <- t(X) %*% W %*% X
opt <- 1/det(I + diag(k)*epsilon)
return(opt)
}
#' Given an initial design and assuming logistic regression, compute the D-optimal objective function
#' (to be used in conjunction with optim)
#'
#' @param t vector of treatments
#' @param z vector of covariate values
#' @param int set to TRUE if treatment-covariate interactions are included in the model
#' @param beta current estimate of regression parameters
#' @param epsilon a small real number used for regularization. If set to zero,
#' no regularization takes place
#' @return value of the D-optimal objective function
#'
#' @export
Dopt.y.t.init <- function(t, z, int, beta, epsilon=0.00001){
if (!is.null(int)){
X <-cbind(1, z, t, z*t)
}else{
X <-cbind(1, z, t)
}
X <- as.matrix(X)
k <- ncol(X)
p <- apply(X, 1, probi, beta=beta) #calculate p(Y=1) for each row of the design matrix
W <- diag(p * (1 - p)) #weights
I <- t(X) %*% W %*% X
opt <- 1/det(I + diag(k)*epsilon)
return(opt)
}
#' Given an information matrix and assuming logistic regression, compute the G-optimal objective function
#' @param I the information matrix
#' @param epsilon a small real number used for regularization. If set to zero,
#' no regularization takes place
#' @return value of the G-optimal objective function
#'
#' @export
#'
calc.y.G <- function(I, epsilon=0.00001){
k <- ncol(I)
I <- as.matrix(I)
X <- as.matrix(cbind(1, expand.grid(rep(list(c(-1,1)),k-1))))
G <- max(diag((X)%*%solve(I + diag(k)*epsilon)%*%t(X)))
return(G)
}
#' Given an information matrix and assuming logistic regression, compute the D-A optimal objective function
#' @param I the information matrix
#' @param A a matrix where each column indicates the linear combination of parameters of interest
#' @param epsilon a small real number used for regularization. If set to zero,
#' no regularization takes place
#' @return value of the D-A optimal objective function
#'
#' @export
calc.y.DA <- function(I, A=t(matrix(c(0, 0, 1))), epsilon=0.00001){
j <- ncol(I)
I <- as.matrix(I)
loss <- sum(diag(A%*%solve(I + epsilon*diag(j))%*%t(A)))
return(loss)
}
#' 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 data generating mechanism
#' @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 lossfunc the objective function to minimize
#' @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 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 beta, probabilities for treatment assignment,
#' proportion of Y=1, proportion tmt=1, values of objective function
#'
#'
#' @export
#'
logit.rand <- function(covar, true.beta, init, int=NULL, lossfunc=calc.y.D, same.start=NULL, rand.start=NULL, bayes=T, u=NULL, true.bvcov=NULL, ...){
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
all.probs <-rep(0, n-init)
if(!is.null(int)){
beta <- rep(0, j+2+j)
}else{
beta <- rep(0, j+2)
}
opt <- yprop.tot <- yprop <- tmtprop <- c()
#starting design
if (!is.null(same.start)) {
D <-same.start
}else if (!is.null(rand.start)) {
D <- cbind(rep(1, init), covar[1:init,], sample(c(-1,1), init, replace=T))
}else if (!is.null(int)) {
D <- logit.coord(covar[1:init,], beta, 2, int=T, lossfunc, ...)
}else{
D <- logit.coord(covar[1:init,], beta, 2, int=NULL, lossfunc, ...)
}
pi <- apply(D, 1, probi, true.beta)
if (!is.null(u)){
y <- ifelse(u[1:init] < pi, 1, 0)
}else{
y <- as.numeric(rbinom(init, 1, pi))
}
#find new estimate of beta by using logistic regression on the first init responses
if(bayes==T){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- beta
for (i in (init+1):n){
new.tmt <- sample(c(-1,1), 1)
#new row of design matrix
if (!is.null(int)){
new.d <- as.numeric(c(1, covar[i, ], new.tmt, covar[i, ]*new.tmt))
} else{
new.d <- as.numeric(c(1, covar[i, ], new.tmt))
}
D <- as.matrix(rbind(D, as.numeric(new.d))) #Add the new row to the design matrix
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){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- rbind(all.beta, beta) #Store all betas
if (!is.null(true.bvcov)){
opt <- c(opt, lossfunc(Imat.beta(D, true.beta), ...))
}else{
opt <- c(opt, lossfunc(Imat.beta(D, beta), ...))
}
yprop.tot <- c(yprop.tot, sum(y==1)/i)
#yprop <- rbind(yprop, table(D[,2],y)[,2] /table(D[,2]))
}
row.names(D) <-NULL
D <- data.frame(D)
results <- list(D=D, y=y, all.beta=all.beta, beta = beta, all.probs = all.probs, #yprop=yprop,
yprop.tot=yprop.tot,
opt=opt
)
return(results)
}
#' Assuming a logistic model for the response, allocate treatment sequentially based on an information matrix-based optimality criterion.
#' Responses are simulated assuming the true parameter values.
#'
#' @param covar a dataframe for the covariates
#' @param true.beta the true parameter values of the data generating mechanism
#' @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 lossfunc the objective function to minimize
#' @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 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 beta, probabilities for treatment assignment,
#' proportion of Y=1, proportion tmt=1, values of objective function
#'
#'
#' @export
#'
logit.des <- function(covar, true.beta, init, int=NULL, lossfunc=calc.y.D, same.start=NULL, rand.start=NULL, stoc=T, bayes=T, u=NULL, true.bvcov=NULL, ...){
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
all.probs <-rep(0, n-init)
if(!is.null(int)){
beta <- rep(0, j+2+j)
}else{
beta <- rep(0, j+2)
}
loss.p <- loss.m <- c()
opt <- yprop.tot <- yprop <- tmtprop <- c()
#starting design
if (!is.null(same.start)) {
D <-same.start
}else if (!is.null(rand.start)) {
D <- cbind(rep(1, init), covar[1:init,], sample(c(-1,1), init, replace=T))
}else if (!is.null(int)) {
D <- logit.coord(covar[1:init,], beta, 2, int=T, lossfunc, ...)
}else{
D <- logit.coord(covar[1:init,], beta, 2, int=NULL, lossfunc, ...)
}
pi <- apply(D, 1, probi, true.beta)
if (!is.null(u)){
y <- ifelse(u[1:init] < pi, 1, 0)
}else{
y <- as.numeric(rbinom(init, 1, pi))
}
#find new estimate of beta by using logistic regression on the first init responses
if(bayes==T){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- beta
for (i in (init+1):n){
#design matrix if tmt 1 is assigned - split into cases including/excluding interactions
if (!is.null(int)){
mat.plus <- rbind(D, c(1, as.numeric(covar[i, ]), 1,as.numeric(covar[i, ]) ))
} else{
mat.plus <- rbind(D, c(1, as.numeric(covar[i, ]), 1))
}
#design matrix if tmt -1 is assigned
if (!is.null(int)){
mat.minus <- rbind(D, c(1, as.numeric(covar[i, ]), -1, -as.numeric(covar[i, ]) ))
} else{
mat.minus <- rbind(D, c(1, as.numeric(covar[i, ]), -1)) #Design matrix with new treatment =-1
}
if (!is.null(true.bvcov)){
d.plus <- lossfunc(Imat.beta(mat.plus, true.beta), ...)
d.minus <- lossfunc(Imat.beta(mat.minus, true.beta), ...)
}else{
d.plus <- lossfunc(Imat.beta(mat.plus, beta), ...)
d.minus <- lossfunc(Imat.beta(mat.minus, beta), ...)
}
#cat(d.plus, d.minus, "\n")
loss.p <- c(loss.p, d.plus)
loss.m <- c(loss.m, d.minus)
probs <- (1/d.minus)/(1/d.plus+1/d.minus)
all.probs[i-init] <- probs
if(stoc==T){
new.tmt <- sample(c(-1,1), 1, prob=c(probs, 1-probs)) #Assign treatments
}else{
if (d.plus > d.minus) {
new.tmt <- -1
} else if (d.plus < d.minus) {
new.tmt <- 1
} else if (d.plus == d.minus) {
new.tmt <- sample(c(-1,1), 1)
}
}
#new row of design matrix
if (!is.null(int)){
new.d <- as.numeric(c(1, covar[i, ], new.tmt, covar[i, ]*new.tmt))
} else{
new.d <- as.numeric(c(1, covar[i, ], new.tmt))
}
D <- as.matrix(rbind(D, as.numeric(new.d))) #Add the new row to the design matrix
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){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- rbind(all.beta, beta) #Store all betas
if (!is.null(true.bvcov)){
opt <- c(opt, lossfunc(Imat.beta(D, true.beta), ...))
}else{
opt <- c(opt, lossfunc(Imat.beta(D, beta), ...))
}
yprop.tot <- c(yprop.tot, sum(y==1)/i)
#yprop <- rbind(yprop, table(D[,2],y)[,2] /table(D[,2]))
}
row.names(D) <-NULL
D <- data.frame(D)
results <- list(D=D, y=y, all.beta=all.beta, beta = beta, all.probs = all.probs, #yprop=yprop,
yprop.tot=yprop.tot,
opt=opt,
loss.p=loss.p, loss.m=loss.m
)
return(results)
}
#' Assuming a logistic model for the response, allocate a continuous treatment sequentially based on an information matrix-based optimality criterion.
#' Responses are simulated assuming the true parameter values.
#'
#' @param covar a dataframe for the covariates
#' @param true.beta the true parameter values of the data generating mechanism
#' @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 lossfunc the objective function to minimize
#' @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 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 beta, probabilities for treatment assignment,
#' proportion of Y=1, proportion tmt=1, values of objective function
#'
#'
#' @export
#'
logit.cont <- function(covar, true.beta, init, int=NULL, lossfunc=Dopt.y.t, same.start=NULL, rand.start=NULL, bayes=T, u=NULL, true.bvcov=NULL, ...){
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
all.probs <-rep(0, n-init)
if(!is.null(int)){
beta <- rep(0, j+2+j)
}else{
beta <- rep(0, j+2)
}
opt <- yprop.tot <- yprop <- tmtprop <- c()
#starting design
if (!is.null(same.start)) {
D <-same.start
}else if (!is.null(rand.start)) {
D <- cbind(rep(1, init), covar[1:init,], runif(init, -1, 1))
}else if (!is.null(int)) {
D <- logit.coord(covar[1:init,], beta, 2, int=T, lossfunc, ...)
}else{
D <- logit.coord(covar[1:init,], beta, 2, int=NULL, lossfunc, ...)
}
pi <- apply(D, 1, probi, true.beta)
if (!is.null(u)){
y <- ifelse(u[1:init] < pi, 1, 0)
}else{
y <- as.numeric(rbinom(init, 1, pi))
}
#find new estimate of beta by using logistic regression on the first init responses
if(bayes==T){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- beta
for (i in (init+1):n){
findopt <- optim(par=runif(1, min=-1, max=1), fn=Dopt.y.t, D=D, z=(covar[i, ]), int=int, beta =beta, method="L-BFGS-B",lower=-1, upper=1)
new.tmt <- findopt$par
if (!is.null(true.bvcov)){
new.opt <- Dopt.y.t(t=as.numeric(new.tmt), D=D, z=(covar[i, ]), int=int, beta =(true.beta))
}else{
new.opt <- findopt$value
}
opt <- c(opt, new.opt)
#new row of design matrix
if (!is.null(int)){
new.d <- as.numeric(c(1, covar[i, ], new.tmt, covar[i, ]*new.tmt))
} else{
new.d <- as.numeric(c(1, covar[i, ], new.tmt))
}
D <- as.matrix(rbind(D, as.numeric(new.d))) #Add the new row to the design matrix
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){
beta <- coef(bayesglm(y~D[,-1], family=binomial(link="logit")))
}else{
beta <- coef(glm(y~D[,-1], family=binomial(link="logit")))
}
all.beta <- rbind(all.beta, beta) #Store all betas
yprop.tot <- c(yprop.tot, sum(y==1)/i)
#yprop <- rbind(yprop, table(D[,2],y)[,2] /table(D[,2]))
}
row.names(D) <-NULL
D <- data.frame(D)
results <- list(D=D, y=y, all.beta=all.beta, beta = beta, #yprop=yprop,
yprop.tot=yprop.tot,
opt=opt)
return(results)
}
#' 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 k an integer for the number of "outer loops"
#' @param int set to T if you allow for treatment-covariate interactions in the model, NULL otherwise
#' @param code set to NULL if (-1,1) coding is used for the treatments. Set to 0 if (0, 1) is used.
#' @param lossfunc the objective function to minimize
#' @param ... further arguments to be passed to <lossfunc>
#' @return design matrix D
#'
#'
#' @export
#'
logit.coord <- function(covar, beta, k, int=NULL, code=NULL, lossfunc=calc.y.D, ... ){
covar <- as.data.frame(covar)
n <- nrow(covar)
j <- ncol(covar) #covar must be a dataframe
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){
if (!is.null(code)){
Dj <- as.matrix(cbind(rep(1, n), covar=covar, tmt=sample(c(0,1), n, replace=T))) #design matrix with random treatment assignment
}else{
Dj <- as.matrix(cbind(rep(1, n), covar=covar, tmt=sample(c(-1,1), n, replace=T))) #design matrix with random treatment assignment
}
if (!is.null(int)){
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
if (!is.null(int)){
D[i, (j+3):(2*j+2)] <- D[i, 2:(j+1)] #add interactino column
}
plusD <- lossfunc(Imat.beta(D, beta), ...) # ... needed in case lossfunc is DA
if (!is.null(code)){
D[i, "tmt"] <- 0
}else{
D[i, "tmt"] <- -1 #calculate criterion for unit i assigned to tmt -1
}
if (!is.null(int)){
D[i, (j+3):(2*j+2)] <- D[i, "tmt"]*D[i, 2:(j+1)]
}
minusD <- lossfunc(Imat.beta(D, beta), ...)
if (plusD < minusD){
D[i,"tmt"] <- 1
if (!is.null(int)){
D[i, (j+3):(2*j+2)] <- D[i, 2:(j+1)]
}
} else {
if(!is.null(code)){
D[i,"tmt"] <- 0
}else{
D[i,"tmt"] <- -1
}
if (!is.null(int)){
D[i, (j+3):(2*j+2)] <- D[i,"tmt"]* D[i, 2:(j+1)]
}
}
}
if (identical(lossfunc(Imat.beta(D, beta), ...), lossfunc(Imat.beta(Dj, beta), ...))) 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, lossfunc, ...)))
D <- Djs[mindet][[1]]
return(D)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.