R/MLM_Exact_Design.R

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

#' rounding algorithm for multinomial logit models
#'
#' @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 the corresponding approximate allocation
#' @param var_names Names for the design factors. Must have the same length asfactor.level. Defaults to "X1", "X2", ...
#' @param det.design the determinant of the approximate design
#' @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 sampled parameter values of beta
#' @param delta2 points with distance less than that will be merged
#' @param L vector: rounding factors
#' @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             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,
#' delta2=1,L=0.5,N=1000,hfunc=hfunc.temp,link=link.temp)
MLM_Exact_Design<-function(J, k.continuous,design_x,design_p,var_names=NULL,det.design,p,ForLion,bvec,bvec_matrix,delta2,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];
    m.design=length(xtemp);
  } else {
    xtemp=design_x[design_p>0,];       # list of design points
    m.design=nrow(xtemp);  # number of design points
  };

  dtemp=as.matrix(stats::dist(xtemp));
  diag(dtemp)=Inf;
  dmin=min(dtemp)
  while ((dmin<delta2)) { #merge closest two neighbors
    #before merging two closest points, save the current state of the design
    x.design_old=xtemp
    p.design_old=ptemp
    m.design_old=m.design

    i1=which.min(apply(dtemp,1,min));
    i2=which.min(dtemp[i1,]);
    if(is.null(dim(design_x))) {
      xnew=(p.design_old[i1]*x.design_old[i1]+p.design_old[i2]*x.design_old[i2])/(p.design_old[i1]+p.design_old[i2]);
      ptemp_mer=c(p.design_old[-c(i1,i2)],p.design_old[i1]+p.design_old[i2]);
      xtemp_mer=c(x.design_old[-c(i1,i2)],xnew);
    } else {
      xnew=(p.design_old[i1]*x.design_old[i1,]+p.design_old[i2]*x.design_old[i2,])/(p.design_old[i1]+p.design_old[i2]);
      ptemp_mer=c(p.design_old[-c(i1,i2)],p.design_old[i1]+p.design_old[i2]);
      xtemp_mer=rbind(x.design_old[-c(i1,i2),],xnew);
    }

    m.design_mer=m.design_old-1

    pnum=p
    X.mat = rep(0,J*pnum*m.design_mer);
    dim(X.mat)=c(J, pnum, m.design_mer)  # initial model matrix X
    for(i in 1:m.design_mer) {
      if(is.null(dim(design_x))){htemp=hfunc(xtemp_mer[i]);}else{htemp=hfunc(xtemp_mer[i,]);}
      X.mat[,,i]=htemp;
    }

    F.mat = 0
    Fi = vector()
    for(i in 1:m.design_mer) {
      if(ForLion==TRUE){
        new_Fi = Fi_MLM_func(X.mat[, ,i], bvec=bvec, link=link)
        Fi = append(Fi, list(new_Fi$F_x))
        F.mat = F.mat + (ptemp_mer[i])*new_Fi$F_x
      } else {
        new_Fi = EW_Fi_MLM_func(X.mat[, ,i], bvec_matrix=bvec_matrix, link)
        Fi = append(Fi, list(new_Fi$F_x))
        F.mat = F.mat + (ptemp_mer[i])*new_Fi$F_x
      };
    }
    eigen_values<-eigen(F.mat)
    min_engenvalue<-min(eigen_values$values)

    if(min_engenvalue<=1e-6){
      xtemp=x.design_old
      ptemp=p.design_old
      m.design=m.design_old
      break
    }else{
      xtemp=xtemp_mer
      ptemp=ptemp_mer
      m.design=m.design_mer
    }
    dtemp=as.matrix(stats::dist(xtemp));
    diag(dtemp)=Inf;
    atemp=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])
        Li<-L[i]
        n2<-length(tconti)
        for(j in 1:n2){
          if (abs((tconti[j]%/%Li)*Li-tconti[j])<Li/2){
            xinew[j]<-(tconti[j]%/%Li)*Li
          }
          else
          {
            xinew[j]<-(tconti[j]%/%Li)*Li+Li
          }
        }
        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==TRUE){
      new_Fi = Fi_MLM_func(X.mat[, ,i], bvec=bvec, link=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=bvec_matrix, link=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,var.names=var_names,rel.efficiency=rel.efficiency);
    class(output) <- "design_output"
    return(output)

}