R/MLM_Exact_Design.R

In ForLion: 'ForLion' Algorithm to Find D-Optimal Designs for Experiments

#' Approximation to exact design algorithm for multinomial logit model
#'
#' @param J number of response levels in the multinomial logit model
#' @param k.continuous number of continuous factors
#' @param design_x the matrix with rows indicating design point which we got from the approximate design
#' @param design_p D-optimal approximate allocation
#' @param det.design the determinant of D-optimal approximate allocation
#' @param p number of parameters
#' @param ForLion TRUE or FALSE, TRUE: this approximate design was generated by ForLion algorithm,
#'                FALSE: this approximate was generated by EW ForLion algorithm
#' @param bvec  If ForLion==TRUE assumed parameter values of model parameters beta, same length of h(y)
#' @param bvec_matrix If ForLion==FALSE the matrix of the bootstrap parameter values of beta
#' @param rel.diff points with distance less than that will be merged
#' @param L rounding factor
#' @param N total number of observations
#' @param hfunc function for obtaining model matrix h(y) for given design point y, y has to follow the same order as n.factor
#' @param link link function, default "continuation", other choices "baseline", "cumulative", and "adjacent"
#'
#' @return x.design matrix with rows indicating design point
#' @return ni.design EW D-optimal or D-optimal exact allocation
#' @return rel.efficiency relative efficiency of the Exact and Approximate Designs
#' @export
#'
#' @examples
#' J=3
#' k.continuous=1
#' design_x<-c(0.0000,103.5451,149.2355)
#' design_p<-c(0.2027, 0.3981, 0.3992)
#' det.design=54016609
#' p=5
#' theta = c(-1.935, -0.02642, 0.0003174, -9.159, 0.06386)
#' hfunc.temp = function(y){
#'    matrix(data=c(1,y,y*y,0,0,0,0,0,1,y,0,0,0,0,0), nrow=3,
#'             ncol=5, byrow=TRUE)
#' }
#' link.temp = "continuation"
#' MLM_Exact_Design(J=J, k.continuous=k.continuous,design_x=design_x,
#' design_p=design_p,det.design=det.design,p=p,ForLion=TRUE,bvec=theta,
#' rel.diff=1,L=0.5,N=1000,hfunc=hfunc.temp,link=link.temp)
MLM_Exact_Design<-function(J, k.continuous,design_x,design_p,det.design,p,ForLion,bvec,bvec_matrix,rel.diff,L,N,hfunc,link){
  ##input
  ##     design_x: the design that want to combine the close design point
  ##     design_p: the proportion for the designx
  ##    xconsist: divide the design

  #####First stage: Combine the close design point
  # Calculate distance matrix
  design_p=design_p/sum(design_p);          # normalize design_p
  ptemp=design_p[design_p>0];        # remove 0 items
  if(is.null(dim(design_x))) {
    xtemp=design_x[design_p>0];
    d.factor=length(design_x);
  } else {
    xtemp=design_x[design_p>0,];       # list of design points
    d.factor=dim(design_x)[2];  # number of factors
  };
  dtemp=as.matrix(stats::dist(xtemp));
  diag(dtemp)=Inf;
  dmin=min(dtemp)
  while ((dmin<rel.diff)) { #merge closest two neighbors
    i1=which.min(apply(dtemp,1,min));
    i2=which.min(dtemp[i1,]);
    if(is.null(dim(design_x))) {
      xnew=(ptemp[i1]*xtemp[i1]+ptemp[i2]*xtemp[i2])/(ptemp[i1]+ptemp[i2]);
      ptemp=c(ptemp[-c(i1,i2)],ptemp[i1]+ptemp[i2]);
      xtemp=c(xtemp[-c(i1,i2)],xnew);
    } else {
      xnew=(ptemp[i1]*xtemp[i1,]+ptemp[i2]*xtemp[i2,])/(ptemp[i1]+ptemp[i2]);
      ptemp=c(ptemp[-c(i1,i2)],ptemp[i1]+ptemp[i2]);
      xtemp=rbind(xtemp[-c(i1,i2),],xnew);
    }
    dtemp=as.matrix(stats::dist(xtemp));
    diag(dtemp)=Inf;
    dmin=min(dtemp);
  };

  combdesign<-xtemp
  combp<-ptemp
  if (k.continuous!=0){
    # ##Second stage: change the value of stack force to 2.5L
    for(i in 1:k.continuous){
      xinew<-NULL
      if(is.null(dim(combdesign))){
        tconti<-as.vector(combdesign)
        n2<-length(tconti)
        for(j in 1:n2){
          if (abs((tconti[j]%/%L)*L-tconti[j])<L/2){
            xinew[j]<-(tconti[j]%/%L)*L
          }else
          {
            xinew[j]<-(tconti[j]%/%L)*L+L
          }
        }
        combdesign<-as.vector(xinew)
      } else{
        tconti<-as.vector(combdesign[ ,i])
        n2<-length(tconti)
        for(j in 1:n2){
          if (abs((tconti[j]%/%L)*L-tconti[j])<L/2){
            xinew[j]<-(tconti[j]%/%L)*L
          }
          else
          {
            xinew[j]<-(tconti[j]%/%L)*L+L
          }
        }
        combdesign[ ,i]<-xinew
      }
    }
  }

  ##Third stage: tranform the approximate design to the exact design
  approximatetoexact.self <- function(n, p, constraints=NULL) {
    # p[1:m] is a real-valued allocation, 0 <= p_i <= constraints[i]/n, sum_i p_i = 1
    # constraints[1:m] n_i <= constraints[i]
    m=length(p);
    if(is.null(constraints)) constraints=rep(n,m);
    ans=floor(n*p);
    k=n-sum(ans);
    if(k>0) {
      stemp=(n*p-ans)*((ans+1) <= constraints);
      otemp=order(-stemp)[1:k];
      ans[otemp]=ans[otemp]+1;
    };
    ans;
  }
  exact<-approximatetoexact.self(N,p=combp,constraints = NULL)
  ##Fourth stage: calculate the determinat of the fisher information matrix and the cost efficieny
  exact_ni<-as.vector(exact)

  pnum=p
  m0=length(exact_ni);     # number of original design points
  X.mat = rep(0,J*pnum*m0);
  dim(X.mat)=c(J, pnum, m0)  # initial model matrix X
  for(i in 1:m0) {
    if(is.null(ncol(combdesign))){htemp=hfunc(combdesign[i]);}else{htemp=hfunc(combdesign[i,]);}
    X.mat[,,i]=htemp;
  }

  F.mat = 0
  Fi = vector()
  for(i in 1:m0) {
    if(ForLion==T){
      new_Fi = Fi_MLM_func(X.mat[, ,i], bvec, link)
      Fi = append(Fi, list(new_Fi$F_x))
      F.mat = F.mat + (exact_ni[i]/N)*new_Fi$F_x
    } else {
      new_Fi = EW_Fi_MLM_func(X.mat[, ,i], bvec_matrix, link)
      Fi = append(Fi, list(new_Fi$F_x))
      F.mat = F.mat + (exact_ni[i]/N)*new_Fi$F_x
    };
  }
    det_exact=det(F.mat)
    rel.efficiency<-(det_exact/det.design)^(1/p)
    #list(x.design=combdesign,ni.design=exact_ni,rel.efficiency=rel.efficiency)
    #define S3 class
    output<-list(x.factor=combdesign,p=combp, ni.design=exact_ni,rel.efficiency=rel.efficiency);
    class(output) <- "design_output"
    return(output)

}