DETMAX for optimal design

Giới Thiệu

Optimal design hay nói chính xác là Optimal experimental design là một nhánh nhỏ của Statistical Design of Experiment. nội dung này xoay quanh ứng dụng của thống kê học trong những nghiên cứu khoa học. Lấy một ví dụ thực nghiệm như sau:

Ta muốn tìm hiểu về liều lượng bức xạ (radiation-dosage) trên khả năng phục hồi của khối u (tumour reduction). Gọi lượng bức xạ đó là \(X\) (treatment), \(X\) sẽ có tất cả 8 mức độ khác nhau (levels), lần lượt là \((x_1,x_2,...,x_3)\). Giả sử có \(N\) bệnh nhân tham gia vào thí nghiệm này, tổng số bệnh nhân sẽ được chia lần lượt thành 8 nhóm, với sỉ số của từng nhóm là \((n_1,n_2,...,n_8)\). N bệnh nhân được đánh số lần lượt là \((x_1,x_2,...,x_N)\). Câu hỏi đặc ra là phải chia số lượng bệnh nhân vào mỗi nhóm như thế nào để kết quả thu được là hiệu quả nhất. ở đây ta có N bệnh nhân, chia thành 8 nhóm, mỗi nhóm có ít nhất 1 người, đây là bài toán tổ hợp có hoàn lại (combination with repetitions). như vậy ta có \[ n_1 +n_2+...+n_8 = N \]

Như vậy có tất cả \(\binom{8+N-1}{N} = \binom{N+7}{N}\), nếu \(N=10\), thì ta có \(\binom{10+7}{10} = \binom{17}{10} = 8436285\) cách khác nhau. mỗi cách gọi là một design. Gọi \(\xi_i\) là design thứ i, ta có sơ đồ như sau

\[ \xi_i =\begin{Bmatrix} x_1 & x_2 & x_3 & x_4 & x_5 & x_6 & x_7 & x_8 \\ n_1 & n_2 & n_3 & n_4 & n_5 & n_6 & n_7 & n_8 \end{Bmatrix}, \]

Câu hỏi đặc ra là chia như thế nào để được design mang lại hiệu quả cao nhất?

models

Chúng ta sẽ sử dụng model bậc hai (quadratic model) để kiểm định. quadratic model co1 dạng tổng quát như sau:

\[ y = \beta_0+\beta_1x+\beta_2x^2+\epsilon \]

Giả sử 8 treatment đó lần lượt là

## [1] -1.0000000 -0.7142857 -0.4285714 -0.1428571  0.1428571  0.4285714
## [7]  0.7142857  1.0000000

Như vậy ta có ma trận design (design matrix) là

##      [,1]       [,2]       [,3]
## [1,]    1 -1.0000000 1.00000000
## [2,]    1 -0.7142857 0.51020408
## [3,]    1 -0.4285714 0.18367347
## [4,]    1 -0.1428571 0.02040816
## [5,]    1  0.1428571 0.02040816
## [6,]    1  0.4285714 0.18367347
## [7,]    1  0.7142857 0.51020408
## [8,]    1  1.0000000 1.00000000

Thuật toán

DETMAX là thuật toán sử dụng trong optimal design để tìm design mang lại hiệu quả nhất theo một tiêu chuẩn \(\Psi\) nào đó. ở đây ta sử dụng \(\Psi_D(.) = |(X^tX)^{-1}|\). Như vậy, ta sẽ phân chia như thế nào để có được một design với tỉ số của hàm số mất mác \(\Psi_D\) ở trên. Dưới đây là Rcode của DETMAX

library(magrittr)
x = seq(-1,1,length.out = 8)
X = vapply(x,.%>%`^`(0:2),numeric(3))%>% t()%>%
{.[sample(nrow(.),nrow(.)),]}

####
detmax = function(X){
  n.point = nrow(X)
  deter = function(X) det(solve(t(X)%*%X))
  outer = TRUE
  iter.out=0
  
  while(outer){
    failset = c(deter(X))
    d = apply(X,1,.%>%{t(.)%*%solve(t(X)%*%X)%*%.})
    option = sample(c('add','delete'),1)
    opt = ifelse(option=='add',which.max(d),which.min(d))
    if(option == "delete"){
      X = X[-opt,]
    }else{
      X = rbind(X,X[opt,])
    }
    obj = deter(X)
    iter.out = iter.out+1
    iter.in = 0
    repeat{
      iter.in = iter.in+1
      ####
      if(nrow(X)<n.point){
        if(sum(obj<failset,na.rm = TRUE)>=1){
          option = 'add'
          failset = rep(NA,length(failset))
          failset = c(failset,obj)
        }else{
          option = 'delete'
          failset = c(failset,obj)
        }
      }else{
        if(sum(obj<failset,na.rm = TRUE)>=1){
          option = 'delete'
          failset = rep(NA,length(failset)+1)
          failset = c(failset,obj)
        }else{
          option = 'add'
          failset = c(failset,obj)
        }
      }
      ######
      if(length(failset)==6 && nrow(X)<n.point) option = "add"
      if(length(failset)==6 && nrow(X)>n.point) option = 'delete'
      #####
      d = apply(X,1,.%>%{t(.)%*%solve(t(X)%*%X)%*%.})
      
      if(option== 'add'){
        choice = which(max(d)-d==0)%>% {.[sample(length(.),1)]}
        X = rbind(X,X[choice,])
      }else{
        choice = which(min(d)-d==0)%>% {.[sample(length(.),1)]}
        X = X[-choice,]
      }
      obj = deter(X)
      if(nrow(X)==n.point) break
    } # inner repeat
    d = X%*%solve(t(X)%*%X)%*%t(X)
    m = vapply(1:nrow(d),function(j){
      vapply(1:ncol(d),function(i){
        d[i,i]-d[j,j]-d[i,i]*d[j,j]+d[i,j]^2
      },numeric(1))
    },numeric(ncol(d)))
    if(max(m)<10^-6){outer = FALSE}
    #print(iter.in)
    #print(iter.out)
  }# outer repeat
  return(X)
}

detmax(X)%>% `[`(,2)%>% table()

## .
##                -1 0.142857142857143                 1 
##                 3                 2                 3

Như vậy ta chia như thế nào để tỉ lệ giữa 3 treatments tại giá trị \((-1,-0.143, 1)\) là \(3:2:3\).