stepwiseCopy<-function(
	formula=NULL, 
	data=NULL, 
	include = NULL, 
	selection = "bidirection", #c("forward", "backward", "bidirection", "score"), 
	select = "SL", #c("AIC", "AICc", "BIC", "CP", "HQ", "HQc", "Rsq", "adjRsq", "SL", "SBC"), 
	sle = 0.15, 
    sls = 0.15, 
	multivarStat = "Pillai", #c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"), 
	weights = NULL, 
    best = NULL
) {
	#convert any logical (boolean) to factor 

	#match args 
    #selection <- match.arg(selection)
    #select <- match.arg(select)
    #multivarStat <- match.arg(multivarStat)
	
	#cant do SL for all possible subsets 
    if (selection == "score" & select == "SL") {
        stop("select = 'SL' is not allowed when specifing selection = 'score'")
    }
	#parse formula 
    if (class(formula) != "formula") {
        stop("class of formula object isn't 'formula'")
    } else {
        termForm <- terms(formula, data = data)
        vars <- as.character(attr(termForm, "variables"))[-1]
        yName <- vars[attr(termForm, "response")]
        xName <- attr(termForm, "term.labels")
        if (attr(termForm, "intercept") == 0) {
            intercept <- "0"
        } else {
            intercept <- "1"
        }
    }
	#test that the vars you want to include are present 
    if (is.character(include)) {
        if (!all(include %in% xName)) {
            stop("variable in include is not included formula or dataset")
        }
        else {
            includeName <- include
            mergeIncName <- paste0(includeName, collapse = " ")
        }
    } else if (is.null(include)) {
        includeName <- NULL
        mergeIncName <- "NULL"
    } else {
        stop("include should be character vector indicating variable to be included in all models")
    }
	#test and apply any weights 
    if (!is.null(weights)) {
        if (length(weights) == nrow(data)) {
            weightData <- data * sqrt(weights)
        }
        else {
            stop("Variable length is different ('(weights)')")
        }
    } else {
        weightData <- data
    }
	
	#make your formula 
    lmFull <- lm(formula, data = weightData)
	#figure class of vars 
    allVarClass <- attr(lmFull$terms, "dataClasses")
	
	#do a table of classes 
    classTable <- as.data.frame(table(allVarClass))
    colnames(classTable) <- c("class", "variable")
	#augment table of classes 
    for (i in names(table(allVarClass))) {
        classTable[names(table(allVarClass)) %in% i, 2] <- paste0(names(allVarClass[allVarClass %in% 
            i]), collapse = " ")
    }
    classTable$class <- paste0(classTable$class, ":")
	#accommodate factors 
    if (any(allVarClass == "factor")) {
        factVar <- names(which(allVarClass == "factor"))
        for (i in factVar) {
            weightData[, i] <- as.factor(as.numeric(weightData[, 
                i]))
        }
    }
	
	#set up data 
    xMatrix <- as.matrix(weightData[, xName])
    qrXList <- qr(xMatrix, tol = 1e-07)
    rank0 <- qrXList$rank
    pivot0 <- qrXList$pivot
    if (rank0 < length(pivot0)) {
        mulcolX <- colnames(qrXList$qr)[pivot0[(rank0 + 1):length(pivot0)]]
        mulcolMergeName <- paste0(mulcolX, collapse = " ")
    } else {
        mulcolX <- NULL
        mulcolMergeName <- "NULL"
    }
    xName <- setdiff(xName, mulcolX)
	#pull off Y
    Y <- as.matrix(lmFull$model[, yName])
	#this is how many models (Y vars) youre fitting 
    nY <- ncol(Y)
	#how many obs 
    nObs <- nrow(data)
	
	#set up stats 
    if (nY == 1) {
        approxF <- "F"
    } else {
        approxF <- multivarStat
        if (any(c(select) == c("BIC", "CP", "Rsq", 
            "adjRsq"))) {
            stop("Can't specify 'BIC','CP','Rsq' or 'adjRsq' when using multivariate multiple regression")
        }
    }
	
	#precheck request 
    if ((select == "CP" | select == "BIC") & lmFull$rank >= 
        nObs) {
        stop("'select' can't specify 'CP' or 'BIC' when variable number is greater than number of observation")
    }
	
	#setup result 
    result <- list()
    ModInf <- matrix(NA, 9, 1)
    ModInf <- cbind(ModInf, matrix(c(yName, mergeIncName, selection, 
        select, sle, sls, approxF, mulcolMergeName, intercept),       #SARA CHANGED:  had sle twice!
        9, 1))
    ModInf <- data.frame(ModInf)
    colnames(ModInf) <- c("", "")
    ModInf[, 1] <- c("Response Variable = ", "Included Variable = ", 
        "Selection Method = ", "Select Criterion = ", 
        "Entry Significance Level(sle) = ", "Stay Significance Level(sls) = ", 
        "Variable significance test = ", "Multicollinearity Terms = ", 
        "Intercept = ")
    if (select == "SL") {
        if (selection == "forward") {
            ModInf <- ModInf[-6, ]
        }
        else if (selection == "backward") {
            ModInf <- ModInf[-5, ]
        }
        else if (selection == "score") {
            ModInf <- ModInf[-c(5:6), ]
        }
    } else {
        ModInf <- ModInf[-c(5:6), ]
    }
    rownames(ModInf) <- 1:nrow(ModInf)
    result$"Basic Information" <- ModInf
    result$"Variable Class" <- classTable
	
	#NOW RUN
	####################################################################
	#first, all possible subsets 
	####################################################################
    if (selection == "score") {
        tempresult <- matrix(NA, 1, 4)
        colnames(tempresult) <- c("NoVariable", "RankModel", 
            select, "VariablesEnteredinModel")
        finalResult <- tempresult
        if (!is.null(includeName)) {
            lmIncForm <- reformulate(c(intercept, includeName), 
                yName)
            lmInc <- lm(lmIncForm, data = weightData)
            tempresult[1, c(1:4)] <- c(length(attr(lmInc$terms, 
                "term.labels")), lmInc$rank, modelFitStat(select, 
                lmInc, "LeastSquare"), paste(c(intercept, 
                includeName), collapse = " "))
            finalResult <- rbind(finalResult, tempresult)
            checkX <- xName[!xName %in% includeName]
        } else {
            checkX <- xName
        }
        for (nv in 1:length(checkX)) {
            comTable <- combn(length(checkX), nv)
            subSet <- NULL
            for (ncom in 1:ncol(comTable)) {
                comVar <- c(intercept, includeName, checkX[comTable[, 
                  ncom]])
                tempFormula <- reformulate(comVar, yName)
                lmresult <- lm(tempFormula, data = weightData)
                tempresult[1, 1:4] <- c(length(attr(lmresult$terms, 
                  "term.labels")), lmresult$rank, modelFitStat(select, 
                  lmresult, "LeastSquare"), paste(comVar, 
                  collapse = " "))
                subSet <- rbind(subSet, tempresult)
            }
            if (is.null(best)) {
                nbest <- nrow(subSet)
            } else {
                if (nrow(subSet) > best) {
                  nbest <- best
                } else {
                  nbest <- nrow(subSet)
                }
            }
            bestSubSet <- as.data.frame(subSet)
            bestSubSet[, 2] <- as.numeric(bestSubSet[, 2])
            if (select == "Rsq" | select == "adjRsq") {
                subResultSort <- bestSubSet[order(bestSubSet[, 
                  2], decreasing = TRUE), ]
            } else {
                subResultSort <- bestSubSet[order(bestSubSet[, 
                  2], decreasing = FALSE), ]
            }
            finalResult <- rbind(finalResult, subResultSort[1:nbest, 
                ])
        }
        finalResult <- finalResult[-1, ]
        rownames(finalResult) <- 1:nrow(finalResult)
        result$Process <- finalResult

	####################################################################
	#stepwise bcs not all possible 
	####################################################################
    } else {	
        subBestPoint <- data.frame(Step = numeric(), EnteredEffect = character(), 
            RemovedEffect = character(), DF = numeric(), NumberEffectIn = numeric(), 
            NumberParmsIn = numeric(), select = numeric())
        colnames(subBestPoint)[7] <- select
        bestPoint <- subBestPoint
		
		#if backwards 
        if (selection == "backward") {
            addIdx <- FALSE
            xModel <- c(intercept, includeName, setdiff(xName, 
                includeName))
            xResidual <- NULL
            if (select == "SL") {
                PIC <- 1
            } else {
                PIC <- modelFitStat(select, lmFull, "LeastSquare")
            }
            bestPoint[1, ] <- c(0, "", "", "", 
                length(attr(lmFull$terms, "term.labels")), 
                lmFull$rank, PIC)
				
		#if forward or bidirectional 
        } else {
			#setup 
            addIdx <- TRUE
            xModel <- c(intercept, includeName)
            xResidual <- setdiff(xName, includeName)
            fmInt <- reformulate(intercept, yName)
            fitInt <- lm(fmInt, data = weightData)
			
			#evaluate "PIC", it depends on "SL" vs other selection 
            if (select == "SL") {
                PIC <- 1
            } else {
                if (intercept == "1") {
                  PIC <- modelFitStat(select, fitInt, "LeastSquare")
                } else {
                  if (select %in% c("Rsq", "adjRsq")) {
                    PIC <- 0
                  } else {
                    PIC <- Inf
                  }
                }
            }
			
			#what to do 
            bestPoint[1, ] <- c(0, intercept, "", fitInt$rank, 
                length(attr(fitInt$terms, "term.labels")), 
                fitInt$rank, PIC)
				
			#force entry of included 
            if (!is.null(includeName)) {
                fmInc <- reformulate(c(intercept, includeName), 
                  yName)
                fitInc <- lm(fmInc, data = weightData)
                if (select == "SL") {
                  PIC <- anova(fitInc, fitInt, test = approxF)[2, 
                    "Pr(>F)"]
                } else {
                  PIC <- modelFitStat(select, fitInc, "LeastSquare")
                }
                subBestPoint[1, ] <- c(0, mergeIncName, "", 
                  anova(fitInc, fitInt, test = approxF)[2, "Df"], 
                  length(attr(fitInt$terms, "term.labels")), 
                  fitInc$rank, PIC)
                bestPoint <- rbind(bestPoint, subBestPoint)
            }
        } #end of backward/forward/bidirectional setup 
		
		#stepwise LOOP 
		kountr<-0 #sara added 
        while (TRUE) {
			kountr<-kountr+1 #sara added 
			dbp(kountr) #sara added 
			
			##################################################################################
			#SARAS ADDED CODE TO AVOID INFINITE LOOP 
			#DID YOU JUST BACK OUT THE VAR YOU ADDED?  IF SO, YOURE DONE 
			##################################################################################
			dbp(bestPoint)
			howManyIter<-nrow(bestPoint)
			
			if(nrow(bestPoint)>1) {
				if (bestPoint[howManyIter,"RemovedEffect"] == bestPoint[howManyIter-1,"EnteredEffect"]) {
					print("HIT SARAS BREAKPOINT: just removed == just entered")
					bestPoint<-bestPoint[1:(nrow(bestPoint)-2),]
					break
				}
			}	
			##################################################################################
			#END OF SARAS ADDED CODE TO AVOID INFINITE LOOP 
			##################################################################################
			
            fm0 <- reformulate(xModel, yName)
            lmAlt <- lm(fm0, data = weightData)
			
			#bifurcate based on T/F addIdx ... was probably initialized as TRUE 
            if (addIdx == TRUE) {
				#if adding, youll check everything not in model 'XResidual' here is vars not in model
                xCheck <- xResidual
                if (length(xCheck) == 0) {
                  break
                }
                xCheckList <- as.list(xCheck)
                names(xCheckList) <- xCheck
                fmX <- lapply(xCheckList, function(x) {
                  reformulate(c(xModel, x), yName)
                })
            } else {
				#if subtracting, youll check everything included  
                xCheck <- setdiff(xModel, c(intercept, includeName))
                if (length(xCheck) == 0) {
                  break
                }
                xCheckList <- as.list(xCheck)
                names(xCheckList) <- xCheck
                fmX <- lapply(xCheckList, function(x) {
                  reformulate(setdiff(xModel, x), yName)
                })
            }
			
			#youre going to fit a reg adding each possible model 
            fitX <- lapply(fmX, function(x) {
                lm(x, data = weightData)
            })
			print("passed fitX, reg with each possible model")

			#here you branch based on selection criteria 
			#--"SL" entails ANOVA 
            if (select == "SL") {
                PICset <- sapply(fitX, function(x) {
                  anova(x, lmAlt, test = approxF)[2, "Pr(>F)"]
                })
                Fset <- sapply(fitX, function(x) {
                  anova(x, lmAlt, test = approxF)[2, "F"]
                })
            } else {
			#--anything else;
			#	-- if youre backing out with only one var to check(?) under certain conditions (?dunno what intercept=="0" means) 
                if (addIdx == FALSE & length(xCheck) == 1 & intercept == "0") {
                  PICset <- Inf
                  names(PICset) <- xCheck
                }
                else {
			#	-- otherwise, fit all the models and keep the stat that you requested as 'select'	
                  PICset <- sapply(fitX, function(x) {
                    modelFitStat(select, x, "LeastSquare")
                  })
                }
            }
			#--IF EITHER: (1) select is Rsq, adjRsq or (2) youre backing out when select is SL 
			#--select your best var 
            if (select == "Rsq" | select == "adjRsq" | (select == "SL" & addIdx == FALSE)) {
                PIC <- max(PICset)
                minmaxVar <- names(which.max(PICset))
                bestLm <- fitX[[minmaxVar]]
            } else {
                PIC <- min(PICset)
                minmaxVar <- names(which.min(PICset))
                bestLm <- fitX[[minmaxVar]]
                if (sum(PICset %in% PIC) > 1 & select == "SL") {
                  Fvalue <- max(Fset)
                  minmaxVar <- names(which.max(Fset))
                  bestLm <- fitX[[minmaxVar]]
                  PIC <- PICset[minmaxVar]
                }
            }
			#--heres a criteria to break 
            if (bestLm$rank == lmAlt$rank & addIdx == TRUE) {
				print("BREAK hit for bestLm$rank == lmAlt$rank & addIdx == TRUE")
                break
            } else {
			#--if you dont break , indicator will be a Boolean
			#	--a) if your select="SL"
                if (select == "SL") {
                  if (addIdx == FALSE) {
                    indicator <- PIC > sls
                  } else {
                    indicator <- PIC < sle
                  }
			#	--b) if your select is "Rsq" or "adjRsq"
                } else if (select == "Rsq" | select == "adjRsq") {
                  indicator <- PIC > as.numeric(bestPoint[nrow(bestPoint),7])
                } else {
			#	--c) if your select is something else
                  indicator <- PIC <= as.numeric(bestPoint[nrow(bestPoint), 7])
                }
				
			#--if "indicator" Boolean is true, you have a model to work with 	
                if (indicator == TRUE) {
                  smr <- summary(bestLm)
				  dbp(smr)
				  #for one Y, extract the pval 
                  if (nY == 1) {
                    f <- smr$fstatistic
					dbp(f)
					if(is.null(f[1])) { #ADDED BY SARA 
                      pval <- NaN					
					} else if (is.nan(f[1])) {
                      pval <- NaN
                    } else {
                      pval <- pf(f[1], f[2], f[3], lower.tail = F)
                    }
                  } else {
				  #for more than one y, youll do something else 
                    for (ny in 1:nY) {
                      f <- smr[[ny]]$fstatistic
                      if (is.nan(f[1])) {
                        pval <- NaN
                      } else {
                        pval <- pf(f[1], f[2], f[3], lower.tail = F)
                      }
                    }
                  }
				  
				  #if you dont have a pval, and its not Rsq or adjRsq, youll stop 
                  if (is.nan(pval) == TRUE & (select != "Rsq" | select != "adjRsq")) {
                    break
                  }
				  
				  #if youre adding, and you didnt stop, add the model into the variable 
                  if (addIdx == TRUE) {
                    xModel <- append(xModel, minmaxVar)
                    xResidual <- setdiff(xResidual, minmaxVar)
                    subBestPoint[1, ] <- c(as.numeric(bestPoint[nrow(bestPoint), 
                      1]) + 1, minmaxVar, "", anova(lmAlt, 
                      bestLm, test = approxF)[2, "Df"], 
                      length(attr(bestLm$terms, "term.labels")), 
                      bestLm$rank, PIC)				  
				  #if youre not adding, then back it out 
                  } else {
                    xResidual <- append(xResidual, minmaxVar)
                    xModel <- setdiff(xModel, minmaxVar)
                    subBestPoint[1, ] <- c(as.numeric(bestPoint[nrow(bestPoint), 
                      1]) + 1, "", minmaxVar, anova(bestLm, 
                      lmAlt, test = approxF)[2, "Df"], 
                      length(attr(bestLm$terms, "term.labels")), 
                      bestLm$rank, PIC)
                  }
				  
				  #append to our "Process" (report 
                  bestPoint <- rbind(bestPoint, subBestPoint)
				  #see what to do next based on if youre bidirection 
                  if (selection == "bidirection") {
                    if (addIdx == FALSE) {
					  print("hit addIdx is FALSE for bidirection")
                      next
                    } else {
					  print("setting addIdx to FALSE for biderction")
                      addIdx <- FALSE
                      next
                    }
                  } else {
					print("because youre not bidirection, youll carry on")
					next
                  }

				  #THIS IS IF INDICATOR IS NOT TRUE 
                } else {
				  print("INDICATOR IS NOT TRUE")
                  if (selection == "bidirection" & addIdx == FALSE) {
					print("setting addIdx to TRUE bcs addIdx is FALSE given bidirection")
                    addIdx <- TRUE
                    next
                  } else {
					print("breaking because NOT (bidrection & addIdx==FALSE")
                    break
                  }
                }	
            }			
        }#this is the end of the "real" iterations 
		
        if (is.null(xModel)) {
            parEst <- NULL
        } else {
            parEst <- summary(lm(reformulate(xModel, yName), 
                data = weightData))
            parEstList <- list()
            if (nY > 1) {
                for (i in names(parEst)) {
                  parEstList[i] <- list(parEst[[i]]$coefficients)
                }
            } else {
                parEstList <- list(parEst$coefficients)
                names(parEstList) <- yName
            }
        }
        bestPoint$DF <- abs(as.numeric(bestPoint$DF))
        bestPoint$DF[is.na(bestPoint$DF)] <- ""
        result$Process <- bestPoint
        result$Variables <- xModel #sara corrected spelling
        result$Coefficients <- parEstList
		
    } #end of stepwise loop 
    return(result)

}