CDsampling (version 0.1.6)

liftone_constrained_MLM: 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,
  label = NULL,
  w00 = NULL,
  link = "cumulative",
  Fi.func = Fi_func_MLM,
  reltol = 1e-05,
  maxit = 500,
  delta = 1e-06,
  epsilon = 1e-08,
  random = TRUE,
  nram = 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"

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()

label

A vector of text strings for subgroups' names, default value NULL

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

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)