R/GLM_Exact_Design.R

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

#' Approximation to exact design algorithm for generalized linear 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 assumed parameter values of model parameters beta, same length of h(y)
#' @param joint_Func_b The prior joint probability distribution of the parameters
#' @param Lowerbounds The lower limit of the prior distribution for each parameter
#' @param Upperbounds The upper limit of the prior distribution for each parameter
#' @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 "logit", other links: "probit", "cloglog", "loglog", "cauchit", "log", "identity"
#'
#' @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
#' k.continuous=1
#' design_x=matrix(c(25, -1, -1,-1, -1 ,
#'                  25, -1, -1, -1, 1,
#'                  25, -1, -1, 1, -1,
#'                  25, -1, -1, 1, 1,
#'                  25, -1, 1, -1, -1,
#'                  25, -1, 1, -1, 1,
#'                  25, -1, 1, 1, -1,
#'                  25, -1, 1, 1, 1,
#'                  25, 1, -1, 1, -1,
#'                  25, 1, 1, -1, -1,
#'                  25, 1, 1, -1, 1,
#'                  25, 1, 1, 1, -1,
#'                  25, 1, 1, 1, 1,
#'                  38.9479, -1, 1, 1, -1,
#'                  34.0229, -1, 1, -1, -1,
#'                  35.4049, -1, 1, -1, 1,
#'                  37.1960, -1, -1, 1, -1,
#'                  33.0884, -1, 1, 1, 1),nrow=18,ncol=5,byrow = TRUE)
#' hfunc.temp = function(y) {c(y,y[4]*y[5],1);};   # y -> h(y)=(y1,y2,y3,y4,y5,y4*y5,1)
#' link.temp="logit"
#'design_p=c(0.0848, 0.0875, 0.0410, 0.0856, 0.0690, 0.0515,
#'           0.0901, 0.0845, 0.0743, 0.0356, 0.0621, 0.0443,
#'           0.0090, 0.0794, 0.0157, 0.0380, 0.0455, 0.0022)
#'det.design=4.552715e-06
#' paras_lowerbound<-c(0.25,1,-0.3,-0.3,0.1,0.35,-8.0)
#' paras_upperbound<-c(0.45,2,-0.1,0.0,0.4,0.45,-7.0)
#'  gjoint_b<- function(x) {
#'  Func_b<-1/(prod(paras_upperbound-paras_lowerbound))
#'  ##the prior distributions are follow uniform distribution
#' return(Func_b)
#' }
#'  GLM_Exact_Design(k.continuous=k.continuous,design_x=design_x,
#'  design_p=design_p,det.design=det.design,p=7,ForLion=FALSE,joint_Func_b=gjoint_b,
#'  Lowerbounds=paras_lowerbound, Upperbounds=paras_upperbound,rel.diff=0,L=1,
#'  N=100,hfunc=hfunc.temp,link=link.temp)

GLM_Exact_Design<-function(k.continuous,design_x,design_p,det.design,p,ForLion,bvec,joint_Func_b, Lowerbounds, Upperbounds,rel.diff,L,N,hfunc,link){

  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
  Xmat.temp0=matrix(0, m0, pnum);
  wvec.temp0=rep(0, m0);
  if(ForLion==T){
  for(i in 1:m0) {
      if(is.null(dim(combdesign))) {
        htemp=Xw_maineffects_self(x=combdesign[i], bvec, link=link, h.func=hfunc);
      } else {
        htemp=Xw_maineffects_self(x=combdesign[i,], bvec, link=link, h.func=hfunc);
      };
      Xmat.temp0[i,]=htemp$X;
      wvec.temp0[i]=htemp$w;
    }
    }else{
      for(i in 1:m0) {
      if(is.null(dim(combdesign))) {
        htemp=EW_Xw_maineffects_self(x=combdesign[i],joint_Func_b, Lowerbounds, Upperbounds, link=link, h.func=hfunc);
      } else {
        htemp=EW_Xw_maineffects_self(x=combdesign[i,],joint_Func_b, Lowerbounds, Upperbounds, link=link, h.func=hfunc);
      };
      Xmat.temp0[i,]=htemp$X;
      wvec.temp0[i]=htemp$E_w;
    }
  };
  det.temp0=det(t(Xmat.temp0 * ((exact_ni/N)*wvec.temp0)) %*% Xmat.temp0);    # |X^T W X|
  rel.efficiency<-(det.temp0/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)

}