# Rbrul version 2.05
# Daniel Ezra Johnson
# April 15 2012

w <- options("warn")
options(warn=2)
a <- try(system("open http://www.danielezrajohnson.com/analytics.html",
	ignore.stderr=T,wait=F),silent=T)
a <- try(system("start http://www.danielezrajohnson.com/analytics.html",
	ignore.stderr=T,wait=F,invisible=T,show.output.on.console=F),silent=T)
options(warn=as.numeric(w))

rbrul <- function(pack=T){
	w <- options("warn")
	options(warn=2)
	pdf <- print.data.frame
	unlockBinding("print.data.frame", asNamespace("base"))
	assign("print.data.frame", rownames.off, envir=asNamespace("base"))
	cntrsts <- options("contrasts")
	rb <- try(rbrul2(pack=pack))
	if ("try-error" %in% class(rb))
		cat(sep="\n","\nSorry, there was an error and Rbrul crashed.
Help is available by emailing rbrul.list@gmail.com.")
	try(rm("addmargins", envir=globalenv()), silent=T)
	assign("print.data.frame", pdf, envir=asNamespace("base"))
	lockBinding("print.data.frame", asNamespace("base"))
	options(contrasts=cntrsts[[1]])
	options(warn=as.numeric(w))
} 

rbrul2 <- function(pack=T){
	if (regexpr("apple",sessionInfo()[[1]]$platform)>0){
		cat("\n","If you are having any trouble getting Rbrul to load,")
		cat("\n","please download and install the latest version of R from http://cran.r-project.org.")
		cat("\n","Then launch R, source Rbrul, and run Rbrul again.")
		cat("\n","Some Windows users may need to run R as an administrator (right-click for this option).")
		cat("\n")}
	if (pack) old <- get.packages() else old <- F
	if (old) {cat(sep="\n","Please go to http://cran.r-project.org/ 
and update to the newest version of R, which is at least 2.14.0."); return()}
	if (!exists("silent")) silent <<- F
	if (!exists("verbose")) verbose <<- T
	if (!exists("gold")) gold <<- T
	if (!exists("centered")) centered <<- T
	if (!exists("blups")) blups <<- T
	if (!exists("thresh")) thresh <<- 0.05
	if (!exists("sim")) sim <<- F
	if (!exists("decimals")) decimals <<- 3
	if (!exists("contr")) contr <<- "contr.sum"
	options(contrasts=c(unordered=contr,ordered="contr.poly"))
	if (!exists("loaded.flag")) loaded.flag <<- F
	if (!exists("fullpath")) fullpath <<- ""
	if (!exists("sepchar")) sepchar <<- ","
	if (!exists("shortpath")) shortpath <<- paste(getwd(),"/",sep="")
	main.pick <<- 99
	while (main.pick != 0){
		cat("\n")
		if (!exists("datafile")) {
			cat(sep="\n","No data loaded.")
			cat("\n")} 
		else {
			cat(sep="\n",paste("Current data file is:",fullpath))
			data.structure()}			
		cat(sep="\n","MAIN MENU")
		if (exists("datafile"))
			cat(sep="\n","1-load/save data 2-adjust data")
		else
			cat(sep="\n","1-load/save data")
		if (exists("datafile"))
			cat(sep="\n","4-crosstabs 5-modeling 6-plotting")
		if (exists("datafile"))
			cat(sep="\n","8-restore data 9-reset 0-exit")
		else
			cat(sep="\n","9-reset 0-exit")
		main.pick <<- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(main.pick)==0) main.pick <<- 99
		if (main.pick == 1) get.data()
		if (main.pick == 2)
			if (exists("datafile")) adjusting() else cat(sep="\n","Please load data first.")
		if (main.pick == 4)
			if (exists("datafile")) crosstab() else cat(sep="\n","Please load data first.")
		if (main.pick == 5)
			if (exists("datafile")) modeling() else cat(sep="\n","Please load data first.")
		if (main.pick == 6)
			if (exists("datafile")) plotting() else cat(sep="\n","Please load data first.")
		if (main.pick == 8)
			if (exists("datafile"))
			{datafile <<- datafile2; loaded.flag <<- T} 
			else cat(sep="\n","No data loaded!")
		if (main.pick == 9) reset()}}

addmargins <- function (A, margin = 1:length(dim(A)), FUN = sum, quiet = FALSE) 
{
    if (is.null(dim(A))) 
        stop("'A' must be an array or table")
    n.sid <- length(margin)
    miss.FUN <- missing(FUN)
    if (length(FUN) == 1 && !is.list(FUN)) {
        fname <- if (!miss.FUN) 
            deparse(substitute(FUN))
        else "total"
        FUN <- list(FUN)
        names(FUN) <- fname
    }
    if (!miss.FUN) {
        add.names <- function(thelist) {
            n <- names(thelist)
            if (is.null(n)) 
                n <- rep("", length(thelist))
            for (i in seq_along(thelist)[-1]) {
                if (!is.call(thelist[[i]])) {
                  if (n[i] == "") 
                    n[i] <- as.character(thelist[[i]])
                }
                else if (as.character(thelist[[i]][[1]]) == "list") 
                  thelist[[i]] <- add.names(thelist[[i]])
            }
            names(thelist) <- n
            thelist
        }
        if (mode(substitute(FUN)) == "call") 
            FUN <- eval(add.names(substitute(FUN)))
        if (is.null(names(FUN))) 
            names(FUN) <- rep("", length(FUN))
    }
    if (length(FUN) != n.sid) {
        if (length(FUN) == 1) 
            FUN <- rep(FUN, n.sid)
        else stop(gettextf("length of FUN, %d,\n does not match the length of the margins, %d", 
            length(FUN), n.sid), domain = NA)
    }
    fnames <- vector("list", n.sid)
    for (i in seq_along(FUN)) {
        fnames[[i]] <- names(FUN)[i]
        if (is.list(FUN[[i]])) {
            topname <- fnames[[i]]
            fnames[[i]] <- names(FUN[[i]])
            blank <- fnames[[i]] == ""
            fnames[[i]][blank] <- seq_along(blank)[blank]
            if (topname == "") {
                fnames[[i]][blank] <- paste("Margin ", margin[i], 
                  ".", fnames[[i]][blank], sep = "")
            }
            else {
                fnames[[i]] <- paste(topname, ".", fnames[[i]], 
                  sep = "")
            }
        }
        else if (fnames[[i]] == "") 
            fnames[[i]] <- paste("Margin", margin[i])
    }
    expand.one <- function(A, margin, FUN, fnames) {
        if (class(FUN) != "list") 
            FUN <- list(FUN)
        d <- dim(A)
        n.dim <- length(d)
        n.mar <- length(FUN)
        newdim <- d
        newdim[margin] <- newdim[margin] + n.mar
        newdimnames <- dimnames(A)
        newdimnames[[margin]] <- c(newdimnames[[margin]], fnames)
        n.new <- prod(newdim)
        skip <- prod(d[1:margin])
        runl <- skip/d[margin]
        apos <- rep(c(rep(TRUE, skip), rep(FALSE, n.mar * runl)), 
            n.new/(skip + n.mar * runl))
        values <- numeric(length(apos))
        values[apos] <- as.vector(A)
        for (i in 1:n.mar) {
            mtab <- if (n.dim > 1) {
                apply(A, (1:n.dim)[-margin], FUN[[i]])
            }
            else FUN[[i]](A)
            select <- rep(FALSE, n.mar)
            select[i] <- TRUE
            mpos <- rep(c(rep(FALSE, skip), rep(select, each = runl)), 
                prod(dim(A))/skip)
            values[mpos] <- as.vector(mtab)
        }
        new.A <- array(values, dim = newdim, dimnames = newdimnames)
        if (inherits(A, "table")) 
            class(new.A) <- c("table", class(new.A))
        new.A
    }
    new.A <- A
    for (i in 1:n.sid) new.A <- expand.one(A = new.A, margin = margin[i], 
        FUN = FUN[[i]], fnames = fnames[[i]])
    if (!quiet && !miss.FUN && n.sid > 1) {
        cat("Margins computed over dimensions\nin the following order:\n")
        for (i in 1:n.sid) cat(paste(i), ": ", names(dimnames(A))[margin[i]], 
            "\n", sep = "")
    }
    new.A
}

adjusting <- function(){
	adjusting.pick <<- 99
	while (!adjusting.pick %in% c(9,0)){
	cat("\n")
	cat(sep="\n","ADJUSTING MENU")
	cat(sep="\n","1-change class 2-rename 3-exclude 4-retain 5-recode")
	cat(sep="\n","6-relevel 7-center/transform 8-count 9-main menu 0-exit")
	cat(sep="\n","10-make full interaction group 11-make partial interaction group")
	adjusting.pick <<- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
	if (length(adjusting.pick)==0) adjusting()
	if (adjusting.pick == 1) reclass()
	if (adjusting.pick == 2) rename()
	if (adjusting.pick == 3) exclude()
	if (adjusting.pick == 4) retain()
	if (adjusting.pick == 5) recode()
	if (adjusting.pick == 6) relev()
	if (adjusting.pick == 7) transform()
	if (adjusting.pick == 8) count()
	if (adjusting.pick == 10) star()
	if (adjusting.pick == 11) partial()
	if (adjusting.pick == 0) main.pick <<- 0}}
	
ben.mer.sim <- function (object, nsim = 1){
	dims <- object@dims
	etasim.fix <- as.vector(object@X %*% fixef(object))
	etasim.reff <- (as(t(object@A) %*% matrix(rnorm(nsim * dims["q"]), nc = nsim),"matrix"))
	if (length(object@muEta) == 0){
		sigma <- lme4:::sigma(object)
		etasim.resid <- matrix(rnorm(nsim * dims["n"]), nc = nsim)
		etasim <- etasim.fix + sigma*(etasim.reff+etasim.resid)
		return(as.vector(etasim))}
	if (length(object@muEta)>0){
		etasim <- etasim.fix+etasim.reff
		musim <- inv.logit(etasim)
		n <- length(musim)
		vsim <- matrix(rbinom(n*nsim,1,prob=musim),nc=nsim)
		return(as.vector(vsim))}}

bin2num <- function(x){
	if (class(x)!="factor") return(x)
	if (length(levels(x))!=2) return(x)
	levels(x)[1] <- 0
	levels(x)[2] <- 1
as.numeric(as.character(x))}

count <- function(){
        cat(sep="\n",paste("Column(s) to count? (press Enter to exit)",
        	get.menu(names(datafile)),sep=""))
        repeat {count.pick <- try(as.numeric(scan(what="character",
		nmax=length(names(datafile)),quiet=T)),silent=T)
                countees <- names(datafile)[count.pick]
                if (!NA %in% countees) break
                cat(sep="\n","Invalid entry, please try again.")}
        if (length(countees)>0)
        for (i in countees){
		count.col <- paste("count",i,sep=".")
		for (j in unique(as.character(datafile[,i])))
			datafile[datafile[,i]==j,count.col] <<- 
				length(datafile[datafile[,i]==j,1])}}

cross <- function(quiet=F,crossees=c(NULL,NULL)){
	if (!quiet){
		cat(sep="\n",paste("Choose two variables to cross? (press Enter to exit)",
        	get.menu(names(datafile)),sep=""))
		repeat {cross.pick <- try(as.numeric(scan(what="character",nmax=2,quiet=T)),silent=T)
			crossees <- names(datafile)[cross.pick]
			if (length(crossees)==2 & !NA %in% crossees) break
        	cat(sep="\n","Invalid entry, please try again.")}}
	x <- datafile[,crossees[1]]
   	y <- datafile[,crossees[2]]
    z <- list("error")
    zz <- vector()
    if (class(x)=="character") x <- as.factor(x)
    if (class(y)=="character") y <- as.factor(y)
    if (class(x)=="factor" & class(y)=="factor"){
		x <- x[drop=T]
		y <- y[drop=T]
		newfac <- 1
		for (i in 1:length(contrasts(x)[1,])){
		for (j in 1:length(contrasts(y)[1,])){
			z[[newfac]] <- get.contrasts(x,i)*get.contrasts(y,j)
			newfac <- newfac + 1}}}
	if (class(x)=="factor" & class(y)%in%c("integer","numeric")|
		class(y)=="factor" & class(x)%in%c("integer","numeric")){
		if (class(y)=="factor") {t <- y; y <- x; x <- t}
		for (i in 1:length(contrasts(x)[1,]))
			z[[i]] <- get.contrasts(x,i)*y}
	if (class(x)%in%c("integer","numeric") & 
		class(y)%in%c("integer","numeric"))
		z <- list(x*y)
	for (i in 1:length(z)){
		if (length(z)==1) datafile[,zz[i]<-paste(crossees[1],".",crossees[2],sep="")] <<- z[[i]]
		else datafile[,zz[i]<-paste(crossees[1],".",crossees[2],".",letters[i],sep="")] <<- z[[i]]}
zz}
                
crosstab <- function(){
	cat(sep="\n",paste("Cross-tab: factor for columns?",get.menu(names(datafile))))
	repeat {columns.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		columns <- names(datafile)[columns.pick]
		if (!NA %in% columns & length(columns.pick)>0) {
			datafile[,columns] <- as.factor(datafile[,columns]); break}
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n",paste("Cross-tab: factor for rows? (",
		get.menu(names(datafile)[!names(datafile)%in%columns],which(!names(datafile)%in%columns),paren=F),
		" Enter-none)",sep=""))
	repeat {rows.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		rows <- names(datafile)[rows.pick]
		if (length(rows.pick)==0) {rows <- 1; break}
		if (!NA %in% rows) {
			datafile[,rows] <- as.factor(datafile[,rows]); break}
		cat(sep="\n","Invalid entry, please try again.")}
	if (rows==1) pages <- 1
	else {
		cat(sep="\n",paste("Cross-tab: factor for 'pages'? (",
		get.menu(names(datafile)[!names(datafile)%in%c(columns,rows)],
			which(!names(datafile)%in%c(columns,rows)),paren=F)," Enter-none)",sep=""))
	repeat {pages.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		pages <- names(datafile)[pages.pick]
		if (length(pages.pick)==0) {pages <- 1; break}
		if (!NA %in% pages) {
			datafile[,pages] <- as.factor(datafile[,pages]); break}
		cat(sep="\n","Invalid entry, please try again.")}}
	if (exists("response")){
		if (response %in% names(datafile)){
			cat(sep="\n",
				paste("Cross-tab: variable for cells? (1-proportion/mean of ",
					response, " 2-change response variable Enter-counts)",sep=""))}
		else{
			cat(sep="\n","Cross-tab: variable for cells? (1-response proportion/mean Enter-counts)")}}
	else{
		cat(sep="\n","Cross-tab: variable for cells? (1-response proportion/mean Enter-counts)")}
	repeat {cells.pick <- try(scan(what="character",nmax=1,quiet=T),silent=T)
		if (length(cells.pick)==0){
			cat(sep="\n","counts")
			result <- addmargins(xtabs(as.formula(paste("~",
				rows,"+",columns,"+",pages,sep="")),datafile))
			return(print(result))}
		if (cells.pick == 1){
			if (!exists("response")) get.response() else
				if (!response%in%names(datafile)) get.response()
				break}
		if (cells.pick == 2){
			get.response()
			break}
		cat(sep="\n","Invalid entry, please try again.")}
	if (class(datafile[,response]) %in% c("integer","numeric"))
		{cat(sep="\n",paste("mean of",response))
		result <- round(addmargins(xtabs(as.formula(paste(response,"~",rows,"+",columns,"+",pages,sep="")),
			datafile))/
			addmargins(xtabs(as.formula(paste("~",rows,"+",columns,"+",pages,sep="")), datafile)),decimals)
		return(print(result))}	
	else {if (!exists("apps")) get.response()
		cat(sep="\n",paste("proportion of",response,"=",apps))
		result <- round(addmargins(xtabs(as.formula(paste("~",
			rows,"+",columns,"+",pages,sep="")),datafile[datafile[,response]%in%apps,]))/
			addmargins(xtabs(as.formula(paste("~",rows,"+",columns,"+",pages)),datafile)),decimals)}
	return(print(result))}

custom.plot <- function(){
	pal <- F
	data.structure()
	cat(sep="\n",paste("Data variables:",get.menu(names(datafile),paren=F)))
	cat("\n")
	if (!exists("model")){
		cat(sep="\n","No model loaded."); combo <- datafile} 
	else {
		mod <- model[[1]]
		if ("glm" %in% class(mod)){
			modelvars <- mod$data
			modelvars$predictions <- mod$fitted.values
			if (mod$family[1]=="gaussian"){
				familia <- "linear"
				modelvars$residuals <- mod$residuals}
			else familia <- "logistic"}
		if ("mer" %in% class(mod)){
			modelvars <- mod@frame
			if (length(mod@muEta)==0){
				modelvars$residuals <- residuals(mod)
				modelvars$predictions <- mod@X %*% fixef(mod)
				familia <- "linear"} 
			else {modelvars$predictions <- inv.logit(mod@X %*% fixef(mod))
				familia <- "logistic"}
			modelvars$predictions.including.random.effects <- fitted(mod)}
		cat(sep="\n",paste("Model variables:",get.menu(names(modelvars),index=
			
(length(names(datafile))+1):(length(names(datafile))+length(names(modelvars))),
			paren=F)))
		if (length(datafile[,1])==length(modelvars[,1]))
			combo <- cbind(datafile,modelvars)
		else {
			cat("\n")
			cat(sep="\n","Warning! Data file and model do not match. Use model variables 
only.")
			combo <- modelvars
			for (i in 1:length(names(datafile)))
				combo <- cbind(dummy=1,combo)}}
	n <- length(combo[,1])
	cat("\n")
	cat(sep="\n","Choose variable for y-axis?")
	repeat{y.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		y <- names(combo)[y.pick]
		if (length(y)==0) next
		if (!NA %in% y & y!="dummy") break
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Choose variable for x-axis?")
	repeat{x.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		x <- names(combo)[x.pick]
		if (length(x)==0) next
		if (!NA %in% x & x!="dummy") break
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Separate (and color) data according to the values of which variable? (press 
Enter to skip)")
	repeat{g.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		g <- names(combo)[g.pick]
		if (length(g)==0) {g <- NULL; gmean.pick <- 0; break}
		if (!NA %in% g & g!="dummy") {
			cat(sep="\n",paste("Also show data (in black) averaged over all values of 
",g,"? (1-yes Enter-no)",sep=""))
			repeat {gmean.pick <- 
try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
				if (length(gmean.pick)==0){
					gmean.pick <- 0
					break}
				if (gmean.pick==1){
					combo2 <- combo
					combo2[,g] <- "XXX"
					combo <- rbind(combo,combo2)
					pal <- T
					break}
				cat(sep="\n","Invalid entry, please try again.")}
			break}
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Split data into horizontal panels according to which variable? (press Enter to 
skip)")
	repeat{cols.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		cols <- names(combo)[cols.pick]
		if (length(cols)==0) {cols <- 1; break}
		if (!NA %in% cols & cols!="dummy") break
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Split data into vertical panels according to which variable? (press Enter to 
skip)")
	repeat{rows.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		rows <- names(combo)[rows.pick]
		if (length(rows)==0) {rows <- 1; break}
		if (!NA %in% rows & rows!="dummy") break
		cat(sep="\n","Invalid entry, please try again.")}
	if (cols==1 & rows==1) strip <- F else strip <- T
	cat(sep="\n","Type of points to plot? (raw points not recommended for binary data)",
		"(0-no points 1-raw points Enter-mean points)")
	repeat{points.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(points.pick)==0) {points.pick <- 0; p <- "pm"; break}
		if (points.pick==0) {p <- NULL; break}
		if (points.pick==1) {p <- "p"; break}
		cat(sep="\n","Invalid entry, please try again.")}
	if ((gmean.pick==1 & points.pick==1)|points.pick==0){
		cat(sep="\n","Scale points according to the number of observations?",
			"Enter size factor between 0.1 and 10 (1 = Enter = default)", 
			"or 0 to not scale points")
		repeat{scale.factor <- 
try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
			if (length(scale.factor)==0) {scale.factor <- 1; break}
			if (scale.factor==0) break
			if (0.1 <= scale.factor & scale.factor <= 10) break
			cat(sep="\n","Invalid entry, please try again.")}}
	cat(sep="\n","Type of lines to plot (raw lines not recommended for binary data)?",
		"0-no lines 1-raw lines Enter-mean lines)")
	repeat{lines.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(lines.pick)==0) {l <- "lm"; break}
		if (lines.pick==0) {l <- NULL; break}
		if (lines.pick==1) {l <- "l"; break}
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Add a reference line? (1-diagonal [y=x] 2-horizontal [y=0] Enter-none)")
	repeat{refline.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(refline.pick)==0) {refline <- ""; break}
		if (refline.pick==1) {refline <- "d"; break}
		if (refline.pick==2) {refline <- "h"; break}
		cat(sep="\n","Invalid entry, please try again.")}
	combo[,y] <- bin2num(combo[,y])
	if (length(g)==0) gg <- "NULL"
		else gg <- paste("combo$",g,sep="")
	form <- as.formula(paste(y,"~",x,"|",cols,"+",rows,sep=""))
	if (strip) strip <- strip.custom(strip.names=c(T,T),strip.levels=c(T,T))
	if (pal){
		superpose.symbol.old <- trellis.par.get("superpose.symbol")
		superpose.line.old <- trellis.par.get("superpose.line")
		superpose.symbol <- superpose.symbol.old
		superpose.line <- superpose.line.old
		superpose.line$col <- c(rainbow(length(unique(combo[,g]))-1),"black")
		superpose.symbol$col <- c(rainbow(length(unique(combo[,g]))-1),"black")
		trellis.par.set(superpose.symbol=superpose.symbol)
		trellis.par.set(superpose.line=superpose.line)}
	go.flag <- T
	if (exists("familia")) {
	if (familia=="logistic" & length(grep("predictions",x))>0){
	go.flag <- F
	print(xyplot(form,groups=eval(parse(text=gg)),type=c(l,p),auto.key=T,
		data=combo,xlim=c(-0.07,1.07),
		panel=function(...){
			if (refline=="d") panel.abline(0,1)
			if (refline=="h") panel.abline(h=0)
			panel.rxyplot(...)},
		strip=strip,scale.factor=scale.factor,
		scales=list(alternating=1)))}}
	if (go.flag) 
		print(xyplot(form,groups=eval(parse(text=gg)),type=c(l,p),auto.key=T,
		data=combo,
		panel=function(...){
			if (refline=="d") panel.abline(0,1)
			if (refline=="h") panel.abline(h=0)
			panel.rxyplot(...)},
		strip=strip,scale.factor=scale.factor,
		scales=list(alternating=1)))
	if (pal){
		trellis.par.set(superpose.symbol=superpose.symbol.old)
		trellis.par.set(superpose.line=superpose.line.old)}}	
# layout for rows and columns
# automatically order factors like in B&P
# prompt for y=x line if y and x are in predicted, response


dan.mer.sim <- function (object,nsim=1,index=NULL){
	fixed <- fixef(object)
	if (length(index)>0) fixed[index] <- 0
	etasim.fix <- as.vector(object@X %*% fixed)
	dims <- object@dims
	etasim.reff <- (as(t(object@A) %*% matrix(rnorm(nsim * dims["q"]), nc = nsim),"matrix"))
	if (length(object@muEta) == 0){
		sigma <- lme4:::sigma(object)
		etasim.resid <- matrix(rnorm(nsim * dims["n"]), nc = nsim)
		etasim <- etasim.fix + sigma*(etasim.reff+etasim.resid)
		return(as.vector(etasim))}
	if (length(object@muEta)>0){
		etasim <- etasim.fix+etasim.reff
		musim <- inv.logit(etasim)
		n <- length(musim)
		vsim <- matrix(rbinom(n*nsim,1,prob=musim),nc=nsim)
		return(as.vector(vsim))}}

data.structure <- function(){
	cat("\n")
	cat(sep="\n","Current data structure:")
	classes <- sapply(datafile,class)
	for (i in 1:length(classes)){
		cat(paste(names(datafile)[i]," (",classes[i]," with ",
			l <- length(u <- unique(datafile[,i])),
				ifelse(l==1, " value): ", " values): "),sep=""))
		cat(sep="\n",paste(paste(u[1:ifelse(l<5,l,5)],
			collapse=" "),ifelse(l<6,"","..."),collapse=""))}
	cat(sep="\n",paste("Total tokens:",length(datafile[,1]),sep=" "))
	cat("\n")}
	
exclude <- function(){
	cat(sep="\n",paste("Factor group to exclude from?",
		get.menu(names(datafile)),sep=""))
	repeat {exclude.group.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		exclude.group <- names(datafile)[exclude.group.pick]
		if (!NA %in% exclude.group) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(exclude.group) == 0) return()
	if (class(datafile[,exclude.group])!="factor"){
		cat(sep="\n",paste("Please change class of",exclude.group,"to 'factor' before excluding."))
		return()}
	fac <- levels(datafile[,exclude.group])
	cat(sep="\n",paste("Factors to exclude from ",exclude.group,"?",get.menu(fac),sep="")) 
	repeat {excluded.factors.pick <- 
		try(as.numeric(scan(what="character",nmax=length(fac),quiet=T)),silent=T)
		excluded.factors <- fac[excluded.factors.pick]
		if (!NA %in% excluded.factors) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(excluded.factors) == 0) return()
	datafile <<- subset(datafile,!datafile[,exclude.group]%in%excluded.factors)
	datafile[,exclude.group] <<- datafile[,exclude.group][drop=T]}

extract <- function(model){
	if ("try-error" %in% class(model)) return(list("Error Message"))
	if ("list" %in% class(model)) {
		cat(sep="\n",model[[2]]); cat("\n")
		model <- model[[1]]}
	extraction <- list()
	group.labels <- vector()
	factor.labels <- vector()
	contrastss <- vector()
	groups <- vector()
	if ("glm" %in% class(model)){
		term.labels <- attr(model$terms,"term.labels")
		factor.slot <- "model$model$"
		intercept <- model$coefficients[1]
		coefficientss <- model$coefficients[-1]} else {
	if ("mer" %in% class(model)){
		term.labels <- c(attr(attr(model@frame,"terms"),"term.labels"),
			names(VarCorr(model)))
		factor.slot <- "model@frame$"
		intercept <- model@fixef[1]
		coefficientss <- c(model@fixef[-1],model@ranef)}}
	if (length(term.labels)>0){
	for (t in term.labels){
		
if (all(split.term(t) %in% continuous)){
	factor.labels <- c(factor.labels,"+1")
	group.labels <- c(group.labels,t)
	contrastss <- c(contrastss,coefficientss[1])
	coefficientss <- coefficientss[-1]} else {
		
if (all(!split.term(t) %in% continuous)) {
	num.coefs <- vector()
	name.coefs <- vector()
	names.coefs <- data.frame()
	for (f in split.term(t)){
		name.coefs <- levels(as.factor(eval(parse(text=paste(factor.slot,f,sep="")))))
		num.coefs <- c(num.coefs,length(name.coefs)-1)
		names.coefs <-
			
rbind(names.coefs,cbind(factorr=f,level=name.coefs,first=c(T,rep(F,length(name.coefs)-1)),
				last=c(rep(F,length(name.coefs)-1),T)))}
	factor.update <- 
levels(as.factor(eval(parse(text=paste(factor.slot,split.term(t),sep="",collapse=":")))))
	factor.labels <- c(factor.labels,factor.update)
	group.labels <- c(group.labels,rep(t,length(factor.update)))
		if (t %in% random){
			contrastss <- c(contrastss,coefficientss[1:length(factor.update)])
			coefficientss <- coefficientss[-c(1:length(factor.update))]} else { 
		if (contr=="contr.treatment")
			for (f in factor.update){
				base.levels <- names.coefs$level[names.coefs$first==T]
				if (any(split.term(f) %in% base.levels))
					contrastss <- c(contrastss,0) else {
						contrastss <- c(contrastss,coefficientss[1])
						coefficientss <- coefficientss[-1]}}
		if (contr=="contr.sum"){
			contrastsss <- vector()
			for (f in factor.update){
				omitted.levels <- names.coefs$level[names.coefs$last==T]
				if (any(split.term(f) %in% omitted.levels))
					contrastsss <- c(contrastsss,NA) else {
						contrastsss <- c(contrastsss,coefficientss[1])
						coefficientss <- coefficientss[-1]}}
			contrastss <- c(contrastss,sumify(contrastsss,
				
length(levels(as.factor(eval(parse(text=paste(factor.slot,split.term(t)[1]))))))))
				}}} else {
				
factor.update <- 
levels(as.factor(eval(parse(text=paste(factor.slot,split.term(t)[!split.term(t)%in%continuous],sep="")))))
	if (split.term(t)[1] %in% continuous)
		factor.labels <- c(factor.labels,paste("+1",factor.update,sep=":")) else
		factor.labels <- c(factor.labels,paste(factor.update,"+1",sep=":"))
	group.labels <- c(group.labels,rep(t,length(factor.update)))
	if (contr=="contr.treatment")
		contrastss <- c(contrastss,0,coefficientss[1:length(factor.update)-1])
	if (contr=="contr.sum")
		contrastss <- c(contrastss,s<-coefficientss[1:length(factor.update)-1],-sum(s))
	coefficientss <- coefficientss[-c(1:length(factor.update)-1)]}}}
	
fram <- data.frame("group"=group.labels,"factor"=factor.labels,"coefs"=contrastss)
	groups <- unique(as.character(fram$group))
	for (i in 1:(length(groups))){
		extraction[[i]] <- fram[fram$group==groups[i],-1]
		names(extraction)[i] <- groups[i]
		names(extraction[[i]]) <- c("factor","coefs")
		if (groups[i] %in% continuous) names(extraction[[i]]) <- 
c("continuous","coefs")
		if (length(split.term(groups[i]))>1){
			names(extraction[[i]]) <- c("factor:factor","coefs")
			if (all(split.term(groups[i]) %in% continuous))
				names(extraction[[i]]) <- c("continuous:continuous","coefs") 
else {
			if (split.term(groups[i])[1] %in% continuous)
				names(extraction[[i]]) <- c("continuous:factor","coefs")
			if (split.term(groups[i])[2] %in% continuous)
				names(extraction[[i]]) <- c("factor:continuous","coefs")}}
		if (groups[i] %in% random) names(extraction[[i]]) <- c("random","coefs") 
		if (names(extraction[[i]])[1] %in% c("factor","factor:factor","random")){
			for (j in 1:length(extraction[[i]][,1])){
				extraction[[i]][j,"tokens"] <- length(which(as.character(
				
eval(parse(text=paste(factor.slot,split.term(groups[i]),sep="",collapse=":")))) ==
				as.character(extraction[[i]][j,1])))
			if (fam == "binomial")
				extraction[[i]][j,proportion] <- 
round(length(which(as.character(
					
eval(parse(text=paste(factor.slot,split.term(groups[i]),sep="",collapse=":")))) ==
					as.character(extraction[[i]][j,1]) &
					
as.character(eval(parse(text=paste(factor.slot,response,sep="")))) == 
					apps)) / extraction[[i]][j,"tokens"],decimals)
			if (fam == "gaussian")
				if (length(split.term(groups[i]))==1)
					extraction[[i]][j,proportion] <- 
round(mean(datafile[,response]						[datafile[,groups[i]] 
== as.character(extraction[[i]][j,1])],na.rm=T),decimals)
				else extraction[[i]][j,proportion] <- 
round(mean(datafile[,response]					
[datafile[,split.term(groups[i])[1]] == split.term(as.character(extraction[[i]][j,1]))[1]&
					datafile[,split.term(groups[i])[2]] == 
split.term(as.character(extraction[[i]][j,1]))[2]]),decimals)}
			if (fam=="binomial" & gold){
				weight <- 
round(get.weights(extraction[[i]]$coefs,extraction[[i]]$tokens),decimals)
				weight[weight==1] <- 
paste(">0.",paste(rep(9,decimals),collapse=""),sep="")
				weight[weight==0] <- 
paste("<0.",paste(rep(0,decimals-1),collapse=""),1,sep="")
				extraction[[i]]$weight <- weight
				if (centered) 
					
names(extraction[[i]])[which(names(extraction[[i]])=="weight")] <- "centered factor weight"
				else
					
names(extraction[[i]])[which(names(extraction[[i]])=="weight")] <- "uncentered factor weight"}}
	extraction[[i]]$coefs <- round(extraction[[i]]$coefs,decimals)
	if (length(extraction[[i]][,1]) > 1) extraction[[i]] <- 
		extraction[[i]][order(extraction[[i]]$coefs,decreasing=T),]
	if (fam == "binomial") names(extraction[[i]])[2] <- "logodds"
	if (fam == "gaussian") names(extraction[[i]])[2] <- "coef"
	if (groups[i] %in% random) {extraction[[i]]$"std dev" <- 
round(attr(eval(parse(text=paste(
			
"VarCorr(model)$",groups[i],sep=""))),"stddev")["(Intercept)"],decimals) 
		if (blups == F)
			extraction[[i]] <- data.frame("random"="","std 
dev"=extraction[[i]][1,"std dev"])}
	}}
	extraction$misc <- data.frame(deviance="",df="",intercept=intercept,mean="")
	if (fam=="binomial" & gold & all(!names(extraction) %in% continuous) & 
!all(names(extraction)%in%c("misc",random))){
		if (centered) extraction$misc$"centered input prob" <- 
round(get.input(extraction),decimals)
		else extraction$misc$"uncentered input prob" <- 
round(get.input(extraction),decimals)} 
	if ("glm" %in% class(model)){
		extraction$misc$deviance <- round(model$deviance,decimals)
		extraction$misc$df <- length(datafile[,response])-model$df.residual} else
	if ("mer" %in% class(model)){
		extraction$misc$deviance <- round(model@deviance[1],decimals)
		extraction$misc$df <- attr(logLik(model),"df")}
	extraction$misc$intercept <- round(mean(intercept),decimals)
	if (fam == "binomial"){
	extraction$misc$mean <- round(length(datafile[datafile[,response]==apps,response])/
		(length(datafile[datafile[,response]==apps,response])+
			length(datafile[datafile[,response]==non.apps,response])),decimals)
	if ("glm" %in% class(model)) {
		n <- length(residuals(model))
		extraction$misc$R2 <- 
			round((1 - exp((model$deviance - model$null.deviance)/n))/
				(1 - exp(-model$null.deviance/n)),decimals)}}
	if (fam == "gaussian"){
	extraction$misc$mean <- round(mean(datafile[,response],na.rm=T),decimals)
	if ("glm" %in% class(model))
		extraction$misc$R2 <- round(1 - model$deviance/model$null.deviance,decimals)}
	names(extraction$misc)[which(names(extraction$misc)=="mean")] <- "grand mean"
	if (fam == "binomial")
		names(extraction$misc)[which(names(extraction$misc)=="R2")] <- "Nagelkerke R2"

extraction}
 
fit.model <- function(data=datafile){
	if (!exists("contr")) contr <<- "contr.sum"
	options(contrasts=c(unordered=contr,ordered="contr.poly"))
		if (length(fixed)==0) fixed <<- 1
		if (length(continuous)==0) continuous <<- 1
		if (length(interactions)==0) interactions <<- 1
		form <- paste(response,
			paste(paste(fixed,collapse="+"),
			paste(continuous,collapse="+"),
			paste(interactions,collapse="+"),
			sep="+"),sep="~")
		if (length(random) > 0){
			form <- paste(form,paste("(1|",random,")",collapse="+",sep=""),sep="+")
			model <- try(glmer(as.formula(form),data,fam),silent=T)}
		else model <- try(glm(as.formula(form),fam,data),silent=T)
model}
	
get.contrasts <- function(x,n){
	cont <- vector()
	for (i in 1:length(x))
		cont[i] <- contrasts(x)[x[i],n]
cont}

get.continuous <- function(){
	continuous <<- vector()
		cat(sep="\n",paste("Are any predictors continuous? (",
			get.menu(predictors,predictors.pick,paren=F)," Enter-none)",sep=""))
		repeat {continuous.pick <<- try(as.numeric(scan(what="character",
			nmax=length(predictors),quiet=T)),silent=T)
			continuous <<- names(datafile)[continuous.pick]
			if (!NA %in% continuous & all(continuous%in%predictors)) break
			cat(sep="\n","Invalid entry, please try again.")}
		if (length(continuous) > 0)
			for (i in 1:length(continuous)){
				contin <- try(as.numeric(as.character(datafile[,continuous[i]])),silent=T)
				if (!"try-error"%in%class(contin)) datafile[,continuous[i]] <<- contin else
					cat(sep="\n",paste("Warning: could not make",continuous[i],"continuous!"))}
		for (j in 1:length(predictors))
			if (!predictors[j]%in%continuous)
				datafile[,predictors[j]] <<- as.factor(datafile[,predictors[j]])[drop=T]}

get.data <- function(){
	if (exists("datafile")) save.data()
	cat("\n")
		if (!exists("datafile")) cat(sep="\n","No data loaded.") else
		cat(sep="\n",paste("Current data file is:",fullpath))
	cat("\n")
	cat(sep="\n","What separates the columns in the data file to open?")
	cat(sep="\n","(c-commas s-semicolons t-tabs tf-token file)")
	cat(sep="\n","Press Enter to exit, keeping current data file, if any.")
	repeat{sep.pick <- try(scan(what="character",nmax=1,quiet=T),silent=T)
		if (length(sep.pick) == 0) return()
		if (sep.pick %in% c("c","s","t","tf")) break}

	if (sep.pick=="c"){
		datafile <<- try(read.table(fullpath <<- file.choose(),
			header=T,sep=",",quote="\"",comment.char=""),silent=T)
 		sepchar <<- ","}
	if (sep.pick=="s"){
		datafile <<- try(read.table(fullpath <<- file.choose(),
            header=T,sep=";",quote="\"",comment.char=""),silent=T)
        sepchar <<- ";"}
	if (sep.pick=="t"){
		datafile <<- try(read.table(fullpath <<- file.choose(),
			header=T,sep="\t",quote="\"",comment.char=""),silent=T)
		sepchar <<- "\t"}
	if (sep.pick=="tf"){
		raw <- try(readLines(fullpath <<- file.choose()),silent=T)
		if ("try-error"%in%class(raw)) datafile <<- raw else {
			datlines <- grep("^\\(",raw)
			dat <- substr(raw[datlines],2,nchar(raw[datlines]))
			spaces <- regexpr(" +",dat)
			comm <- substr(dat,spaces+attr(spaces,"match.length"),nchar(dat))
			rump <- substr(dat,1,ifelse(spaces>-1,spaces-1,nchar(dat)))
			minw <- min(nchar(rump))
			maxw <- max(nchar(rump))
			if (minw!=maxw){
				cat(sep="\n","Warning! Tokens not all the same length. Omit short token(s) or drop final column(s)?")
				cat(sep="\n","(o-omit tokens d-drop columns)")
				repeat{drop.pick <- try(scan(what="character",nmax=1,quiet=T),silent=T)
					if (length(drop.pick) > 0)
					if (drop.pick %in% c("d","o")) break}
			if (drop.pick=="o"){
				rump <- rump[which(nchar(rump)==maxw)]
				minw <- maxw}}
		datafile <<- data.frame(dummy=1:length(rump))
		for (i in 1:minw)
			datafile[,i+1] <<- as.factor(substr(rump,i,i))
		datafile <<- datafile[,-1]
		sepchar <<- ","
		names(datafile) <<- 1:minw
		cat("\n")
		print(datafile[1:5,])
		cat("\n")
		repeat {cat(sep="\n","Enter labels for the columns (press Enter twice when done)")
			names <- try(scan(what="character",nmax=minw,quiet=T))
			if (!"try-error"%in%names & length(names)>0) break
			cat(sep="\n","Enter at least two legitimate column names!")}
		names(datafile) <<- names
		datafile <<- datafile[,1:length(names(datafile)[!is.na(names(datafile))])]
		if (max(spaces)>-1)
			datafile$comments <<- ifelse(spaces>-1,comm,"")}}
	
	if ("try-error"%in%class(datafile)){
		cat("\n")
		cat(sep="\n",paste("Error in loading file."))
		cat(sep="\n",datafile[1])
		rm(datafile,envir=globalenv())
		return()}

	lastslash <- gregexpr("/",fullpath)[[1]][length(gregexpr("/",fullpath)[[1]])]
	if (lastslash == -1)
	lastslash <- gregexpr("\\\\",fullpath)[[1]][length(gregexpr("\\\\",fullpath)[[1]])]
	try(setwd(shortpath <<- substr(fullpath,1,lastslash)),silent=T)

	droprow <- vector()
	for (i in 1:length(datafile[,1]))
		if (all(is.na(datafile[i,]))) droprow <- c(droprow,i)
	if (length(droprow)>0) datafile <<- datafile[-droprow,]
	dropcolumn <- vector()
	for (j in 1:length(datafile[1,]))
		if (all(is.na(datafile[,j]))) dropcolumn <- c(dropcolumn,j)
	if (length(dropcolumn)>0) datafile <<- datafile[,-dropcolumn]
	
	loaded.flag <<- T
	datafile2 <<- datafile
}


get.input <- function(extraction){
	tally <- 0
	if (length(extraction)>1){
		for (i in 1:(length(extraction)-1)){
			if (names(extraction[[i]])[1] == "factor") {
				if (centered) tally <- 
				tally + sum(extraction[[i]]$logodds)/length(extraction[[i]]$logodds) else
				tally <- tally + sum(extraction[[i]]$logodds*extraction[[i]]$tokens)/
					sum(extraction[[i]]$tokens)}}}
inv.logit(extraction$misc$intercept+tally)}

get.interactions <- function(){
	interactions <<- vector()
	if (length(predictors)==1) return()
	interactions.update <- vector()
	for (i in 1:sum(1:length(predictors)-1)){
	cat(sep="\n",paste("Consider a pairwise interaction? Choose two predictors at a time. (",
		get.menu(predictors,predictors.pick,paren=F)," Enter-done)",sep=""))
    repeat {interactions.pick <-
    		try(as.numeric(scan(what="character",nmax=2,quiet=T)),silent=T)
    		interactions.update <- paste(names(datafile)[interactions.pick],collapse=":")
    		if (any(split.term(interactions.update) %in% random)){
    	cat(sep="\n","Interactions involving grouping factors (random slopes) are not yet supported. Try again? (Enter-done)")
    			next}
    		if (length(interactions.pick)==2 & !"NA" %in% split.term(interactions.update)
    			& !interactions.update %in% interactions)
    			if (!identical(split.term(interactions.update)[1],
    				split.term(interactions.update)[2])) break
    		if (length(interactions.pick)==0) break
    		cat(sep="\n","Invalid entry, please try again.")}
	if (length(interactions.pick)==0) break
	interactions <<- c(interactions,interactions.update)}}

get.menu <- function(variable,index=1:length(variable),paren=T){
	if (paren)
		menu <- paste(" (",paste(index,"-",variable,sep="",collapse=" "),")",sep="")
	if (!paren)
		menu <- paste(index,"-",variable,sep="",collapse=" ")
menu}
		
get.p <- function(model1, model2, trim=F){
if ("try-error" %in% class(model1) & "try-error" %in% class(model2)){
	p.flag <<- "error code ="
	p <- ifelse(model2=="interaction", -888, -999)
	return(p)}
if ("try-error" %in% class(model1) & !"try-error" %in% class(model2)){
	mod2 <- extract(model2)
	p.flag <<- "log-likelihood ="
	p <- -mod2$misc$deviance/2
	return(p)}
if (!"try-error" %in% class(model1) & "try-error" %in% class(model2)){
	p.flag <<- "error code ="
	p <- ifelse(model2=="interaction", 888, 999)
	return(p)}
if (!"try-error" %in% class(model1) & !"try-error" %in% class(model2)){
	p.flag <<- "p ="
	if ("glm" %in% class(model1) & "glm" %in% class(model2)){
		t1 <- attr(attr(terms(model1),"factors"),"dimnames")[[2]]
		t2 <- attr(attr(terms(model2),"factors"),"dimnames")[[2]]
		if (!(all(t1%in%t2)|(all(t2%in%t1)))) return("error")
		if (model1$df.residual == model2$df.residual) return(1)
		if (length(t1) < length(t2)) {m1 <- model1; m2 <- model2}
		if (length(t1) > length(t2)) {m1 <- model2; m2 <- model1}
		if (length(t1) == length(t2)){
				if (trim == F) return(1)
				if (trim == T) {
					if (model1$df.residual > model2$df.residual) {m1 <- model1; m2 <- model2}
					if (model1$df.residual < model2$df.residual) {m1 <- model2; m2 <- model1}}}		
		if (m1$family$family=="binomial") 
			p <- pchisq(m1$deviance-m2$deviance,m1$df.residual-m2$df.residual,lower.tail=FALSE)
		if (m1$family$family=="gaussian")
			p <- pf(((m1$deviance-m2$deviance)/(m1$df.residual-m2$df.residual))/(m2$deviance/m2$df.residual),
				m1$df.residual-m2$df.residual,m2$df.residual,lower.tail=FALSE)
		return(p)}
	if ("mer" %in% class(model1)){
		if ("mer" %in% class(model2))
			terms <- attr(attr(model2@frame,"terms"),"term.labels")
		else
			terms <- attr(attr(model2$model,"terms"),"term.labels")
		if (trim == F & identical(attr(attr(model1@frame,"terms"),"term.labels"),terms)){
			increase <- abs(mean(attr(logLik(model2),"df")-attr(logLik(model1),"df")))
			if (increase==0) return(1)
			p <- mean(pchisq(abs(-2*(logLik(model2)-logLik(model1))),
				c(increase-1,increase),lower.tail=FALSE))
			return(p)}
		if (identical(sapply(VarCorr(model1),length),sapply(VarCorr(model2),length))){
			if (sim)
				p <- sim.fz(model1,model2,threshold,min=100,max=500,confidence=.01)
			if (!sim)
				p <- pchisq(abs(mean(model2@deviance[1]-model1@deviance[1])),
					abs(mean(attr(logLik(model2),"df")-attr(logLik(model1),"df"))),
					lower.tail=FALSE)
			return(p)}}}
"error"}

# fix the abs() stuff

get.packages <- function(){
w <- options("warn")
sess <- sessionInfo()
if ((as.numeric(sess$R.version$major) == 2 & as.numeric(sess$R.version$minor) < 12.0) | 
        as.numeric(sess$R.version$major) < 2) {
   	          options(warn=as.numeric(w))
	          return(T)}
if (!exists("inv.logit")){
cat("\n")
try1 <- try(suppressPackageStartupMessages(require(boot)),silent=T)
if (try1==F|"try-error"%in%class(try1)){
cat(sep="\n","Installing boot package...")
cat("\n")
install.packages("boot")
try1 <- try(suppressPackageStartupMessages(require(boot)),silent=T)}}
if (!exists("binconf")){
cat("\n")
try4 <- try(suppressPackageStartupMessages(require(Hmisc)),silent=T)
if (try4==F|"try-error"%in%class(try4)){
cat(sep="\n","Installing Hmisc package...")
cat("\n")
install.packages("Hmisc")
try4 <- try(suppressPackageStartupMessages(require(Hmisc)),silent=T)}}
if (!exists("trellis.focus")){
cat("\n")
try3 <- try(suppressPackageStartupMessages(require(lattice)),silent=T)
if (try3==F|"try-error"%in%class(try3)){
cat(sep="\n","Installing lattice package...")
cat("\n")
install.packages("lattice")
try3 <- try(suppressPackageStartupMessages(require(lattice)),silent=T)}}
if (!exists("glmer")){
cat("\n")
try2 <- try(suppressPackageStartupMessages(require(lme4)),silent=T)
if (try2==F|"try-error"%in%class(try2)){
cat(sep="\n","Installing lme4 package...")
cat("\n")
install.packages("lme4")}}
if (installed.packages()["lme4","Version"]<"0.999375-35"){
cat("\n")
cat(sep="\n","Updating lme4 package...")
cat("\n")
detach(package:lme4)
try(update.packages(ask=F,oldPkgs="lme4"),silent=T)}
try2 <- try(suppressPackageStartupMessages(require(lme4)),silent=T)
options(warn=as.numeric(w))
F}

get.predictors <- function(){
	if (exists("predictors") & !loaded.flag) old.predictors <- predictors
	singletons <- vector()
	for (i in names(datafile)[-response.pick])
		if (length(levels(datafile[,i]))==1) singletons <- c(singletons,i)
	cat(sep="\n",paste("Choose predictors (independent variables) by number",
		if (exists("old.predictors"))
			if (!response %in% old.predictors) 
				paste(", or Enter to keep ",paste(old.predictors,collapse=" & "),sep=""),
		 get.menu(names(datafile)[!names(datafile)%in%c(response,singletons)],
			which(!names(datafile)%in%c(response,singletons))),sep=""))
	repeat {predictors.pick <<- 
		try(as.numeric(scan(what="character",nmax=length(names(datafile))-length(c(response,singletons)),
			quiet=T)),silent=T)
		if (length(predictors.pick)==0 & exists("old.predictors")){
			if (!response %in% old.predictors){
				predictors <<- old.predictors
				predictors.pick <<- which(names(datafile)%in%predictors)
				return(T)}}
		predictors <<- names(datafile)[predictors.pick]
		if (!NA %in% predictors & all(!predictors%in%c(response,singletons)) & length(predictors)>0) break
		if (length(predictors)==0) cat(sep="\n","Please choose at least one predictor.")
		else cat(sep="\n","Invalid entry, please try again.")}
	for (i in 1:length(predictors))
		if (NA %in% datafile[,predictors[i]])
			cat(sep="\n",paste("Danger! NA's exist in ",predictors[i],"!",sep=""))
loaded.flag <<- F
F}

get.random <- function(){
random <<- vector()
categorical <- predictors[!predictors%in%continuous]
if (length(categorical)>0){
	cat(sep="\n",paste("Any grouping factors (random effects)? (",
		get.menu(categorical,predictors.pick[!predictors.pick%in%continuous.pick],paren=F)," Enter-none)",sep=""))
	repeat {random.pick <- try(as.numeric(scan(what="character",
		nmax=length(categorical),quiet=T)),silent=T)
		random <<- names(datafile)[random.pick]
		if (!NA %in% random & all(random%in%categorical)) break
		cat(sep="\n","Invalid entry, please try again.")}}}

get.response <- function(){
	if (exists("response") & !loaded.flag) old.response <- response
	singletons <- vector()
	names(datafile)[which(names(datafile)=="response")] <<- "response."
	for (i in names(datafile))
		if (length(levels(datafile[,i]))==1) singletons <- c(singletons,i)
	cat(sep="\n",paste("Choose response (dependent variable) by number",
		if (exists("old.response")) paste(", or Enter to keep ",old.response,sep=""),
			get.menu(names(datafile)[!names(datafile)%in%singletons],
			which(!names(datafile)%in%singletons)),sep=""))
	repeat {response.pick <<- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(response.pick) == 0 & exists("old.response")) {
			response <<- old.response
			response.pick <<- which(names(datafile)==response)
			return(list(response,datafile[,response]))}
		response <<- names(datafile)[response.pick]
		if (!NA %in% response & length(response)>0) break
		if (length(response)==0) cat(sep="\n","Please choose a response variable.")
		else cat(sep="\n","Invalid entry, please try again.")}
		cat(sep="\n","Type of response? (1-continuous Enter-binary)")
		repeat {fam.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
			if (length(fam.pick)==0) {fam <<- "binomial"; break} 
			if (fam.pick == 1) {fam <<- "gaussian"; break}
			cat(sep="\n","Invalid entry, please try again.")}
		if (fam == "binomial"){
			if (exists("old.binary") & exists("old.response"))
				datafile[,old.response] <<- old.binary
			if (!exists("old.binary"))
				old.binary <<- datafile[,response]
			datafile[,response] <<- as.factor(datafile[,response])[drop=T]
			fac <- levels(datafile[,response])
			cat(sep="\n",paste("Choose application value(s) by number?",get.menu(fac)))
			repeat {application.values.pick <- 
				try(as.numeric(scan(what="character",nmax=length(fac)-1,quiet=T)),silent=T)
				application.values <- fac[application.values.pick]
				if (length(application.values) > 0)
					if (!NA %in% application.values)
						break
				cat(sep="\n","Invalid entry, please try again.")}
			apps <<- paste(application.values,collapse="+")
			for (i in 1:length(fac))
				if (fac[i] %in% application.values) fac[i] <- apps
			if (length(fac[fac!=apps]) == 1) non.application.values <- fac[fac!=apps] else {
				cat(sep="\n",paste("Choose non-application value(s) by number, or press Enter for all remaining factors?"))
				cat(sep="\n",get.menu(fac[fac!=apps],which(fac!=apps)))
				repeat {non.application.values.pick <- try(as.numeric(scan(what="character",
					nmax=length(fac)-length(application.values),quiet=T)),silent=T)
					non.application.values <- fac[non.application.values.pick]
					if (!NA %in% non.application.values & all(non.application.values%in%fac[fac!=apps]))
						break
					cat(sep="\n","Invalid entry, please try again.")}
				if (length(non.application.values) == 0)
					non.application.values <- fac[fac!=apps]}
			non.apps <<- paste(non.application.values,collapse="+")
			proportion <<- paste(apps,"/",apps,"+",non.apps,sep="")
			for (i in 1:length(fac[fac!=apps]))
				if (fac[fac!=apps][i] %in% non.application.values) fac[fac!=apps][i] <- non.apps
			levels(datafile[,response]) <<- fac
			datafile[,response] <<- relevel(datafile[,response],non.apps)
			datafile <<- subset(datafile,datafile[,response] %in% c(non.apps,apps))
			datafile[,response] <<- datafile[,response][drop=T]}
		if (fam == "gaussian"){
			proportion <<- "mean"
			resp <- try(as.numeric(as.character(datafile[,response])),silent=T)
			if (!"try-error"%in% class(resp)) datafile[,response] <<- resp
			else {cat(sep="\n",paste("Could not make",response,"continuous!")); get.response()}}
	if (NA %in% datafile[,response]) cat(sep="\n","Warning! NA's exist in response!")
list(response,datafile[,response])}
	
get.settings <- function(){
	cat(sep="\n","Run silently? (1-yes Enter-no)")
	repeat {silent.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(silent.pick) == 0) {silent <<- F; break}
		if (silent.pick == 1) {silent <<- T; break}
		cat(sep="\n","Invalid entry, please try again.")}
	if (!silent){
		cat(sep="\n","Hide intermediate models' details? (1-yes Enter-no)")
		repeat {verbose.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
			if (length(verbose.pick) == 0) {verbose <<- T; break}
			if (verbose.pick == 1) {verbose <<- F; break}
			cat(sep="\n","Invalid entry, please try again.")}
		cat(sep="\n","Hide factor weights? (1-yes Enter-no)")
		repeat {gold.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
			if (length(gold.pick) == 0) {gold <<- T; break}
			if (gold.pick == 1) {gold <<- F; centered <<- F; break}
			cat(sep="\n","Invalid entry, please try again.")}
		if (gold){
			cat(sep="\n","Center factors [recommended]? (1-no Enter-yes)")
			repeat {centered.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
				if (length(centered.pick) == 0) {centered <<- T; break}
				if (centered.pick == 1) {centered <<- F; break}
				cat(sep="\n","Invalid entry, please try again.")}}
		cat(sep="\n","Show random effect estimates (BLUPs)? (1-no Enter-yes)")
		repeat {blups.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
			if (length(blups.pick) == 0) {blups <<- T; break}
			if (blups.pick == 1) {blups <<- F; break}
			cat(sep="\n","Invalid entry, please try again.")}}
	cat(sep="\n","Threshold p-value for fixed effect significance? (1-automatic, Enter-0.05, or type another value)")
	repeat {thresh.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(thresh.pick) == 0) {thresh <<- 0.05; break}
		if (0 < thresh.pick & thresh.pick < 1) {thresh <<- thresh.pick; break}
		if (thresh.pick == 1) {thresh <<- "A"; break}
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Use slow but more accurate [?] simulation method for fixed effect significance? (1-yes Enter-no)")
	repeat {sim.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(sim.pick) == 0) {sim <<- F; break}
		if (sim.pick == 1) {sim <<- T; break}
		cat(sep="\n","Invalid entry, please try again.")}
	cat(sep="\n","Number of significant digits/decimal places to display? (Enter-3)")
	repeat {decimals.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(decimals.pick) == 0) {decimals <<- 3; break}
		if (!"try-error" %in% class(decimals.pick)) {decimals <<- decimals.pick; break}}
	cat(sep="\n","Type of contrasts to use for factors? (1-treatment Enter-sum)")
	repeat {contr.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(contr.pick) == 0) {contr <<- "contr.sum"; break}
		if (contr.pick == 1) {contr <<- "contr.treatment"; break}
		cat(sep="\n","Invalid entry, please try again.")}
	options(contrasts=c(unordered=contr,ordered="contr.poly"))}

get.variables <- function(){
	r <- get.response(); response <<- r[[1]]; datafile[,response] <<- r[[2]]
	g <- get.predictors()
	if (!g) {get.continuous(); get.random(); get.interactions()}	
	set.variables()
	load.variables()
variables}

get.weights <- function(logodds,tokens=rep(1,length(logodds))){
if (length(which(!is.na(logodds)))<2) weights <- NA else {
	real.logodds <- logodds[which(!is.na(logodds))]
	if (real.logodds[1] != 0) real.logodds <- real.logodds-real.logodds[1]
	real.tokens <- tokens[which(!is.na(logodds))]
	num.weights <- length(real.logodds)
	right <- real.logodds
	if (centered) left <- rep(1/num.weights,num.weights) else
		left <- real.tokens/num.weights
	for (i in 1:(num.weights-1)){
		left <- rbind(left,c(-1,rep(0,i-1),1,rep(0,num.weights-1-i)))}
	real.weights <- inv.logit(solve(as.matrix(left),right))
	weights <- rep(NA,length(logodds))
	weights[which(!is.na(logodds))] <- real.weights}
weights}

load.variables <- function(){
	response <<- variables[[1]]
	if (names(variables)[1]=="response.binary") fam <<- "binomial"
	if (names(variables)[1]=="response.continuous") fam <<- "gaussian"
	fixed <<- vector()
	continuous <<- vector()
	interactions <<- vector()
	random <<- vector()
	for (i in 2:length(variables)){
		if (names(variables)[i]=="fixed.factor")
			fixed <<- c(fixed,variables[[i]])
		if (names(variables)[i]=="fixed.continuous")
			continuous <<- c(continuous,variables[[i]])
		if (names(variables)[i]=="fixed.interaction")
			interactions <<- c(interactions,variables[[i]])
		if (names(variables)[i]=="random.intercept")
			random <<- c(random,variables[[i]])}}
			
modeling <- function(){
modeling.pick <<- 99
while (!modeling.pick %in% c(9,0)){
	cat("\n")
	if (!exists("variables")) {cat(sep="\n","No variables chosen."); cat("\n")} else {
		cat(sep="\n","Current variables are:")
		for (i in 1:length(variables)){
			cat(paste(names(variables)[i],": ",
				paste(variables[[i]],collapse=" "),sep=""))
			if (names(variables)[i]=="response.binary")
				cat(paste(" (",apps," vs. ",non.apps,")",sep=""))
			cat("\n")}
	cat("\n")}
	cat(sep="\n","MODELING MENU")
	cat(sep="\n","1-choose variables 2-one-level 3-step-up 4-step-down 5-step-up/step-down")
	cat(sep="\n","6-trim 7-plotting 8-settings 9-main menu 0-exit")
	cat(sep="\n","10-chi-square test")
	modeling.pick <<- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
	if (length(modeling.pick)==0) modeling()
	if (modeling.pick == 1) get.variables()
	if (modeling.pick == 2) step.one()
	if (modeling.pick == 3) step.up()
	if (modeling.pick == 4) step.down()
	if (modeling.pick == 5) step.updown()
	if (modeling.pick == 6) trim()
	if (modeling.pick == 7) plotting()
	if (modeling.pick == 8) get.settings()
	if (modeling.pick == 10) test.chisq()
	if (modeling.pick == 9) plotting.pick <<- 9
	if (modeling.pick == 0) {plotting.pick <<- 0; main.pick <<- 0}}}

panel.rsuperpose <- function (x, y = NULL, subscripts, groups, panel.groups = panel.rxyplot, 
    col = NA, col.line = superpose.line$col, col.symbol = superpose.symbol$col, 
    pch = superpose.symbol$pch, cex = superpose.symbol$cex, fill = superpose.symbol$fill, 
    font = superpose.symbol$font, fontface = superpose.symbol$fontface, 
    fontfamily = superpose.symbol$fontfamily, lty = superpose.line$lty, 
    lwd = superpose.line$lwd, alpha = superpose.symbol$alpha, 
    type = "p", ..., distribute.type = FALSE, scale.factor = 0) 
{
    if (distribute.type) {
        type <- as.list(type)
    }
    else {
        type <- unique(type)
        wg <- match("g", type, nomatch = NA)
        if (!is.na(wg)) {
            panel.grid(h = -1, v = -1)
            type <- type[-wg]
        }
        type <- list(type)
    }
    x <- as.numeric(x)
    if (!is.null(y)) 
        y <- as.numeric(y)
    if (length(x) > 0) {
        if (!missing(col)) {
            if (missing(col.line)) 
                col.line <- col
            if (missing(col.symbol)) 
                col.symbol <- col
        }
        superpose.symbol <- trellis.par.get("superpose.symbol")
        superpose.line <- trellis.par.get("superpose.line")
        vals <- if (is.factor(groups)) 
            levels(groups)
        else sort(unique(groups))
        nvals <- length(vals)
        col <- rep(col, length = nvals)
        col.line <- rep(col.line, length = nvals)
        col.symbol <- rep(col.symbol, length = nvals)
        pch <- rep(pch, length = nvals)
        fill <- rep(fill, length = nvals)
        lty <- rep(lty, length = nvals)
        lwd <- rep(lwd, length = nvals)
        alpha <- rep(alpha, length = nvals)
        cex <- rep(cex, length = nvals)
        font <- rep(font, length = nvals)
        fontface <- rep(fontface, length = nvals)
        fontfamily <- rep(fontfamily, length = nvals)
        type <- rep(type, length = nvals)
        panel.groups <- if (is.function(panel.groups)) 
            panel.groups
        else if (is.character(panel.groups)) 
            get(panel.groups)
        else eval(panel.groups)
        subg <- groups[subscripts]
        ok <- !is.na(subg)
        for (i in seq_along(vals)) {
            id <- ok & (subg == vals[i])
            typehack <- type[[i]]
            if (vals[i]=="XXX"){
            	typehack[which(typehack=="p")] <- "pm"
            	typehack[which(typehack=="l")] <- "lm"}
            if (any(id)) {
                args <- list(x = x[id], subscripts = subscripts[id], 
                  pch = pch[i], cex = cex[i], font = font[i], 
                  fontface = fontface[i], fontfamily = fontfamily[i], 
                  col = col[i], col.line = col.line[i], col.symbol = col.symbol[i], 
                  fill = fill[i], lty = lty[i], lwd = lwd[i], 
                  alpha = alpha[i], type = typehack, group.number = i, scale.factor = scale.factor,
                  ...)
                if (!is.null(y)) 
                  args$y <- y[id]
                do.call(panel.groups, args)
            }
        }
    }
}

panel.rxyplot <- function (x, y, type = "p", groups = NULL, pch = if (is.null(groups)) 
plot.symbol$pch else superpose.symbol$pch, 
    col, col.line = if (is.null(groups)) plot.line$col else superpose.line$col, 
    col.symbol = if (is.null(groups)) plot.symbol$col else superpose.symbol$col, 
    font = if (is.null(groups)) plot.symbol$font else superpose.symbol$font, 
    fontfamily = if (is.null(groups)) plot.symbol$fontfamily else superpose.symbol$fontfamily, 
    fontface = if (is.null(groups)) plot.symbol$fontface else superpose.symbol$fontface, 
    lty = if (is.null(groups)) plot.line$lty else superpose.line$lty, 
    cex = if (is.null(groups)) plot.symbol$cex else superpose.symbol$cex, 
    fill = if (is.null(groups)) plot.symbol$fill else superpose.symbol$fill, 
    lwd = if (is.null(groups)) plot.line$lwd else superpose.line$lwd, 
    horizontal = FALSE, ..., jitter.x = FALSE, jitter.y = FALSE, 
    factor = 0.5, amount = NULL, scale.factor = 0) 
{
    if (all(is.na(x) | is.na(y))) 
        return()
    x <- as.numeric(x)
    y <- as.numeric(y)

				xx <- x[!is.na(y)]
				yy <- y[!is.na(y)]
				y.avg <- tapply(yy,xx,mean)
				xvals <- try(as.numeric(names(y.avg)),silent=T)
				if ("try-error"%in%class(xvals)) xvals <- as.factor(names(y.avg))
				ord <- order(xvals)
				ym <- y.avg[ord]
				xm <- xvals[ord]
				if (scale.factor==0) {nn <- 1; s <- 1} else
				{nn <- tapply(yy,xx,length); s <- scale.factor}

    plot.symbol <- trellis.par.get("plot.symbol")
    plot.line <- trellis.par.get("plot.line")
    superpose.symbol <- trellis.par.get("superpose.symbol")
    superpose.line <- trellis.par.get("superpose.line")
    if (!missing(col)) {
        if (missing(col.line)) 
            col.line <- col
        if (missing(col.symbol)) 
            col.symbol <- col
    }
    if (!is.null(groups)) 
        panel.rsuperpose(x, y, type = type, groups = groups, pch = pch, 
            col.line = col.line, col.symbol = col.symbol, font = font, 
            fontfamily = fontfamily, fontface = fontface, lty = lty, 
            cex = cex, fill = fill, lwd = lwd, horizontal = horizontal, 
            panel.groups = panel.rxyplot, jitter.x = jitter.x, 
            jitter.y = jitter.y, factor = factor, amount = amount, scale.factor = scale.factor,
            ...)
    else {
        if ("o" %in% type || "b" %in% type) 
            type <- c(type, "p", "l")
        if ("g" %in% type) 
            panel.grid(h = -1, v = -1)
        if ("p" %in% type) 
            panel.points(x = if (jitter.x) 
                jitter(x, factor = factor, amount = amount)
            else x, y = if (jitter.y) 
                jitter(y, factor = factor, amount = amount)
            else y, cex = cex, fill = fill, font = font, fontfamily = fontfamily, 
                fontface = fontface, col = col.symbol, pch = pch, 
                ...)
		if ("pm" %in% type)
            panel.points(x = if (jitter.x) 
                jitter(xm, factor = factor, amount = amount)
            else xm, y = if (jitter.y) 
                jitter(ym, factor = factor, amount = amount)
            else ym, cex = cex*sqrt(nn)*s, fill = fill, font = font, fontfamily = fontfamily, 
                fontface = fontface, col = col.symbol, pch = pch, 
                ...)
        if ("l" %in% type) 
            panel.lines(x = x, y = y, lty = lty, col = col.line, 
                lwd = lwd, ...)
		if ("lm" %in% type) 
            panel.lines(x = xm, y = ym, lty = lty, col = col.line, 
                lwd = lwd, ...)
        if ("h" %in% type) {
            if (horizontal) 
                panel.lines(x = x, y = y, type = "H", lty = lty, 
                  col = col.line, lwd = lwd, ...)
            else panel.lines(x = x, y = y, type = "h", lty = lty, 
                col = col.line, lwd = lwd, ...)
        }
        if ("s" %in% type) {
            ord <- if (horizontal) 
                sort.list(y)
            else sort.list(x)
            n <- length(x)
            xx <- numeric(2 * n - 1)
            yy <- numeric(2 * n - 1)
            xx[2 * 1:n - 1] <- x[ord]
            yy[2 * 1:n - 1] <- y[ord]
            xx[2 * 1:(n - 1)] <- x[ord][-1]
            yy[2 * 1:(n - 1)] <- y[ord][-n]
            panel.lines(x = xx, y = yy, lty = lty, col = col.line, 
                lwd = lwd, ...)
        }
        if ("S" %in% type) {
            ord <- if (horizontal) 
                sort.list(y)
            else sort.list(x)
            n <- length(x)
            xx <- numeric(2 * n - 1)
            yy <- numeric(2 * n - 1)
            xx[2 * 1:n - 1] <- x[ord]
            yy[2 * 1:n - 1] <- y[ord]
            xx[2 * 1:(n - 1)] <- x[ord][-n]
            yy[2 * 1:(n - 1)] <- y[ord][-1]
            panel.lines(x = xx, y = yy, lty = lty, col = col.line, 
                lwd = lwd, ...)
        }
        if ("r" %in% type) 
            panel.lmline(x, y, col = col.line, lty = lty, lwd = lwd, 
                ...)
        if ("smooth" %in% type) 
            panel.loess(x, y, horizontal = horizontal, col = col.line, 
                lty = lty, lwd = lwd, ...)
        if ("a" %in% type) 
            panel.linejoin(x, y, horizontal = horizontal, lwd = lwd, 
                lty = lty, col.line = col.line, ...)
    }
}

partial <- function(quiet=F,partialrees=c(NULL,NULL)){
	cat(sep="\n",paste("Choose two factors to make an interaction group, to use 
*along with* the original two?\n(press Enter to exit)",
	get.menu(names(datafile)),sep=""))
    repeat {partial.pick <- try(as.numeric(scan(what="character",nmax=2,quiet=T)),silent=T)
    partialrees <- names(datafile)[partial.pick]
    if (length(partialrees)==2 & !NA %in% partialrees)
    		if (!identical(partialrees[1],partialrees[2])) break
    cat(sep="\n","Invalid entry, please try again.")}
    x <- datafile[,partialrees[1]]
    y <- datafile[,partialrees[2]]
    if (!(class(x)%in%c("factor","character") & class(y)%in%c("factor","character"))){
    	cat(sep="\n","Those are not both factors!")
    	return()}
    fac.x <- levels(x)
	cat(sep="\n",paste("Factor(s) from ",partialrees[1]," to cross with factor(s) from ",partialrees[2],"?",get.menu(fac.x),sep="")) 
	repeat {x.pick <- 
		try(as.numeric(scan(what="character",nmax=length(fac.x),quiet=T)),silent=T)
		x.cross <- fac.x[x.pick]
		if (!NA %in% x.cross) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(x.cross) == 0) return()
	fac.y <- levels(y)
	cat(sep="\n",paste("Factor(s) from ",partialrees[2]," to cross with factor(s) from ",partialrees[1],"?",get.menu(fac.y),sep="")) 
	repeat {y.pick <- 
		try(as.numeric(scan(what="character",nmax=length(fac.y),quiet=T)),silent=T)
		y.cross <- fac.y[y.pick]
		if (!NA %in% y.cross) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(y.cross) == 0) return()
	cat(sep="\n", paste("Name for this interaction, or press Enter to use",
		paste(paste(x.cross, collapse=""), paste(y.cross, collapse=""), sep="_"),"?"))
	interaction.pick <- try(scan(what="character", nmax=1, quiet=T), silent=T)
		if (length(interaction.pick) == 0)
			interaction.name <- paste(paste(x.cross,collapse=""), paste(y.cross,collapse=""), sep="_")
		else
			interaction.name <- interaction.pick
	cat(sep="\n", paste("Name for the other levels, or press Enter to use 'other'?"))
	other.pick <- try(scan(what="character", nmax=1, quiet=T), silent=T)
		if (length(other.pick) == 0)
			other.name <- "other"
		else
			other.name <- other.pick
	datafile[,paste(partialrees,collapse="_")] <<-
		ifelse(x %in% x.cross & y %in% y.cross, interaction.name, other.name)}
	
	
plotting <- function(){
plotting.pick <<- 99
while (!plotting.pick %in% c(9,0)){
	cat("\n")
	cat(sep="\n","PLOTTING MENU")
	cat(sep="\n","1-custom scatterplot")
	if (exists("datafile"))
	cat(sep="\n","5-modeling 7-save plot 9-main menu 0-exit")
	else
	cat(sep="\n","7-save plot 9-main menu 0-exit")
	plotting.pick <<- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
	if (length(plotting.pick)==0) plotting.pick <<- 0
	if (plotting.pick == 1) custom.plot()
	if (plotting.pick == 5) modeling()
	if (plotting.pick == 7) cat(sep="\n","Ha! Not yet implemented!")
	if (plotting.pick == 9) modeling.pick <<- 9
	if (plotting.pick == 0) {modeling.pick <<- 0; main.pick <<- 0}}}
	
## HOW TO GET INTERACTIONS OUT
read.variables <- function(){
	model <- model[[1]]
	if ("glm" %in% class(model)){
		fam <<- as.character(model$family$family)
		response <<- names(model$model)[1]
		if (fam == "binomial"){
			non.apps <<- levels(model$model[,1])[1]
			apps <<- levels(model$model[,1])[2]}
		predictors <<- names(model$model)[-1]
		continuous <<- predictors[which(sapply(model$model,class)[-1]%in%c("integer","numeric"))]
		random <<- character()}
	if ("mer" %in% class(model)){
		fam <<- ifelse(length(model@muEta)>0,"binomial","gaussian")
		response <<- names(model@frame)[1]
		if (fam == "binomial"){
			non.apps <<- levels(model@frame[,1])[1]
			apps <<- levels(model@frame[,1])[2]}
		predictors <<- names(model@frame)[-1]
		continuous <<- predictors[which(sapply(model@frame,class)[-1]%in%c("integer","numeric"))]
		random <<- names(model@flist)}}

reclass <- function(){
	data.structure()
	cat(sep="\n",paste("Change class of which variable?",
		get.menu(names(datafile)),sep=""))
	repeat {reclassee.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		reclassee <- names(datafile)[reclassee.pick]
		if (!NA %in% reclassee) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(reclassee)==0) return()
	repeat{
		cat(sep="\n",paste("Change",reclassee,
			"to which class? (f-factor c-continuous [integer/numeric])"))
		reclass.pick <- try(scan(what="character",nmax=1,quiet=T),silent=T)
		if (length(reclass.pick)==0) return()
		if (reclass.pick=="f"|reclass.pick=="F"){
			target <- try(as.factor(datafile[,reclassee]),silent=T)
			if (!"try-error"%in%class(target)) datafile[,reclassee] <<- target
				else cat(sep="\n","Failed to change class!")
			return()}
		if (reclass.pick=="c"|reclass.pick=="C"){
			target <- try(as.numeric(as.character(datafile[,reclassee])),silent=T)
			if (!"try-error"%in%class(target)) datafile[,reclassee] <<- target
				else cat(sep="\n","Failed to change class!")
			return()}
		cat(sep="\n","Invalid entry, please try again.")}}
	
recode <- function(){
	cat(sep="\n",paste("Factor group to recode? (press Enter to exit)",
		get.menu(names(datafile)),sep=""))
	repeat {recode.group.pick <-
		try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		recode.group <- names(datafile)[recode.group.pick]
		if (!NA %in% recode.group) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(recode.group) == 0) return()
	if (class(datafile[,recode.group])!="factor"){
		cat(sep="\n",paste("Please change class of",recode.group,"to 'factor' before recoding."))
		return()}
	fac <- levels(datafile[,recode.group])
	repeat{
	cat(sep="\n",paste("Factor(s) of ",recode.group," to recode together? (",
		get.menu(fac,paren=F)," Enter-done)",sep=""))
	repeat {source.factors.pick <- 
		try(as.numeric(scan(what="character",nmax=length(fac),quiet=T)),silent=T)
		source.factors <- fac[source.factors.pick]
		if (!NA %in% source.factors) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(source.factors) == 0) break
	sf <- paste(source.factors,sep="",collapse=" ")
	cat(sep="\n",paste("Recode ",sf," as what?",sep=""))
	target.factor <- scan(what="character",nmax=1,quiet=T)
	if (length(target.factor) == 0) break
	for (i in 1:length(fac))
		if (fac[i] %in% source.factors) fac[i] <- target.factor}
	cat(sep="\n","Recode to new column? (Yes-type new column name No-press Enter)")
	repeat {new.pick <- try(scan(what="character",nmax=1,quiet=T),silent=T)
		if (length(new.pick) == 0) {levels(datafile[,recode.group]) <<- fac; break}
		a <- try(datafile[,new.pick] <<- datafile[,recode.group],silent=T)
		if (!"try-error"%in%class(a)) {levels(datafile[,new.pick]) <<- fac; break} 
		cat(sep="\n","Invalid column name.")}}

relev <- function(){
	cat(sep="\n",paste("Factor group to relevel? (press Enter to exit)",
		get.menu(names(datafile)),sep=""))
	repeat {relev.group.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		relev.group <- names(datafile)[relev.group.pick]
		if (!NA %in% relev.group) break
		if (class(datafile[,relev.group])!="factor"){
			cat(sep="\n",paste("Please change class of",relev.group,"to 'factor' before releveling."))
			break}
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(relev.group) == 0) return()
	fac <- levels(datafile[,relev.group])
	cat(sep="\n",paste("For",relev.group,"current baseline is",fac[1]))
	cat(sep="\n",paste("Choose new baseline (by number)?",get.menu(fac)))
	repeat {relev.pick <- 
		try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		baseline <- fac[relev.pick]
		if (!NA %in% baseline) break
		cat ("Invalid entry, please try again.")}
	if (length(baseline) == 0) return()
	datafile[,relev.group] <<- relevel(datafile[,relev.group],baseline)}

rename <- function(){
	cat(sep="\n",paste("Variable(s) to rename? (press Enter to exit)",
		get.menu(names(datafile)),sep=""))
	repeat {rename.group.pick <- try(as.numeric(scan(what="character",nmax=length(names(datafile)),quiet=T)),silent=T)
		rename.group <- names(datafile)[rename.group.pick]
		if (!NA %in% rename.group) break
		cat(sep="\n","Invalid entry, please try again.")}
		if (length(rename.group) == 0) return()
		for (group in rename.group){
			cat(sep="\n",paste("Type new name for ",group,"?",sep=""))
			repeat{
				new.name <- scan(what="character",nmax=1,quiet=T)
				t <- try((names(datafile)[which(names(datafile)==group)] <<- new.name),silent=T)
				if ("try-error" %in% class(t))
					cat(sep="\n","Invalid entry, please try again.")
					else break}}}

reset <- function(){
	vars <- c("silent",
		"verbose",
		"gold",
		"centered",
		"blups",
		"thresh",
		"sim",
		"decimals",
		"interactions",
		"contr",
		"main.pick",
		"datafile",
		"fullpath",
		"model",
		"continuous",
		"continuous.pick",	
		"sepchar",
		"p.flag",
		"predictors",
		"predictors.pick",
		"random",
		"response",
		"response.pick",
		"fam",
		"apps",	
		"non.apps",	
		"proportion",
		"fixed",
		"adjusting.pick",
		"modeling.pick",
		"variables",
		"threshold",
		"up",
		"up.preds",
		"up.model",
		"down",	
		"down.preds",	
		"down.model",
		"one",
		"one.preds",
		"one.model",
		"loaded.flag",
		"combo",
		"old.binary")
	for (v in vars)
		try(rm(list=eval(v),envir=globalenv()),silent=T)
	source("http://www.danielezrajohnson.com/Rbrul.R")
	rbrul()
	main.pick <<- 0}

retain <- function(){
	cat(sep="\n",paste("Factor group to retain using?",
		get.menu(names(datafile)),sep=""))
	repeat {retain.group.pick <- 
try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		retain.group <- names(datafile)[retain.group.pick]
		if (!NA %in% retain.group) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(retain.group) == 0) return()
	if (class(datafile[,retain.group])!="factor"){
		cat(sep="\n",paste("Please change class 
of",retain.group,"to 'factor' before retaining."))
		return()}
	fac <- levels(datafile[,retain.group])
	cat(sep="\n",paste("Factors to retain from 
",retain.group,"?",get.menu(fac),sep="")) 
	repeat {retained.factors.pick <- 
		try(as.numeric(scan(what="character",nmax=length(fac),quiet=T)),silent=T)
		retained.factors <- fac[retained.factors.pick]
		if (!NA %in% retained.factors) break
		cat(sep="\n","Invalid entry, please try again.")}
	if (length(retained.factors) == 0) return()
	datafile <<- subset(datafile,datafile[,retain.group]%in%retained.factors)
	datafile[,retain.group] <<- datafile[,retain.group][drop=T]}

rownames.off <- function (x, ..., digits = NULL, quote = FALSE, right = TRUE, 
    row.names = FALSE) 
{
    n <- length(row.names(x))
    if (length(x) == 0L) {
        cat(gettextf("data frame with 0 columns and %d rows\n", 
            n))
    }
    else if (n == 0L) {
        print.default(names(x), quote = FALSE)
        cat(gettext("<0 rows> (or 0-length row.names)\n"))
    }
    else {
        m <- as.matrix(format.data.frame(x, digits = digits, 
            na.encode = FALSE))
        if (!isTRUE(row.names)) 
            dimnames(m)[[1L]] <- if (identical(row.names, FALSE)) 
                rep.int("", n)
            else row.names
        print(m, ..., quote = quote, right = right)
    }
    invisible(x)
}

save.data <- function(){
	cat("\n")
        if (!exists("datafile")) {
		cat(sep="\n","No data loaded.")
                cat("\n")}
	else cat(sep="\n",paste("Current data file is:",fullpath))
	data.structure()
	cat(sep="\n","Save current data to a file? (Yes-type filename No-press Enter)")
	repeat {save.pick <- try(scan(what="character",nmax=1,quiet=T),silent=T)
		if (length(save.pick) == 0) break
		a <- try(write.table(datafile,save.pick,row.names=F,sep=sepchar,quote=F),silent=T)
		if (!"try-error"%in%class(a)){
			if (length(grep("/",save.pick))>0)
			fullpath <<- save.pick
			else
			fullpath <<- paste(shortpath,save.pick,sep="")
			break}
		cat(sep="\n","Error in saving file.")}}

set.variables <- function(){
	type <- "continuous"
	if (length(fam)>0)
		if (fam=="binomial") type <- "binary"
		variables <<- list(response=response,
			fixed.factor=predictors[!predictors%in%continuous & !predictors%in%random],
			fixed.continuous=continuous,
			fixed.interaction=interactions,
			random.intercept=random)
		names(variables)[1] <<- paste("response",type,sep=".")
		blanks <- vector()
		for (j in 1:length(variables))
			if (length(variables[[j]])==0) blanks <- c(blanks,j)
		if (length(blanks)>0) variables <<- variables[-blanks]}

sim.fz <- function(mod0,mod1,threshold=.05,min=1,max=1000,confidence=.01,mix=F){
	t0 <- attr(attr(terms(mod0),"factors"),"dimnames")[[2]]
	t1 <- attr(attr(terms(mod1),"factors"),"dimnames")[[2]]
	if (length(t0) < length(t1)){
		m0 <- mod0; m1 <- mod1 
		term <- which(!t1%in%t0)
		index <- which(attr(m1@X,"assign")==term)}
	if (length(t0) > length(t1)){
		m0 <- mod1; m1 <- mod0
		term <- which(!t0%in%t1)
		index <- which(attr(m1@X,"assign")==term)}
	if (length(t0) == length(t1)) return(1)
	if (any(length(m0@muEta)))
		fz.obs <- as.vector(fixef(m1)[index] %*% 
			solve(vcov(m1)[index,index]) %*% 
			fixef(m1)[index]) else { 
		ss <- sum((.Call(lme4:::mer_update_projection, m1)[[2]][index])^2)
        df <- length(index)
        ms <- ss/df
        fz.obs <- ms/(lme4:::sigma(m1)^2)}
	fz <- vector()
	run <- 1
	repeat{
		if (mix & contr=="contr.sum") d <- dan.mer.sim(m1,index=index)
		else d <- dan.mer.sim(m0)
		s1 <- try(refit(m1,d),silent=T)
		if ("try-error"%in%class(s1)) next
		if (any(length(m0@muEta)))
			fz[run] <- as.vector(fixef(s1)[index] %*% 
				solve(vcov(s1)[index,index]) %*%
				fixef(s1)[index]) else {
			ss <- sum((.Call(lme4:::mer_update_projection, s1)[[2]][index])^2)
        	df <- length(index)
        	ms <- ss/df
        	fz[run] <- ms/(lme4:::sigma(s1)^2)}
		a <- mean(fz>=fz.obs)
		cat(sep="\n",paste("p =",round(a,decimals),"in",run,"runs"))
		if (a > threshold) 
			b <- ifelse(binconf(length(fz[fz>=fz.obs]),run,confidence)[,"Lower"]>threshold,T,F)
		if (a < threshold) 
			b <- ifelse(binconf(length(fz[fz>=fz.obs]),run,confidence)[,"Upper"]<threshold,T,F)
		if (a == threshold) b <- F
		if ((b & run>=min) | run==max) return(a)
		run <- run + 1}}

split.term <- function(interaction){
	colon <- regexpr(":",interaction)
	result <- vector()
	index <- 1
	result[index] <- substr(interaction,1,colon-1)
	if (nchar(result[1])>0) index <- 2
	result[index] <- substr(interaction,colon+1,nchar(interaction))
result}

star <- function(){
	cat(sep="\n",paste("Choose two factors to make an interaction group, to use 
*instead* of the original two?\n(press Enter to exit)",
	get.menu(names(datafile)),sep=""))
    repeat {star.pick <- try(as.numeric(scan(what="character",nmax=2,quiet=T)),silent=T)
    starrees <- names(datafile)[star.pick]
    if (length(starrees)==2 & !NA %in% starrees)
    		if (!identical(starrees[1],starrees[2])) break
    cat(sep="\n","Invalid entry, please try again.")}
    x <- datafile[,starrees[1]]
    y <- datafile[,starrees[2]]
    if (class(x)%in%c("factor","character") & class(y)%in%c("factor","character"))
    	datafile[,paste(starrees,collapse=".")] <<- as.factor(paste(x,y,sep=".")) 
    else
    	cat(sep="\n","Those are not both factors!")}

step.down <- function(){
	if (!exists("variables")) return(cat(sep="\n","Please choose variables first."))
	load.variables()
	fix <- fixed
	cont <- continuous
	int <- interactions
	in.preds <- c(fixed,continuous,interactions)
	keep.going <- T
	steppe <- 0
	run <- 1
	best.model <- fit.model()
	if ("try-error" %in% class(best.model)) error.count <- 1 else error.count <- 0
	best.run <- 0
	if (!silent){
	cat("\n")
	cat(sep="\n","STEPPING DOWN...")
	cat("\n")
	cat(sep="\n",paste("STEP 0 - Run 0 - full model with",paste(in.preds,collapse=" + ")))
	cat("\n")
	if (verbose) print(extract(best.model),row.names=F)}
	while (keep.going & length(in.preds)>0){
		steppe <- steppe + 1
		if (thresh == "A") threshold <<- 0.05/length(in.preds) else threshold <<- thresh
		compare <- vector()
		mod <- list()
		mod.run <- vector()
		if (!silent) {cat("\n")
		cat(sep="\n",paste("STEP",steppe,"- dropping from Run",best.run,"model with",paste(in.preds,collapse=" + ")))}
		if (length(in.preds)>0)
		for (i in 1:length(in.preds)){
			if (in.preds[i]%in%sapply(in.preds[in.preds%in%int],split.term)){
				mod[[i]] <- "interaction"; class(mod[[i]]) <- "try-error"
				compare[i] <- get.p(best.model,mod[[i]])
				next}
			cat("\n")
			cat(sep="\n",paste("Trying without ",in.preds[i],"...",sep=""))
			selection <- in.preds[-i]
			continuous <<- cont[cont%in%selection]
			fixed <<- fix[fix%in%selection]
			interactions <<- int[int%in%selection]
			mod[[i]] <- fit.model()
			mod.run[i] <- run
			compare[i] <- get.p(best.model,mod[[i]])
			if (!silent) {cat("\n")
			if (verbose) print(extract(mod[[i]]),row.names=F)
			above.preds <- ifelse(length(in.preds[-i])==0,
				"no fixed predictors",
				paste(in.preds[-i],collapse=" + "))
			cat(sep="\n",paste("Run ",run," (above) with ",above.preds,
				" is worse than Run ",best.run,
				ifelse(above.preds=="no fixed predictors",""," also")," including ",in.preds[i],
				", ",p.flag," ",signif(compare[i],decimals),sep=""))}
			run <- run + 1}
		repeat{
		if (all(compare %in% c(-888, -999))){
			if(any(in.preds%in%interactions))
				best.model <-	mod[[xxx <- (which(in.preds%in%interactions)[1])]]
			else
				best.model <- mod[[xxx <- 1]]
			best.run <- mod.run[xxx]
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No successful run last step or this step, dropping first free term ",
					in.preds[xxx],"...",sep=""))}
			in.preds <- in.preds[-xxx]
			error.count <- error.count + 1
			break}
		if (all(compare %in% c(888, 999))){
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No successful run this step, best model from last step is Run ",best.run,sep=""))}
			keep.going <- F
			error.count <- error.count + 1
			break}
		if (max(compare)<0){
			best.model <- mod[[yyy <- which(compare==max(compare))[1]]]
			best.run <- mod.run[yyy]
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No successful run last step, best model this step is Run ",best.run,", dropping ",
					in.preds[yyy],"...",sep=""))}
			in.preds <- in.preds[-yyy]
			error.count <- error.count + 1
			break}
		if ((zzz <- max(compare[compare<888])) >= threshold){
				best.model <- mod[[yyy <- which(compare==zzz)[1]]]
				best.run <- mod.run[yyy]
				if (!silent) {
					cat("\n")
					cat(sep="\n",paste("Dropping ",in.preds[yyy],"...",sep=""))}
				in.preds <- in.preds[-yyy]
				break} 
		if (zzz < threshold){
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("All remaining predictors are significant, best model from last step is Run ",best.run,sep=""))}
			keep.going <- F	
			break}
		error <- c(success,ful)}}
	if (length(in.preds)==0){
		in.preds <- "no fixed predictors"
		compare <- 1234}
	down0 <- paste(ifelse(length(random)>0,paste(random," [random] and ",sep="",collapse=""),""),
		paste(in.preds[order(compare)]," (",
		ifelse(compare==1234,"",signif(compare[order(compare)],decimals)),")",
		sep="",collapse=" + "),"\n[p-values dropping from full model]",sep="")
	down1 <- gsub("\\(999\\)", "[error]", down0)
	down2 <- gsub("\\(-999\\)", "[error]", down1)
	down3 <- gsub("\\(888\\)", "[main effect]", down2)
	down4 <- gsub("\\(-888\\)", "[main effect]", down3)
	down <<- down4
	down.preds <<- in.preds
	fixed <<- fix
	continuous <<- cont
	interactions <<- int
model <<- list(best.model,down,contr)
down.model <<- best.model
if (!silent) {
	cat("\n")
	cat(sep="\n",paste("BEST STEP-DOWN MODEL IS WITH",down))
	if (verbose) {
		cat("\n"); print(extract(down.model),row.names=F)}
	if (error.count > 0)
		cat(sep="\n",paste("THERE WERE ERRORS ON",error.count,"OF",steppe+1,"STEPS!"))}
down.model}

step.one <- function(){
	if (!exists("variables")) return(cat(sep="\n","Please choose variables first."))	
	load.variables()
	fix <- fixed
	cont <- continuous
	int <- interactions
	one.preds <<- c(fixed,continuous,interactions)
	in.preds <- c(fixed,continuous,interactions)
	one.model <<- fit.model()
	compare <- vector()
	mod <- list()
	if (length(in.preds)>0)
	for (i in 1:length(in.preds)){
		if (in.preds[i]%in%sapply(in.preds[in.preds%in%int],split.term)){
			mod[[i]] <- "interaction"; class(mod[[i]]) <- "try-error"
			compare[i] <- get.p(one.model,mod[[i]])
			next}
		selection <- in.preds[-i]
		continuous <<- cont[cont%in%selection]
		fixed <<- fix[fix%in%selection]
		interactions <<- int[int%in%selection]
		mod[[i]] <- fit.model()
		compare[i] <- get.p(one.model,mod[[i]])}
	one0 <- paste(ifelse(length(random)>0,paste(random," [random] and ",sep="",collapse=""),""),
               	paste(in.preds[order(compare)]," (",signif(compare[order(compare)],decimals),")",
                sep="",collapse=" + "),sep="")
    one1 <- gsub("\\(999\\)", "[error]", one0)
	one2 <- gsub("\\(-999\\)", "[error]", one1)
	one3 <- gsub("\\(888\\)", "[main effect]", one2)
	one4 <- gsub("\\(-888\\)", "[main effect]", one3)
	one <<- one4
	fixed <<- fix
	continuous <<- cont
	interactions <<- int
	if (!silent){
	cat("\n")
	cat(sep="\n",paste("ONE-LEVEL ANALYSIS WITH",one))
	if (verbose) print(extract(one.model),row.names=F)}
model <<- list(one.model,one,contr)
one.model}

step.up <- function(){
	if (!exists("variables")) return(cat(sep="\n","Please choose variables first."))
	load.variables()
	fix <- fixed
	cont <- continuous
	int <- interactions
	out.preds <- c(fixed,continuous,interactions)
	in.preds <- "no fixed predictors"
	p.add <- vector()
	keep.going <- T
	steppe <- 0
	run <- 1
	fixed <<- 1
	continuous <<- 1
	interactions <<- 1
	best.model <- fit.model()
	if ("try-error" %in% class(best.model)) error.count <- 1 else error.count <- 0
	best.run <- 0
	if (!silent){
	cat("\n")
	cat(sep="\n","STEPPING UP...")
	cat("\n")
	cat(sep="\n","STEP 0 - Run 0 - model with no fixed predictors")
	cat("\n")
	if (verbose) print(extract(best.model),row.names=F)}
	while (keep.going & length(out.preds)>0){
		steppe <- steppe + 1
		if (thresh == "A")
			threshold <<- 0.05/length(out.preds) else threshold <<- thresh
		compare <- vector()
		mod <- list()
		mod.run <- vector()
		if (!silent){
		cat("\n")
		cat(sep="\n",paste("STEP",steppe,"- adding to Run",best.run,"model with",paste(in.preds,collapse=" + ")))}
		for (i in 1:length(out.preds)){
			if (out.preds[i]%in%int & !all(split.term(out.preds[i])%in%in.preds)){
				mod[[i]] <- "interaction"; class(mod[[i]]) <- "try-error"
				compare[i] <- get.p(best.model,mod[[i]])
				next}
			cat("\n")
			cat(sep="\n",paste("Trying with ",out.preds[i],"...",sep=""))
			selection <- c(in.preds,out.preds[i])
			continuous <<- cont[cont%in%selection]
			fixed <<- fix[fix%in%selection]
			interactions <<- int[int%in%selection]
			mod[[i]] <- fit.model()
			mod.run[i] <- run
			compare[i] <- get.p(best.model,mod[[i]])
			if (!silent) {cat("\n")
			if (verbose) print(extract(mod[[i]]),row.names=F)
			above.preds <- ifelse(all(in.preds=="no fixed predictors"),
				out.preds[i],
				paste(paste(in.preds,collapse=" + "),out.preds[i],sep=" + "))
			cat(sep="\n",paste("Run ",run," (above) with ",above.preds,
				" is better than Run ",best.run," without ",out.preds[i],", ",
				p.flag," ",signif(compare[i],decimals),sep=""))}
			run <- run + 1}
	repeat{
		if (all(compare %in% c(-888, -999))){
			if (any(!out.preds%in%int))
				best.model <- mod[[xxx <- (which(!out.preds%in%int)[1])]]
			else
				best.model <- mod[[xxx <- 1]]
			best.run <- mod.run[xxx]
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No successful run last step or this step, adding first free term ",
					out.preds[xxx],"...",sep=""))}
			in.preds <- in.preds[in.preds!="no fixed predictors"]
			in.preds <- c(in.preds,out.preds[xxx])
			out.preds <- out.preds[!out.preds%in%in.preds]
			error.count <- error.count + 1
			break}
		if (all(compare %in% c(888, 999))){
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No successful run this step, best model from last step is Run 
",best.run,sep=""))}
			keep.going <- F
			error.count <- error.count + 1
			break}
		if (max(compare)<0){
			best.model <- mod[[yyy <- which(compare==max(compare))[1]]]
			best.run <- mod.run[yyy]
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No successful run last step, best model this step is Run ",best.run,", adding ",
					out.preds[yyy],"...",sep=""))}
			in.preds <- c(in.preds,out.preds[yyy])
			out.preds <- out.preds[!out.preds%in%in.preds]
			error.count <- error.count + 1
			break}			
		if (min(compare)<threshold){
			best.model <- mod[[zzz <- (which(compare==min(compare))[1])]]
			best.run <- mod.run[zzz]
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("Adding ",out.preds[zzz],"...",sep=""))}
			in.preds <- in.preds[in.preds!="no fixed predictors"]
			in.preds <- c(in.preds,out.preds[zzz])
			out.preds <- out.preds[!out.preds%in%in.preds]
			p.add <- c(p.add,signif(min(compare),decimals))
			break}
		if (min(compare)>=threshold){
			if (!silent) {
				cat("\n")
				cat(sep="\n",paste("No significant improvement this step, best model from last step is Run ",best.run,sep=""))}
			keep.going <- F
			break}
		error <- c(success,ful)}}	
	if (length(p.add)==0) p.add <- 1234
	up0 <- paste(ifelse(length(random)>0,
			paste(random," [random] and ",sep="",collapse=""),""),
			paste(in.preds," (",ifelse(p.add==1234,"",signif(p.add,decimals)),")",
			sep="",collapse=" + "),"\n[p-values building from null model]",sep="")
	up1 <- gsub("\\(999\\)", "[error]", up0)
	up2 <- gsub("\\(-999\\)", "[error]", up1)
	up3 <- gsub("\\(888\\)", "[main effect]", up2)
	up4 <- gsub("\\(-888\\)", "[main effect]", up3)
	up <<- up4
	fixed <<- fix
	continuous <<- cont
	interactions <<- int
	up.preds <<- in.preds
model <<- list(best.model,up,contr)
up.model <<- best.model
if (!silent) {cat("\n")
cat(sep="\n",paste("BEST STEP-UP MODEL IS WITH",up))
if (verbose) {cat("\n"); print(extract(up.model),row.names=F)}
if (error.count > 0)
	cat(sep="\n",paste("THERE WERE ERRORS ON",error.count,"OF",steppe+1,"STEPS!"))}
up.model}

step.updown <- function(){
	if (!exists("variables")) return(cat(sep="\n","Please choose variables first."))
	step.up()
	step.down()
	cat("\n")
	if (!silent){
		cat(sep="\n", paste("BEST STEP-UP MODEL IS WITH", up))
		cat(sep="\n")
		cat(sep="\n", paste("BEST STEP-DOWN MODEL IS WITH", down))
		cat("\n")
	if (all(up.preds%in%down.preds) & all(down.preds%in%up.preds))
		cat(sep="\n","STEP-UP AND STEP-DOWN MATCH!\n") else 
		cat(sep="\n","STEP-UP AND STEP-DOWN MISMATCH!!\n")}}

sumify <- function(data,ncol){
	data <- matrix(data,ncol=ncol)
	for (r in 1:(nrow(data)-1))
		data[r,which(is.na(data[r,]))] <- -sum(data[r,],na.rm=T)
	if (nrow(data)>1) data[nrow(data),] <- -colSums(data,na.rm=T)
as.vector(data)}

test.chisq <- function(){
	cat(sep="\n",paste("Enter first deviance or log likelihood."))
	repeat {
		d1 <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (!"try-error"%in%class(d1)) break
		cat(sep="\n",paste("Invalid entry, please try again."))}
	cat(sep="\n",paste("Enter second deviance or log likelihood."))
	repeat {
		d2 <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (!"try-error"%in%class(d2)) break
		cat(sep="\n",paste("Invalid entry, please try again."))}
	cat(sep="\n",paste("If these were log likelihood values, press 1. Press Enter if they were deviances."))
	repeat {
		ll.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (length(ll.pick)==0) break
		if (ll.pick==1) {d1<-2*d1; d2<-2*d2; break}
		cat(sep="\n",paste("Invalid entry, please try again."))}
	cat(sep="\n",paste("Enter difference in degrees of freedom."))
	repeat {
		df <- try(as.integer(scan(what="character",nmax=1,quiet=T)),silent=T)
		if (!"try-error"%in%class(df)) break
		cat(sep="\n",paste("Invalid entry, please try again."))}
	xx <- abs(d1-d2)
	cat(sep="\n",paste("Chi-square = ",xx,", df = ",df,", p = ",signif(pchisq(xx,df,lower=F),decimals),sep=""))}
	
	
transform <- function(){
        cat(sep="\n",paste("Continuous variable(s) to center or transform? (press Enter to exit)",
        	get.menu(names(datafile)),sep=""))
        repeat {transform.pick <- try(as.numeric(scan(what="character",
		nmax=length(names(datafile)),quiet=T)),silent=T)
                transformees <- names(datafile)[transform.pick]
                if (!NA %in% transformees) break
                cat(sep="\n","Invalid entry, please try again.")}
        if (length(transformees)>0)
        for (i in transformees){
        		if (!class(datafile[,i])%in%c("integer","numeric")){
        			cat(sep="\n",paste("Please change class of",i,"to 'continuous' before transforming."))
        			next}
                cat(sep="\n",paste("Do what to ",i,
                	"? (1-log 2-logit 3-center (rezero) 4-z-score Enter-nothing)",sep=""))
                repeat {trans.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
                        if (length(trans.pick)==0) break
                        if (trans.pick %in% 1:6) break
                        cat(sep="\n","Invalid entry, please try again.")}
                if (length(trans.pick)==0) break
                if (trans.pick == 1){
                        trans.col <- paste("log",i,sep=".")
                        datafile[,trans.col] <<- log(datafile[,i])}
                if (trans.pick == 2){
                        trans.col <- paste("logit",i,sep=".")
                        if (max(datafile[,i],na.rm=T)<=1 & min(datafile[,i],na.rm=T)>=0){
                        datafile[,trans.col] <<- logit(datafile[,i])} else
                        datafile[,trans.col] <<- NA}
                if (trans.pick == 3){
                        trans.col <- paste("ctr",i,sep=".")
                        cat(sep="\n",paste("Center (rezero) on what value? (mean is ",signif(mean(datafile[,i],na.rm=T),5),")",sep=""))
                        repeat {center.pick <- try(as.numeric(scan(what="character",nmax=1,quiet=T)),silent=T)
                                if (is.real(center.pick)&length(center.pick)>0) break
                               	cat(sep="\n","Invalid entry, please try again.")}
                        datafile[,trans.col] <<- datafile[,i]-center.pick}
                if (trans.pick == 4){
                        trans.col <- paste("z",i,sep=".")
                        datafile[,trans.col] <<- (datafile[,i]-mean(datafile[,i],na.rm=T))/
				sd(datafile[,i],na.rm=T)}}}
				
trim <- function(){
	if (!exists("model")) return(cat(sep="\n","No model loaded."))
	orig.datafile <- datafile
	orig.model <- model
	ff <- vector()
	if (length(fixed) == 0)
		fixed3 <- character() 
	if (length(fixed) > 0)
		for (f in fixed){
			ff <- c(ff, length(levels(datafile[,f])))
			fixed3 <- fixed[which(ff > 2)]}
	if (length(fixed3) == 0){
		cat(sep="\n")
		cat(sep="\n", "This model has no fixed-effect predictors (with 3 or more levels) to trim.")
		return()}
	big.model <- step.one()
	if (length(fixed3) == 1){
		fixeds <- fixed3}
	if (length(fixed3) > 1){
		cat(sep="\n")
		cat(sep="\n",
			paste("Choose predictors to combine levels within; press Enter for all",
			get.menu(fixed3), ".", sep=""))
	repeat {fixed.pick <- try(as.numeric(scan(what="character",nmax=length(fixed3),quiet=T)),silent=T)
		if (length(fixed.pick) == 0) {fixeds <- fixed3; break}
		if (!NA %in% fixed3[fixed.pick]) {fixeds <- fixed3[fixed.pick]; break}
		cat(sep="\n","Invalid entry, please try again.")}}
	chooselevels <- vector()
	for (f in fixeds){
		for (l in levels(datafile[,f])){
			chooselevels <- c(chooselevels, paste(f, l, sep="."))}}
	cat(sep="\n")
	cat(sep="\n",
		paste("Choose levels to consider combining; press Enter for all",
		get.menu(chooselevels), ".", sep=""))
	repeat {levels.pick <- 
		try(as.numeric(scan(what="character",nmax=length(chooselevels),quiet=T)),silent=T)
		if (length(levels.pick) == 0) break
		if (length(levels.pick) == 1) {cat(sep="\n", "Choose at least two levels."); next} 	
		if (!NA %in% chooselevels[levels.pick]) {chooselevels <- chooselevels[levels.pick]; break}
		cat(sep="\n", "Invalid entry, please try again.")}
	for (loop in 1:1000){
		ff <- vector()
		for (f in fixeds)
			ff <- c(ff, length(levels(datafile[,f])))
		fixeds <- fixeds[which(ff > 2)]
		if (length(fixeds) == 0){
			cat(sep="\n")
			cat(sep="\n", "This model has no fixed-effect predictors (with 3 or more levels) to trim.")
			break}
		pairframe <- data.frame()
		choices <- data.frame()
		for (f in fixeds){
			l <- levels(datafile[,f])
			choices <- rbind(choices, cbind(left=combn(l, 2)[1,], right=combn(l, 2)[2,]))
			for (ch in 1:length(choices[,1])){
				if (paste(f, choices$left[ch], sep=".") %in% chooselevels &
					paste(f, choices$right[ch], sep=".") %in% chooselevels){
					pairframe <- rbind(pairframe,
						cbind(factor=f, left=as.character(choices$left[ch]),
							right=as.character(choices$right[ch])))}}}
		if (length(pairframe) == 0) break
		for (i in 1:length(pairframe[,1])){
			little.datafile <- datafile
			levs <- levels(little.datafile[, as.character(pairframe$factor[i])])
			leftnum <- which(levs == pairframe$left[i])
			rightnum <- which(levs == pairframe$right[i])
			levs[c(leftnum, rightnum)] <- 
				paste(pairframe$left[i], pairframe$right[i], sep="+")
			levels(little.datafile[, as.character(pairframe$factor[i])]) <- levs
			pairframe$p[i] <- round(get.p(big.model, fit.model(little.datafile), trim=T), decimals)}
		pairframe <- cbind(row = 1:length(pairframe[,1]), pairframe)
		pairframe2 <- cbind(pairframe, "     ", pairframe[,-1], pairframe[,1])
		names(pairframe2) <- c(names(pairframe), " ",
			paste(" ", c(names(pairframe)[-1], names(pairframe[1])), sep=""))
		names(pairframe2)[c(2, 7)] <- "predictor"
		for (pf in c(7, 8, 9, 11, 10))
			pairframe2[, pf] <- pairframe2[order(pairframe2[,10]), pf]
		cat(sep="\n")
		print(pairframe2)
		cat(sep="\n")
		cat(sep="\n", "Type a row number to combine a pair of levels ('left' and 'right'). Press Enter to stop.")
		repeat {combine.pick <- 
			try(as.numeric(scan(what="character", nmax=1, quiet=T)), silent=T)
			if (length(combine.pick) == 0) break
			if (combine.pick %in% pairframe$row) break
			cat(sep="\n","Invalid entry, please try again.")}
		if (length(combine.pick) == 0) break
		pair <- pairframe[combine.pick,]
		levs <- levels(datafile[, as.character(pair$factor)])
		leftnum <- which(levs == pair$left)
		rightnum <- which(levs == pair$right)
		cat(sep="\n")
		cat(sep="\n", paste("Type new name for combination or press Enter to use ",
			(lr <- paste(pair$left,pair$right,sep="+")),"?", sep=""))
		repeat {combo.pick <-
			try(scan(what="character", nmax=1, quiet=T), silent=T)
			if (length(combo.pick) == 0){
				levs[c(leftnum, rightnum)] <- lr
				chooselevels <- c(chooselevels, paste(pair$factor, lr, sep="."))
				break}
			else {tlr <- try(levs[c(leftnum, rightnum)] <- combo.pick)
				if (!"try-error" %in% class(tlr)){
					chooselevels <- c(chooselevels, paste(pair$factor,combo.pick,sep="."))
					break}}
			cat(sep="\n", "Invalid entry, please try again.")}
		levels(datafile[, as.character(pair$factor)]) <<- levs
		big.model <- step.one()}
	cat(sep="\n")
	cat(sep="\n", "Press 1 to keep trimmed data file and model, press Enter to keep originals.")
	repeat {keep.pick <-
		try(as.numeric(scan(what="character", nmax=1, quiet=T)), silent=T)
		if (length(keep.pick) == 0){
			datafile <<- orig.datafile 
			model <<- orig.model
			break}
		if (keep.pick == 1){
			break}
		cat(sep="\n", "Invalid entry, please try again.")}}