R/GLM_Exact_Design.R

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

#' rounding algorithm for generalized linear models
#'
#' @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 Integral_based TRUE or FALSE, whether or not integral-based EW D-optimality is used, FALSE indicates sample-based EW D-optimality is used.
#' @param b_matrix     matrix of bootstrapped or simulated parameter values.
#' @param joint_Func_b prior distribution function of model parameters
#' @param Lowerbounds  vector of lower ends of ranges of prior distribution for model parameters.
#' @param Upperbounds  vector of upper ends of ranges of prior distribution for model parameters.
#' @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 "logit", other links: "probit", "cloglog", "loglog", "cauchit", "log", "identity"
#'
#' @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
#' 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,Integral_based=TRUE,
#'  joint_Func_b=gjoint_b,Lowerbounds=paras_lowerbound, Upperbounds=paras_upperbound,
#'  delta2=0,L=1,N=100,hfunc=hfunc.temp,link=link.temp)

GLM_Exact_Design<-function(k.continuous,design_x,design_p,var_names=NULL,det.design,p,ForLion,bvec,Integral_based,
                           b_matrix,joint_Func_b, Lowerbounds, Upperbounds,delta2,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];
    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
    Xmat.temp0=matrix(0, m.design_mer, pnum);
    wvec.temp0=rep(0, m.design_mer);
    if(ForLion==TRUE){
      for(i in 1:m.design_mer) {
        if(is.null(dim(design_x))) {
          htemp=Xw_maineffects_self(x=xtemp_mer[i], b=bvec, link=link, h.func=hfunc);
        } else {
          htemp=Xw_maineffects_self(x=xtemp_mer[i,],b=bvec, link=link, h.func=hfunc);
        };
        Xmat.temp0[i,]=htemp$X;
        wvec.temp0[i]=htemp$w;
      }
    }else{
      for(i in 1:m.design_mer) {
        if(is.null(dim(design_x))) {
          if(Integral_based==TRUE)
          {htemp=EW_Xw_maineffects_self(x=xtemp_mer[i],Integral_based=Integral_based, joint_Func_b=joint_Func_b, Lowerbounds=Lowerbounds, Upperbounds=Upperbounds, link=link, h.func=hfunc);}
          else
          {htemp=EW_Xw_maineffects_self(x=xtemp_mer[i],Integral_based=Integral_based, b_matrix=b_matrix, link=link, h.func=hfunc)}
        } else {
          if(Integral_based==TRUE)
          {htemp=EW_Xw_maineffects_self(x=xtemp_mer[i,],Integral_based=Integral_based,joint_Func_b=joint_Func_b, Lowerbounds=Lowerbounds, Upperbounds=Upperbounds, link=link, h.func=hfunc);}
          else
          {htemp=EW_Xw_maineffects_self(x=xtemp_mer[i,],Integral_based=Integral_based,b_matrix=b_matrix, link=link, h.func=hfunc);}
        };
        Xmat.temp0[i,]=htemp$X;
        wvec.temp0[i]=htemp$E_w;
      }
    };
    F.mat=t(Xmat.temp0 * (ptemp_mer*wvec.temp0)) %*% Xmat.temp0;

    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;
    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])
        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
  Xmat.temp0=matrix(0, m0, pnum);
  wvec.temp0=rep(0, m0);
  if(ForLion==TRUE){
  for(i in 1:m0) {
      if(is.null(dim(combdesign))) {
        htemp=Xw_maineffects_self(x=combdesign[i], b=bvec, link=link, h.func=hfunc);
      } else {
        htemp=Xw_maineffects_self(x=combdesign[i,], b=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))) {
        if(Integral_based==TRUE)
          {htemp=EW_Xw_maineffects_self(x=combdesign[i],Integral_based=Integral_based, joint_Func_b=joint_Func_b, Lowerbounds=Lowerbounds, Upperbounds=Upperbounds, link=link, h.func=hfunc);}
        else
          {htemp=EW_Xw_maineffects_self(x=combdesign[i],Integral_based=Integral_based, b_matrix=b_matrix, link=link, h.func=hfunc)}
      } else {
        if(Integral_based==TRUE)
          {htemp=EW_Xw_maineffects_self(x=combdesign[i,],Integral_based=Integral_based,joint_Func_b=joint_Func_b, Lowerbounds=Lowerbounds, Upperbounds=Upperbounds, link=link, h.func=hfunc);}
        else
          {htemp=EW_Xw_maineffects_self(x=combdesign[i,],Integral_based=Integral_based,b_matrix=b_matrix, 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,var.names=var_names,rel.efficiency=rel.efficiency);
  class(output) <- "design_output"
  return(output)

}