#################################################################################### # # # SELC: An R Function for the Sequential Elimination of Level Combinations # # Algorithm # # # # Inputs: # # initialdesign : matrix of information that contains the initial # # design and the corresponding responses. # # Levels of each factor must be positive integers. # # forbiddenarray : matrix of information that contains the forbidden # # array design points. # # Levels of each factor must be positive integers. # # If no forbiddenarray exists, then specify # # forbiddenarray=cbind(0,0,...,0) # # strength : Strength of the forbidden array # # order : Order of the forbidden array # # level_true : Number of levels of each factor. Must have the form: # # level_true = c(1,2,3,...) # # direction : direction of desired response # # 1 = maximum is desired # # 0 = minimum is desired # # number_of_offspring : Number of desired offspring to generate # # seed : Seed value for random operations. Default value # # is clock time # # # # Output: # # An object is output from the function that contains 2 items: # # newruns : A matrix of design points that contains the new # # offspring # # currentforbiddenarray : A matrix that contains the design points of the # # forbidden array plus the design points specified # # by strength # # # # Example code to call SELC function: # # >initialdesign = as.matrix(read.table("c:/initialdesign.txt")) # # >forbiddenarray = as.matrix(read.table("c:/forbiddenarray.txt")) # # >source("c:/SELC.R") # # >args(SELC) # # >next_gen = SELC(initialdesign,forbiddenarray,2,2,c(5,34,241),1,5) # # # # References: # # Mandal, A., Wu, C.F.J., and Johnson, K. (2006). SELC: Sequential elimination # # of level combinations by means of modified genetic algorithms. # # Technometrics, 48(2), 273-283. # # # # Mandal, A., Johnson, K., Wu, C.F.J., and Bornemeier, D. (2007). Identifying # # promising compounds in drug discovery: Genetic algorithms and some new # # statistical techniques. Journal of Chemical Information and Modeling, 47(3), # # 981-988. (DOI: 10.1021/ci600556v) Highlighted in Drug Discovery News # # Online: http://www.drugdiscoverynews.com/index.php?newsarticle=1318 # # # # Written by Kjell Johnson, Abhyuday Mandal and Tan Ding # # March 2008 # # Version 1.0 # #################################################################################### #################################################################################### # Functions needed for SELC # #################################################################################### # Given model matrix D, modelmatrix generates the design matrix X # Design matrix contains main effects, quadratic effects and 2fi's (linear by linear only). # Notes : The model matrix does not contain an intercept # Linear/quadratic system is used modelmatrix=function(D) { k = ncol(D) l = nrow(D) X = c() # First the main effects - standardized to have norm 1 p=2 for (i in 1:k) { X = cbind(X,D[,i]/(sum(abs(D[,i])^p)^(1/p))) } # Then the quadratic effects X = cbind(X,X^2) # Then the interaction effects c = 1 #Counter f = factorial(k)/(2*factorial(k-2)) Y = matrix(0,l,f) for (i in 1:(k-1)) { for (j in (i+1):k) { Y[,c]=X[,i]*X[,j] c=c+1 } } X=cbind(X,Y) return(X) } # Function to evaluate the statistical significance of each factor in the model evalsig=function(x,y) { x = modelmatrix(x) if (qr(x)\$rank < ncol(x)) { print("More runs required in initial design") sigfactor=matrix(0,1,(ncol(x)-1)) }else { g = lm(y ~ x) sg = summary(g) sgcoef = sg\$coef sgcoef = as.matrix(sgcoef) pvalues = sgcoef[,4] n = length(pvalues) sigfactor=matrix(0,1,n) for (i in 1:n) { if (pvalues[i]<0.05) sigfactor[i] = 1 } } return(sigfactor) } # Function to choose the elements of x with weighted probability choose_with_prob=function(x,number) { k = ncol(x) l = nrow(x) x = matrix(x,k*l,1) x = x*(x > 0) prob = x/sum(x) cumulative_prob = cumsum(prob) candidates = c() while (NROW(candidates) < number) { u = runif(1) temp = (u >= cumulative_prob)* matrix(1:nrow(prob)) candidates = rbind(candidates, max(temp)+1) candidates = unique(candidates) } return(candidates) } # Function to perform crossover crossover=function(p1,p2) { crossover_location = floor(runif(1)*length(p1)+1) p1new = c( p1[1:crossover_location-1], p2[crossover_location:NROW(p1)]) p2new = c( p2[1:crossover_location-1], p1[crossover_location:NROW(p1)]) newoff = rbind(p1new, p2new) return(newoff) } # Function to evaluate the significance of the interaction terms # if interactiton_number defines the number of significant interactions # The nonzero entries of the (k-1)th column and (k)th row of the matrix # facmat corresponds to the 2fi's of factor k is_there_imp_2fi=function(k,sigfactor,num_main_effects) { n=length(sigfactor) sigfactor = matrix(sigfactor,n,1) facmat = matrix(0,num_main_effects-1,num_main_effects-1) count = 1 for ( i in 1:(num_main_effects-1)) { for ( j in i:(num_main_effects-1)) { facmat[i,j] = count count = count + 1 } } if (k == 1) { values = facmat[1,] } else if (k==num_main_effects) { values = facmat[,k-1] } else { values = c(facmat[,k-1], facmat[k,]) } values = values[values>0] + 2*num_main_effects r=length(values) values = matrix(values,r,1) interaction_number = max(values*sigfactor[values]) return(interaction_number) } # Function to identify the parent factors identify=function(i,sigfactor,num_main_effects) { # i = effect number # z = main effect/quadratic effect/parents value = i # this is the indicator of presence of significant 2fi. if ( i <= num_main_effects) { index = i value = is_there_imp_2fi(index,sigfactor,num_main_effects) } else if ( i <= 2*num_main_effects) { index = i-num_main_effects value = is_there_imp_2fi(i,sigfactor,num_main_effects) } if (value == 0) { z=index } else { i = value int_num=i-2*num_main_effects count=0 count1=int_num num_ints=0 while (count1>0) { count=count+1 count1=count1-(num_main_effects-count) } count1=0 for (j in 1:(count-1) ) { count1=count1+(num_main_effects-j) } count1=int_num-count1 count1=count+count1 z = cbind(count, count1) } return(z) } # Function to calculate the mutation probability ith factor prob_main_effect=function(x,y,i,level_true) { num_level_true = level_true[i] level_observed = sort(unique(x[,i])) num_level_observed =length(level_observed) ymean = c(1:num_level_true*0) for (j in 1:num_level_observed) { ymean[level_observed[j]] = mean(y[x[,i]==level_observed[j]]) } ymean = ymean * (ymean > 0) ymean = ymean/sum(ymean) uniform_prob = 1/num_level_true * c(1:num_level_true*0+1) alpha = (num_level_observed / num_level_true) prob = alpha*ymean + (1-alpha)*uniform_prob return(prob) } # Function to calculate mutation for factor i with weight scheme given by the vector prob one_factor_mutation=function(p1,x,y,i,level_true) { prob1 = prob_main_effect(x,y,i,level_true) cumulative_prob1 = cumsum(prob1) u = runif(1) temp = (u >= cumulative_prob1) * (1:length(prob1)) p1[i] = max(temp) + 1 return(p1) } # Function to calculate the mutation probability for 2fis prob_2fi=function(x,y,i,j,level_true) { yprob = matrix(0,level_true[i],level_true[j]) for (ik in 1:level_true[i]) { for (jk in 1:level_true[j]) { if ((length(y[x[,i]==ik & x[,j]==jk])) > 0) { yprob[ik,jk] = mean(y[x[,i]==ik & x[,j]==jk]) } } } yprob = yprob*(yprob>0) yprob = yprob/sum(sum(yprob)) uniform_prob = 1/(level_true[i]*level_true[j]) * matrix(1,level_true[i],level_true[j]) alpha = sum(sum(yprob > 0))/(level_true[i]*level_true[j]) prob = alpha*yprob + (1-alpha)*uniform_prob return(prob) } # Function to calcuate the mutation probability for 2fis with weight scheme given by the vector prob. two_factor_mutation=function(p1,x,y,i,j,level_true) { prob = prob_2fi(x,y,i,j,level_true) cumulative_prob = cumsum(prob) u = runif(1) temp = (u < cumulative_prob) * (1:length(prob)) index = min(temp[temp>0]) number_of_rows = level_true[i] number_of_col = level_true[j] row = index - (ceiling(index/number_of_rows)-1)*number_of_rows; col = ceiling(index/number_of_rows); p1[i] = row p1[j] = col return(p1) } # Function to perform mutation mutation=function(x,y,sigfactor,p1,num_main_effects,level_true) { mutation_location = floor(runif(1)*length(p1)+1) if (sigfactor[mutation_location] == 1) { i = identify(mutation_location,sigfactor,num_main_effects) if (length(i)==1) { p1 = one_factor_mutation(p1,x,y,i,level_true) } else { p1 = two_factor_mutation(p1,x,y,i[1],i[2],level_true) } } else { p1[mutation_location] = floor(runif(1)*level_true[mutation_location]+1) } return (p1) } # Function to check element of one matrix with another check = function(oldmatrix,newrow,order) { oldmatrix=as.matrix(oldmatrix) newrow=as.matrix(newrow) n = NROW(oldmatrix) check = 1 for (i in 1:n) { if ( sum(newrow == oldmatrix[i,]) >= order ) check = 0 break } return(check) } # Function to check the input arguments for feasibility feasibility = function(initialdesign,forbiddenarray,strength,order,level_true,direction) { proceed = 1; a = nrow(initialdesign) b = ncol(initialdesign) c = nrow(forbiddenarray) d = ncol(forbiddenarray) f = ncol(t(as.matrix(level_true))) if (b != (d+1)) { print("Number of factors in the forbidden array do not match the number of factors in the initial design") proceed = 0 }else if (f != d) { print("Number of factors in level_true must be the same as the number of levels in the forbidden array") proceed = 0 }else { m1=apply(initialdesign[,1:b-1],2,max) m2=apply(forbiddenarray,2,max) m3=apply(rbind(m1,m2),2,max) if (sum(level_true >= m3) != d) { print("There exits levels in initialdesign or forbiddenarray that are greater than the levels specified in level_true") proceed = 0 } } if (strength > (a-2)) { print("Strength must be less than or equal to:") print(a-2) proceed = 0 } if (order < 1) { print("Order must be at least 1") proceed = 0 } if (order > (b-1)) { print("Order cannot be greater than:") print(b-1) proceed = 0 } if (direction!=0 & direction!=1) { print("Direction must be 0 or 1") proceed = 0 } return(proceed) } #################################################################################### # # # Main function to implement SELC # # # #################################################################################### SELC=function(initialdesign, forbiddenarray, strength, order, level_true, direction, number_of_offspring, seed_num=-1) { proceed = feasibility(initialdesign,forbiddenarray,strength,order,level_true,direction) if (proceed == 1) { if (seed_num>0) { set.seed(seed_num) } m = ncol(initialdesign) num_main_effects = m-1 x = as.matrix(initialdesign[,-m]) y = as.matrix(initialdesign[,m]) if (direction==1) { y = (y - min(y))/(max(y) - min(y)); }else if (direction==0) { y = (max(y) - y)/(max(y) - min(y)); }else { print("Check the direction") return } yc = y-mean(y) if (length(level_true)!=num_main_effects) { print("Check the number of factors") newrums=matrix() new=matrix() return } for (i in 1:num_main_effects) { if (length(unique(x[,i]))>level_true[i]) { print("Check the number of factors") newruns=matrix() newforbiddenarray=matrix() return } } #identify significant factors sigfactor = evalsig(x,yc) newruns = c() # Add new runs to the forbidden array based on the initial design yy = unique(sort(y)) index = matrix(0*1:length(y)) ijk = 1 while ( sum(index) < strength) { index = index + (y==yy[ijk]) ijk = ijk+1 } index = index*matrix(1*1:length(y)) index = index[index>0] newforbiddenarray = x[index,] currentforbiddenarray = rbind(forbiddenarray,newforbiddenarray) #Reproduction while (length(newruns) < (m-1)*number_of_offspring) { population_size = 2 candidates = choose_with_prob(y,population_size) p1 = x[candidates[1],] p2 = x[candidates[2],] # Crossover newoff = crossover(p1,p2) p1new = newoff[1,] p2new = newoff[2,] # Mutation p1new = mutation(x,y,sigfactor,p1new,num_main_effects,level_true) p2new = mutation(x,y,sigfactor,p2new,num_main_effects,level_true) # Check on eligibility indicator11 = check(newforbiddenarray,p1new,order) indicator12 = check(forbiddenarray,p1new,num_main_effects) indicator21 = check(newforbiddenarray,p2new,order) indicator22 = check(forbiddenarray,p2new,num_main_effects) if ((indicator11+indicator12)==2) {indicator1=1} else {indicator1=0} if ((indicator21+indicator22)==2) {indicator2=1} else {indicator2=0} if (indicator1 == 1) { indicator3 = check(x,p1new,num_main_effects) if (indicator3 == 1) { newruns = rbind(newruns, p1new) } } if (indicator2 == 1) { indicator4 = check(x,p2new,num_main_effects) if (indicator4 == 1) { newruns = rbind(newruns, p2new) } } } newruns = newruns[1:number_of_offspring,] return(list(newruns=newruns,currentforbiddenarray=currentforbiddenarray)) }else { print("WARNING: FUNCTION STOPPED BECAUSE OF ERRORS") } }