Package 'CDsampling'

Title: 'CDsampling': Constraint Sampling in Paid Research Studies
Description: In the context of paid research studies and clinical trials, budget considerations and patient sampling from available populations are subject to inherent constraints. We introduce the 'CDsampling' package, which integrates optimal design theories within the framework of constrained sampling. This package offers the possibility to find both D-optimal approximate and exact allocations for samplings with or without constraints. Additionally, it provides functions to find constrained uniform sampling as a robust sampling strategy with limited model information. Our package offers functions for the computation of the Fisher information matrix under generalized linear models (including regular linear regression model) and multinomial logistic models.To demonstrate the applications, we also provide a simulated dataset and a real dataset embedded in the package. Yifei Huang, Liping Tong, and Jie Yang (2025)<doi:10.5705/ss.202022.0414>.
Authors: Yifei Huang [aut, cre], Liping Tong [aut], Jie Yang [aut]
Maintainer: Yifei Huang <[email protected]>
License: MIT + file LICENSE
Version: 0.1.2
Built: 2024-11-20 07:33:06 UTC
Source: CRAN

Help Index

Convert the approximate allocation (proportion) to exact allocation (integer) with bounded constraint (ni <= Ni)

Description

Convert the approximate allocation (proportion) to exact allocation (integer) with bounded constraint (ni <= Ni)

Usage

approxtoexact_constrained_func(
  n,
  w,
  m,
  beta = NULL,
  link = NULL,
  X = NULL,
  Fdet_func = Fdet_func_GLM,
  iset_func = NULL
)

Arguments

n

Sample size, must be a positive integer
w

Approximate allocation/proportion, must be a real-valued vector, can get from running liftone_constrained_GLM or liftone_constrained_MLM
m

The number of sampling groups
beta

Model parameter coefficients, default to be NULL for use in constrained uniform sampling
link

Link function of GLM or MLM, if used for GLM model (GLM_T is T), options are "identity", "logit", "probit", "cloglog", "loglog". If used for MLM (GLM_T is F), options are "continuation", "cumulative", "adjacent", and "baseline"
X

Design matrix of the model for GLM or MLM, default to be NULL for use in constrained uniform sampling
Fdet_func

determinant of Fisher information matrix function, Fdet_func can be self-defined, or use "Fdet_func_GLM", "Fdet_func_MLM" in the package, default is Fdet_func_GLM
iset_func

self-defined function for checking which index of sampling group fall within constraint if add 1 more subject (I set, see Algorithm 2 in Huang, Tong, Yang (2023)), two example functions are provided in the package, iset_func_trial and iset_func_trauma

Value

allocation is the exact allocation or integer value of the number of subjects sampled from the group

allocation.real is the proportion or the approximate allocation of the number of subjects sampled from the group

det.maximum is the maximum of |F| from the current exact allocation

Examples

beta = c(0, 3, 3, 3) #main effect model beta_0, beta_1, beta_21, beta_22
X.liftone=matrix(data=c(1,0,0,0,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,0,1,1,0,1), ncol=4, byrow=TRUE)
exact_design = approxtoexact_constrained_func(n=200, w=c(0.25, 0.20, 0.05, 0.50, 0.00, 0.00),
m=6, beta=beta, link='logit', X=X.liftone, Fdet_func=Fdet_func_GLM, iset_func=iset_func_trial)

Convert the approximate allocation (proportion) to exact allocation (integer) without constraint

Description

Convert the approximate allocation (proportion) to exact allocation (integer) without constraint

Usage

approxtoexact_func(n, w)

Arguments

n

Sample size, must be a positive integer
w

Approximate allocation/proportion, must be a real-valued vector, can get from running liftone_constrained_GLM or liftone_constrained_MLM

Value

allocation is the exact allocation or integer value of the number of subjects sampled from the group

Examples

exact_design = approxtoexact_func(n=600, w=c(0.2593526, 0.0000000, 0.0000000,
0.1565024, 0.2891565, 0.0000000, 0.0000000, 0.2949885))

Find (constrained) uniform exact allocation of the study for bounded design

Description

Find (constrained) uniform exact allocation of the study for bounded design

Usage

bounded_uniform(Ni, nsample)

Arguments

Ni

a vector with size m, upper bound for exact design of each category/stratification group, if unconstrained, use Inf vector, the function will return unbounded uniform allocation
nsample

a number, the sample size

Value

n is the constrained/unconstrained uniform exact allocation

Examples

bounded_uniform(Ni=c(50, 40, 10, 200, 150, 50), nsample=200)

Fisher information matrix of generalized linear model (GLM)

Description

Fisher information matrix of generalized linear model (GLM)

Usage

F_func_GLM(w, beta, X, link = "logit")

Arguments

w

allocation (can be exact or approximate)
beta

GLM model covariate coefficient
X

model matrix
link

link function, default"logit", choose from "logit", "cloglog", "loglog", "probit", and "identity"(for regular linear regression)

Value

object of class "matrix_output", Fisher information matrix given X and model parameter beta

Examples

w = c(1/3,1/3, 1/3)
beta = c(0.5, 0.5, 0.5)
X = matrix(data=c(1,-1,-1,1,-1,1,1,1,-1), byrow=TRUE, nrow=3)
F_func_GLM(w=w, beta=beta, X=X, link='logit')

The Fisher information matrix of multinomial logistic model (MLM)

Description

The Fisher information matrix of multinomial logistic model (MLM)

Usage

F_func_MLM(w, beta, X, link)

Arguments

w

allocation (can be exact or approximate)
beta

MLM model covariate coefficient
X

MLM model matrix
link

link function of Multinomial logistic regression model, options are "baseline", "cumulative", "adjacent", or "continuation"

Value

The Fisher information matrix of MLM model

Examples

w = rep(1/8, 8)
Xi=rep(0,5*12*8) #response levels * num of parameters * num of design points
dim(Xi)=c(5,12,8)
#design matrix
Xi[,,1] = rbind(c( 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,2] = rbind(c( 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,3] = rbind(c( 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,4] = rbind(c( 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,5] = rbind(c( 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,6] = rbind(c( 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,7] = rbind(c( 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,8] = rbind(c( 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
bvec_temp = c(-4.047, -2.225, -0.302, 1.386, 4.214, 3.519, 2.420,
1.284, -0.131, -0.376, -0.237, -0.120)
link_temp = "cumulative"

F_func_MLM(w=w, beta=bvec_temp, X=Xi, link=link_temp)

Determinant of Fisher information matrix for GLM

Description

Determinant of Fisher information matrix for GLM

Usage

Fdet_func_GLM(w, beta, X, link = "logit")

Arguments

w

allocation (can be exact or approximate)
beta

GLM model covariate coefficient
X

model matrix
link

link function, default"logit", choose from "logit", "cloglog", "loglog", "probit", and "identity"(for regular linear regression)

Value

the determinant of Fisher information matrix given X and model parameter beta

Examples

w = c(1/3,1/3, 1/3)
beta = c(0.5, 0.5, 0.5)
X = matrix(data=c(1,-1,-1,1,-1,1,1,1,-1), byrow=TRUE, nrow=3)
Fdet_func_GLM(w=w, beta=beta, X=X, link='logit')

Determinant of Fisher information matrix of multinomial logistic model (MLM)

Description

Determinant of Fisher information matrix of multinomial logistic model (MLM)

Usage

Fdet_func_MLM(w, beta, X, link)

Arguments

w

allocation (can be exact or approximate)
beta

MLM model covariate coefficient
X

MLM model matrix
link

link function of Multinomial logistic regression model, options are "baseline", "cumulative", "adjacent", or "continuation"

Value

Determinant of the Fisher information matrix of MLM model

Examples

w = rep(1/8, 8)
Xi=rep(0,5*12*8) #response levels * num of parameters * num of design points
dim(Xi)=c(5,12,8)
#design matrix
Xi[,,1] = rbind(c( 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,2] = rbind(c( 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,3] = rbind(c( 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,4] = rbind(c( 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,5] = rbind(c( 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,6] = rbind(c( 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,7] = rbind(c( 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,8] = rbind(c( 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
bvec_temp = c(-4.047, -2.225, -0.302, 1.386, 4.214, 3.519,
2.420, 1.284, -0.131, -0.376, -0.237, -0.120)
link_temp = "cumulative"

Fdet_func_MLM(w=w, beta=bvec_temp, X=Xi, link=link_temp)

Determinant function to be used for finding constrained uniform samplings

Description

Determinant function to be used for finding constrained uniform samplings

Usage

Fdet_func_unif(w, beta = NULL, X = NULL, link = NULL)

Arguments

w

allocation (can be exact or approximate)
beta

use NULL (default to be NULL)
X

use NULL (default to be NULL)
link

use NULL (default to be NULL)

Value

product of all allocation

Examples

Fdet_func_unif(w=c(0.2,0.2,0.2,0.2,0.2))

Generate Fisher information matrix F_x at a design point x_i for Multinomial logistic regression model

Description

Generate Fisher information matrix F_x at a design point x_i for Multinomial logistic regression model

Usage

Fi_func_MLM(X_x, beta, link)

Arguments

X_x

model matrix given design point x_i (for example, X_x = h.func(x_i), where h.func transforms a design point to a model matrix)
beta

parameter coefficients in the Multinomial logistic regression model, the order of coefficients in bvec and the order of design points in X_x should be consistent
link

link function of Multinomial logistic regression model, options are "baseline", "cumulative", "adjacent", or "continuation"

Value

F_x is the Fisher information matrix at design point x_i (with model matrix X_x);

U_x is a middle step matrix for calculation of F_x, details see Corollary 3.1 in Bu, X., Majumdar, D., & Yang, J. (2020). D-optimal designs for multinomial logistic models

Examples

X_x_temp = rbind(c( 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
bvec_temp = c(-4.047, -2.225, -0.302, 1.386, 4.214, 3.519,
2.420, 1.284, -0.131, -0.376, -0.237, -0.120)
link_temp = "cumulative"
Fi_func_MLM(X_x=X_x_temp, beta=bvec_temp, link=link_temp)

trauma_data example (see Huang, Tong, Yang (2023)) specific function for finding index set that if allocation of that index add "1", the new allocation still falls within the constraint Used in approxtoexact_constrained_func()

Description

trauma_data example (see Huang, Tong, Yang (2023)) specific function for finding index set that if allocation of that index add "1", the new allocation still falls within the constraint Used in approxtoexact_constrained_func()

Usage

iset_func_trauma(allocation)

Arguments

allocation

the exact allocation

Value

list of TRUE and FALSE, if TRUE, it means the allocation of this index will fall out of the constraint with more subject; if TURE, it means the allocation of this index can add more subjects

Examples

iset_func_trauma(allocation=c(50,30,10,10,100,100,200,10))

trial_data example (see Huang, Tong, Yang (2023)) specific function for finding index set that if allocation of that index add "1", the new allocation still falls within the constraint Used in approxtoexact_constrained_func()

Description

trial_data example (see Huang, Tong, Yang (2023)) specific function for finding index set that if allocation of that index add "1", the new allocation still falls within the constraint Used in approxtoexact_constrained_func()

Usage

iset_func_trial(allocation)

Arguments

allocation

the exact allocation

Value

list of TRUE and FALSE, if TRUE, it means the allocation of this index will fall out of the constraint with more subject; if TURE, it means the allocation of this index can add more subjects

Examples

iset_func_trial(allocation=c(50,30,10,100,100,40))

Find constrained D-optimal approximate design for generalized linear models (GLM)

Description

Find constrained D-optimal approximate design for generalized linear models (GLM)

Usage

liftone_constrained_GLM(
  X,
  W,
  g.con,
  g.dir,
  g.rhs,
  lower.bound,
  upper.bound,
  reltol = 1e-05,
  maxit = 500,
  random = TRUE,
  nram = 3,
  w00 = NULL,
  epsilon = 1e-12
)

Arguments

X

Model matrix, with nrow = num of design points and ncol = num of parameters
W

Diagonal of W matrix in Fisher information matrix, can be calculated from W_func_GLM in package
g.con

A matrix of numeric constraint coefficients, one row per constraint, on column per variable (to be used in as const.mat lp() and mat in Rglpk_solve_LP())
g.dir

Vector of character strings giving the direction of the constraint: each value should be one of "<," "<=," "=," "==," ">," or ">=". (In each pair the two values are identical.) to be used as const.dir in lp() and dir in Rglpk_solve_LP()
g.rhs

Vector of numeric values for the right-hand sides of the constraints. to be used as const.rhs in lp() and rhs in Rglpk_solve_LP().
lower.bound

A function to determine lower bound r_i1 in Step 3 of Constrained lift-one algorithm from Yifei, H., Liping, T., Yang, J. (2023) Constrained D-optimal design for paid research study
upper.bound

A function to determine upper bound r_i2 in Step 3 of Constrained lift-one algorithm from Yifei, H., Liping, T., Yang, J. (2023) Constrained D-optimal design for paid research study
reltol

The relative convergence tolerance, default value 1e-5
maxit

The maximum number of iterations, default value 500
random

TRUE or FALSE, if TRUE then the function will run with additional "nram" number of random initial points, default to be TRUE
nram

When random == TRUE, the function will generate nram number of initial points, default is 3
w00

Specified initial design proportion; default to be NULL, this will generate a random initial design
epsilon

A very small number, for comparison of >0, <0, ==0, to reduce errors, default 1e-12

Value

w is the approximate D-optimal design

w0 is the initial design used to get optimal design w

maximum is the maximized |F| value

itmax is the number of iterations

convergence is TRUE or FALSE, if TRUE means the reported design is converged

deriv.ans is the derivative from step 6 of constrained lift-one algorithm

gmax is the maximum g function in step 8 of constrained lift-one algorithm

reason is the lift-one loops break reason, either "all derivatives <=0" or "gmax <=0"

Examples

#Example 6 in Section 3.4 of Yifei, H., Liping, T., Yang, J. (2025)
#Constrained D-optimal design for paid research study

#main effect model beta_0, beta_1, beta_21, beta_22
beta = c(0, -0.1, -0.5, -2)

#gives the 6 categories (0,0,0), (0,1,0),(0,0,1),(1,0,0),(1,1,0),(1,0,1)
X.liftone=matrix(data=c(1,0,0,0,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,0,1,1,0,1),
ncol=4, byrow=TRUE)

#calculate W matrix based on beta's under logit link
W_matrix=W_func_GLM(X= X.liftone, b=beta)

m=6 #number of categories
nsample = 200
rc = c(50, 40, 10, 200, 150, 50)/nsample
g.con = matrix(0,nrow=(2*m+1), ncol=m)
g.con[1,] = rep(1, m)
g.con[2:(m+1),] = diag(m)
g.con[(m+2):(2*m+1), ] = diag(m)
g.dir = c("==", rep("<=", m), rep(">=", m))
g.rhs = c(1, rc, rep(0, m))

lower.bound=function(i, w){
  nsample = 200
  rc = c(50, 40, 10, 200, 150, 50)/nsample
  m=length(w) #num of categories
  temp = rep(0,m)
  temp[w>0]=1-pmin(1,rc[w>0])*(1-w[i])/w[w>0];
  temp[i]=0;
  max(0,temp);
}
upper.bound=function(i, w){
  nsample = 200
  rc = c(50, 40, 10, 200, 150, 50)/nsample
  m=length(w) #num of categories
  rc[i];
  min(1,rc[i]);
}

approximate_design = liftone_constrained_GLM(X=X.liftone, W=W_matrix,
g.con=g.con, g.dir=g.dir, g.rhs=g.rhs, lower.bound=lower.bound,
upper.bound=upper.bound, reltol=1e-10, maxit=100, random=TRUE, nram=4,
w00=NULL, epsilon = 1e-8)

Find constrained D-optimal designs for Multinomial Logit Models (MLM)

Description

Find constrained D-optimal designs for Multinomial Logit Models (MLM)

Usage

liftone_constrained_MLM(
  m,
  p,
  Xi,
  J,
  beta,
  lower.bound,
  upper.bound,
  g.con,
  g.dir,
  g.rhs,
  w00 = NULL,
  link = "cumulative",
  Fi.func = Fi_func_MLM,
  reltol = 1e-05,
  maxit = 500,
  delta = 1e-06,
  epsilon = 1e-08,
  random = TRUE,
  nram = 3
)

Arguments

m

The number of design points; it is usually the number of combinations of all the stratification factors
p

The number of parameters in the MLM model
Xi

Model matrix, a J by p by m 3D array of predictors for separate response category at all design points(input to determine ppo,npo,po)
J

The number of response levels
beta

A p*1 vector, parameter coefficients for MLM, the order of beta should be consistent with Xi
lower.bound

A function to determine lower bound r_i1 in Step 3 of Constrained lift-one algorithm from Yifei, H., Liping, T., Yang, J. (2023) Constrained D-optimal design for paid research study
upper.bound

A function to determine upper bound r_i2 in Step 3 of Constrained lift-one algorithm from Yifei, H., Liping, T., Yang, J. (2023) Constrained D-optimal design for paid research study
g.con

A matrix of numeric constraint coefficients, one row per constraint, on column per variable (to be used in as const.mat lp() and mat in Rglpk_solve_LP())
g.dir

Vector of character strings giving the direction of the constraint: each value should be one of "<," "<=," "=," "==," ">," or ">=". (In each pair the two values are identical.) to be used as const.dir in lp() and dir in Rglpk_solve_LP()
g.rhs

Vector of numeric values for the right-hand sides of the constraints. to be used as const.rhs in lp() and rhs in Rglpk_solve_LP()
w00

Specified initial design proportion; default to be NULL, this will generate a random initial design
link

Link function of MLM, default to be "cumulative", options from "continuation", "cumulative", "adjacent", and "baseline"
Fi.func

A function for calculating Fisher information at a specific design point, default to be Fi_func_MLM function in the package
reltol

The relative convergence tolerance, default value 1e-5
maxit

The maximum number of iterations, default value 500
delta

A very small number, used in alpha_star calculation, default to be 1e-6.
epsilon

A very small number, for comparison of >0, <0, ==0, to reduce errors, default 1e-12
random

TRUE or FALSE, if TRUE then the function will run with additional "nram" number of random initial points, default to be TRUE
nram

When random == TRUE, the function will generate nram number of initial points, default is 3

Value

w is the approximate D-optimal design

w0 is the initial design used to get optimal design w

Maximum is the maximized |F| value

itmax is the number of iterations

convergence is TRUE or FALSE, if TRUE means the reported design is converged

deriv.ans is the derivative from step 6 of constrained lift-one algorithm

gmax is the maximum g function in step 8 of constrained lift-one algorithm

reason is the lift-one loops break reason, either "all derivatives <=0" or "gmax <=0"

Examples

#Example 8 of Trauma data example in Yifei, H., Liping, T., Yang, J. (2025)
#Constrained D-optimal design for paid research study

J = 5    # number of categories,  >= 3
p = 12    # number of parameters
m = 8    # number of design points
nsample=600 #collect 600 samples finally from the 802 subjects
lower.bound <- function(i, w0){
n = 600
constraint = c(392,410)
if(i <= 4){
  a.lower <- (sum(w0[5:8])-(constraint[2]/n)*(1-w0[i]))/(sum(w0[5:8]))
}
else{
  a.lower <- (sum(w0[1:4])-(constraint[1]/n)*(1-w0[i]))/(sum(w0[1:4]))
}
a.lower
}
upper.bound <- function(i, w0){
  n = 600
  constraint = c(392,410)
  if(i <= 4){
    b.upper <- ((constraint[1]/n)*(1-w0[i]) - (sum(w0[1:4])-w0[i]))/(1-sum(w0[1:4]))
  }
  else{
    b.upper <- ((constraint[2]/n)*(1-w0[i]) - (sum(w0[5:8])-w0[i]))/(1-sum(w0[5:8]))
  }
  b.upper
}


constraint = c(392,410)
g.con = matrix(0,nrow=length(constraint)+1+m, ncol=m)
g.con[2:3,] = matrix(data=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1), ncol = m, byrow=TRUE)
g.con[1,] = rep(1, m)
g.con[4:(length(constraint)+1+m), ] = diag(1, nrow=m)
g.dir = c("==", "<=","<=", rep(">=",m))
g.rhs = c(1, ifelse((constraint/nsample<1),constraint/nsample,1), rep(0, m))
Xi=rep(0,J*p*m)
dim(Xi)=c(J,p,m)
Xi[,,1] = rbind(c( 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,2] = rbind(c( 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,3] = rbind(c( 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,4] = rbind(c( 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 0, 0, 0, 0, 0, 0, 0),
              c( 0, 0, 0, 0, 0, 0, 1, 4, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,5] = rbind(c( 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,6] = rbind(c( 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,7] = rbind(c( 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Xi[,,8] = rbind(c( 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
thetavec = c(-4.3050, -0.0744,  4.3053, -2.3334, -0.3290, 3.4773,
-0.1675, -0.3609, 2.7358, 1.2935, -0.1612, 1.4899)
set.seed(123)
liftone_constrained_MLM(m=m, p=p, Xi=Xi, J=J, beta=thetavec, lower.bound=lower.bound,
upper.bound=upper.bound, g.con=g.con,g.dir=g.dir, g.rhs=g.rhs, w00=NULL,
link='cumulative', Fi.func = Fi_func_MLM, reltol=1e-5, maxit=500,
delta = 1e-6, epsilon=1e-8, random=TRUE, nram=1)

Unconstrained lift-one algorithm to find D-optimal allocations for GLM

Description

Unconstrained lift-one algorithm to find D-optimal allocations for GLM

Usage

liftone_GLM(
  X,
  W,
  reltol = 1e-05,
  maxit = 500,
  random = TRUE,
  nram = 3,
  w00 = NULL
)

Arguments

X

Model matrix, with nrow = num of design points and ncol = num of parameters
W

Diagonal of W matrix in Fisher information matrix, can be calculated from W_func_GLM in package
reltol

The relative convergence tolerance, default value 1e-5
maxit

The maximum number of iterations, default value 500
random

TRUE or FALSE, if TRUE then the function will run with additional "nram" number of random initial points, default to be TRUE
nram

When random == TRUE, the function will generate nram number of initial points, default is 3
w00

Specified initial design proportion; default to be NULL, this will generate a random initial design

Value

w is the approximate D-optimal design

w0 is the initial design used to get optimal design w

Maximum is the maximized |F| value

itmax is the number of iterations

convergence is TRUE or FALSE, if TRUE means the reported design is converged

Examples

beta = c(0.5, 0.5, 0.5)
X = matrix(data=c(1,-1,-1,1,-1,1,1,1,-1), byrow=TRUE, nrow=3)
W_matrix = W_func_GLM(X=X, beta=beta)
w00 = c(1/6, 1/6, 2/3)
approximate_design = liftone_GLM(X=X, W=W_matrix, reltol=1e-10, maxit=100,
random=FALSE, nram=3, w00=w00)

Unconstrained lift-one algorithm to find D-optimal allocations for MLM

Description

Unconstrained lift-one algorithm to find D-optimal allocations for MLM

Usage

liftone_MLM(
  m,
  p,
  Xi,
  J,
  beta,
  link = "continuation",
  Fi.func = Fi_func_MLM,
  reltol = 1e-05,
  maxit = 500,
  w00 = NULL,
  random = TRUE,
  nram = 3
)

Arguments

m

The number of design points; it is usually the number of combinations of all the stratification factors
p

The number of parameters in the MLM model
Xi

Model matrix, a J by p by m 3D array of predictors for separate response category at all design points(input to determine ppo,npo,po)
J

The number of response levels
beta

A p*1 vector, parameter coefficients for MLM, the order of beta should be consistent with Xi
link

Link function of MLM, default to be "cumulative", options from "continuation", "cumulative", "adjacent", and "baseline"
Fi.func

A function for calculating Fisher information at a specific design point, default to be Fi_func_MLM function in the package
reltol

The relative convergence tolerance, default value 1e-5
maxit

The maximum number of iterations, default value 500
w00

Specified initial design proportion; default to be NULL, this will generate a random initial design
random

TRUE or FALSE, if TRUE then the function will run with additional "nram" number of random initial points, default to be TRUE
nram

When random == TRUE, the function will generate nram number of initial points, default is 3

Value

w is the approximate D-optimal design

w0 is the initial design used to get optimal design

Maximum is the maximized |F| value

itmax is the number of iterations

convergence is TRUE or FALSE, if TRUE means the reported design is converged

Examples

J = 5    # number of categories,  >= 3
p = 12    # number of parameters
m = 8    # number of design points
Xi=rep(0,J*p*m) #J*p*m=5*12*8
dim(Xi)=c(J,p,m)
#design matrix
Xi[,,1] = rbind(c( 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,2] = rbind(c( 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,3] = rbind(c( 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,4] = rbind(c( 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,5] = rbind(c( 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,6] = rbind(c( 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,7] = rbind(c( 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1),
               c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

Xi[,,8] = rbind(c( 1, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 1, 4, 1, 0, 0, 0, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 1, 4, 1, 0, 0, 0),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 1),
                c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

thetavec = c(-4.047, -0.131, 4.214, -2.225, -0.376, 3.519,
-0.302, -0.237,  2.420, 1.386,  -0.120,  1.284)

liftone_MLM(m=m, p=p, Xi=Xi, J=J, beta=thetavec, link = "cumulative",
Fi.func=Fi_func_MLM, reltol=1e-5, maxit=5000, w00=NULL, random=TRUE, nram=3)

Print Method for list_output Objects

Description

Custom print method for objects of class 'list_output'.

Usage

## S3 method for class 'list_output'
print(x, ...)

Arguments

x

An object of class 'list_output'.
...

Additional arguments (ignored).

Value

Invisibly returns 'x' (the input object).

Print Method for matrix_list Objects

Description

Custom print method for objects of class 'matrix_list'.

Usage

## S3 method for class 'matrix_list'
print(x, ...)

Arguments

x

An object of class 'matrix_list'.
...

Additional arguments (ignored).

Value

Invisibly returns 'x' (the input object).

Print Method for matrix_output Objects

Description

Custom print method for objects of class 'matrix_output'.

Usage

## S3 method for class 'matrix_output'
print(x, ...)

Arguments

x

An object of class 'matrix_output'.
...

Additional arguments (ignored).

Value

Invisibly returns 'x' (the input object).

Trauma data with multinomial response

Description

The data frame saves data from the trauma trial data from Chuang-Stein and Agresti (1997).

Usage

trauma_data

Format

A data frame with 802 rows and 5 variables:

Severity

severity of the trauma symptoms, mild or moderate/severe

Dose

dose levels applied to the patients, 4 levels, placebo, low, medium and high

Label

stratification group in terms of severity and dose

Outcome

treatment outcome, 5 levels, death, vegetative state, major disability, minor disability and good recovery

ID

patient ID, 1-802

Source

Chuang-Stein and Agresti (1997)

Examples

data(trauma_data) #lazy loading

Generated clinical trial data with binary response

Description

Generated with logistic regression model:

logit(P(Yij=1∣genderi,agei1,agei2))=3∗genderi+3∗agei1+3∗agei2logit(P(Y_{ij} = 1 | gender_i, age_{i1}, age_{i2})) = 3*gender_i +3*age_{i1} +3*age_{i2}

The data frame can be used to run GLM clinical trial example in Huang, Tong, Yang (2023)

Usage

trial_data

Format

A data frame with 500 rows and 6 variables:

gender

gender of the patients

age_1

1 or 0, whether or not the patient belongs to 18-25 age group

age_2

1 or 0, whether or not the patient belongs to 26-64 age group

label

stratification group in terms of gender and age, 1 to 6

Y

treatment effective or not, Y=1 means treatment is effective to the patient

ID

patient ID, 1-500

Source

Generated pseudo clinical trial data to serve as an example.

Examples

data(trial_data) #lazy loading

Calculate the diagonal elements nu of Fisher information matrix

Description

Calculate the diagonal elements nu of Fisher information matrix

Usage

W_func_GLM(X, beta, link = "logit")

Arguments

X

Model matrix
beta

Parameters of GLM model
link

GLM link function, default is "logit", options are "logit", "probit", "cloglog", "loglog", "identity", identity is the same as ordinary linear regression

Value

the diagonal element nu of GLM Fisher information matrix, can be used as w in liftone_constrained_GLM

Examples

beta = c(0, 3, 3, 3) #main effect model beta_0, beta_1, beta_21, beta_22
#gives the 6 categories (0,0,0), (0,1,0),(0,0,1),(1,0,0),(1,1,0),(1,0,1)
X.liftone=matrix(data=c(1,0,0,0,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,0,1,1,0,1), ncol=4, byrow=TRUE)
#calculate diagonal elements of W based on beta's under logit link
W=W_func_GLM(X= X.liftone, beta=beta, link="logit")