---
title: "Introduction to CDsampling"
author: "Y.H."
output: 
  rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to CDsampling}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```



```{css, echo=FALSE}
.note-box {
  border: 1px solid #ccc;
  background: white;
  padding: 10px;
  margin: 10px 0;
  border-radius: 5px;
}
```

# Table of Contents

- [Computation of Fisher information matrix](#Computation_Fisher)
  - [Example 1: GLM Fisher information matrix](#example1_GLM_Fisher)
  - [Example 2: MLM Fisher information matrix](#example2_MLM_Fisher)
- [Applications in paid research studies](#Applications)
  - [Example 3: trial_data & constrained sampling with GLM](#example3_trial_data)
  - [Example 4: trauma_data & constrained sampling with MLM](#example4_trauma_data)
- [References](#References)

```{r setup}
library(CDsampling)
```

In the context of paid research studies and clinical trials, budget considerations and patient sampling from available populations are subject to inherent constraints. **CDsampling** 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 sampling with or without constraints. 

* Additionally, it provides functions to find constrained uniform sampling as a robust sampling strategy with limited model information. 

* It also provides tool for computation of Fisher information matrix of the generalized linear models (GLMs) including regular linear regression model and the multinomial logistic models (MLMs).

* Two datasets are embedded in the package for application examples.

<a id="Computation_Fisher"></a>

# Computation of Fisher information matrix

<a id="example1_GLM_Fisher"></a> 

## Example 1: GLM Fisher information matrix 
Consider a research study with a simple logistic regression model $$\log(\frac{\mu_i}{1-\mu_i}) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}$$ where $\mu_i = E(Y_i\mid {\mathbf x}_i)$, ${\mathbf x}_i = (x_{i1}, x_{i2})^\top \in \{(-1, -1), (-1, +1), (+1, -1)\}$ and parameters $\boldsymbol \beta = (\beta_0, \beta_1, \beta_2) = (0.5, 0.5, 0.5)$. In this example, we have $m=3$ design points $(x_{i1}, x_{i2})^\top \in \{(-1, -1), (-1, +1), (+1, -1)\}$, the design matrix $\mathbf X$ (with the first column for the intercept): 
$$\begin{bmatrix} 
1 & -1 & -1 \\
1 & -1 & 1 \\
1 & 1 & -1\\
\end{bmatrix}$$.

To calculate Fisher information matrix of the design with GLM, we can use *F_func_GLM( )* in the package with input of approximate allocation $w$, coefficients $\boldsymbol \beta$, and design matrix $\mathbf X$. 

```{r}
beta = c(0.5, 0.5, 0.5) #coefficients
X = matrix(data=c(1,-1,-1,1,-1,1,1,1,-1), byrow=TRUE, nrow=3) #design matrix
w = c(1/3,1/3,1/3) #approximate allocation
CDsampling::F_func_GLM(w=w, beta=beta, X=X, link='logit')
```

<a id="example2_MLM_Fisher"></a>

## Example 2: MLM Fisher information matrix
Consider a research study with cumulative non-proportional odds multinomial logit model with $J=5$ response levels and covariates $(x_{i1}, x_{i2})=\{(1,0),(2,0),(3,0), (4,0),(1,1),(2,1),(3,1),(4,1)\}$. The model can be written as:$$\log(\frac{\pi_{i1}}{\pi_{i2}+\dots+\pi_{i5}}) = \beta_{11}+\beta_{12}x_{i1}+\beta_{13}x_{i2}$$
    $$\log(\frac{\pi_{i1}+\pi_{i2}}{\pi_{i3}+\pi_{i4}+\pi_{i5}}) = \beta_{21}+\beta_{22}x_{i1}+\beta_{23}x_{i2}$$
    $$\log(\frac{\pi_{i1}+\pi_{i2}+\pi_{i3}}{\pi_{i4}+\pi_{i5}}) = \beta_{31}+\beta_{32}x_{i1}+\beta_{33}x_{i2}$$
    $$\log(\frac{\pi_{i1}+\dots+\pi_{i4}}{\pi_{i5}}) = \beta_{41}+\beta_{42}x_{i1}+\beta_{43}x_{i2}$$
    where $i=1,\dots,8$. We have $m=8$ design points and $p=12$ parameters. We assume the parameters $\boldsymbol \beta = (\beta_{11},  \beta_{12},  \beta_{13}, \beta_{21},  \beta_{22},  \beta_{23},  \beta_{31},  \beta_{32},  \beta_{33}, \beta_{41}, \beta_{42},  \beta_{43})^\top = (-4.047, -0.131, 4.214, -2.225, -0.376,$  $3.519, -0.302, -0.237,  2.420, 1.386,  -0.120, 1.284)^\top$. The approximate allocation for the eight design points is $\mathbf w = (1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8)^\top$. 
    
To calculate the Fisher information matrix with MLM, we can use the *F_func_MLM( )* with input of approximate allocation $w$, covariate coefficients $\boldsymbol \beta$, design matrix $\mathbf X$, and multinomial logit model (cumulative for this example). 

The design matrix incorporates all $8$ design points of covariates $(x_{i1}, x_{i2})$ specified by the cumulative logit model's four equations. For example, when $x_{i1}=1$, $x_{i2}=0$, the design matrix takes the following format: 

$$\begin{bmatrix} 
1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{bmatrix}$$.

```{r}
J=5; p=12; m=8; #response levels; num of parameters; num of design points
beta = 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) 
Xi=rep(0,J*p*m) #design matrix   
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))
 alloc = rep(1/8,m) #approximate allocation
 CDsampling::F_func_MLM(w=alloc, beta=beta, X=Xi, link='cumulative')
```
    
<a id="Applications"></a>

# Applications in paid research studies

<a id="example3_trial_data"></a>

## Example 3: trial_data & constrained sampling with GLM

The *trial_data* is a simulated dataset containing information for $N=500$ volunteers with gender, age, and final efficacy information. The covariates considered in this example are: 

1). **Gender:** 

- $0$: for female 

- $1$ for male


2). **Age group:** 

- $age\_1 = 0$ and $age\_2 = 0$: ages $18\sim25$

- $age\_1 = 1$ and $age\_2 = 0$: ages $26\sim64$

- $age\_1 = 0$ and $age\_2 = 1$: ages $65$ and above. 

There are $m=6$ design points, that is, the number of combinations of gender and age groups $(x_{gender\_i}, x_{age\_1i}, x_{age\_2i})$: 

1). $(0,0,0)$: Female, 18-25 $N_1=50$)

2). $(0,1,0)$: Female, 26-64 ($N_2=40$)  

3). $(0,0,1)$: Female, 65+ ($N_3=10$)

4). $(1,0,0)$: Male, 18-25 ($N_4=200$)

5). $(1,1,0)$: Male, 26-64 ($N_5=150$) 

6). $(1,0,1)$: Male, 65+ ($N_6=50$) 


Suppose that a sample of $n=200$ participants is required due to budget limit.Our goal is to find the constrained D-optimal allocation $(w_1, w_2, \dots, w_6)$ with feasible allocation $$S = \{(w_1, \ldots, w_m)^T \in S_0 \mid n w_i \leq N_i, i=1, \ldots, m\}.$$


We use constrained lift-one algorithm *liftone_constrained_GLM( )* to find the locally D-optimal approximate sampling allocations with the input of design matrix $X$, $\mathbf W$ matrix which is the result returned from *W_func_GLM( )*, constraints setup (g.con, g.dir, and g.rhs), and boundaries in searching for lift-one weight (step 3 in the constrained lift-one algorithm, see reference). 


We consider the logistic regression model for $j=1,\dots,m$, $i=1,\dots,n_j$ with $\boldsymbol \beta=(\beta_{1}, \beta_{21}, \beta_{22})= (0,3,3,3)$: 

\begin{equation}\label{eq:trial_logistic_model}
{\rm logit} \{P(Y_{ij}=1 \mid x_{gender\_i}, x_{age\_1i}, x_{age\_2i})\} = \beta_0 + \beta_1 x_{gender\_i} + \beta_{21} x_{age\_1i} + \beta_{22} x_{age\_2i} 
\end{equation}

Use the following R codes to define the coefficients, sample size, and design matrix: 

```{r}
beta = c(0, 3, 3, 3) #coefficients 
#design matrix X
X=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) 
nsample=200 #sample size
```

To run the *liftone_constrained_GLM( )* function, we also need to the $\mathbf W$ matrix from the calculation of Fisher information matrix, we can use the *W_func_GLM( )* function in the package: 

```{r}
W_matrix=CDsampling::W_func_GLM(X=X, b=beta, link="logit") #W matrix
```

Lastly, we also need to define the constraints (number of patients from different gender and age group) and boundaries for constrained sampling (please see the reference for details of lower bound $r_{i1}$ and upper bound $r_{i2}$): 

```{r}
rc = c(50, 40, 10, 200, 150, 50)/nsample #constraints for each subgroup
m = 6
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])
}
```

We can define an optional subgroups label for the constrained approximate allocations: 
```{r}
label = c("F, 18-25", "F, 26-64", "F, >=65", "M, 18-25", "M, 26-64", "M, >=65")
```


Now, we can run the constrained lift-one algorithm to find optimal approximate allocations for the defined subgroups: 

```{r}
set.seed(2025)
approximate_design = CDsampling::liftone_constrained_GLM(X=X, W=W_matrix, g.con=g.con, g.dir=g.dir, 
                                               g.rhs=g.rhs, lower.bound=lower.bound, 
                                               upper.bound=upper.bound, label=label, reltol=1e-10, 
                                               maxit=100, random=TRUE, nram=4, w00=NULL, 
                                               epsilon=1e-8)

print(approximate_design)
```

The output contains several key components: 

- $w$: the converged D-optimal approximate 

- $w_0$: random initial weights used for the optimization  

- *maximum*: the achieved maximum determinant of the Fisher information matrix

- *reason*: criteria for lift-one loop termination with either "all derivative <=0", and "gmax <=0"  (see reference).

To find the exact allocation (integer value of allocation), we can use the *approxtoexact_constrained_func( )* with input of sample size $n$, approximate allocation found by constrained lift-one algorithm, number of design point $m$, coefficients *beta*, link type of GLM ("logit" in this example), **Fdet_func=Fdet_func_GLM**, and design matrix $X$: 
```{r}
exact_design = CDsampling::approxtoexact_constrained_func(n=200, w=approximate_design$w, m=6, 
                                                           beta=beta, link='logit', X=X, 
                                                           Fdet_func=Fdet_func_GLM, 
                                                          iset_func=iset_func_trial, label=label) 

print(exact_design)
```

The *allocation* represents the final D-optimal exact allocation of the constrained sampling, while allocation.real is the approximate allocation found in the previous step. 

<div class="note-box">
**Note:**
When prior estimates of beta coefficients are unavailable, consider using **EW D-optimal allocations** as an alternative option. EW D-optimal method requires prior distributions of beta coefficients, which substitutes the coefficients value with expectation of prior distributions in the computaion. (Please see reference for details). 
</div>

If there is no knowledge of beta coefficients, we recommend using constrained uniform design with *approxtoexact_constrained_func( )*. We adjust **Fdet_func =Fdet_func_unif** and the approximate application $w$ to be $1$ subject in each subgroup. In this example, we have $n=200$ sample thus $w=(1/200, 1/200, 1/200, 1/200, 1/200, 1/200)$: 

```{r}
w00 = rep(1/200, 6)
unif_design = CDsampling::approxtoexact_constrained_func(n=200, w=w00, m=6, beta=NULL, 
  link=NULL, X=NULL, Fdet_func=Fdet_func_unif, iset_func=iset_func_trial)

print(unif_design)
```

The *allocation* represents the constrained uniform exact allocation, while allocation.real is the initial allocation we provided. 


<a id="example4_trauma_data"></a>

## Example 4: trauma_data & constrained sampling with MLM
 
The *trauma_data* contains information on $N=802$ trauma patients with five ordered treatment outcomes: 

1). Death ($1$) 

2). Vegetative state ($2$)

3). Major disability ($3$)

4). Minor disability ($4$)

5). Good recovery ($5$).


We consider two key covariates:

1). **Treatment dose levels**: 

  - $x_{i1} = 1$: Placebo
  
  - $x_{i1} = 2$: Low dose 
  
  - $x_{i1} = 3$: Medium dose 
  
  - $x_{i1} = 4$: High dose. 
  

2). **Trauma symptoms**:  

  - $x_{i2} = 0$: mild symptoms ($392$ mild symptoms patients)
  
  - $x_{i2} = 1$: moderate /severe symptoms ($410$ moderate/severe symptoms patients),


This creates $m=8$ design points from the combination of all covariates levels. 
  
  In this study, we plan to enroll $n=600$ patients and the collection of all feasible approximate allocation is: $$S=\{(w_1, \ldots, w_8)^\top \in S_0\mid n(w_1+w_2+w_3+w_4) \leq 392, n(w_5+w_6+w_7+w_8) \leq 410\}.$$ 
  
  The cumulative non-proportional logit model is specified as:$$\log(\frac{\pi_{i1}}{\pi_{i2}+\dots+\pi_{i5}}) = \beta_{11}+\beta_{12}x_{i1}+\beta_{13}x_{i2}$$
    $$\log(\frac{\pi_{i1}+\pi_{i2}}{\pi_{i3}+\pi_{i4}+\pi_{i5}}) = \beta_{21}+\beta_{22}x_{i1}+\beta_{23}x_{i2}$$
    $$\log(\frac{\pi_{i1}+\pi_{i2}+\pi_{i3}}{\pi_{i4}+\pi_{i5}}) = \beta_{31}+\beta_{32}x_{i1}+\beta_{33}x_{i2}$$
    $$\log(\frac{\pi_{i1}+\dots+\pi_{i4}}{\pi_{i5}}) = \beta_{41}+\beta_{42}x_{i1}+\beta_{43}x_{i2}$$

with the assumed parameter vector $\boldsymbol\beta = (\hat\beta_{11}, \hat\beta_{12}, \hat\beta_{13}, \hat\beta_{21}, \hat\beta_{22}, \hat\beta_{23}, \hat\beta_{31}, \hat\beta_{32}, \hat\beta_{33}, \hat\beta_{41}, \hat\beta_{42}, \hat\beta_{43})^\top  = (-4.047, -0.131, 4.214, -2.225, -0.376,\\ 3.519, -0.302, -0.237,  2.420, 1.386,  -0.120,  1.284)^\top$. 


For this example, we can reuse the set-up of design matrix and beta coefficients from [Example 2: MLM Fisher information matrix](#example2_MLM_Fisher): 

```{r}
J=5; p=12; m=8; #response levels; num of parameters; num of design points
beta = 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)
Xi=rep(0,J*p*m) #design matrix
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))
```


Similarly, we define the sample size, constraints, and boundaries using the following R codes: 

```{r}
nsample=600 #sample size
constraint = c(392, 410)  #mild:severe

#boundaries r_i1 and r_i2 in constrained lift-one
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
}

#constraints of the example
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))

```


Define an optional subgroups label for the constrained approximate allocations: 
```{r}
label = c("Placebo-Mild", "Low-Mild", "Medium-Mild", "High-Mild", "Placebo-Severe", "Low-Severe", "Medium-Severe", "High-Severe")
```

Then, we can use the constrained lift-one algorithm to find the locally D-optimal approximate sampling with *liftone_constrained_MLM( )* function in the package: 

```{r}
set.seed(123)
approx_design = CDsampling::liftone_constrained_MLM(m=m, p=p, Xi=Xi, J=J, beta=beta, 
                                       lower.bound=lower.bound, upper.bound=upper.bound, 
                                        g.con=g.con, g.dir=g.dir, g.rhs=g.rhs, label=label, w00=NULL, 
                                        link='cumulative', Fi.func = Fi_func_MLM, 
                                        reltol=1e-5, maxit=500, delta = 1e-6, 
                                        epsilon=1e-8, random=TRUE, nram=1)

print(approx_design)
```

Similar to [Example 3: trial_data & constrained sampling with GLM](#example3_trial_data), $w$ reports the D-optimal approximate allocation, while $w_0$ denotes the random initial weights used in the optimization procedure. *maximum* reports the determinant of Fisher information matrix, and *reason* is the termination criteria for the lift-one algorithm. 

To convert the approximate allocation to exact allocation, we use the function *approxtoexact_constrained_func( )* in the package, with **Fdet_func=Fdet_func_MLM**: 

```{r}
exact_design = CDsampling::approxtoexact_constrained_func(n=600, w=approx_design$w, m=8, 
 beta=beta, link='cumulative', X=Xi, label=label, Fdet_func=Fdet_func_MLM, 
 iset_func=iset_func_trauma)

print(exact_design)
```


If there is no coefficients information, we can find the constrained uniform sampling with *approxtoexact_constrained_func( )* using **Fdet_func=Fdet_func_unif**: 

```{r}
iset_func_trauma <- function(allocation){
  iset = rep(1,8)
  if(sum(allocation[1:4])>=392){iset[1:4]=0}
  if(sum(allocation[5:8])>=410){iset[5:8]=0}
  return(iset)
 }

unif_design = CDsampling::approxtoexact_constrained_func(n=600, w=rep(1/600,8), m=8, beta=NULL, 
 link=NULL, X=NULL, label=label, Fdet_func=Fdet_func_unif, iset_func=iset_func_trauma)

unif_design
```



<div class="note-box">
**Note:**  
The reported determinant of the Fisher information matrix (`det.maximum`) is *model-specific* and depends on the choice of `Fdet_func`. For meaningful comparisons:   
- Use `Fdet_func_GLM( )` for GLMs with the derived allocations  
- Use `Fdet_func_MLM( )` for MLMs with the derived allocations 

**Do not compare** values across functions.  
</div>

<a id="References"></a>

# References
- Huang Y., Tong L., Yang J. (2025). *Constrained D-optimal Design for Paid Research Study*, Statistica Sinica, DOI: 10.5705/ss.202022.0414.



